octave_packages/m/general/gradient.m

   1 ## Copyright (C) 2000-2012 Kai Habel
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn  {Function File} {@var{dx} =} gradient (@var{m})
  21 ## @deftypefnx {Function File} {[@var{dx}, @var{dy}, @var{dz}, @dots{}] =} gradient (@var{m})
  22 ## @deftypefnx {Function File} {[@dots{}] =} gradient (@var{m}, @var{s})
  23 ## @deftypefnx {Function File} {[@dots{}] =} gradient (@var{m}, @var{x}, @var{y}, @var{z}, @dots{})
  24 ## @deftypefnx {Function File} {[@dots{}] =} gradient (@var{f}, @var{x0})
  25 ## @deftypefnx {Function File} {[@dots{}] =} gradient (@var{f}, @var{x0}, @var{s})
  26 ## @deftypefnx {Function File} {[@dots{}] =} gradient (@var{f}, @var{x0}, @var{x}, @var{y}, @dots{})
  27 ##
  28 ## Calculate the gradient of sampled data or a function.  If @var{m}
  29 ## is a vector, calculate the one-dimensional gradient of @var{m}.  If
  30 ## @var{m} is a matrix the gradient is calculated for each dimension.
  31 ##
  32 ## @code{[@var{dx}, @var{dy}] = gradient (@var{m})} calculates the one
  33 ## dimensional gradient for @var{x} and @var{y} direction if @var{m} is a
  34 ## matrix.  Additional return arguments can be use for multi-dimensional
  35 ## matrices.
  36 ##
  37 ## A constant spacing between two points can be provided by the
  38 ## @var{s} parameter.  If @var{s} is a scalar, it is assumed to be the spacing
  39 ## for all dimensions.
  40 ## Otherwise, separate values of the spacing can be supplied by
  41 ## the @var{x}, @dots{} arguments.  Scalar values specify an equidistant
  42 ## spacing.
  43 ## Vector values for the @var{x}, @dots{} arguments specify the coordinate for
  44 ## that
  45 ## dimension.  The length must match their respective dimension of @var{m}.
  46 ##
  47 ## At boundary points a linear extrapolation is applied.  Interior points
  48 ## are calculated with the first approximation of the numerical gradient
  49 ##
  50 ## @example
  51 ## y'(i) = 1/(x(i+1)-x(i-1)) * (y(i-1)-y(i+1)).
  52 ## @end example
  53 ##
  54 ## If the first argument @var{f} is a function handle, the gradient of the
  55 ## function at the points in @var{x0} is approximated using central
  56 ## difference.  For example, @code{gradient (@@cos, 0)} approximates the
  57 ## gradient of the cosine function in the point @math{x0 = 0}.  As with
  58 ## sampled data, the spacing values between the points from which the
  59 ## gradient is estimated can be set via the @var{s} or @var{dx},
  60 ## @var{dy}, @dots{} arguments.  By default a spacing of 1 is used.
  61 ## @seealso{diff, del2}
  62 ## @end deftypefn
  63
  64 ## Author:  Kai Habel <kai.habel@gmx.de>
  65 ## Modified: David Bateman <dbateman@free.fr> Added NDArray support
  66
  67 function varargout = gradient (m, varargin)
  68
  69   if (nargin < 1)
  70     print_usage ();
  71   endif
  72
  73   nargout_with_ans = max(1,nargout);
  74   if (ismatrix (m))
  75     [varargout{1:nargout_with_ans}] = matrix_gradient (m, varargin{:});
  76   elseif (isa (m, "function_handle"))
  77     [varargout{1:nargout_with_ans}] = handle_gradient (m, varargin{:});
  78   elseif (ischar(m))
  79     [varargout{1:nargout_with_ans}] = handle_gradient (str2func (m), varargin{:});
  80   else
  81     error ("gradient: first input must be an array or a function");
  82   endif
  83
  84 endfunction
  85
  86 function varargout = matrix_gradient (m, varargin)
  87   transposed = false;
  88   if (isvector (m))
  89     ## make a row vector.
  90     transposed = (size (m, 2) == 1);
  91     m = m(:).';
  92   endif
  93
  94   nd = ndims (m);
  95   sz = size (m);
  96   if (length(sz) > 1)
  97     tmp = sz(1); sz(1) = sz(2); sz(2) = tmp;
  98   endif
  99
 100   if (nargin > 2 && nargin != nd + 1)
 101     print_usage ();
 102   endif
 103
 104   ## cell d stores a spacing vector for each dimension
 105   d = cell (1, nd);
 106   if (nargin == 1)
 107     ## no spacing given - assume 1.0 for all dimensions
 108     for i = 1:nd
 109       d{i} = ones (sz(i) - 1, 1);
 110     endfor
 111   elseif (nargin == 2)
 112     if (isscalar (varargin{1}))
 113       ## single scalar value for all dimensions
 114       for i = 1:nd
 115         d{i} = varargin{1} * ones (sz(i) - 1, 1);
 116       endfor
 117     else
 118       ## vector for one-dimensional derivative
 119       d{1} = diff (varargin{1}(:));
 120     endif
 121   else
 122     ## have spacing value for each dimension
 123     if (length(varargin) != nd)
 124       error ("gradient: dimensions and number of spacing values do not match");
 125     endif
 126     for i = 1:nd
 127       if (isscalar (varargin{i}))
 128         d{i} = varargin{i} * ones (sz(i) - 1, 1);
 129       else
 130         d{i} = diff (varargin{i}(:));
 131       endif
 132     endfor
 133   endif
 134
 135   m = shiftdim (m, 1);
 136   for i = 1:min (nd, nargout)
 137     mr = rows (m);
 138     mc = numel (m) / mr;
 139     Y = zeros (size (m), class (m));
 140
 141     if (mr > 1)
 142       ## Top and bottom boundary.
 143       Y(1,:) = diff (m(1:2, :)) / d{i}(1);
 144       Y(mr,:) = diff (m(mr-1:mr, :) / d{i}(mr - 1));
 145     endif
 146
 147     if (mr > 2)
 148       ## Interior points.
 149       Y(2:mr-1,:) = ((m(3:mr,:) - m(1:mr-2,:))
 150           ./ kron (d{i}(1:mr-2) + d{i}(2:mr-1), ones (1, mc)));
 151     endif
 152
 153     ## turn multi-dimensional matrix in a way, that gradient
 154     ## along x-direction is calculated first then y, z, ...
 155
 156     if (i == 1)
 157       varargout{i} = shiftdim (Y, nd - 1);
 158       m = shiftdim (m, nd - 1);
 159     elseif (i == 2)
 160       varargout{i} = Y;
 161       m = shiftdim (m, 2);
 162     else
 163       varargout{i} = shiftdim (Y, nd - i + 1);
 164       m = shiftdim (m, 1);
 165     endif
 166   endfor
 167
 168   if (transposed)
 169     varargout{1} = varargout{1}.';
 170   endif
 171 endfunction
 172
 173 function varargout = handle_gradient (f, p0, varargin)
 174   ## Input checking
 175   p0_size = size (p0);
 176
 177   if (numel (p0_size) != 2)
 178     error ("gradient: the second input argument should either be a vector or a matrix");
 179   endif
 180
 181   if (any (p0_size == 1))
 182     p0 = p0 (:);
 183     dim = 1;
 184     num_points = numel (p0);
 185   else
 186     num_points = p0_size (1);
 187     dim = p0_size (2);
 188   endif
 189
 190   if (length (varargin) == 0)
 191     delta = 1;
 192   elseif (length (varargin) == 1 || length (varargin) == dim)
 193     try
 194       delta = [varargin{:}];
 195     catch
 196       error ("gradient: spacing parameters must be scalars or a vector");
 197     end_try_catch
 198   else
 199     error ("gradient: incorrect number of spacing parameters");
 200   endif
 201
 202   if (isscalar (delta))
 203     delta = repmat (delta, 1, dim);
 204   elseif (!isvector (delta))
 205     error ("gradient: spacing values must be scalars or a vector");
 206   endif
 207
 208   ## Calculate the gradient
 209   p0 = mat2cell (p0, num_points, ones (1, dim));
 210   varargout = cell (1, dim);
 211   for d = 1:dim
 212     s = delta (d);
 213     df_dx = (f (p0{1:d-1}, p0{d}+s, p0{d+1:end})
 214            - f (p0{1:d-1}, p0{d}-s, p0{d+1:end})) ./ (2*s);
 215     if (dim == 1)
 216       varargout{d} = reshape (df_dx, p0_size);
 217     else
 218       varargout{d} = df_dx;
 219     endif
 220   endfor
 221 endfunction
 222
 223 %!test
 224 %! data = [1, 2, 4, 2];
 225 %! dx = gradient (data);
 226 %! dx2 = gradient (data, 0.25);
 227 %! dx3 = gradient (data, [0.25, 0.5, 1, 3]);
 228 %! assert (dx, [1, 3/2, 0, -2]);
 229 %! assert (dx2, [4, 6, 0, -8]);
 230 %! assert (dx3, [4, 4, 0, -1]);
 231 %! assert (size_equal(data, dx));
 232
 233 %!test
 234 %! [Y,X,Z,U] = ndgrid (2:2:8,1:5,4:4:12,3:5:30);
 235 %! [dX,dY,dZ,dU] = gradient (X);
 236 %! assert (all(dX(:)==1));
 237 %! assert (all(dY(:)==0));
 238 %! assert (all(dZ(:)==0));
 239 %! assert (all(dU(:)==0));
 240 %! [dX,dY,dZ,dU] = gradient (Y);
 241 %! assert (all(dX(:)==0));
 242 %! assert (all(dY(:)==2));
 243 %! assert (all(dZ(:)==0));
 244 %! assert (all(dU(:)==0));
 245 %! [dX,dY,dZ,dU] = gradient (Z);
 246 %! assert (all(dX(:)==0));
 247 %! assert (all(dY(:)==0));
 248 %! assert (all(dZ(:)==4));
 249 %! assert (all(dU(:)==0));
 250 %! [dX,dY,dZ,dU] = gradient (U);
 251 %! assert (all(dX(:)==0));
 252 %! assert (all(dY(:)==0));
 253 %! assert (all(dZ(:)==0));
 254 %! assert (all(dU(:)==5));
 255 %! assert (size_equal(dX, dY, dZ, dU, X, Y, Z, U));
 256 %! [dX,dY,dZ,dU] = gradient (U, 5.0);
 257 %! assert (all(dU(:)==1));
 258 %! [dX,dY,dZ,dU] = gradient (U, 1.0, 2.0, 3.0, 2.5);
 259 %! assert (all(dU(:)==2));
 260
 261 %!test
 262 %! [Y,X,Z,U] = ndgrid (2:2:8,1:5,4:4:12,3:5:30);
 263 %! [dX,dY,dZ,dU] = gradient (X+j*X);
 264 %! assert (all(dX(:)==1+1j));
 265 %! assert (all(dY(:)==0));
 266 %! assert (all(dZ(:)==0));
 267 %! assert (all(dU(:)==0));
 268 %! [dX,dY,dZ,dU] = gradient (Y-j*Y);
 269 %! assert (all(dX(:)==0));
 270 %! assert (all(dY(:)==2-j*2));
 271 %! assert (all(dZ(:)==0));
 272 %! assert (all(dU(:)==0));
 273 %! [dX,dY,dZ,dU] = gradient (Z+j*1);
 274 %! assert (all(dX(:)==0));
 275 %! assert (all(dY(:)==0));
 276 %! assert (all(dZ(:)==4));
 277 %! assert (all(dU(:)==0));
 278 %! [dX,dY,dZ,dU] = gradient (U-j*1);
 279 %! assert (all(dX(:)==0));
 280 %! assert (all(dY(:)==0));
 281 %! assert (all(dZ(:)==0));
 282 %! assert (all(dU(:)==5));
 283 %! assert (size_equal(dX, dY, dZ, dU, X, Y, Z, U));
 284 %! [dX,dY,dZ,dU] = gradient (U, 5.0);
 285 %! assert (all(dU(:)==1));
 286 %! [dX,dY,dZ,dU] = gradient (U, 1.0, 2.0, 3.0, 2.5);
 287 %! assert (all(dU(:)==2));
 288
 289 %!test
 290 %! x = 0:10;
 291 %! f = @cos;
 292 %! df_dx = @(x) -sin (x);
 293 %! assert (gradient (f, x), df_dx (x), 0.2);
 294 %! assert (gradient (f, x, 0.5), df_dx (x), 0.1);
 295
 296 %!test
 297 %! xy = reshape (1:10, 5, 2);
 298 %! f = @(x,y) sin (x) .* cos (y);
 299 %! df_dx = @(x, y) cos (x) .* cos (y);
 300 %! df_dy = @(x, y) -sin (x) .* sin (y);
 301 %! [dx, dy] = gradient (f, xy);
 302 %! assert (dx, df_dx (xy (:, 1), xy (:, 2)), 0.1)
 303 %! assert (dy, df_dy (xy (:, 1), xy (:, 2)), 0.1)
 304