Negative values from DFT/IDFT calculation

Question

I'm creating a low-pass filter in frequency domain with 0.015 cutoff frequency and the corresponding high-pass filter like below:

auto getGaussianFilter(int32_t size) {
    const auto sigmaf = 0.015f*size; // cutoff freq.

    cv::Mat kernel = cv::Mat::zeros(size, size, CV_32FC1);

    for (auto fy = -size/2; fy < size/2; ++fy)
        for (auto fx = -size/2; fx < size/2; ++fx) 
            kernel.at<float>(fy + size/2, fx + size/2) = std::exp(-(fy*fy+fx*fx)/(2*sigmaf*sigmaf));
    return kernel;
} 

...

lowpassFilter  = getGaussianFilter(256);
highpassFilter = 1 - lowpassFilter; 

, and these are applied to images like following snippet.

std::vector<cv::Mat> planes { src, cv::Mat::zeros(src.size(), src.type()) }; 
std::vector<cv::Mat> fplanes { filter, filter };
cv::Mat complex, complexFilter;
cv::merge(planes, complex);
cv::merge(fplanes, complexFilter);

cv::dft(complex, complex, cv::DFT_COMPLEX_OUTPUT);
shift(complex); // swapping quadrants
cv::mulSpectrums(complex, complexFilter, complex, cv::DFT_ROWS);
shift(complex); 
cv::idft(complex, complex, cv::DFT_SCALE | cv::DFT_REAL_OUTPUT); 

Input images are 256x256 and in logarithmic scale, and the low-pass filter is applied to the square of the high-pass filtered input images.

And then I want to take a square root of the outcome of the final idft, but it contains negative values; thus, several NaN values appear.

applyFilter(src, dst, highpassFilter);
cv::pow(dst, 2, dst);
applyFilter(dst, dst, lowpassFilter);
cv::sqrt(dst, dst);  // NaN !

Why are there negative values, and how can I deal with them in order to take the square root ?

---

EDIT: added the code of shift

void shift(cv::Mat& src) {
    src = src(cv::Rect(0, 0, src.cols & -2, src.rows & -2)); 
    const auto cy = src.rows/2, cx = src.cols/2;
    cv::Mat q0(src, cv::Rect(0, 0, cx, cy));
    cv::Mat q1(src, cv::Rect(cx, 0, cx, cy));
    cv::Mat q2(src, cv::Rect(0, cy, cx, cy));
    cv::Mat q3(src, cv::Rect(cx, cy, cx, cy));
    cv::Mat tmp;
    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);

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

Answer 1

I didn't catch this one until the third or fourth read through your code...

std::vector<cv::Mat> fplanes { filter, filter };
cv::merge(fplanes, complexFilter);

You create a complex-valued filter here, which is not conjugate symmetric. Applying this filter creates a frequency spectrum which is not conjugate symmetric, and hence you get a non-real inverse transform.

What you want to do is have a purely real filter:

std::vector<cv::Mat> fplanes { filter, cv::Mat::zeros(filter(), filter()) };
cv::merge(fplanes, complexFilter);

Think of frequency-domain filters in terms of a magnitude and a phase (rather than a real and imaginary component). The magnitude changes how you change the amplitude of each frequency component in your image, the phase shifts each frequency component. In general (except in very specific cases) you do not want to shift frequency components in your image, as it will mess up your image. Thus, keep the phase of your frequency-domain filters at 0. A zero phase means that each value is non-negative and real.

That said, your sequence of shifting the frequency spectrum of the image, applying the filter, and shifting back:

cv::dft(complex, complex, cv::DFT_COMPLEX_OUTPUT);
shift(complex); // swapping quadrants
cv::mulSpectrums(complex, complexFilter, complex, cv::DFT_ROWS);
shift(complex); 
cv::idft(complex, complex, cv::DFT_SCALE | cv::DFT_REAL_OUTPUT); 

is identical to simply shifting the filter:

cv::dft(complex, complex, cv::DFT_COMPLEX_OUTPUT);
shift(complexFilter); // swapping quadrants
cv::mulSpectrums(complex, complexFilter, complex, 0);
cv::idft(complex, complex, cv::DFT_SCALE | cv::DFT_REAL_OUTPUT); 

Note also that I replaced cv::DFT_ROWS with 0. According to the documentation, that flag indicates that each image row is an independent 1D transform. My guess is that this only matters for CCS-packed matrices, but it's better to leave it out anyway.