idivergence.m

3D shape reconstruction matlab code. It used shape from defocus technique with divergence. You can r

% Depth from Defocus and Image Restoration 
% by Minimizing I-Divergence
%
% Copyright 2006 Paolo Favaro ([e-mail])
% 
% School of Engineering and Physical Sciences
% Heriot-Watt University, Edinburgh, UK
% 
% Last revision: May 2006
%
% This program can be used only for research purposes.
% This program is distributed WITHOUT ANY WARRANTY; 
% without even the implied warranty of MERCHANTABILITY 
% or FITNESS FOR A PARTICULAR PURPOSE.

clear all
close all

No = 0;
Yes = 1;

% set the maximum number of iterations
MaxIterations = 100;

fprintf('Image restoration and shape estimation\n');
fprintf('I-divergence approach\n\n');
% shape type (wave, slope, box, sin)
notValid = Yes;
while notValid
    shape = input(['What test do you want to perform?'...
        '\nChoose wave, slope, box, sin, plane, realdata: '],'s'); 
    notValid = ~strcmp(shape,'wave')&...
        ~strcmp(shape,'slope')&...
        ~strcmp(shape,'box')&...
        ~strcmp(shape,'sin')&...
        ~strcmp(shape,'plane')&...
        ~strcmp(shape,'realdata');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Data set initialization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% choose between loading or generating defocused images
LoadImage = strcmp(shape,'realdata');
if LoadImage == No
    % generate synthetic image
    fprintf('Generating synthetic dataset...');
    % set image size
    M = 101; % use odd numbers
    N = 101; % use odd numbers
    % lens parameters
    F = 12e-3; % focal length
    Fnumb = 2; % F number
    D = F/Fnumb; % aperture (effective lens diameter)
    gamma = 2e4; % calibration parameter
    p = 1./(1/F-1./[.52 .85]); % distance CCD-lens
    K = length(p); % number of focus settings
    % generate depth map
    [X,Y] = meshgrid([-(N-1)/2:(N-1)/2],[-(M-1)/2:(M-1)/2]);
    switch shape
        case 'wave'
            DepthMap = (1+cos(-(X.^2+Y.^2)/200))./...
                sqrt(X.^2/8600+Y.^2/8600+1)/2*.33+.52;
        case 'slope'
            % along X
            DepthMap = (X-min(X(:)))/(max(X(:))-min(X(:)))*.33+.52;
        case 'box'
            block1 = floor(M/5:M/5*3);
            block2 = floor(M/5*2:M/5*4);
            block3 = floor(N/10:N/2);
            block4 = floor(N/10*4:N/10*9);
            DepthMap = ones(M,N)*.52;
            DepthMap(block1,block3) = ...
                DepthMap(block1,block3)+.2;
            DepthMap(block2,block4) = ...
                DepthMap(block2,block4)+.13;
        case 'sin'
            DepthMap = (sin(X/3)+1)/2*.33+.52;
        otherwise % plane
            DepthMap = ones(M,N)*.33/2+.52;
    end
    % compute corresponding blur maps
    for k=1:K
        BlurMap(:,:,k) = gamma*D/2*abs(1-p(k).*(1/F-1./DepthMap));
    end
    % generate random texture (with different intensities)
    Radiance = [50*rand(floor(M/3),N);...
        120*rand(floor(M/3),N);...
        255*rand(M-2*floor(M/3),N)]+1; % skip zero
    % add space-varying blurring to challenge algorithm with
    % different textures (high-freq to low-freq)
    TextureDepth = [zeros(M,floor(N/3))...
        .4*ones(M,floor(N/3))...
        .6*ones(M,N-2*floor(N/3))];
    Radiance = defocusApprox(Radiance,TextureDepth);
    % generate defocused images
    I = defocusApprox(Radiance,BlurMap);
    % 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 = 2e-1; % regularization term step
    maxDelta = 1; % max amount of relative blur
    % regularization constant 
    % to avoid division by zero in the preconditioner
    epsilon = max(I(:))/255;
    fprintf('done\n');
