I am actually want to convert this blur detection into C++. As a beginner in OpenCV, I am actually following this for conversion, But maybe I am getting it wrong. Here is my approach. I have to use DFT instead of FFT in C++.
(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)
I converted this by
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);
Then, in the next part
# 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)
I converted this in this way
// 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);
Then the tutorial stated
# 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)
And I converted this in this way
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];
I know it is a mess. I want to get the blur value like the tutorial says.
I think this comes close to the original Python code
#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;
}