푸리에 도메인의 디 블러 링 알고리즘

Question

푸리에 변환 도메인에서 디콘 볼 루션 알고리즘을 구현하려고합니다. MATLAB에서 직접 구현을 해봤고 OpenCV 라이브러리를 사용하여 C++로 변환하려고합니다.

기본적으로 내가하는 일은 입력 이미지에서 그라디언트를 추출하는 것입니다. 변환 된 도메인에서 일부 작업을 수행 한 다음 다시 공간 영역으로 돌아옵니다.

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)은 x 방향의 그래디언트 커널의 FFT입니다 (예 : DFT([1,-1])). 흐린 이미지의 크기에 제로 패딩됩니다. mux는 디컨 볼 루션을 수행하기위한 보조 변수입니다.

역 DFT를 수행하기 전에 변환 된 도메인에서 크기 및 위상별로 나누기를 결정했습니다.

가우스 커널에서, 내 변수의 변환 역/어쩌면 전방에, 어쩌면 부문에서, 내 코드에서 오류입니다 몰라요 ...

누군가가 나를 도울 수 있다면, 나는 고맙게 생각한다.

여기에 코드의 중요한 부분이 있습니다 (게시하기 전에 간단히 표시 했으므로 시도 할 때 좋은 흐려짐 결과가 기대되지 않습니다. 기본적으로 기대되는 것은 시각적으로 즐거운 출력 이미지입니다) :

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); 

는

Answer 1

문제는 필터의 푸리에 변환에 있었다 감사드립니다. 변환하기 전에 필터 커널을 이동해야합니다. 이것은 psl2otf 함수가 Matlab에서하는 것과 같습니다.

float sigma=1.0; 
short int ksize=13; // always odd 

Mat ker1D=getGaussianKernel(ksize,sigma,CV_32F); //1D gaussian kernel 
fkernel.create(myImage.size(),CV_32F); // 

// zero-padding 
fkernel.setTo(Scalar::all(0)); 
temp=ker1D*ker1D.t(); // 2D gaussian kernel, (Gaussian filter is separable) 

int r=(ksize-1)/2; //radius 

// shifting 

temp(Rect(r,r,r+1,r+1)).copyTo(fkernel(Rect(0,0,r+1,r+1))); // q1 
temp(Rect(r,0,r+1,r)).copyTo(fkernel(Rect(0,fkernel.cols-r,r+1,r))); // q2 
temp(Rect(0,r,r,r+1)).copyTo(fkernel(Rect(fkernel.rows-r,0,r,r+1))); // q3 
temp(Rect(0,0,r,r)).copyTo(fkernel(Rect(fkernel.rows-r,fkernel.cols-r,r,r))); // q4 
// DFT 
Mat planes[] = {Mat_<float>(fkernel), Mat::zeros(fkernel.size(), CV_32F)}; 
Mat ffkernel; 
merge(planes, 2, ffkernel); 
dft(ffkernel, ffkernel,DFT_COMPLEX_OUTPUT);