Ported Octave-Forge Function (Was: Irregularly gridded Discrete Laplacian Operator)

[David Bateman]

Find attached a patch for the current CVS of Octave that ports the del2
function from octave-forge. It made a number of changes to make it
matlab compatible, including

* Support of NDArrays
* Use linear extrapolation for the boundaries rather than a nearest
neighbor approximation
* Change the treatment of irregular gridding to give the same result as
matlab

The patch also adds del2 to the arith.txi file to document its existence..

Regards
David

-- 
David Bateman                                address@hidden
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob) 
91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax) 

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary

*** ./doc/interpreter/arith.txi.orig47  2007-07-19 14:17:12.770858865 +0200
--- ./doc/interpreter/arith.txi 2007-07-19 14:18:26.246127311 +0200
***************
*** 30,35 ****
--- 30,37 ----
  
  @DOCSTRING(cplxpair)
  
+ @DOCSTRING(del2)
+ 
  @DOCSTRING(exp)
  
  @DOCSTRING(factor)
*** ./scripts/general/del2.m.orig47     2007-07-19 14:15:57.623673128 +0200
--- ./scripts/general/del2.m    2007-07-19 14:18:38.765491291 +0200
***************
*** 0 ****
--- 1,156 ----
+ ## Copyright (C) 2000  Kai Habel
+ ## Copyright (C) 2007  David Bateman
+ ##
+ ## This file is part of Octave.
+ ##
+ ## Octave is free software; you can redistribute it and/or modify it
+ ## under the terms of the GNU General Public License as published by
+ ## the Free Software Foundation; either version 2, or (at your option)
+ ## any later version.
+ ##
+ ## Octave is distributed in the hope that it will be useful, but
+ ## WITHOUT ANY WARRANTY; without even the implied warranty of
+ ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ ## General Public License for more details.
+ ##
+ ## You should have received a copy of the GNU General Public License
+ ## along with Octave; see the file COPYING.  If not, write to the Free
+ ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+ ## 02110-1301, USA.
+ 
+ ## -*- texinfo -*-
+ ## @deftypefn {Function File} address@hidden =} del2 (@var{m})
+ ## @deftypefnx {Function File} address@hidden =} del2 (@var{m}, @var{h})
+ ## @deftypefnx {Function File} address@hidden =} del2 (@var{m}, @var{dx}, 
@var{dy}, @dots{})
+ ##
+ ## Calculates the discrete Laplace operator. If @var{m} is a matrix this is
+ ## defined as
+ ##
+ ## @iftex
+ ## @tex
+ ## $d = {1 \over 4} \left( {d^2 \over dx^2} M(x,y) + {d^2 \over dy^2} M(x,y) 
\right)$
+ ## @end tex
+ ## @end iftex
+ ## @ifnottex
+ ## @example
+ ## @group
+ ##       1    / d^2            d^2         \
+ ## D  = --- * | ---  M(x,y) +  ---  M(x,y) | 
+ ##       4    \ dx^2           dy^2        /
+ ## @end group
+ ## @end example
+ ## @end ifnottex
+ ##
+ ## The above to continued to N-dimensional arrays calculating the second
+ ## derivative over the higher dimensions.
+ ##
+ ## The spacing between evaluation points may be defined by @var{h}, which is a
+ ## scalar defining the spacing in all dimensions. Or alternative, the spacing
+ ## in each dimension may be defined separately by @var{dx}, @var{dy}, etc. 
+ ## Scalar spacing values give equidistant spacing, whereas vector spacing 
+ ## values can be used to specify variable spacing. The length of the vectors
+ ## must match the respective dimension of @var{m}. The default spacing value
+ ## is 1.
+ ##
+ ## You need at least 3 data points for each dimension. Boundary points are
+ ## calculated as the linear extrapolation of the interior points.
+ ##
+ ## @seealso{gradient, diff}
+ ## @end deftypefn
+ 
+ ## Author:  Kai Habel <address@hidden>
+ 
+ function D = del2 (M, varargin)
+   
+   if (nargin < 1)
+     print_usage ();
+   endif
+ 
+   nd = ndims (M);
+   sz = size (M);
+   dx = cell (1, nd);
+   if (nargin == 2 || nargin == 1)
+     if (nargin == 1)
+       h = 1;
+     else
+       h = varargin{1}
+     endif
+     for i = 1 : nd
+       if (isscalar (h))
+       dx{i} = h * ones (sz (i), 1);
+       else
+       if (length (h) == sz (i))
+         dx{i} = diff (h)(:);
+       else
+         error ("dimensionality mismatch in %d-th spacing vector", i);
+       endif
+       endif
+     endfor
+   elseif (nargin - 1 == nd)
+     ## Reverse dx{1} and dx{2} as the X-dim is the 2nd dim of the ND array
+     tmp = varargin{1};
+     varargin{1} = varargin{2};
+     varargin{2} = tmp;
+ 
+     for i = 1 : nd
+       if (isscalar (varargin{i}))
+       dx{i} = varargin{i} * ones (sz (i), 1);
+       else
+       if (length (varargin{i}) == sz (i))
+         dx{i} = diff (varargin{i})(:);
+       else
+         error ("dimensionality mismatch in %d-th spacing vector", i);
+       endif
+       endif
+     endfor
+   else
+     print_usage ();
+   endif
+ 
+   idx = cell (1, nd);
+   for i = 1: nd
+     idx{i} = ":";
+   endfor
+ 
+   D = zeros (sz);
+   for i = 1: nd
+     if (sz(i) >= 3)
+       DD = zeros (sz);
+       idx1 = idx2 = idx3 = idx;
+ 
+       ## interior points
+       idx1{i} = 1 : sz(i) - 2;
+       idx2{i} = 2 : sz(i) - 1;
+       idx3{i} = 3 : sz(i);
+       szi = sz;
+       szi (i) = 1;
+ 
+       h1 = repmat (shiftdim (dx{i}(1 : sz(i) - 2), 1 - i), szi);
+       h2 = repmat (shiftdim (dx{i}(2 : sz(i) - 1), 1 - i), szi);
+       DD(idx2{:}) = ((M(idx1{:}) - M(idx2{:})) ./ h1 + ...
+                    (M(idx3{:}) - M(idx2{:})) ./ h2) ./ (h1 + h2);
+ 
+       ## left and right boundary
+       if (sz(i) == 3)
+       DD(idx1{:}) = DD(idx3{:}) = DD(idx2{:});
+       else
+       idx1{i} = 1;
+       idx2{i} = 2;
+       idx3{i} = 3;
+       DD(idx1{:}) = (dx{i}(1) + dx{i}(2)) / dx{i}(2) * DD (idx2{:}) - ...
+           dx{i}(1) / dx{i}(2) * DD (idx3{:});
+ 
+       idx1{i} = sz(i);
+       idx2{i} = sz(i) - 1;
+       idx3{i} = sz(i) - 2;
+       DD(idx1{:}) =  (dx{i}(sz(i) - 1) + dx{i}(sz(i) - 2)) / ...
+           dx{i}(sz(i) - 2) * DD (idx2{:}) - ...
+           dx{i}(sz(i) - 1) / dx{i}(sz(i) - 2) * DD (idx3{:});
+       endif
+ 
+       D += DD;
+     endif
+   endfor
+ 
+   D = D ./ nd;
+ endfunction
*** ./scripts/general/Makefile.in.orig47        2007-07-19 14:17:44.681238505 
+0200
--- ./scripts/general/Makefile.in       2007-07-19 14:17:51.190907905 +0200
***************
*** 22,29 ****
  
  SOURCES = __isequal__.m __splinen__.m accumarray.m arrayfun.m bicubic.m \
    bitcmp.m bitget.m bitset.m blkdiag.m cart2pol.m cart2sph.m cell2mat.m \
