defocusunif.m

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

function [I] = defocusUnif(Radiance,sigma)
% Generate defocused images (uniform blur)
%
% 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.

sigmamax = max(sigma(:));
win = ceil(sigmamax*3);

% extend domain
Radiance_ext = ones(size(Radiance,1)+...
    2*win,size(Radiance,2)+...
    2*win,size(Radiance,3));
Radiance_ext(win+1:end-win,win+1:end-win,:) = Radiance;
if win>0
    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,:);
end

[M,N,K] = size(Radiance_ext);
K = length(sigma); % number of focus settings
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]);
% allocate space for the defocused images
I = ones(M,N,K);

% compute a defocused image for each focus setting
for k=1:K
    % integrate over pixel area
    if (sigma(k)>tolerance)
        wsi = wi/sqrt(2)/sigma(k);
        wsj = wj/sqrt(2)/sigma(k);
        wdx = .5/sqrt(2)/sigma(k);
        derf = erf(wsi+wdx)-erf(wsi-wdx);
        PSF = .25*derf.*derf';
    else
        PSF = (wi==0).*(wj==0); % Dirac delta
    end
    I(:,:,k) = conv2(Radiance_ext(:,:,k),PSF,'same');
end
% restore original domain
I = I(win+1:M-win,win+1:N-win,:);

return