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)