grad

 grad - gradient, forward differences

   [gx,gy] = grad(M, options);
 or
   g = grad(M, options);

   options.bound = 'per' or 'sym'
   options.order = 1 (backward differences)
                 = 2 (centered differences)

   Works also for 3D array.
   Assme that the function is evenly sampled with sampling step 1.

   See also: div.

   Copyright (c) Gabriel Peyre

0001 function [fx,fy,fz] = grad(M, options)
0002 
0003 % grad - gradient, forward differences
0004 %
0005 %   [gx,gy] = grad(M, options);
0006 % or
0007 %   g = grad(M, options);
0008 %
0009 %   options.bound = 'per' or 'sym'
0010 %   options.order = 1 (backward differences)
0011 %                 = 2 (centered differences)
0012 %
0013 %   Works also for 3D array.
0014 %   Assme that the function is evenly sampled with sampling step 1.
0015 %
0016 %   See also: div.
0017 %
0018 %   Copyright (c) Gabriel Peyre
0019 
0020 
0021 options.null = 0;
0022 bound = getoptions(options, 'bound', 'sym');
0023 order = getoptions(options, 'order', 1);
0024 
0025 
0026 % retrieve number of dimensions
0027 nbdims = 2;
0028 if size(M,1)==1 || size(M,2)==1
0029     nbdims = 1;
0030 end
0031 if size(M,1)>1 && size(M,2)>1 && size(M,3)>1
0032     nbdims = 3;
0033 end
0034 
0035 
0036 if strcmp(bound, 'sym')    
0037     if order==1
0038         fx = M([2:end end],:,:)-M;
0039     else
0040         fx = ( M([2:end end],:,:)-M([1 1:end-1],:,:) )/2;
0041         % boundary
0042         fx(1,:,:) = M(2,:,:)-M(1,:,:);
0043         fx(end,:,:) = M(end,:,:)-M(end-1,:,:);
0044     end
0045     if nbdims>=2
0046         if order==1
0047             fy = M(:,[2:end end],:)-M;
0048         else
0049             fy = ( M(:,[2:end end],:)-M(:,[1 1:end-1],:) )/2;
0050             % boundary
0051             fy(:,1,:) = M(:,2,:)-M(:,1,:);
0052             fy(:,end,:) = M(:,end,:)-M(:,end-1,:);
0053         end
0054     end
0055     if nbdims>=3
0056         if order==1
0057             fz = M(:,:,[2:end end])-M;
0058         else
0059             fz = ( M(:,:,[2:end end])-M(:,:,[1 1:end-1]) )/2;
0060             % boundary
0061             fz(:,:,1) = M(:,:,2)-M(:,:,1);
0062             fz(:,:,end) = M(:,:,end)-M(:,:,end-1);
0063         end
0064     end
0065 else
0066     if order==1
0067         fx = M([2:end 1],:,:)-M;
0068     else
0069         fx = ( M([2:end 1],:,:)-M([end 1:end-1],:,:) )/2;
0070     end
0071     if nbdims>=2
0072         if order==1
0073             fy = M(:,[2:end 1],:)-M;
0074         else
0075             fy = ( M(:,[2:end 1],:)-M(:,[end 1:end-1],:) )/2;
0076         end
0077     end
0078     if nbdims>=3
0079         if order==1
0080             fz = M(:,:,[2:end 1])-M;
0081         else
0082             fz = ( M(:,:,[2:end 1])-M(:,:,[end 1:end-1]) )/2;
0083         end
0084     end
0085 end
0086 
0087 if nargout==1
0088     if nbdims==2
0089         fx = cat(3,fx,fy);
0090     elseif nbdims==3
0091         fx = cat(4,fx,fy,fz);
0092     end
0093 end