fftshift c++ implemetation for openCV

Question

I have already looked in this question fftshift/ifftshift C/C++ source code

I'm trying to implement fftshift from matlab

this is the code from the matlab function for 1D array

numDims = ndims(x);
    idx = cell(1, numDims);
    for k = 1:numDims
        m = size(x, k);
        p = ceil(m/2);
        idx{k} = [p+1:m 1:p];
    end
y = x(idx{:});

my c++/openCV code is, what fftshift basically does is swap the values from a certain pivot place.
since I can't seem to understand how is the matrix built in opencv for complex numbers.
it says here
http://docs.opencv.org/modules/core/doc/operations_on_arrays.html#dft
CCS (complex-conjugate-symmetrical
I thought it will be easier to split the complex numbers into real and imaginary and swap them. and then merge back to one matrix.

cv::vector<float> distanceF (f.size());

//ff = fftshift(ff);
cv::Mat ff;
cv::dft(distanceF, ff, cv::DFT_COMPLEX_OUTPUT);

//Make place for both the complex and the real values
cv::Mat planes[] = {cv::Mat::zeros(distanceF.size(),1, CV_32F), cv::Mat::zeros(distanceF.size(),1, CV_32F)};
cv::split(ff, planes);                   // planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))

int numDims = ff.dims;
for (int i = 0; i < numDims; i++)
{
    int m = ff.rows;
    int p = ceil(m/2);

}

my problem is that because of my input to the DFT is a vector<float> I can't seem to be able to create planes mat in order to split the complex numbers?

Can you think how a better way to make the swap of the values inside the cv::mat data struct?

Answer 1

There are no implementations in earlier answers that work correctly for odd-sized images.

fftshift moves the origin from the top-left to the center (at size/2). ifftshift moves the origin from the center to the top-left. These two actions are identical for even sizes, but differ for odd-sizes.

For an odd size, fftshift swaps the first (size+1)/2 pixels with the remaining size/2 pixels, which moves the pixel at index 0 to size/2. ifftshift does the reverse, swapping the first size/2 pixels with the remaining (size+1)/2 pixels. This code is the most simple implementation of both these actions that I can come up with. (Note that (size+1)/2 == size/2 if size is even.)

bool forward = true; // true for fftshift, false for ifftshift
cv::Mat img = ...;   // the image to process

// input sizes
int sx = img.cols;
int sy = img.rows;

// size of top-left quadrant
int cx = forward ? (sx + 1) / 2 : sx / 2;
int cy = forward ? (sy + 1) / 2 : sy / 2;

// split the quadrants
cv::Mat top_left(img, cv::Rect(0, 0, cx, cy));
cv::Mat top_right(img, cv::Rect(cx, 0, sx - cx, cy));
cv::Mat bottom_left(img, cv::Rect(0, cy, cx, sy - cy));
cv::Mat bottom_right(img, cv::Rect(cx, cy, sx - cx, sy - cy));

// merge the quadrants in right order
cv::Mat tmp1, tmp2;
cv::hconcat(bottom_right, bottom_left, tmp1);
cv::hconcat(top_right, top_left, tmp2);
cv::vconcat(tmp1, tmp2, img);

This code makes a copy of the full image twice, but it is easy and quick to implement. A more performant implementation would swap values in-place. This answer has correct code to do so on a single line, it would have to be applied to each column and each row of the image.

Answer 2

Here is what I do (quick and dirty, can be optimized):

// taken from the opencv DFT example (see opencv/samples/cpp/dft.cpp within opencv v440 sourcecode package)
cv::Mat fftshift(const cv::Mat& mat){
    
    // create copy to not mess up the original matrix (ret is only a "window" over the provided matrix)
    cv::Mat cpy;
    mat.copyTo(cpy);

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

    // rearrange the quadrants of Fourier image so that the origin is at the image center
    int cx = ret.cols/2;
    int cy = ret.rows/2;
    cv::Mat q0(ret, cv::Rect(0, 0, cx, cy));   // Top-Left - Create a ROI per quadrant
    cv::Mat q1(ret, cv::Rect(cx, 0, cx, cy));  // Top-Right
    cv::Mat q2(ret, cv::Rect(0, cy, cx, cy));  // Bottom-Left
    cv::Mat q3(ret, cv::Rect(cx, cy, cx, cy)); // Bottom-Right

    cv::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);

    return ret;
}

