DEMO 2: Shape from Defocus via I-Divergence

Shape from Defocus and Image Restoration

Demonstration of Shape from Defocus via I-Divegence Minimization method

clear all
close all

Data set initialization

load input images and initialize parameters

load('DataSet.mat');
% first resize...
scalingFactor = 1/2;
I1 = imresize(I1,scalingFactor); % limit number of computations
I2 = imresize(I2,scalingFactor); % limit number of computations
% ...then select region of interest
I(:,:,2) = I1(1:end-10,60:end-50);
I(:,:,1) = I2(1:end-10,60:end-50);
I = I+(I<1e-3);% make sure images are strictly positive
[M,N,K] = size(I);
% set the maximum number of iterations
MaxIterations = 100;
% lens parameters
F = 12e-3; % focal length
Fnumb = 2; % F number
D = F/Fnumb; % aperture (effective lens diameter)
gamma = 3e4; % calibration parameter
p = 1./(1/F-1./[.52 .85]); % distance CCD-lens
% select the domain where the data term is considered
win = 3;
mask = ones(M,N);
mask(1:win,:) = 0;
mask(end-win+1:end,:) = 0;
mask(:,1:win) = 0;
mask(:,end-win+1:end) = 0;
dt_data = 2e-1; % data term step
dt_reg = 1e0; % regularization term step
dx = 1e-3; % discrete gradient step
maxDelta = 1; % max amount of relative blur
% regularization constant
% to avoid division by zero in the preconditioner
epsilon = 40*max(I(:))/255;
K = length(p); % number of focus settings
% preconditioning threshold
thresh = 6e-2;

Initialize unknown shape and radiance

initialize depth map with a plane and the radiance with one input image

DepthMap_e = sum(p)*F/(sum(p)-2*F)*ones(M,N);
Radiance_e = I(:,:,1);

Start deblurring and shape estimation

integrate the gradient-flow for MaxIterations iterations

for i=1:MaxIterations
    % Compute approximate gradients
    for k=1:K
        % compute blur maps
        BlurMap_e(:,:,k) = gamma*D/2*abs(1-p(k).*...
            (1/F-1./DepthMap_e));
        BlurMap_dx(:,:,k) = gamma*D/2*abs(1-p(k).*...
            (1/F-1./(DepthMap_e+dx)));
        % limit the amount of blur (to prevent out
        % of memory in defocusApprox)
        BlurMap_e(:,:,k) = BlurMap_e(:,:,k).*...
            (BlurMap_e(:,:,k)<2)+2*(BlurMap_e(:,:,k)>=2);
        BlurMap_dx(:,:,k) = BlurMap_dx(:,:,k).*...
            (BlurMap_dx(:,:,k)<2)+2*(BlurMap_dx(:,:,k)>=2);
    end
    % compute synthetic defocused images and gradients
    [I_e,F_e,I_dx] = defocusApprox(Radiance_e,BlurMap_e,I,BlurMap_dx);
    % update radiance (Lucy-Richardson iteration)
    Radiance_e = Radiance_e.*F_e;
    % compute gradient of image wrt depth map
    dIds = (I_dx-I_e)/dx;
    % compute gradient of i-divergence wrt depth map
    dEds = sum((1-I./I_e).*dIds,3);
    % compute preconditioning (positive definite kernel)
    precond = (epsilon+sum(abs(dIds),3)).^2./(epsilon+sum(I_e,3));
    dEds = dEds./precond;
    % precondition gradient
    dEds = (abs(dEds) <  thresh).*dEds+...
        (abs(dEds) >= thresh).*thresh.*sign(dEds);
    % at the boundary use only regularization
    dDatads = dEds.*mask;
    % combine data term and regularization term
    % make sure the iterations are an even number
    nIt = 2*(floor(dt_reg/1e-1)+1);
    beta_data = dt_data/nIt;
    beta_regularization = dt_reg/nIt;
    for tt=1:nIt
        % compute gradient of the cost functional
        % (data term + regularization) wrt the depth map
        if tt-floor(tt/2)*2==0
            dDepth = beta_data*dDatads - beta_regularization*...
                laplacian(DepthMap_e,'forward');
        else
            dDepth = beta_data*dDatads - beta_regularization*...
                laplacian(DepthMap_e,'backward');
        end
        % update depth map
        DepthMap_e = DepthMap_e - dDepth;
    end
end

Results visualization

visualize input images, the estimated radiance and the estimated depth map (as a mesh and in gray-scale)

figure(1)
subplot(211)
image([I(:,:,1) 200*ones(size(I,1),10) ...
    I(:,:,2) 200*ones(size(I,1),10) ...
    Radiance_e])
axis equal
axis off
set(gca,'XTick',[],'YTick',[]);
set(gca,'Box','off');
colormap gray(256)
title('Input #1/Input #2/Estimated Radiance')
drawnow;
step = ceil(min([size(DepthMap_e,1) ...
    size(DepthMap_e,2)])/50);
pix2mm = 1/600;
[X,Y] = meshgrid(pix2mm*[1:step:size(DepthMap_e,2)],...
    pix2mm*[1:step:size(DepthMap_e,1)]);
axisZ = 1./(1/F-1./p);
subplot(223)
h = surf(X,-DepthMap_e(1:step:end,1:step:end),-Y,DepthMap_e(1:step:end,1:step:end));
shading interp
hold off
axis([0 N -axisZ(2) -axisZ(1) -M 0])
set(gca,'XTick',[],'YTick',[],'ZTick',[]);
colormap gray(256)
title('Estimated Depth Map')
view([-1 -2 1]);
axis equal
axis off
set(h,'EdgeColor','r');
drawnow;
subplot(224)
imagesc(DepthMap_e);
colormap gray(256)
axis equal
set(gca,'XTick',[],'YTick',[]);
set(gca,'Box','off');
title('Estimated Depth Map')
drawnow;