Getting a doubled, mirrored image on the output of IDFT-CodePudding

I would appreciate some help understanding the output of a small educational program I put together. I am new to OpenCV and don't have much C experience.

The goal of the script is to perform the following:

Load an image
Perform DFT on the image
Apply a circular binary mask to the spectrum, where the radius of the circle can be increased by hitting a key on the keyboard (essentially applying a very crude filter to the image)
Display the result of the inverse transform of the spectrum after the mask was applied

I have the basic functionality working: I can load the image, perform DFT, view the output image and increase the radius of the circle (advancing through a for-loop with the circle radius following i), and see the result of the inverse transform of the modified spectrum.

However I do not understand why the output is showing a vertically flipped copy of the input superimposed on the image (see example below with Lena.png). This is not my intended result. When I imshow() the inverse DFT result without applying the mask, I get a normal, non-doubled image. But after applying the mask, the IDFT output looks like this:

I am not looking for a quick solution: I would really appreciate if someone more experienced could ask leading questions to help me understand this result so that I can try to fix it myself.

My code:

#include <opencv2/core/core.hpp>
#include <opencv2/highgui.hpp>
#include <opencv2/imgcodecs.hpp>
#include <opencv2/imgproc.hpp>
#include <iostream>

using namespace cv;

void expand_img_to_optimal(Mat &padded, Mat &img);
void fourier_transform(Mat &img);

int main(int argc, char **argv)
{
    Mat input_img;
    input_img = imread("Lena.png" , IMREAD_GRAYSCALE);

    if (input_img.empty())
    {
        fprintf(stderr, "Could not Open image\n\n");
        return -1;
    }

    fourier_transform(input_img);
    return 0;
}

void fourier_transform(Mat &img)
{
    Mat padded;
    expand_img_to_optimal(padded, img);

    Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
    Mat complexI;
    merge(planes, 2, complexI);

    dft(complexI, complexI, DFT_COMPLEX_OUTPUT);
    
    // For-loop to increase mask circle radius incrementally
    for (int i=0; i<400; i =10) {
        Mat mask = Mat::ones(complexI.size(), complexI.type());
        mask.convertTo(mask, CV_8U);

        Mat dest = Mat::ones(complexI.size(), complexI.type());

        circle(mask, Point(mask.cols/2, mask.rows/2), i, 0, -1, 8, 0);

        complexI.copyTo(dest, mask);    

        Mat inverseTransform;
        idft(dest, inverseTransform, DFT_INVERSE|DFT_REAL_OUTPUT);
        normalize(inverseTransform, inverseTransform, 0, 1, NORM_MINMAX);
        imshow("Reconstructed", inverseTransform);
        waitKey(0);
    }
}

void expand_img_to_optimal(Mat &padded, Mat &img) {
    int row = getOptimalDFTSize(img.rows);
    int col = getOptimalDFTSize(img.cols);
    copyMakeBorder(img, padded, 0, row - img.rows, 0, col - img.cols, BORDER_CONSTANT, Scalar::all(0));
}

CodePudding user response：

This happens because you are inverse-transforming a frequency spectrum that is not conjugate-symmetric around the origin.

The origin of the frequency domain image is the top-left pixel. Your disk mask must be centered there. The frequency domain is periodic, so that the part of the mask that extends to the left of the image wraps around and comes in to the right edge, same with top and bottom edges.

The easiest way to generate a proper mask is to

create it with the origin at (mask.cols/2, mask.rows/2), like you already do, and then
apply the ifftshift operation.

OpenCV doesn’t have a ifftshift function, this Q&A has code that implements the fftshift, but is correct only for even-sized arrays. This is ifftshift, moving the pixel at mask.cols/2 to 0:

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

// input origin
int cx = sx / 2;
int cy = sy / 2;

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

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

It would be more efficient to do in-place, this makes two copies, but it is quick to implement.

CodePudding user response：

Firstly I'd like to thank @Cris Luengo for his helpful input on implementing the ifftshift in OpenCV.

In the end, the problem with my code was in this line:

Mat mask = Mat::ones(complexI.size(), complexI.type());

Instead of using the type of complexI, it looks like I should have used the type of img:

Mat mask = Mat::ones(complexI.size(), img.type());

Why? I'm not sure yet. Still trying to understand. Here is my complete code that is working how I intended:

#include <opencv2/core/core.hpp>
#include <opencv2/highgui.hpp>
#include <opencv2/imgcodecs.hpp>
#include <opencv2/imgproc.hpp>
#include <iostream>

using namespace cv;

void expand_img_to_optimal(Mat &padded, Mat &img);
void fourier_transform(Mat &img);
void ifft_shift(Mat &mask);                

int main(int argc, char **argv)
{
    Mat input_img;
    input_img = imread("Lena.png" , IMREAD_GRAYSCALE);

    if (input_img.empty())
    {
        fprintf(stderr, "Could not Open image\n\n");
        return -1;
    }

    fourier_transform(input_img);
    return 0;
}

void fourier_transform(Mat &img)
{
    Mat padded;
    expand_img_to_optimal(padded, img);

    Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
    Mat complexI;
    merge(planes, 2, complexI);

    dft(complexI, complexI, DFT_COMPLEX_OUTPUT);

    for (float i=0; i<4000; i =2) {
        // Create disk mask matrix
        Mat mask = Mat::ones(complexI.size(), CV_8U);
        circle(mask, Point(mask.cols/2, mask.rows/2), i, 0, -1, 8, 0);

        // Perform ifft shift
        ifft_shift(mask);

        // Destination matrix for masked spectrum
        Mat dest;
        complexI.copyTo(dest, mask);

        // Perform inverse DFT
        Mat inverseTransform;
        idft(dest, inverseTransform, DFT_INVERSE|DFT_REAL_OUTPUT);
        normalize(inverseTransform, inverseTransform, 0, 1, NORM_MINMAX);
        imshow("Reconstructed", inverseTransform);
        waitKey(0);
    }
}

void expand_img_to_optimal(Mat &padded, Mat &img) {
    int row = getOptimalDFTSize(img.rows);
    int col = getOptimalDFTSize(img.cols);
    copyMakeBorder(img, padded, 0, row - img.rows, 0, col - img.cols, BORDER_CONSTANT, Scalar::all(0));
}

void ifft_shift(Mat &mask) {
    // input sizes
    int sx = mask.cols;
    int sy = mask.rows;

    // input origin
    int cx = sx / 2;
    int cy = sy / 2;

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

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