% 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 2004function 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