Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > sci.image.processing > #4202 > unrolled thread

problem about replace opencv dft() function with fftw api

Back to article view | Back to sci.image.processing

Contents

  problem about replace opencv dft() function with fftw api Zhi Sun <sunzhimailbox@gmail.com> - 2017-01-10 05:39 -0800
    Re: problem about replace opencv dft() function with fftw api Martin Brown <|||newspam|||@nezumi.demon.co.uk> - 2017-01-12 08:22 +0000

#4202 — problem about replace opencv dft() function with fftw api

Hi guy,

we are encountering an issue about programing fftw api to replace cv:dft to speed-up image/video process,but the result output of fftw looks like not correctly.

the background is that,we have a video analysis program to handle video frame one by one, saying that a serial of images (640x480), and we originally call cv::dft,that is slow. the video is 20 fps.

we have use opencv(2.4.13) sample code about dft, and implement a function (fftd) to replace the cv::dft as below.


is there anything we missed? ,maybe there are some misunderstanding about the FFT output, any suggestion would be appreciated!

--------------------------------------

cv::Mat fftd( cv::Mat img, bool backwards)
{

#if defined(USE_FFTW)	//not workable

	if (img.channels() == 1){
		cv::Mat planes[] = {cv::Mat_<float> (img), cv::Mat_<float>::zeros(img.size())};
		cv::merge(planes, 2, img);
	}
		
	if ( backwards){
		unsigned t[2];
		t[0] = Microseconds();	
		cv::dft(img, img, backwards ? (cv::DFT_INVERSE | cv::DFT_SCALE) : 0 );
		t[1] = Microseconds();
		printf( "cv:dft,  %6d usecs, img:%dx%d, channels:%d,backward:%d\n",t[1]-t[0],img.cols,img.rows,img.channels(),backwards);
		return img;
    }

	unsigned t[2];
	unsigned t1[2];
	unsigned t2[2];
	
	cv::Mat img1;
    cv::Mat img2;
	uchar *img1_data;
	uchar *img2_data;
	
	fftwf_complex    *data_in;
    fftwf_complex    *fft;           
    fftwf_plan       plan_f;
	
	int             width, height, step;
    int             i, j, k;
	
	int type = img.type();
	CV_Assert( type == CV_32FC1 || type == CV_32FC2 || type == CV_64FC1 || type == CV_64FC2 );
	
	t[0] = Microseconds();	
		
	img1 = img;
	img2 = img;
	
	width  	  = img.size().width;
    height 	  = img.size().height;
    step	  = img.step;
    img1_data = ( uchar* ) img1.data;
	img2_data = ( uchar* ) img2.data;
	//printf("width=%d,height=%d,step=%d\n",width,height,step);
	//width=34,height=34,step=272
	
	data_in = fftwf_alloc_complex(width * height);
    fft     = fftwf_alloc_complex(width * height);
	
	t1[0] = Microseconds();
	plan_f = fftwf_plan_dft_1d( width * height, data_in, fft, backwards ? FFTW_BACKWARD : FFTW_FORWARD,  FFTW_ESTIMATE );
	t1[1] = Microseconds();
	
	for( i = 0, k = 0 ; i < height ; i++ ) {
        for( j = 0 ; j < width ; j++ ) {
			float *p = (float*)&img1_data[i * step + j * sizeof(float)];
			data_in[k][0] = p[0];
			data_in[k][1] = p[1];
            //data_in[k][0] = ( float )img1_data[i * step + j];
            //data_in[k][1] = 0.0;
            k++;
        }
    }

	t2[0] = Microseconds();	
	fftwf_execute( plan_f );
	t2[1] = Microseconds();
	
	for( i = 0, k = 0 ; i < height ; i++ ) {
		for( j = 0 ; j < width ; j++ ) {
			float *p = (float*)&img2_data[i * step + j * sizeof(float)];
			p[0] = fft[k][0];
			p[1] = fft[k][1];
			k++;
		}
	}
	
	if (backwards) img2 *= 1.d / (float) (width * height);
	
    fftwf_destroy_plan( plan_f );
    fftwf_free( data_in );
    fftwf_free( fft );
	
	t[1] = Microseconds();
	//printf( "DFT,fftw   %6d usecs, plan:%6d usecs, exec:%6d usecs, img:%dx%d, channels:%d,backwards:%d\n",t[1]-t[0],t1[1]-t1[0],t2[1]-t2[0],img.cols,img.rows,img.channels(),backwards);
	
	img = img2;

#else

    if (img.channels() == 1)
    {
        cv::Mat planes[] = {cv::Mat_<float> (img), cv::Mat_<float>::zeros(img.size())};
        //cv::Mat planes[] = {cv::Mat_<double> (img), cv::Mat_<double>::zeros(img.size())};
        cv::merge(planes, 2, img);
    }
    cv::dft(img, img, backwards ? (cv::DFT_INVERSE | cv::DFT_SCALE) : 0 );

#endif

#if 1
	static int cnt = 0;
	printf("\nNo.%d\t",++cnt);
	for (int i = 0; i < img.rows; i++){
		for (int j = 0; j < img.cols; j++){
			if(i == 0 && j < 10)
				printf("(%.2f,%.2f) ",img.at<cv::Vec2d > (i, j)[0],img.at<cv::Vec2d > (i, j)[1]);
		}
	}
	printf("\n");
#endif

    return img;
}



