Алгоритм деблринга в области Фурье

Я пытаюсь реализовать алгоритм деконволюции в преобразованной области Фурье. У меня есть рабочая реализация в Matlab, и я хочу перевести ее на C ++, используя библиотеку OpenCV.

По сути, я занимаюсь извлечением градиентов из входного изображения, делаю некоторые вещи в преобразованной области, а затем возвращаюсь в область пространства.

Проблематичной частью для меня является выполнение этого (поэлементного) разделения:

DFT(im) = ( conj( DFT(f) ) * DFT(image) + L2 * conj( DFT( gradKernel-x ) ) * DFT(mux) )+ ... ) / ( norm( DFT(f) )^2 + L2 * norm(gradKernel-x)^2 + ... )

f является гауссовым ядром, которое определено в коде.
DFT( gradKernel-x ) это БПФ ядра градиента в направлении х, т.е. DFT([1,-1]) дополненный нулями до размера размытого изображения. mux — это вспомогательная переменная для деконволюции.

Я решил выполнить деление по амплитуде и фазе отдельно в преобразованной области перед выполнением обратного ДПФ.

Я не знаю, где ошибка в моем коде, может быть, в делении, может быть в прямом / обратном преобразовании моих переменных, в гауссовом ядре …

Если кто-то может мне помочь, я был бы очень благодарен.

Вот критическая часть кода (обратите внимание, что я упростил его перед публикацией, так что не ожидайте хорошего результата размытия, если вы попробуете, в основном то, что я ожидаю от этого, это визуально приятное изображение на выходе):

imH00=imread("Cameraman256.png",0);
if(!imH00.data)
{
std::cout<<  "Could not open or find the image" << std::endl ;
return -1;
}

imH00.convertTo(imH00,CV_32F,1.0/255.0,0.0);

// Gaussian Kernel
Mat ker1D=getGaussianKernel(ksize,sigma,CV_32F);
fkernel.create(imH00.size(),CV_32F);
// zero-padding
fkernel.setTo(Scalar::all(0));
temp=ker1D*ker1D.t();
temp.copyTo(fkernel(Rect(0,0,temp.rows,temp.cols)));

// Fourier transform
Mat planes[] = {Mat_<float>(fkernel), Mat::zeros(fkernel.size(), CV_32F)};
Mat ffkernel;
merge(planes, 2, ffkernel);
dft(ffkernel, ffkernel,DFT_COMPLEX_OUTPUT);

// Gradient filter in frequency domain, trying to do something similar to psf2otf([1;-1],size(imH00)); in Matlab.
dx=Mat::zeros(imH00.size(),CV_32F);
dx.at<float>(0,0)=1;
dx.at<float>(0,1)=-1;
Mat dxplanes[] = {Mat_<float>(dx), Mat::zeros(dx.size(), CV_32F)};
Mat fdx;
merge(dxplanes, 2, fdx);
dft(fdx, fdx,DFT_COMPLEX_OUTPUT);

dy=Mat::zeros(imH00.size(),CV_32F);
dy.at<float>(0,0)=1;
dy.at<float>(1,0)=-1;
Mat dyplanes[] = {Mat_<float>(dy), Mat::zeros(dy.size(), CV_32F)};
Mat fdy;
merge(dyplanes, 2, fdy);
dft(fdy, fdy,DFT_COMPLEX_OUTPUT);

// Denominators

Mat den1;
Mat den2;
Mat den21;
Mat den22;

// ||fdx||^2 and ||fdy||^2
mulSpectrums(fdx,fdx,den21,DFT_COMPLEX_OUTPUT,true);
mulSpectrums(fdy,fdy,den22,DFT_COMPLEX_OUTPUT,true);
add(den21,den22,den2);

mulSpectrums(ffkernel,ffkernel,den1,0,true);

imHk=imH00.clone();

mux=Mat::zeros(imH00.size(),CV_32F);
muy=Mat::zeros(imH00.size(),CV_32F);

// FFT imH00
Mat fHktplanes[] = {Mat_<float>(imH00), Mat::zeros(imH00.size(), CV_32F)};
Mat fHkt;
merge(fHktplanes, 2, fHkt);
dft(fHkt, fHkt,DFT_COMPLEX_OUTPUT);

std::cout<<"starting deconvolution"<<std::endl;
for (int j=0; j<4; j++)
{
// Deconvolution

// Gradients

Mat ddx(1,2,CV_32F,Scalar(0));
ddx.at<float>(0,0)=1;
ddx.at<float>(0,1)=-1;
filter2D(imHk,dHx,CV_32F,ddx);

Mat ddy(2,1,CV_32F,Scalar(0));
ddy.at<float>(0,0)=1;
ddy.at<float>(1,0)=-1;
filter2D(imHk,dHy,CV_32F,ddy);mux=Scalar((float)-0.5*L1/L2);
add(mux,dHx,mux);

muy=Scalar((float)-0.5*L1/L2);
add(muy,dHy,muy);

// FFT mux, muy
Mat fmuxplanes[] = {Mat_<float>(mux), Mat::zeros(mux.size(), CV_32F)};
Mat fmux;
merge(fmuxplanes, 2, fmux);
dft(fmux, fmux,DFT_COMPLEX_OUTPUT);

Mat fmuyplanes[] = {Mat_<float>(muy), Mat::zeros(muy.size(), CV_32F)};
Mat fmuy;
merge(fmuyplanes, 2, fmuy);
dft(fmuy, fmuy,DFT_COMPLEX_OUTPUT);

Mat num1,num2,num3,num,den;

// Spectrums multiplication, complex conjugate
mulSpectrums(fHkt,ffkernel,num1,DFT_COMPLEX_OUTPUT,true);
mulSpectrums(fmux,fdx,num2,DFT_COMPLEX_OUTPUT,true);
mulSpectrums(fmuy,fdy,num3,DFT_COMPLEX_OUTPUT,true);

add(num2,num3,num2);
add(num1,L2*num2,num);
add(den1,L2*den2,den);

// Division in polar coordinates

Mat auxnum[2];
split(num,auxnum);
Mat auxden[2];
split(den,auxden);

Mat numMag,numPhase;
magnitude(auxnum[0],auxnum[1],numMag);
phase(auxnum[0],auxnum[1],numPhase);

Mat denMag,denPhase;
magnitude(auxden[0],auxden[1],denMag);
phase(auxden[0],auxden[1],denPhase);

Mat division[2];
divide(numMag,denMag,division[0]);
division[1]=numPhase-denPhase;

polarToCart(division[0],division[1],division[0],division[1]);
Mat fHk;
merge(division,2,fHk);

Mat imHkaux;
Mat planesfHk[2];
dft(fHk, imHkaux, DFT_INVERSE+DFT_SCALE);
split(imHkaux,planesfHk);
imHk=planesfHk[0]; // imHk is the Real part
}
imHk.convertTo(imHk,CV_8U,255);
imshow( "Deblurred image", imHk);

Спасибо