• problem about replace opencv dft() function with fftw api

    From Zhi Sun@21:1/5 to All on Tue Jan 10 05:39:12 2017
    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

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From Martin Brown@21:1/5 to Zhi Sun on Thu Jan 12 08:22:34 2017
    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

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)