% SHAPELETSURF - reconstructs surface from surface normals using shapelets
%
% Function reconstructs an estimate of a surface from its surface normals by
% correlating the surface normals with that those of a bank of shapelet
% basis functions. The correlation results are summed to produce the
% reconstruction.  The sumation of shapelet basis functions results in an
% implicit integration of the surface while enforcing surface continuity.
%
% Note that the reconstruction is only valid up to a scale factor. However
% the reconstruction process is very robust to noise and to missing data
% values.  Reconstructions (up to positive/negative shape ambiguity) are
% possible where there is an ambiguity of pi in tilt values.  Low quality
% reconstructions are also possible with just slant, or just tilt data
% alone.
% 
%
% Usage:
%  recsurf = shapletsurf(slant, tilt, nscales, minradius, mult, opt)
%                                       6        1        2
% Arguments:
%            slant     - 2D array of surface slant values across image.
%            tilt      - 2D array of surface tilt values.
%            nscales   - number of shapelet scales to use.
%            minsigma  - sigma of smallest scale Gaussian shapelet.
%            mult      - scaling factor between successive shapelets.
%
%   opt can be the string:
%           'slanttilt' - reconstruct using both slant and tilt (default).
%           'tiltamb'   - reconstruct assuming tilt ambiguity of pi.
%           'slant'     - reconstruct with slant only.
%           'tilt'      - reconstruct with tilt only.
%  
% Returns:
%           recsurf     - reconstructed surface.
%
% Remember when viewing the surface you should use 
%   >> axis ij 
% So that the surface corresponds to the image slant and tilt axes
%
% References:
%
% Peter Kovesi, "Surface Normals to Surfaces via Shapelets"
% Proceedings Australia-Japan Advanced Workshop on Computer Vision
% Adelaide, 9-11 September 2003
%
% Peter Kovesi, "Shapelets Correlated with Surface Normals Produce
% Surfaces". Technical Report 03-003, October 2003.
% http://www.csse.uwa.edu.au/~pk/research/pkpapers/shapelets-03-002.pdf

% Copyright (c) 2003-2005 Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% http://www.csse.uwa.edu.au/
% 
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
% 
% The above copyright notice and this permission notice shall be included in 
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.

% July      2003 - Original version.
% September 2003 - Correction to reconstruction with tilt ambiguity.
% October   2003 - Changed to use Gaussian shapelets.
% February  2004 - Convolutions done via fft for speed.
% March     2004 - Padding of slant and tilt data to accommodate large filters.