int main(int argc, char ** argv)
{
    help(argv[0]);

    const char* filename = argc >=2 ? argv[1] : "lena.jpg";

    Mat I = imread(filename, CV_LOAD_IMAGE_GRAYSCALE);
    if( I.empty())
        return -1;

    Mat padded;                            //expand input image to optimal size
    int m = getOptimalDFTSize( I.rows );
    int n = getOptimalDFTSize( I.cols ); // on the border add zero values
    copyMakeBorder(I, padded, 0, m - I.rows, 0, n - I.cols, BORDER_CONSTANT, Scalar::all(0));

    Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
    Mat complexI;
    merge(planes, 2, complexI);         // Add to the expanded another plane with zeros

	#if 1
	complexI = fftd(complexI);			//the function we are testing!!!
	#else
    dft(complexI, complexI);            // this way the result may fit in the source matrix
	#endif

    // compute the magnitude and switch to logarithmic scale
    // => log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
    split(complexI, planes);                   // planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))
    magnitude(planes[0], planes[1], planes[0]);// planes[0] = magnitude
    Mat magI = planes[0];

    magI += Scalar::all(1);                    // switch to logarithmic scale
    log(magI, magI);

    // crop the spectrum, if it has an odd number of rows or columns
    magI = magI(Rect(0, 0, magI.cols & -2, magI.rows & -2));

    // rearrange the quadrants of Fourier image  so that the origin is at the image center
    int cx = magI.cols/2;
    int cy = magI.rows/2;

    Mat q0(magI, Rect(0, 0, cx, cy));   // Top-Left - Create a ROI per quadrant
    Mat q1(magI, Rect(cx, 0, cx, cy));  // Top-Right
    Mat q2(magI, Rect(0, cy, cx, cy));  // Bottom-Left
    Mat q3(magI, Rect(cx, cy, cx, cy)); // Bottom-Right

    Mat tmp;                           // swap quadrants (Top-Left with Bottom-Right)
    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);

    q1.copyTo(tmp);                    // swap quadrant (Top-Right with Bottom-Left)
    q2.copyTo(q1);
    tmp.copyTo(q2);

    normalize(magI, magI, 0, 1, CV_MINMAX); // Transform the matrix with float values into a
                                            // viewable image form (float between values 0 and 1).

    imshow("Input Image"       , I   );    // Show the result
    imshow("spectrum magnitude", magI);
    waitKey();

    return 0;
}



---------------------------------------

thanks
-zhi

[toc] | [next] | [standalone]

#4205

On 10/01/2017 13:39, Zhi Sun wrote:
> Hi guy,
>
> we are encountering an issue about programing fftw api to replace
> cv:dft to speed-up image/video process,but the result output of fftw
> looks like not correctly.
>

A couple of quick points from the code.

It looks like you are asking FFTW to do a 1d transform of size 640x480 
rather than a 2D transform 640,480. Is this what you intend?

Also when loading the complex array the image needs to go in the real 
components and the imaginary components identically zero. The complex 
array needs to be N+1 long to accommodate the alternating component.

> the background is that,we have a video analysis program to handle
> video frame one by one, saying that a serial of images (640x480), and
> we originally call cv::dft,that is slow. the video is 20 fps.

It will be a little bit harder to get going but using the Real to 
complex conjugate symmetric form of the FFTW library will give you 
transforms that can be done in place.

Also comment out the tranform and check that the conversion to and from 
the real array is actually working correctly. That code isn't pretty.

> we have use opencv(2.4.13) sample code about dft, and implement a
> function (fftd) to replace the cv::dft as below.
>
>
> is there anything we missed? ,maybe there are some misunderstanding
> about the FFT output, any suggestion would be appreciated!

I suspect FFTW is doing what you asked it to do - give or take any 
indexing errors that might have crept in elsewhere. I'd suggest testing 
the 1D transform on the first line of the image as a starting point now.
(after you have checked that the convertion to and from reals is OK)

Then alter it to use a 2D transform and you stand a chance.


> cv::Mat fftd( cv::Mat img, bool backwards) {
>
>  plan_f = fftwf_plan_dft_1d( width * height, data_in, fft, backwards ? FFTW_BACKWARD : FFTW_FORWARD,  FFTW_ESTIMATE );

> t1[1] = Microseconds();  for( i = 0, k = 0 ; i < height ; i++ ) {
> for( j = 0 ; j < width ; j++ ) { float *p = (float*)&img1_data[i *
> step + j * sizeof(float)]; data_in[k][0] = p[0]; data_in[k][1] =
> p[1]; //data_in[k][0] = ( float )img1_data[i * step + j];
> //data_in[k][1] = 0.0; k++; } }
>
> t2[0] = Microseconds(); fftwf_execute( plan_f ); t2[1] =

Also you will eventually want to move the planning and destroy plan 
outside the main FFT routine since you will presumably be transforming a 
lot of images of exactly the same size. Hope this helps.

-- 
Regards,
Martin Brown

[toc] | [prev] | [standalone]

Back to top | Article view | sci.image.processing

csiph-web