No BSD License  

Highlights from
Inverse Gradient

Inverse Gradient

by

 

23 Jun 2005 (Updated )

The inverse function of gradient().

invgradient(dFx,dFy,varargin)
%INVGRADIENT Basically the opposite of gradient(), aka Inverse Gradient.
%   F_bar = INVGRADIENT(dFx,dFy) reverses GRADIENT(F), where unit  
%   spacing is assumed.
%   F_bar = INVGRADIENT(dFx,dFy,H), where H is scalar, uses H as the
%   spacing between points.
%   F_bar = INVGRADIENT(dFx,dFy,Hx,Hy) is basically the same, except that
%   both x- and y-spacings are specified.
%   Obviously, the DC offset of the integrated function is arbitrary. Also,
%   INVGRADIENT is limited to 2-dimensions at the moment, ie no
%   3-dimensional matrices, or vectors. If dFx and dFy aren't
%   self-consistent then no exact solution can be integrated. In this case
%   the optimal solution, in the least squares sense, will be found. So 
%   if you do F_bar = INVGRADIENT(dFx,dFy), F_bar will be found such that 
%   doing [dFx_ dFy_] = GRADIENT(F_bar) will give you dFx_ and dFy_ that 
%   are as close as possible to the original dFx and dFy.
%  
%   Patrick Lu
%   April 5, 2005
%   If this doesn't work right, email patlu@nospam.stanford.edu, w/o the nospam.

function F_bar = invgradient(dFx,dFy,varargin)
[rows cols] = size(dFx);
if(length(varargin)==0)
    xspacing = 1;
    yspacing = 1;
elseif(length(varargin)==1)
    xspacing = varargin{1};
    yspacing = varargin{1};
else
    xspacing = varargin{1};
    yspacing = varargin{2};
end

if (rows<3 | cols < 3)
    fprintf('Please keep the input size at least 3x3\n');    
end
%create Gx, which takes the x-derivative such that F*Gx = dFx
Gx = [-1; 1; zeros(cols-2,1)];

middle_column = [-.5;0;.5;zeros(cols-3,1)];
for k = 1:cols-2

    Gx = [Gx,middle_column]; 
    %permute
    middle_column=circshift(middle_column,1);
end
Gx = [Gx, [zeros(cols-2,1); -1; 1]];
Gx = Gx/xspacing;


%create Gy, which takes the y-derivative such that Gy*F = dFy
Gy = [-1, 1, zeros(1,rows-2)];
middle_row = [-.5,0,.5,zeros(1,rows-3)];
for k=1:rows-2
    Gy = [Gy;middle_row];
    %permute
    middle_row=circshift(middle_row,[0 1]);
end
Gy = [Gy; [zeros(1,rows-2),-1,1]];
Gy = Gy/yspacing;


%put them into G2
G2 = [];
for k = 1:cols
    G2_row = [Gx(:,k) zeros(cols,rows-1)];
    G2_row = reshape(G2_row',1,rows*cols);
    for l=1:rows
        G2 = [G2;G2_row];
        G2_row = circshift(G2_row,[0 1]);
    end
end

G2_2 = [];

for k = 1:cols
    G2_2 = blkdiag(G2_2,Gy);
end

G2 = [G2;G2_2];

F_bar = pinv(G2)*[reshape(dFx,rows*cols,1);reshape(dFy,rows*cols,1)];
F_bar = reshape(F_bar,rows,cols);

Contact us