% FRANKOTCHELLAPPA  - Generates integrable surface from gradients
%
% An implementation of Frankot and Chellappa'a algorithm for constructing
% an integrable surface from gradient information.
%
% Usage:      z = frankotchellappa(dzdx,dzdy)
%
% Arguments:  dzdx,  - 2D matrices specifying a grid of gradients of z
%             dzdy     with respect to x and y.
%
% Returns:    z      - Inferred surface heights.

% Reference:
%
% Robert T. Frankot and Rama Chellappa
% A Method for Enforcing Integrability in Shape from Shading
% IEEE PAMI Vol 10, No 4 July 1988. pp 439-451
%
% Note this code just implements the surface integration component of the
% paper (Equation 21 in the paper).  It does not implement their shape from
% shading algorithm.

% Copyright (c) 2004 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.

% October 2004

function z = frankotchellappa(dzdx,dzdy)
    
    if ~all(size(dzdx) == size(dzdy))
      error('Gradient matrices must match');
    end

    [rows,cols] = size(dzdx);
    
    % The following sets up matrices specifying frequencies in the x and y
    % directions corresponding to the Fourier transforms of the gradient
    % data.  They range from -0.5 cycles/pixel to + 0.5 cycles/pixel. The
    % fiddly bits in the line below give the appropriate result depending on
    % whether there are an even or odd number of rows and columns
    
    [wx, wy] = meshgrid(([1:cols]-(fix(cols/2)+1))/(cols-mod(cols,2)), ...
			([1:rows]-(fix(rows/2)+1))/(rows-mod(rows,2)));
    
    % Quadrant shift to put zero frequency at the appropriate edge
    wx = ifftshift(wx); wy = ifftshift(wy);

    DZDX = fft2(dzdx);   % Fourier transforms of gradients
    DZDY = fft2(dzdy);

    % Integrate in the frequency domain by phase shifting by pi/2 and
    % weighting the Fourier coefficients by their frequencies in x and y and
    % then dividing by the squared frequency.  eps is added to the
    % denominator to avoid division by 0.
    
    Z = (-j*wx.*DZDX -j*wy.*DZDY)./(wx.^2 + wy.^2 + eps);  % Equation 21
    
    z = real(ifft2(Z));  % Reconstruction