I'm trying to implement a deconvolution algorithm in Fourier transformed domain. I've a working implementation by myself in Matlab, and I want to translate it to C++ using OpenCV library.
Basically what I do is to extract the gradients from an input image, I do some stuff in transformed domain and then I come back to space domain.
The problematic part to me is to perform this (element wise) division:
DFT(im) = ( conj( DFT(f) ) * DFT(image) + L2 * conj( DFT( gradKernel-x ) ) * DFT(mux) )+ ... ) / ( norm( DFT(f) )^2 + L2 * norm(gradKernel-x)^2 + ... )
f is a gaussian kernel that is defined in the code.
DFT( gradKernel-x ) is the FFT of the gradient kernel in x direction, i.e. DFT([1,-1]) zero-padded to the size of the blurred image. mux is an auxiliary variable to perform the deconvolution.
I decided to perform the division by magnitude and phase separately in transformed domain before performing an inverse DFT.
I don't know where is the error in my code, maybe in the division, maybe in the forward/inverse transform of my variables, in the gaussian kernel...
If somebody can help me, I'd very grateful.
Here is the critical part of the code (note that I simplified it before posting, so don't expect a good deblurring result if you try it, basically what I expect from this is a visually pleasant output image):
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);
Thank you
