kuwahara filter (фильтр кувахары)

[mrgloom 242]

заимплементил фильтр кувахары(для одноканальных изображений),на выходе вокруг изображения черная рамка, убирать лень.

есть предположение что можно считать его быстрее, но пока не знаю как.

http://rsbweb.nih.gov/ij/plugins/kuwahara.html

пример работы

IplImage* kuwahara_grey(IplImage* img,int kernel_size)

{

	int nchannels= img->nChannels;

	int w=img->width;

	int h=img->height;

	IplImage* out= cvCreateImage(cvSize(w,h),IPL_DEPTH_8U,nchannels);


	int dt=kernel_size/2;


	//TODO: расширить изображение-убрать черные края.


	if(nchannels==1)  //только для одноканальных

	{

		for( int y=2*dt; y <h-2*dt ; ++y ) 

		{           

			uchar* ptr = (uchar*) (img->imageData + y*img->widthStep);

			uchar* p_out = (uchar*) (out->imageData + y*out->widthStep);

			for( int x=2*dt; x <w-2*dt; ++x )

			{   

				 int gr=0;

				 int sum=0;

				 double mean1=0;

				 double mean2=0;

				 double mean3=0;

				 double mean4=0;

				 double disp1=0;

				 double disp2=0;

				 double disp3=0;

				 double disp4=0;


				 int k=0;

				 //1

				 sum=0;

				 mean1=0;

				 uchar* pd=ptr;

				 for (int k1=0;k1<=dt;++k1)

				 {

					for(int k2=-dt;k2<=0;++k2)

					{

						sum+=pd[x+k2];

						++k;

					}

					pd+= img->widthStep;

				 }

				 mean1=sum/((dt+1)*(dt+1));


				 pd=ptr;

				 disp1=0;

				 for (int k1=0;k1<=dt;++k1)

				 {

					for(int k2=-dt;k2<=0;++k2)

					{

						disp1+=(pd[x+k2]-mean1)*(pd[x+k2]-mean1);

					}

					pd+= img->widthStep;

				 }


				 //2

				 sum=0;

				 mean2=0;

				 pd=ptr+dt*img->widthStep;

				 for (int k1=0;k1<=dt;++k1) 

				 {

					for(int k2=-dt;k2<=0;++k2)

					{

						sum+=pd[x+k2];

					}

					pd+= img->widthStep;

				 }

				 mean2=sum/((dt+1)*(dt+1));


				 pd=ptr+dt*img->widthStep;

				 disp2=0;

				 for (int k1=0;k1<=dt;++k1)

				 {

					for(int k2=-dt;k2<=0;++k2)

					{

						disp2+=(pd[x+k2]-mean2)*(pd[x+k2]-mean2);

					}

					pd+= img->widthStep;

				 }


				 //3

				 sum=0;

				 mean3=0;

				 pd=ptr;

				 for (int k1=0;k1<=dt;++k1)

				 {

					for(int k2=0;k2<=dt;++k2)

					{

						sum+=pd[x+k2];

					}

					pd+= img->widthStep;

				 }

				 mean3=sum/((dt+1)*(dt+1));


				 pd=ptr;

				 disp3=0;

				 for (int k1=0;k1<=dt;++k1)

				 {

					for(int k2=0;k2<=dt;++k2)

					{

						disp3+=(pd[x+k2]-mean3)*(pd[x+k2]-mean3);

					}

					pd+= img->widthStep;

				 }


				 //4

				 sum=0;

				 mean4=0;

				 pd=ptr+dt*img->widthStep;

				 for (int k1=0;k1<=dt;++k1)

				 {

					for(int k2=0;k2<=dt;++k2)

					{

						sum+=pd[x+k2];

					}

					pd+= img->widthStep;

				 }

				 mean4=sum/((dt+1)*(dt+1));


				 pd=ptr+dt*img->widthStep;

				 disp4=0;

				 for (int k1=0;k1<=dt;++k1)

				 {

					for(int k2=0;k2<=dt;++k2)

					{

						disp4+=(pd[x+k2]-mean4)*(pd[x+k2]-mean4);

					}

					pd+= img->widthStep;

				 }


				gr=0;

				if((disp1<=disp2)&&(disp1<=disp3)&&(disp1<=disp4))

					gr=mean1;

				if((disp2<=disp1)&&(disp2<=disp3)&&(disp2<=disp4))

					gr=mean2;

				if((disp3<=disp1)&&(disp3<=disp2)&&(disp3<=disp4))

					gr=mean3;	

				if((disp4<=disp1)&&(disp4<=disp2)&&(disp4<=disp3))

					gr=mean4;


				p_out[x]= gr;

			}

		}

	}

	return out;

}

[ngry 0]

есть предположение что можно считать его быстрее, но пока не знаю как.

Вот так быстрее считать дисперсию. Ёщё можно погуглить по запросу "Integral Image":

double m1 = 0.0, m2 = 0.0;

numType x;

for (i = 0; i < N; i++) {

   x = randNum();

   m1 += x; 

   m2 += x*x;

}


m1 /= N;

m2 /= N;


printf("       Average: %f\n", m1);

printf("Std. deviation: %f\n", sqrt(m2-m1*m1));

[Smorodov 579]

Интересно как в сравнении со встроенным OpenCV -шным:

meanstddev

[mrgloom 242]

так std.dev. же не x^2-m^2, а (x-m)^2.

или у интегрального изображения что то меняется?

[Smorodov 579]

Я встречал метод быстрого вычисления стандартного отклонения и среднего значения при помощи интегрального изображения, правда не помню где.

По поводу связи среднего и стандартного отклонения с моментами:

http://mathworld.wolfram.com/Moment.html

[agmt 0]

так std.dev. же не x^2-m^2, а (x-m)^2.

Так m - особое число и данные формулы равны: Sum(2*m*x)=2*m*Sum(x)=2*m*N*m=2*N*m**2

[mrgloom 242]

http://ru.wikipedia.org/wiki/%D0%94%D0%B8%D1%81%D0%BF%D0%B5%D1%80%D1%81%D0%B8%D1%8F_%D1%81%D0%BB%D1%83%D1%87%D0%B0%D0%B9%D0%BD%D0%BE%D0%B9_%D0%B2%D0%B5%D0%BB%D0%B8%D1%87%D0%B8%D0%BD%D1%8B

http://ru.wikipedia.org/wiki/%D0%9C%D0%B0%D1%82%D0%B5%D0%BC%D0%B0%D1%82%D0%B8%D1%87%D0%B5%D1%81%D0%BA%D0%BE%D0%B5_%D0%BE%D0%B6%D0%B8%D0%B4%D0%B0%D0%BD%D0%B8%D0%B5

ну матожидание линейно, но что то я всё равно не понял как они раскрыли квадрат модуля.

[Smorodov 579]

Из Бишопа (Pattern Recognition and Machine Learning. Bishop).

[mrgloom 242]

ну так я и не понял как делается expanding the square, вроде как должно 3 члена получится же?

[Smorodov 579]

Вот как оно делается:

[mrgloom 242]

хмм, да всё верно.

п.с. можно юзать латех

http://www.codecogs.com/latex/eqneditor.php

[Smorodov 579]

Да у меня есть чем редактировать, просто лень иногда