来自 DFT/IDFT 计算的负值

Question

我正在创建一个截止频率为 0.015 的频域低通滤波器和相应的高通滤波器，如下所示：

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; 

，这些应用于像下面的sn-p这样的图像。

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

输入图像为 256x256 且采用对数刻度，低通滤波器应用于高通滤波输入图像的平方。

然后我想取最终 idft 的结果的平方根， 但它包含负值；因此，出现了几个 NaN 值。

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

为什么会有负值，我该如何处理它们才能取平方根？

---

EDIT：添加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

直到第三次或第四次阅读您的代码时，我才发现这个...

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

您在这里创建一个复值滤波器，它不是共轭对称的。应用此滤波器会创建一个非共轭对称的频谱，因此您会得到一个非实数逆变换。

你想做的是有一个纯粹的过滤器：

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

从幅度和相位（而不是实部和虚部）的角度考虑频域滤波器。幅度改变了您如何改变图像中每个频率分量的幅度，相移每个频率分量。一般来说（除了非常特殊的情况）你不想改变图像中的频率分量，因为它会弄乱你的图像。因此，将频域滤波器的相位保持在 0。零相位意味着每个值都是非负的和实数。

也就是说，您移动图像频谱、应用过滤器和向后移动的顺序：

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

等同于简单地移动过滤器：

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

还请注意，我将cv::DFT_ROWS 替换为 0。根据the documentation，该标志表示每个图像行都是独立的一维变换。我的猜测是，这只对 CCS 压缩矩阵很重要，但最好还是不要使用它。