Выполнение свертки в частотной области вручную, но с неправильным выводом изображения в CPP / Opencv

Question

Я выполнил следующие шаги: - 1. Вычислил дфт изображения 2. Вычислил дфт ядра (но 1-й добавил его к размеру изображения) 3. Умножил действительные и мнимые части обоих дфтов по отдельности 4. Вычислил обратный дфт я пыталсядля отображения изображений на каждом промежуточном этапе, но окончательное изображение получается почти черным, за исключением углов. Выходное преобразование Фурье изображения после умножения и его обратный выходной ДФТ

входное изображение

enter code here

#include <iostream>
#include <stdlib.h>
#include <opencv2/opencv.hpp>
#include <stdio.h>
int r=100;
#define SIGMA_CLIP 6.0f
using namespace cv;
using namespace std;

void updateResult(Mat complex)
{
Mat work;
idft(complex, work);
Mat planes[] = {Mat::zeros(complex.size(), CV_32F), Mat::zeros(complex.size(), CV_32F)};
split(work, planes);                // planes[0] = Re(DFT(I)), planes[1] = Im(DFT(I))

magnitude(planes[0], planes[1], work);    // === sqrt(Re(DFT(I))^2 + Im(DFT(I))^2)
normalize(work, work, 0, 1, NORM_MINMAX);
imshow("result", work);
}
void shift(Mat magI) {

// crop if it has an odd number of rows or columns
magI = magI(Rect(0, 0, magI.cols & -2, magI.rows & -2));

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);
}
Mat updateMag(Mat complex )
{

Mat magI;
Mat planes[] = {Mat::zeros(complex.size(), CV_32F), Mat::zeros(complex.size(), CV_32F)};
split(complex, planes);                // planes[0] = Re(DFT(I)), planes[1] = Im(DFT(I))

magnitude(planes[0], planes[1], magI);    // sqrt(Re(DFT(I))^2 + Im(DFT(I))^2)

// switch to logarithmic scale: log(1 + magnitude)
magI += Scalar::all(1);
log(magI, magI);

shift(magI);
normalize(magI, magI, 1, 0, NORM_INF); // Transform the matrix with float values into a
         return magI;                                 // viewable image form (float between values 0 and 1).
//imshow("spectrum", magI);
 }
 Mat createGausFilterMask(Size imsize, int radius) {

// call openCV gaussian kernel generator
double sigma = (r/SIGMA_CLIP+0.5f);
Mat kernelX = getGaussianKernel(2*radius+1, sigma, CV_32F);
Mat kernelY = getGaussianKernel(2*radius+1, sigma, CV_32F);
// create 2d gaus
Mat kernel = kernelX * kernelY.t();

int w = imsize.width-kernel.cols;
int h = imsize.height-kernel.rows;

int r = w/2;
int l = imsize.width-kernel.cols -r;

int b = h/2;
int t = imsize.height-kernel.rows -b;

Mat ret;
copyMakeBorder(kernel,ret,t,b,l,r,BORDER_CONSTANT,Scalar::all(0));

return ret;

}