// reverse the swapping of fftshift. (-> reverse the quadrant swapping)
cv::Mat ifftshift(const cv::Mat& mat){

    // create copy to not mess up the original matrix (ret is only a "window" over the provided matrix)
    cv::Mat cpy;
    mat.copyTo(cpy);

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

    // rearrange the quadrants of Fourier image so that the origin is at the image center
    int cx = ret.cols/2;
    int cy = ret.rows/2;
    cv::Mat q0(ret, cv::Rect(0, 0, cx, cy));   // Top-Left - Create a ROI per quadrant
    cv::Mat q1(ret, cv::Rect(cx, 0, cx, cy));  // Top-Right
    cv::Mat q2(ret, cv::Rect(0, cy, cx, cy));  // Bottom-Left
    cv::Mat q3(ret, cv::Rect(cx, cy, cx, cy)); // Bottom-Right

    cv::Mat tmp; // swap quadrants (Bottom-Right with Top-Left)
    q3.copyTo(tmp);
    q0.copyTo(q3);
    tmp.copyTo(q0);
    q2.copyTo(tmp); // swap quadrant (Bottom-Left with Top-Right)
    q1.copyTo(q2);
    tmp.copyTo(q1);

    return ret;
}

Answer 3

I know, this is quite an old thread, but I found it today while looking for a solution to shift the fft-result. and maybe the little function I wrote with the help of this site and other sources, could be helpful for future readers searching the net and ending up here too.

bool FftShift(const Mat& src, Mat& dst)
{
  if(src.empty()) return true;

  const uint h=src.rows, w=src.cols;        // height and width of src-image
  const uint qh=h>>1, qw=w>>1;              // height and width of the quadrants

  Mat qTL(src, Rect(   0,    0, qw, qh));   // define the quadrants in respect to
  Mat qTR(src, Rect(w-qw,    0, qw, qh));   // the outer dimensions of the matrix.
  Mat qBL(src, Rect(   0, h-qh, qw, qh));   // thus, with odd sizes, the center
  Mat qBR(src, Rect(w-qw, h-qh, qw, qh));   // line(s) get(s) omitted.

  Mat tmp;
  hconcat(qBR, qBL, dst);                   // build destination matrix with switched
  hconcat(qTR, qTL, tmp);                   // quadrants 0 & 2 and 1 & 3 from source
  vconcat(dst, tmp, dst);

  return false;
}

Answer 4

I have been implementing it myself based on this post, I used Fabian implementation which is working fine. But there is a problem when there is an odd number of row or column, the shift is then not correct.

You need then to padd your matrix and after to get rid of the extra row or column.

   {bool flag_row = false;
    bool flag_col = false;

    if( (inputMatrix.rows % 2)>0)
    {
        cv::Mat row = cv::Mat::zeros(1,inputMatrix.cols, CV_64F);  
        inputMatrix.push_back(row);
        flag_row =true;
    }

    if( (inputMatrix.cols % 2)>0)
    {
        cv::Mat col = cv::Mat::zeros(1,inputMatrix.rows, CV_64F);  
        cv::Mat tmp;
        inputMatrix.copyTo(tmp);
        tmp=tmp.t();
        tmp.push_back(col);
        tmp=tmp.t();
        tmp.copyTo(inputMatrix);

        flag_col = true;
    }

    int cx = inputMatrix.cols/2.0;
    int cy = inputMatrix.rows/2.0;

    cv::Mat outputMatrix;
    inputMatrix.copyTo(outputMatrix);

    // rearrange the quadrants of Fourier image
    // so that the origin is at the image center
    cv::Mat tmp;
    cv::Mat q0(outputMatrix, cv::Rect(0, 0, cx, cy));
    cv::Mat q1(outputMatrix, cv::Rect(cx, 0, cx, cy));
    cv::Mat q2(outputMatrix, cv::Rect(0, cy, cx, cy));
    cv::Mat q3(outputMatrix, cv::Rect(cx, cy, cx, cy));

    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);

    q1.copyTo(tmp);
    q2.copyTo(q1);
    tmp.copyTo(q2);

    int row = inputMatrix.rows;
    int col = inputMatrix.cols;
    if(flag_row)
    {
        outputMatrix = Tools::removerow(outputMatrix,row/2-1);
    }
    if(flag_col)
    {
        outputMatrix = Tools::removecol(outputMatrix,col/2-1);
    }

    return outputMatrix;

