Code covered by the BSD License

### John D'Errico (view profile)

23 Jan 2006 (Updated )

Vectorized & fully sparse 1-d, 2-d, & 3-d integrated gradients

```function fhat = intgrad2(fx,fy,dx,dy,f11)
%
% arguments: (input)
%  fx,fy - (ny by nx) arrays, as gradient would have produced. fx and
%          fy must both be the same size. Note that x is assumed to
%          be the column dimension of f, in the meshgrid convention.
%
%          nx and ny must both be at least 2.
%
%          fx and fy will be assumed to contain consistent gradient
%          information. If they are inconsistent, then the generated
%          gradient will be solved for in a least squares sense.
%
%          Central differences will be used where possible.
%
%     dx - (OPTIONAL) scalar or vector - denotes the spacing in x
%          if dx is a scalar, then spacing in x (the column index
%          of fx and fy) will be assumed to be constant = dx.
%          if dx is a vector, it denotes the actual coordinates
%          of the points in x (i.e., the column dimension of fx
%          and fy.) length(dx) == nx
%
%          DEFAULT: dx = 1
%
%     dy - (OPTIONAL) scalar or vector - denotes the spacing in y
%          if dy is a scalar, then the spacing in x (the row index
%          of fx and fy) will be assumed to be constant = dy.
%          if dy is a vector, it denotes the actual coordinates
%          of the points in y (i.e., the row dimension of fx
%          and fy.) length(dy) == ny
%
%          DEFAULT: dy = 1
%
%    f11 - (OPTIONAL) scalar - defines the (1,1) eleemnt of fhat
%          after integration. This is just the constant of integration.
%
%          DEFAULT: f11 = 0
%
% arguments: (output)
%   fhat - (nx by ny) array containing the integrated gradient
%
% Example usage 1: (Note x is uniform in spacing, y is not.)
%  xp = 0:.1:1;
%  yp = [0 .1 .2 .4 .8 1];
%  [x,y]=meshgrid(xp,yp);
%  f = exp(x+y) + sin((x-2*y)*3);
%
%  Time required was 0.06 seconds
%
% Example usage 2: Large grid, 101x101
%  xp = 0:.01:1;
%  yp = 0:.01:1;
%  [x,y]=meshgrid(xp,yp);
%  f = exp(x+y) + sin((x-2*y)*3);
%
%  Time required was 4 seconds

% Author; John D'Errico
% Current release: 2
% Date of release: 1/27/06

% size
if (length(size(fx))>2) || (length(size(fy))>2)
error 'fx and fy must be 2d arrays'
end
[ny,nx] = size(fx);
if (nx~=size(fy,2)) || (ny~=size(fy,1))
error 'fx and fy must be the same sizes.'
end
if (nx<2) || (ny<2)
error 'fx and fy must be at least 2x2 arrays'
end

% supply defaults if needed
if (nargin<3) || isempty(dx)
% default x spacing is 1
dx = 1;
end
if (nargin<4) || isempty(dy)
% default y spacing is 1
dy = 1;
end
if (nargin<5) || isempty(f11)
% default integration constant is 0
f11 = 0;
end

% if scalar spacings, expand them to be vectors
dx=dx(:);
if length(dx) == 1
dx = repmat(dx,nx-1,1);
elseif length(dx)==nx
% dx was a vector, use diff to get the spacing
dx = diff(dx);
else
error 'dx is not a scalar or of length == nx'
end
dy=dy(:);
if length(dy) == 1
dy = repmat(dy,ny-1,1);
elseif length(dy)==ny
% dy was a vector, use diff to get the spacing
dy = diff(dy);
else
error 'dy is not a scalar or of length == ny'
end

if (length(f11) > 1) || ~isnumeric(f11) || isnan(f11) || ~isfinite(f11)
error 'f11 must be a finite scalar numeric variable.'
end

% build gradient design matrix, sparsely. Use a central difference
% in the body of the array, and forward/backward differences along
% the edges.

% A will be the final design matrix. it will be sparse.
% The unrolling of F will be with row index running most rapidly.
rhs = zeros(2*nx*ny,1);
% but build the array elements in Af
Af = zeros(2*nx*ny,6);
L = 0;

% do the leading edge in x, forward difference
indx = 1;
indy = (1:ny)';
ind = indy + (indx-1)*ny;
rind = repmat(L+(1:ny)',1,2);
cind = [ind,ind+ny];
dfdx = repmat([-1 1]./dx(1),ny,1);
Af(L+(1:ny),:) = [rind,cind,dfdx];
rhs(L+(1:ny)) = fx(:,1);
L = L+ny;

% interior partials in x, central difference
if nx>2
[indx,indy] = meshgrid(2:(nx-1),1:ny);
indx = indx(:);
indy = indy(:);
ind = indy + (indx-1)*ny;
m = ny*(nx-2);

rind = repmat(L+(1:m)',1,2);
cind = [ind-ny,ind+ny];

dfdx = 1./(dx(indx-1)+dx(indx));
dfdx = dfdx*[-1 1];

Af(L+(1:m),:) = [rind,cind,dfdx];

rhs(L+(1:m)) = fx(ind);

L = L+m;
end

% do the trailing edge in x, backward difference
indx = nx;
indy = (1:ny)';
ind = indy + (indx-1)*ny;
rind = repmat(L+(1:ny)',1,2);
cind = [ind-ny,ind];
dfdx = repmat([-1 1]./dx(end),ny,1);
Af(L+(1:ny),:) = [rind,cind,dfdx];
rhs(L+(1:ny)) = fx(:,end);
L = L+ny;

% do the leading edge in y, forward difference
indx = (1:nx)';
indy = 1;
ind = indy + (indx-1)*ny;
rind = repmat(L+(1:nx)',1,2);
cind = [ind,ind+1];
dfdy = repmat([-1 1]./dy(1),nx,1);
Af(L+(1:nx),:) = [rind,cind,dfdy];
rhs(L+(1:nx)) = fy(1,:).';
L = L+nx;

% interior partials in y, use a central difference
if ny>2
[indx,indy] = meshgrid(1:nx,2:(ny-1));
indx = indx(:);
indy = indy(:);
ind = indy + (indx-1)*ny;
m = nx*(ny-2);

rind = repmat(L+(1:m)',1,2);
cind = [ind-1,ind+1];

dfdy = 1./(dy(indy-1)+dy(indy));
dfdy = dfdy*[-1 1];

Af(L+(1:m),:) = [rind,cind,dfdy];

rhs(L+(1:m)) = fy(ind);

L = L+m;
end

% do the trailing edge in y, backward diffeence
indx = (1:nx)';
indy = ny;
ind = indy + (indx-1)*ny;
rind = repmat(L+(1:nx)',1,2);
cind = [ind-1,ind];
dfdy = repmat([-1 1]./dy(end),nx,1);
Af(L+(1:nx),:) = [rind,cind,dfdy];
rhs(L+(1:nx)) = fy(end,:).';

% finally, we can build the rest of A itself, in its sparse form.
A = sparse(Af(:,1:2),Af(:,3:4),Af(:,5:6),2*nx*ny,nx*ny);

% Finish up with f11, the constant of integration.
% eliminate the first unknown, as f11 is given.
rhs = rhs - A(:,1)*f11;

% Solve the final system of equations. They will be of
% full rank, due to the explicit integration constant.
% Just use sparse \
fhat = A(:,2:end)\rhs;
fhat = reshape([f11;fhat],ny,nx);

```