function recsurf = shapeletsurf(varargin)
    
    [slant, tilt, nscales, minsigma, mult, opt] = checkargs(varargin(:));
    
    if strcmp(opt,'tiltamb')  % If we have an ambiguity of pi in the tilt 
        tilt = tilt*2;        % work with doubled angles.
    end

    [rows,cols] = size(slant);

    % Check size of largest filter and, if necessary, pad slant and tilt
    % data with zeros so that the largest filter can be accommodated.
    % If this is not done wraparound in the convolutions via the FFT produce
    % artifacts in the reconstruction.
    % Treat max size as +- 3 sigma
    maxsize = ceil(6*minsigma*mult^(nscales-1)); 
    paddingapplied = 0;

    if rows < maxsize | cols < maxsize  % padding needed
      paddingapplied = 1;
      rowpad = max(0,round(maxsize-rows));
      colpad = max(0,round(maxsize-cols));
      fprintf('Warning: To accommodate the largest filter size the\n');
      fprintf('slant and tilt data is being padded from %dx%d to %dx%d \n', ...
                      rows,cols, rows+rowpad, cols+colpad);

      slant = [slant zeros(rows, colpad)
               zeros(rowpad, cols+colpad)];
      tilt  = [tilt  zeros(rows, colpad)
               zeros(rowpad, cols+colpad)];
      origrows = rows; origcols = cols;    % Remember original size.
      [rows,cols] = size(slant);           % Update current size.
    end
    
    % Precompute some values for speed of execution.  Note that because we
    % generally want to use shapelets at quite large scales relative to the
    % size of the image correlations are done in the frequency domain for
    % speed. (Note that in the code below the conjugate of the fft is used
    % because we want the correlation, not convolution.)
    surfgrad = tan(slant);  SURFGRAD= fft2(surfgrad);
    sintilt = sin(tilt);    SINTILT = fft2(sintilt);
    costilt = cos(tilt);    COSTILT = fft2(costilt);
    SURFGRADSINTILT = fft2(surfgrad.*sintilt);  
    SURFGRADCOSTILT = fft2(surfgrad.*costilt);
    
    for s = 1:nscales
        fprintf('scale %d out of %d\r',s,nscales);
        
        % Use a Gaussian filter shape as the shapelet basis function
        % as the phase distortion in the reconstruction should be zero.
        f = gaussianf(minsigma*mult^(s-1), rows, cols);

        [fdx,fdy] = gradient(f);                  % filter gradients
        [fslant,ftilt] = grad2slanttilt(fdx,fdy); % filter slants and tilts
        
        if strcmp(opt,'tiltamb')
            ftilt = ftilt*2;
        end    
        
        % Now perform the correlations (via the fft) as required depending
        % on the options selected. 
        sinftilt = sin(ftilt);   
        cosftilt = cos(ftilt);   
        filtgrad = tan(fslant);   
        
        filtgradsintilt = filtgrad.*sinftilt;
        filtgradcostilt = filtgrad.*cosftilt;
        
        if strcmp(opt,'slanttilt')     % both slant and tilt data available
            FILTGRADSINTILT = fft2(filtgrad.*sinftilt);
            FILTGRADCOSTILT = fft2(filtgrad.*cosftilt);
            fim{s} = real(ifft2(conj(FILTGRADCOSTILT) .* SURFGRADCOSTILT)) + ...
                     real(ifft2(conj(FILTGRADSINTILT) .* SURFGRADSINTILT));
            
        elseif strcmp(opt,'tiltamb')   % assume tilt ambiguity of pi
            FILTGRADSINTILT = fft2(filtgrad.*sinftilt);
            FILTGRADCOSTILT = fft2(filtgrad.*cosftilt);
            FILTGRAD = fft2(filtgrad);
            fim{s} = (real(ifft2(conj(FILTGRADCOSTILT) .* SURFGRADCOSTILT)) + ...
                      real(ifft2(conj(FILTGRADSINTILT) .* SURFGRADSINTILT)) + ...
                      real(ifft2(conj(FILTGRAD) .* SURFGRAD)))/2;    
            
        elseif strcmp(opt,'tilt');     % tilt only reconstruction
            SINFTILT = fft2(sinftilt);
            COSFTILT = fft2(cosftilt);
            fim{s} = real(ifft2(conj(COSFTILT) .* COSTILT)) + ...
                     real(ifft2(conj(SINFTILT) .* SINTILT));      
            
        elseif strcmp(opt,'slant');    % just use slant and ignore tilt
            FILTGRAD = fft2(filtgrad);
            fim{s} = real(ifft2(conj(FILTGRAD) .* SURFGRAD));  
        end
        
    end
    fprintf('\n');
    
    % Reconstruct by adding filtered outputs
    
    recsurf = zeros(size(slant));
    for s = 1:nscales
        recsurf = recsurf + fim{s};
    end


    if paddingapplied  % result is padded - extract the bit we want.
	recsurf = recsurf(1:origrows, 1:origcols);
    end
    
%-------------------------------------------------------------------------
% Function to generate a Gaussian filter for use as a shapelet.
% Usage:
%      f = gaussian(sigma, rows, cols)
%
% Arguments:
%      sigma      - standard deviation of Gaussian
%      rows, cols - size of filter to create

function f = gaussianf(sigma, rows, cols)

    [x,y] = meshgrid( [1:cols]-(fix(cols/2)+1), [1:rows]-(fix(rows/2)+1));
    r = sqrt(x.^2 + y.^2);
    f = fftshift(exp(-r.^2/(2*sigma^2)));

%--------------------------------------------------------------------------
% Function to check argument values and set defaults

function [slant, tilt, nscales, minradius, mult, opt] = checkargs(arg);
    
    if length(arg)<5 
        error('too few arguments');
    elseif length(arg)>6
        error('too many arguments');
    end
    
    slant     = arg{1};
    tilt      = arg{2};
    nscales   = arg{3};
    minradius = arg{4};
    mult      = arg{5};
    
    if length(arg) == 5
        opt = 'slanttilt'; % Default is to assume both slant and tilt values
                           % are valid.
    else
        opt = arg{6};
    end
    
    if ~all(size(slant)==size(tilt))
        error('slant and tilt matrices must match');
    end
    
    if nscales < 1
        error('number of scales must be 1 or greater');
    end
    
    if minradius < .1
        error('minimum radius of shapelet must be greater than .1');
    end
    
    if mult < 1
        error('scaling factor between successive filters should be greater than 1');
    end