!   circshift.m common_size.m cplxpair.m cumtrapz.m deal.m diff.m flipdim.m \
!   fliplr.m flipud.m gradient.m ind2sub.m int2str.m interp1.m \
    interp2.m interp3.m interpn.m interpft.m is_duplicate_entry.m isa.m \
    isdefinite.m isdir.m isequal.m isequalwithequalnans.m isscalar.m \
    issquare.m issymmetric.m isvector.m logical.m logspace.m lookup.m mod.m \
--- 22,29 ----
  
  SOURCES = __isequal__.m __splinen__.m accumarray.m arrayfun.m bicubic.m \
    bitcmp.m bitget.m bitset.m blkdiag.m cart2pol.m cart2sph.m cell2mat.m \
!   circshift.m common_size.m cplxpair.m cumtrapz.m deal.m del2.m diff.m \
!   flipdim.m fliplr.m flipud.m gradient.m ind2sub.m int2str.m interp1.m \
    interp2.m interp3.m interpn.m interpft.m is_duplicate_entry.m isa.m \
    isdefinite.m isdir.m isequal.m isequalwithequalnans.m isscalar.m \
    issquare.m issymmetric.m isvector.m logical.m logspace.m lookup.m mod.m \

2007-07-19  David Bateman  <address@hidden>

        * general/del2.m: New function for discrete laplacian operator.
        * Makefile.in: Include del2.m in SOURCES.

2007-07-19  David Bateman  <address@hidden>

        * interpreter/arith.txi: Document del2.