//code reference https://docs.opencv.org/2.4/doc/tutorials/core/discrete_fourier_transform/discrete_fourier_transform.html
int main( int argc, char** argv )
{ 

String file;
file = "lena.png";

Mat image = imread(file, CV_LOAD_IMAGE_GRAYSCALE);

Mat padded;                             

int m = getOptimalDFTSize( image.rows );
int n = getOptimalDFTSize( image.cols );  
copyMakeBorder(image, padded, 0, m - image.rows, 0, n -image.cols, BORDER_CONSTANT, Scalar::all(0));//expand input image to optimal size , on the border add zero values

Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
Mat complexI;
merge(planes, 2, complexI);  
dft(complexI, complexI); //computing dft
split(complexI, planes); //image converted to complex and real dft here

Mat mask = createGausFilterMask(padded.size(),r );  // Forming the gaussian filter
Mat mplane[] = {Mat_<float>(mask), Mat::zeros(mask.size(), CV_32F)};
Mat kernelcomplex;
merge(mplane, 2, kernelcomplex); 

dft(kernelcomplex, kernelcomplex);

split(kernelcomplex, mplane);// splitting the dft of kernel to real and complex 
mplane[1]=mplane[0]; //overwriting imaginary values with real values of kernel dft
Mat kernel_spec;
merge(mplane, 2, kernel_spec);
mulSpectrums(complexI, kernel_spec, complexI, DFT_ROWS);
Mat magI=updateMag(complexI);

namedWindow( "image fourier", CV_WINDOW_AUTOSIZE );

imshow("spectrum magnitude", magI);


updateResult(complexI); //converting to viewable form, computing idft

waitKey(0); 

return 0;
}

Какой шаг идет не так?Или мне не хватает какой-то концепции?

---

Отредактировал код с помощью Cris, и теперь он отлично работает.

Answer 1

Есть два сразу очевидных вопроса:

1. Гауссово является действительным и симметричным.Его преобразование Фурье тоже должно быть.Если ДПФ вашего ядра имеет ненулевой мнимый компонент, вы делаете что-то не так.

   Вероятно, вы ошибаетесь в том, что ваше ядро ​​имеет источник в середине изображения, а нев верхнем левом образце.Это та же проблема, что и в этот другой вопрос .Решение заключается в использовании эквивалента MATLAB ifftshift, реализация которого показана в документации OpenCV («шаг 6, Обрезка и перестановка») .

2. Чтобы применить свертку, вам нужно умножить два ДПФ вместе, а не на реальные части и мнимые части ДПФ.Умножение двух комплексных чисел a+ib и c+id приводит к ac-bd+iad+ibc, а не ac+ibd.

   Но поскольку ДПФ вашего ядра должны быть только действительными, вы можете просто умножить вещественную составляющуюядро с действительными и мнимыми компонентами изображения: (a+ib)c = ac+ibc.

---

Кажется очень окольным, что вы делаете с комплексными изображениями.Почему бы не позволить OpenCV обрабатывать все это для вас?Вероятно, вы можете * просто сделать что-то вроде этого:

Mat image = imread(file, CV_LOAD_IMAGE_GRAYSCALE);

// Expand input image to optimal size, on the border add zero values
Mat padded;                             
int m = getOptimalDFTSize(image.rows);
int n = getOptimalDFTSize(image.cols);  
copyMakeBorder(image, padded, 0, m - image.rows, 0, n -image.cols, BORDER_CONSTANT, Scalar::all(0));

// Computing DFT
Mat DFTimage;
dft(padded, DFTimage); 

// Forming the Gaussian filter
Mat kernel = createGausFilterMask(padded.size(), r);
shift(kernel);
Mat DFTkernel;
dft(kernel, DFTkernel);

// Convolution
mulSpectrums(DFTimage, DFTkernel, DFTimage, DFT_ROWS);

// Display Fourier-domain result
Mat magI = updateMag(DFTimage);
imshow("spectrum magnitude", magI);

// IDFT
Mat work;
idft(complex, work); // <- NOTE! Don't inverse transform log-transformed magnitude image!

Обратите внимание, что результат в области Фурье на самом деле является специальным представлением комплексно-сопряженного симметричного ДПФ, предназначенного для экономии места и вычислений.Чтобы вычислить полный комплексный вывод, добавьте DFT_COMPLEX_OUTPUT к вызову dft и DFT_REAL_OUTPUT к вызову idft (этот последний затем принимает симметрию и выдает реальный вывод, избавляя вас от хлопотвычисления величины).

^{* Я говорю, вероятно, потому что я не скомпилировал ничего из этого ... Если что-то не так, пожалуйста, дайте мне знать, или отредактируйте ответ и исправьте его.}