function outDepth = zDFD(inImgRGB1, inImgRGB2, inKernel1, inKernel2, inSigma, inRange, inFlag);
% written by Changyin Zhou @ Columbia U, September 2009

%   This function is to do DFD from two defocused images
%   Input:
%    1. inImgRGB1, inImgRGB2: two RGB defocused image
%    2. inKernel1, inKernel2: two aperture patterns
%    3. inSigma: noise level
%    4. inRange: the possible depth range
%    5. inFlag: 0 = simply minimum energy; 1 = do 3 order polynomial fitting
%   Output:
%    1. Depth map

if nargin<7
    inFlag = 0;
end;

[hei, wid, ch] = size(inImgRGB1);

% Generate the prior matrix
global AStar;
AStar = eMakePrior(hei, wid)+0.00000000001;
lambda = inSigma^2*hei*wid./AStar;

inKernelSize = size(inKernel1, 1);
throughout1 = sum(inKernel1(:));
throughout2 = sum(inKernel2(:));
   
allConfidence = zeros(hei, wid, length(inRange));
count = 0;
for dd = inRange %Try different depths
    count = count+1;
    if abs(dd)<1 % |blur size<=1| ==> in focus
        if dd < 0
            dd = -1;
        else
            dd = 1;
        end;            
    end;
    %Scale the defocus kernel
    [kernel1, kernel2] = eScaleKernel(hei, wid, inKernel1, inKernel2, dd, throughout1, throughout2);
    for ch = 1:3 %Evaluate the three channels
        inImg1 = fft2(inImgRGB1(:, :, ch));
        inImg2 = fft2(inImgRGB2(:, :, ch));

        F0_hat = (inImg1.*conj(kernel1)+inImg2.*conj(kernel2))./(kernel1.*conj(kernel1)+kernel2.*conj(kernel2)+lambda+0.001); 
        error = abs(ifft2(F0_hat.*kernel1 - inImg1)) + abs(ifft2(F0_hat.*kernel2 - inImg2));
        allConfidence(:, :, count) = allConfidence(:, :, count)+ error;
    end;
end;
clear inImgRGB1 inImgRGB2 F0_hat error kernel1 kernel2 smoothing lambda

if inFlag==0 %Simply find the depth yielding minimum energy
    smoothing = fft2(fftshift(fspecial('gaussian', [hei, wid],  0.2*abs(dd)+2)));
    count = 0;
    for dd = inRange
        count =  count + 1;
        allConfidence(:, :, count) = abs(ifft2((fft2(allConfidence(:, :, count)).*smoothing)));
    end;

    [temp, outDepth] = min(allConfidence, [], 3);
    outDepth = inRange(outDepth);
    outDepth_GC = 0;
    return;
end;

if inFlag ==1 % polynomial depth interpolation
    smoothing = fft2(fftshift(fspecial('gaussian', [hei, wid],  0.1*abs(dd)+1.5)));
    count = 0;
    for dd = inRange
        count =  count + 1;
        allConfidence(:, :, count) = abs(ifft2((fft2(allConfidence(:, :, count)).*smoothing)));
    end;

    [temp, outDepthNaive] = min(allConfidence, [], 3);
    outDepthNaive(outDepthNaive==1|outDepthNaive==2) =3;
    outDepthNaive(outDepthNaive==count|outDepthNaive==count-1) =count-2;
    
    ind = [1:hei*wid]'; layer = hei*wid;
    five = zeros(5, layer);
    five(1, :) = allConfidence((outDepthNaive(:)-3)*layer+ind);
    five(2, :) = allConfidence((outDepthNaive(:)-2)*layer+ind);
    five(3, :) = allConfidence((outDepthNaive(:)-1)*layer+ind);
    five(4, :) = allConfidence((outDepthNaive(:)-0)*layer+ind);
    five(5, :) = allConfidence((outDepthNaive(:)+1)*layer+ind);
    clear allConfidence

    A = [(-2)^3, (-2)^2, (-2)^1, 1;...
        (-1)^3, (-1)^2, (-1)^1, 1;...
        (-0)^3, (-0)^2, (-0)^1, 1;...
        (1)^3, (1)^2, (1)^1, 1;...
        (2)^3, (2)^2, (2)^1, 1;...
        ];
    A = inv(A'*A)*A';
    five = A*five; %[a b c d]'

    temp = (4*five(2, :).^2-12*five(1, :).*five(3, :));
    select = temp<0;
    temp = sqrt(temp);
    offset1 = (-2*five(2, :)+temp)./(6*five(1, :));
    offset2 = (-2*five(2, :)-temp)./(6*five(1, :));
    offset1(select) = 0;
    offset2(select) = 0;
    select  = abs(offset1)>abs(offset2);
    offset1(select) = offset2(select);
    offset1(offset1>1) = 1;
    offset1(offset1<-1) = -1;
    offset1 = reshape(offset1, [hei, wid]);

    step = inRange(15)-inRange(14);
    outDepth = inRange(outDepthNaive) + offset1*step;
    outDepth_GC = 0;
    return;
end;
