defocusapprox.m

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

function [I,F,I_dx] = defocusApprox(Radiance,BlurMap,I0,dBlurMap)
% Generate defocused images and deblurring gradient using 
% approximate-but-fast kernel
%
% Check the approximation by running the test code in the very end!
%
% 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.

if nargin<4
    sigmamax = max(BlurMap(:));
else
    sigmamax = max([BlurMap(:);dBlurMap(:)]);
end
win = ceil(sigmamax*3);

% extend domain
Radiance_ext = ones(size(Radiance)+2*win);
Radiance_ext(win+1:end-win,win+1:end-win) = Radiance;
Radiance_ext(1:win,:) = [Radiance(win+1:-1:2,win+1:-1:2) ...
    Radiance(win+1:-1:2,:) ...
    Radiance(win+1:-1:2,end-1:-1:end-win)];
Radiance_ext(end-win+1:end,:) = [...
    Radiance(end-1:-1:end-win,end-1:-1:end-win) ...
    Radiance(end-1:-1:end-win,:) ...
    Radiance(end-1:-1:end-win,end-1:-1:end-win)];
Radiance_ext(win+1:end-win,1:win) = Radiance(:,win+1:-1:2);
Radiance_ext(win+1:end-win,end-win+1:end) = ...
    Radiance(:,end-1:-1:end-win);

[M,N] = size(Radiance_ext);
K = size(BlurMap,3); % number of focus settings

BlurMap_ext = ones(M,N,K);
BlurMap_ext(win+1:end-win,win+1:end-win,:) = BlurMap;
BlurMap_ext(1:win,:,:) = [BlurMap(win+1:-1:2,win+1:-1:2,:) ...
    BlurMap(win+1:-1:2,:,:) ...
    BlurMap(win+1:-1:2,end-1:-1:end-win,:)];
BlurMap_ext(end-win+1:end,:,:) = [...
    BlurMap(end-1:-1:end-win,end-1:-1:end-win,:) ...
    BlurMap(end-1:-1:end-win,:,:) ...
    BlurMap(end-1:-1:end-win,end-1:-1:end-win,:)];
BlurMap_ext(win+1:end-win,1:win,:) = BlurMap(:,win+1:-1:2,:);
BlurMap_ext(win+1:end-win,end-win+1:end,:) = ...
    BlurMap(:,end-1:-1:end-win,:);

tolerance = 1e-5; % threshold for Dirac delta approx. of kernel

W = win*2+1; % max size of the kernel window
% compute window coordinates
[wi,wj] = meshgrid([-win:win],[-win:win]);
% precompute integrating grid
wip4 = (wi+.5)/sqrt(2)*1.694;
wim4 = (wi-.5)/sqrt(2)*1.694;
% compute pixel coordinates
[ti,tj] = meshgrid([win+1:N-win],[win+1:M-win]);
ti = ti(:);
tj = tj(:);
% allocate memory for the kernel
H = zeros((M-W+1)*(N-W+1),W*W);
% allocate space for the defocused images
I = ones(M,N,K);
% allocate space for the kernel
Kernel = zeros(M,N,W,W);
if nargin>2
    I0_ext = ones(M,N,K);
    I0_ext(win+1:end-win,win+1:end-win,:) = I0;
    I0_ext(1:win,:,:) = [I0(win+1:-1:2,win+1:-1:2,:) ...
        I0(win+1:-1:2,:,:) ...
        I0(win+1:-1:2,end-1:-1:end-win,:)];
    I0_ext(end-win+1:end,:,:) = [...
        I0(end-1:-1:end-win,end-1:-1:end-win,:) ...
        I0(end-1:-1:end-win,:,:) ...
        I0(end-1:-1:end-win,end-1:-1:end-win,:)];
    I0_ext(win+1:end-win,1:win,:) = I0(:,win+1:-1:2,:);
    I0_ext(win+1:end-win,end-win+1:end,:) = ...
        I0(:,end-1:-1:end-win,:);
    % initialize Lucy-Richardson variables
    H0 = zeros(M,N,K);
    F0 = zeros(M,N,K);
    F = ones(M,N);
    Ratio = zeros(M,N,K);
