想从图像中检测模糊，但无法正确处理

Posted 2023-02-23

【中文标题】想从图像中检测模糊，但无法正确处理

【英文标题】：Want to detect blur from image, but couldn't get it right

【发布时间】：2020-06-27 15:11:15

【问题描述】：

我实际上想将此blur detection 转换为C++。作为OpenCV的初学者，我实际上是在关注this进行转换，但也许我弄错了。这是我的方法。我必须在 C++ 中使用 DFT 而不是 FFT。

(h, w) = image.shape
(cX, cY) = (int(w / 2.0), int(h / 2.0))
# compute the FFT to find the frequency transform, then shift
# the zero frequency component (i.e., DC component located at
# the top-left corner) to the center where it will be more
# easy to analyze
fft = np.fft.fft2(image)
fftShift = np.fft.fftshift(fft)

我把这个转换成

Mat I = imread( samples::findFile( filename ), IMREAD_GRAYSCALE);
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
dft(complexI, complexI, DFT_COMPLEX_OUTPUT);  // this way the result may fit in the source matrix

#For DFT shift as python code
// 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);

那么，下一部分

# zero-out the center of the FFT shift (i.e., remove low
# frequencies), apply the inverse shift such that the DC
# component once again becomes the top-left, and then apply
# the inverse FFT
fftShift[cY - size:cY + size, cX - size:cX + size] = 0
fftShift = np.fft.ifftshift(fftShift)
recon = np.fft.ifft2(fftShift)

我是这样转换的

// construct a Mat object to zero out of the center, here size = 60
Mat H;
Mat H(complexI.size(), CV_32F, Scalar(1));
float D = 0, D0 = 60;

for (int u = 0; u < H.rows; u++)

    for (int  v = 0; v < H.cols; v++)
    
        D = sqrt((u - scr.rows / 2)*(u - scr.rows / 2) + (v - scr.cols / 2)*(v - scr.cols / 2));
        if (D < D0)
        
            H.at<float>(u, v) = 0;
        
    

Mat planesH[] =  Mat_<float>(H.clone()), Mat_<float>(H.clone()) ;

Mat planes_dft[] =  complexI, Mat::zeros(complexI.size(), CV_32F) ;
split(complexI, planes_dft);

Mat planes_out[] =  Mat::zeros(complexI.size(), CV_32F), Mat::zeros(complexI.size(), CV_32F) ;
planes_out[0] = planesH[0].mul(planes_dft[0]);
planes_out[1] = planesH[1].mul(planes_dft[1]);

merge(planes_out, 2, complexIH);

#for Dft shift
Mat p0(complexIH, Rect(0, 0, cx, cy));   // Top-Left - Create a ROI per quadrant
Mat p1(complexIH, Rect(cx, 0, cx, cy));  // Top-Right
Mat p2(complexIH, Rect(0, cy, cx, cy));  // Bottom-Left
Mat p3(complexIH, Rect(cx, cy, cx, cy)); // Bottom-Right
p0.copyTo(tmp);
p3.copyTo(p0);
tmp.copyTo(p3);
p1.copyTo(tmp);                    // swap quadrant (Top-Right with Bottom-Left)
p2.copyTo(p1);
tmp.copyTo(p2);

Mat recon;
dft(complexIH, recon, DFT_INVERSE);

然后教程说明

# compute the magnitude spectrum of the reconstructed image,
# then compute the mean of the magnitude values
magnitude = 20 * np.log(np.abs(recon))
mean = np.mean(magnitude)
# the image will be considered "blurry" if the mean value of the
# magnitudes is less than the threshold value
return (mean, mean <= thresh)

我是这样转换的

Mat planes2[] = Mat_<float>(complexIH), Mat::zeros(complexIH.size(), CV_32F);
// compute the magnitude and switch to logarithmic scale
// => log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
split(recon, planes2);                   // planes2[0] = Re(DFT(I), planes2[1] = Im(DFT(I))
magnitude(planes2[0], planes2[1], planes2[0]);// planes2[0] = magnitude
Mat output = planes2[0];

output += Scalar::all(1);                    // switch to logarithmic scale
log(output, output);
float avg = mean(magI)[0];

