mediansmoothing.m

3D shape reconstruction matlab code. It used shape from defocus technique with least squares. You ca

function F = mediansmoothing(f,visualize)
% function F = mediansmoothing(f,visualize)
%
% This function performs median filtering and then
% Gaussian smoothing
%
% 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.

On = 1;
Off = 0;

if nargin<2
    % turn visualization on or off
    visualize = Off;
elseif (visualize ~= On) & (visualize ~= Off)
    % turn visualization on or off
    visualize = Off;
end

% window of median filtering and Gaussian smoothing
win = 7;

% median filtering
k = 0;
F = zeros(win*win,(size(f,1)-win+1)*(size(f,2)-win+1));
for i=0:win-1
    for j=0:win-1
        k = k+1;
        F(k,:) = reshape(f(i+[1:size(f,1)-win+1],...
                 j+[1:size(f,2)-win+1]),...
                 1,(size(f,1)-win+1)*(size(f,2)-win+1));
    end
end
F = median(F,1);
Temp = reshape(F,size(f,1)-win+1,size(f,2)-win+1);
F = median(f(:))*ones(size(f));
F(floor((win-1)/2)+[1:size(f,1)-win+1],...
  floor((win-1)/2)+[1:size(f,2)-win+1]) = Temp;

% smoothing
[x,y] = meshgrid([-(win-1)/2:(win-1)/2]/(win-1)*5,...
                 [-(win-1)/2:(win-1)/2]/(win-1)*5);
weights = exp(-x.^2-y.^2);
weights = weights/sum(weights(:));
F = conv2(F-mean(F(:)),weights,'same')+mean(F(:));

if visualize
   figure(2)
   subplot(121)
   imagesc(f);
   title('Estimated depth map');
   colormap gray(256)
   axis equal
   axis off
   subplot(122)
   imagesc(F);
   title('Estimated depth map (smoothed)');
   colormap gray(256)
   axis equal
   axis off
   drawnow;
end
return