end
if nargin>3
    dBlurMap_ext = ones(M,N,K);
    dBlurMap_ext(win+1:end-win,win+1:end-win,:) = dBlurMap;
    dBlurMap_ext(1:win,:,:) = [...
        dBlurMap(win+1:-1:2,win+1:-1:2,:) ...
        dBlurMap(win+1:-1:2,:,:) ...
        dBlurMap(win+1:-1:2,end-1:-1:end-win,:)];
    dBlurMap_ext(end-win+1:end,:,:) = [...
        dBlurMap(end-1:-1:end-win,end-1:-1:end-win,:) ...
        dBlurMap(end-1:-1:end-win,:,:) ...
        dBlurMap(end-1:-1:end-win,end-1:-1:end-win,:)];
    dBlurMap_ext(win+1:end-win,1:win,:) = ...
        dBlurMap(:,win+1:-1:2,:);
    dBlurMap_ext(win+1:end-win,end-win+1:end,:) = ...
        dBlurMap(:,end-1:-1:end-win,:);
    % initialize gradient wrt shape
    I_dx = ones(M,N,K);
    H_dx = zeros((M-W+1)*(N-W+1),W*W);
end

transpose = reshape([1:W^2],W,W)';

i = ti(1:(M-W+1)*(N-W+1));
j = tj(1:(M-W+1)*(N-W+1));
ind = j+(i-1)*size(BlurMap_ext,1)-...
      size(BlurMap_ext,1)*size(BlurMap_ext,2);

% rearrange radiance
R = zeros((M-W+1)*(N-W+1),W*W);
index1 = j+(i-1)*size(Radiance_ext,1);
[u,v] = meshgrid([-win:win],[-win:win]);
index2 = u(transpose(:))'+v(transpose(:))'*size(Radiance_ext,1);
index = index1*ones(1,(2*win+1)*(2*win+1))+...
        ones((M-W+1)*(N-W+1),1)*index2;
R = Radiance_ext(index);

% compute a defocused image for each focus setting
for k=1:K
    
    ind = ind+size(BlurMap_ext,1)*size(BlurMap_ext,2);
    sigma = BlurMap_ext(ind);
    dirac = sigma<=tolerance;
    isigma = (1-dirac)./(sigma+dirac);
    wsip = isigma(:)*wip4(:)';
    wsim = isigma(:)*wim4(:)';
    derf = (1-exp(-abs(wsip))).*sign(wsip)-...
           (1-exp(-abs(wsim))).*sign(wsim);
    H = .25*derf.*derf(:,transpose(:))+...
        dirac(:)*((wi(:)==0).*(wj(:)==0))';

    if nargin>3
        sigma_dx = dBlurMap_ext(ind);
        dirac = sigma_dx<=tolerance;
        isigma_dx = (1-dirac)./(sigma_dx+dirac);
        wsip = isigma_dx(:)*wip4(:)';
        wsim = isigma_dx(:)*wim4(:)';
        derf = (1-exp(-abs(wsip))).*sign(wsip)-...
               (1-exp(-abs(wsim))).*sign(wsim);
        H_dx = .25*derf.*derf(:,transpose(:))+...
               dirac(:)*((wi(:)==0).*(wj(:)==0))';
    end
    
    I(win+1:M-win,win+1:N-win,k) = ...
        reshape(sum(H.*R,2),M-W+1,N-W+1);
    if nargin>3
        I_dx(win+1:M-win,win+1:N-win,k) = ...
            reshape(sum(H_dx.*R,2),M-W+1,N-W+1);
    end
    Kernel(win+1:M-win,win+1:N-win,1:W,1:W) = ...
        reshape(H,M-W+1,N-W+1,W,W);

    if nargin>2
        % compute Lucy-Richardson multiplying factor
        Ratio(win+1:M-win,win+1:N-win) = ...
            I0_ext(win+1:M-win,win+1:N-win,k)./...
            I(win+1:M-win,win+1:N-win,k);
        for w1=1:W
            for w2=1:W
                i = [win+1:N-win]-(w1-win-1);
                j = [win+1:M-win]-(w2-win-1);
                Ker = Kernel(j,i,w1,w2);
                H0(win+1:M-win,win+1:N-win,k) = ...
                    H0(win+1:M-win,win+1:N-win,k)+Ker;
                F0(win+1:M-win,win+1:N-win,k) = ...
                    F0(win+1:M-win,win+1:N-win,k)+Ker.*Ratio(j,i);
            end
        end
    end

end
if nargin>2
    % compute Lucy-Richardson multiplying factor
    HH = sum(H0(win+1:M-win,win+1:N-win,:),3);
    F = sum(F0(win+1:M-win,win+1:N-win,:),3)./HH;
end
if nargin>3
    % restore original domain
    I_dx = I_dx(win+1:M-win,win+1:N-win,:);
end
% restore original domain
I = I(win+1:M-win,win+1:N-win,:);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% testing the approximation of erf
% via the exponential
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x = [-5:.1:5];
Original = erf(x);
Approximate = (1-exp(-abs(x)*1.694)).*sign(x);

figure(1)
plot(x,Original,'b')
hold on
plot(x,Approximate,'r');
hold off
legend('Original','Approximate');
drawnow