else
    % load images
    fprintf('Loading dataset...');
    load('DataSet.mat');
    % first resize
    scalingFactor = 1/2;
    I1 = imresize(I1,scalingFactor); % limit number of computations
    I2 = imresize(I2,scalingFactor); % limit number of computations
    % 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);
    % 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
    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
    fprintf('done\n')
end

% discrete gradient step
dx = 1e-3;

% preconditioning threshold
thresh = 6e-2;

% choose whether to use the fast initialization or not
% init type (Y,N)
notValid = Yes;
while notValid
    init = input(['Do you want to use '...
                  'the fast initialization (Y/N)?'],'s'); 
    notValid = ~strcmp(init,'Y')&...
        ~strcmp(init,'N');
end

if strcmp(init,'Y')
    % initialize depth map and radiance
    [DepthMap_e,Radiance_e] = initialize(I,maxDelta,win,...
                              gamma,F,D,p,'I-Divergence');
else
    % initialize depth map with a plane
    DepthMap_e = sum(p)*F/(sum(p)-2*F)*ones(M,N);
    % initialize radiance with one image
    Radiance_e = I(:,:,1);
end

% start deblurring iterations
fprintf(['Image restoration and shape estimation'...
    ' (%3i iterations)'],MaxIterations);
for i=1:MaxIterations
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % visualize estimated radiance and depth map
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    figure(1)
    set(gcf,'name',sprintf('Iteration %3i out of %3i',i,MaxIterations));
    if LoadImage == No
        subplot(221)
        image([I(:,:,1) 200*ones(size(I,1),10) I(:,:,2)])
        axis equal
        axis off
        set(gca,'XTick',[],'YTick',[]);
        set(gca,'Box','off');
        colormap gray(256)
        title('Input #1/#2')
        drawnow;
        subplot(222)
        imagesc([Radiance_e 200*ones(size(I,1),10) Radiance])
        axis equal
        axis off
        set(gca,'XTick',[],'YTick',[]);
        set(gca,'Box','off');
        colormap gray(256)
        title('Estimated/True Radiance')
        drawnow;
        step = ceil(min([size(DepthMap_e,1) ...
               size(DepthMap_e,2)])/50);
        pix2mm = 1/200;
        [X,Y] = meshgrid(pix2mm*[1:step:size(DepthMap_e,2)],...
                pix2mm*[1:step:size(DepthMap_e,1)]);
        [X2,Y2] = meshgrid(pix2mm*(N+[1:step:size(DepthMap_e,2)]),...
                pix2mm*([1:step:size(DepthMap_e,1)]));
        axisZ = 1./(1/F-1./p);
        subplot(223)
        h1 = surf(X2,-DepthMap(1:step:end,1:step:end),-Y2,DepthMap(1:step:end,1:step:end));
        shading interp
        hold on
        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 2*N -axisZ(2) -axisZ(1) -M 0])
        set(gca,'XTick',[],'YTick',[],'ZTick',[]);
        colormap gray(256)
        title('Estimated(red)/True(blue) Depth Map')
        %legend('True Depth','Estimated Depth');
        view([-1 -2 1]);
        axis equal
        axis off
        set(h1,'EdgeColor','b');
        set(h,'EdgeColor','r');
        drawnow;
        subplot(224)
        imagesc([DepthMap_e ...
                 max(DepthMap(:))*ones(size(I,1),10) ...
                 DepthMap]);
        colormap gray(256)
        %axis off
        axis equal
        set(gca,'XTick',[],'YTick',[]);
        set(gca,'Box','off');
        title('Estimated/True Depth Map')
        drawnow;
    else
        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')
        %legend('True Depth','Estimated Depth');
        view([-1 -2 1]);
        axis equal
        axis off
        set(h,'EdgeColor','r');
        drawnow;
        subplot(224)
        imagesc(DepthMap_e);
        colormap gray(256)
        %axis off
        axis equal
        set(gca,'XTick',[],'YTick',[]);
        set(gca,'Box','off');
        title('Estimated Depth Map')
        drawnow;
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % end visualization
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % compute I-div gradient
    for k=1:K
        % compute corresponding 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 idiv 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
fprintf('\ndone\n');

return