### Patrick Lu (view profile)

23 Jun 2005 (Updated )

%   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.

[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);