Answer 5

Ok, this thread is may be out of date in the meantime but maybe for other users.. Take a look at the samples:

opencv/samples/cpp/dft.cpp (line 66 - 80)

int cx = mag.cols/2;
int cy = mag.rows/2;

// rearrange the quadrants of Fourier image
// so that the origin is at the image center
Mat tmp;
Mat q0(mag, Rect(0, 0, cx, cy));
Mat q1(mag, Rect(cx, 0, cx, cy));
Mat q2(mag, Rect(0, cy, cx, cy));
Mat q3(mag, Rect(cx, cy, cx, cy));

q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);

q1.copyTo(tmp);
q2.copyTo(q1);
tmp.copyTo(q2);

I think that's a short and clean way for different dimensions.

Answer 6

How about using adjustROI and copyTo instead of .at()? It would certainly be more efficient:

Something in the lines of (for your 1D case):

Mat shifted(ff.size(),ff.type());
pivot = ff.cols / 2;
ff(Range::all(),Range(pivot + 1, ff.cols)).copyTo(shifted(Range::all(),Range(0,pivot)));
ff(Range::all(),Range(0,pivot+1)).copyTo(shifted(Range::all(),Range(pivot,ff.cols)));

For the 2D case, two more lines should be added, and the rows ranges modified...

Answer 7

In Matlab's implementation, the main code are the two lines:

idx{k} = [p+1:m 1:p];
y = x(idx{:});

The first one obtains the correct index order against the original one; then the second one assigns the output array according to the index order. Therefore, if you want to re-write Matlab's implementation without data swapping, you need to allocate a new array and assign the array.

Answer 8

this is for future reference: been tested and is bit accurate for 1D

    cv::Mat ff;
    cv::dft(distanceF, ff, cv::DFT_ROWS|cv::DFT_COMPLEX_OUTPUT);

    //Make place for both the complex and the real values
    cv::Mat planes[] = {cv::Mat::zeros(distanceF.size(),1, CV_32F), cv::Mat::zeros(distanceF.size(),1, CV_32F)};
    cv::split(ff, planes);                   // planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))

    int m = planes[0].cols;
    int pivot = ceil(m/2);
    //duplicate FFT results with Complex conjugate in order to get exact matlab results

    for (int i = pivot + 1, k = pivot; i < planes[1].cols; i++, k--)
    {
         planes[1].at<float>(i) = planes[1].at<float>(k) * -1; 
         planes[0].at<float>(i) = planes[0].at<float>(k);
    }   

    //TODO maybe we need to see what happens for pair and odd ??
    float im  = planes[1].at<float>(0);
    float re = planes[0].at<float>(0);

    for (int i = 0; i < pivot; i++)
    {   
        //IM
        planes[1].at<float>(i) = planes[1].at<float>(pivot + i +1); 
        planes[1].at<float>(pivot +i +1) = planes[1].at<float>(i +1);

        //Real
        planes[0].at<float>(i) = planes[0].at<float>(pivot + i +1); 
        planes[0].at<float>(pivot +i +1) = planes[0].at<float>(i +1);
    }
    planes[1].at<float>(pivot) = im;
    planes[0].at<float>(pivot) = re;