我知道这是一团糟。我想像教程说的那样得到模糊值。

【问题讨论】：

你试过这个blog 吗？这看起来很简单，您也将获得该方法的 C++ 实现。

那个使用拉普拉斯算子，但我想使用上面提到的FFT。你能指导我吗？

“我必须在 C++ 中使用 DFT 而不是 FFT。” cv::DFT 函数实现 FFT。 FFT 是一系列高效计算 DFT 的算法。

对不起，我只是 OpenCV 的初学者。不知道。

【参考方案1】：

我认为这接近于原始 Python 代码

#include <iostream>
#include <opencv2/opencv.hpp>

using namespace cv;
using namespace std;

int main(int argc, char **argv) 
    if (argc <= 1) 
        fprintf(stderr, "Error: missing image file\n");
        return 1;
    
    string image_file = argv[1];

    cout << "Processing " << image_file << std::endl;
    Mat frame = imread(image_file, IMREAD_GRAYSCALE);

    // Go float
    Mat fImage;
    frame.convertTo(fImage, CV_32F);

    // FFT
    cout << "Direct transform...\n";
    Mat fourierTransform;
    dft(fImage, fourierTransform, DFT_SCALE|DFT_COMPLEX_OUTPUT);

    int Wd = frame.cols;
    int Ht = frame.rows;
    int cx = Wd/2;
    int cy = Ht/2;
    int Sw = 60;
    int Sh = 60;

    //center low frequencies in the middle
    //by shuffling the quadrants.
    Mat q0(fourierTransform, Rect(0, 0, cx, cy));       // Top-Left - Create a ROI per quadrant
    Mat q1(fourierTransform, Rect(cx, 0, cx, cy));      // Top-Right
    Mat q2(fourierTransform, Rect(0, cy, cx, cy));      // Bottom-Left
    Mat q3(fourierTransform, 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);

    // Block the low frequencies
    fourierTransform(Rect(cx-Sw,cy-Sh,2*Sw,2*Sh)).setTo(0);

    //shuffle the quadrants to their original position
    Mat orgFFT;
    fourierTransform.copyTo(orgFFT);
    Mat p0(orgFFT, Rect(0, 0, cx, cy));       // Top-Left - Create a ROI per quadrant
    Mat p1(orgFFT, Rect(cx, 0, cx, cy));      // Top-Right
    Mat p2(orgFFT, Rect(0, cy, cx, cy));      // Bottom-Left
    Mat p3(orgFFT, Rect(cx, cy, cx, cy));     // Bottom-Right

    p0.copyTo(tmp);
    p3.copyTo(p0);
    tmp.copyTo(p3);

    p1.copyTo(tmp);                                     // swap quadrant (Top-Right with Bottom-Left)
    p2.copyTo(p1);
    tmp.copyTo(p2);

    // IFFT
    cout << "Inverse transform...\n";
    Mat invFFT;
    Mat logFFT;
    double minVal,maxVal;

    dft(orgFFT, invFFT, DFT_INVERSE|DFT_REAL_OUTPUT);

    //img_fft = 20*numpy.log(numpy.abs(img_fft))
    invFFT = cv::abs(invFFT);
    cv::minMaxLoc(invFFT,&minVal,&maxVal,NULL,NULL);

    //check for impossible values
    if(maxVal<=0.0)
        cerr << "No information, complete black image!\n";
        return 1;
    

    cv::log(invFFT,logFFT);
    logFFT *= 20;

    //result = numpy.mean(img_fft)
    cv::Scalar result= cv::mean(logFFT);
    cout << "Result : "<< result.val[0] << endl;

    // Back to 8-bits
    Mat finalImage;
    logFFT.convertTo(finalImage, CV_8U);

    // show if you like
    imshow("Input", frame);
    imshow("Result", finalImage);
    cv::waitKey();

    return 0;

【讨论】：

这个答案只有一个错误。 cv::log 不在位，因此需要将其结果分配给 invFFT 矩阵。

修复了cv::abs 错误。它的返回值没有被使用。它可能是负面的，随后让cv::log 崩溃。