Patch to add quiver3 and surfnorm

[David Bateman]

Here is a patch I worked on over the weekend that adds the quiver3
function. Matlab has a strange arrow head in the xy plane for the arrows
of quiver3 and so if you are looking at the plot in the xz or yz plane
then the arrowhead tend to disappear. This is the same as in matlab
itself and so I don't see that we can do much..

As for the surfnorm function, I noticed this existed as the test
function for quiver3 on the matlab website and so implemented it to test
quiver3.. However, for me it makes more sense to have to normal vectors
at the center of each quadrilateral patch defined by the meshgrid than
at the vertices. The reason is that the definition of the normal vector
is very clear, and its something like this normal vector that I've
always used in the past for integral equation formulations of
electromagnetic fields. This is NOT what surnorm does it appears. What
surfnorm does is estimate the normal vectors at the vertices from the
normal vectors at the centers of the quadrilaterals defined by the
meshgrid.

There is an issue with this estimation in that the vertices at the edge
of the meshgrid must be extrapolated. This means that closed surfaces,
like those returned by the "sphere" function, will have funny
extrapolation effects. Matlab seems not to suffer from this problem too
much and so I suspect that use a cubic extrapolation, whereas at the
moment I've limited myself to a linear extrapolation. This is not a real
issue to me however. Does anyone else care?

D.

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

*** ./scripts/plot/Makefile.in.orig25   2007-11-26 12:08:44.624230974 +0100
--- ./scripts/plot/Makefile.in  2007-11-26 12:09:03.728254999 +0100
***************
*** 113,118 ****
--- 113,119 ----
    polar.m \
    print.m \
    quiver.m \
+   quiver3.m \
    replot.m \
    ribbon.m \
    scatter.m \
***************
*** 132,137 ****
--- 133,139 ----
    surf.m \
    surface.m \
    surfc.m \
+   surfnorm.m \
    text.m \
    title.m \
    view.m \
*** ./scripts/plot/quiver3.m.orig25     2007-11-26 12:08:31.071923140 +0100
--- ./scripts/plot/quiver3.m    2007-11-26 12:20:34.704115468 +0100
***************
*** 0 ****
--- 1,97 ----
+ ## 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 3 of the License, 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, see
+ ## <http://www.gnu.org/licenses/>.
+ 
+ ## -*- texinfo -*-
+ ## @deftypefn {Function File} {} quiver3 (@var{u}, @var{v}, @var{w})
+ ## @deftypefnx {Function File} {} quiver3 (@var{x}, @var{y}, @var{z}, 
@var{u}, @var{v}, @var{w})
+ ## @deftypefnx {Function File} {} quiver3 (@dots{}, @var{s})
+ ## @deftypefnx {Function File} {} quiver3 (@dots{}, @var{style})
+ ## @deftypefnx {Function File} {} quiver3 (@dots{}, 'filled')
+ ## @deftypefnx {Function File} {} quiver3 (@var{h}, @dots{})
+ ## @deftypefnx {Function File} address@hidden =} quiver3 (@dots{})
+ ##
+ ## Plot the @code{(@var{u}, @var{v}, @var{w})} components of a vector field 
in 
+ ## an @code{(@var{x}, @var{y}), @var{z}} meshgrid. If the grid is uniform, 
you 
+ ## can specify @var{x}, @var{y} @var{z} as vectors.
+ ##
+ ## If @var{x}, @var{y} and @var{z} are undefined they are assumed to be
+ ## @code{(1:@var{m}, 1:@var{n}, 1:@var{p})} where @address@hidden, @var{n}] = 
+ ## size(@var{u})} and @address@hidden = max (size (@var{w}))}.
+ ##
+ ## The variable @var{s} is a scalar defining a scaling factor to use for
+ ##  the arrows of the field relative to the mesh spacing. A value of 0 
+ ## disables all scaling. The default value is 1.
+ ##
+ ## The style to use for the plot can be defined with a line style @var{style}
+ ## in a similar manner to the line styles used with the @code{plot} command.
+ ## If a marker is specified then markers at the grid points of the vectors are
+ ## printed rather than arrows. If the argument 'filled' is given then the
+ ## markers as filled.
+ ##
+ ## The optional return value @var{h} provides a list of handles to the 
+ ## the parts of the vector field (body, arrow and marker).
+ ##
+ ## @example
+ ## @group
+ ## [x, y, z] = peaks (25);
+ ## surf (x, y, z);
+ ## hold on;
+ ## [u, v, w] = surfnorm (x, y, z / 10);
+ ## quiver3 (x, y, z, u, v, w);
+ ## @end group
+ ## @end example
+ ##
+ ## @seealso{plot}
+ ## @end deftypefn
+ 
+ function retval = quiver3 (varargin)
+ 
+   if (nargin < 2)
+     print_usage ();
+   elseif (isscalar (varargin{1}) && ishandle (varargin{1}))
+     h = varargin {1};
+     if (! strcmp (get (h, "type"), "axes"))
+       error ("quiver: expecting first argument to be an axes object");
+     endif
+     oldh = gca ();
+     unwind_protect
+       axes (h);
+       newplot ();
+       tmp = __quiver__ (h, varargin{2:end});
+     unwind_protect_cleanup
+       axes (oldh);
+     end_unwind_protect
+   else
+     newplot ();
+     tmp = __quiver__ (gca (), 1, varargin{:});
+   endif
+ 
+   if (nargout > 0)
+     retval = tmp;
+   endif
+ 
+ endfunction
+ 
+ %!demo
+ %! [x,y]=meshgrid (-1:0.1:1); 
+ %! z=sin(2*pi*sqrt(x.^2+y.^2)); 
+ %! theta=2*pi*sqrt(x.^2+y.^2)+pi/2;
+ %! quiver3(x,y,z,sin(theta),cos(theta),ones(size(z)));
+ %! hold on; 
+ %! mesh(x,y,z); 
+ %! hold off;
*** ./scripts/plot/surfnorm.m.orig25    2007-11-26 12:09:23.765231036 +0100
--- ./scripts/plot/surfnorm.m   2007-11-26 12:19:44.998634645 +0100
***************
*** 0 ****
--- 1,155 ----
+ ## 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 3 of the License, 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, see
+ ## <http://www.gnu.org/licenses/>.
+ 
+ ## -*- texinfo -*-
+ ## @deftypefn {Function File} {} surfnorm (@var{x}, @var{y}, @var{z})
+ ## @deftypefnx {Function File} {} surfnorm (@var{z})
+ ## @deftypefnx {Function File} address@hidden, @var{ny}, @var{nz}] =} 
surfnorm (@dots{})
+ ## @deftypefnx {Function File} {} surfnorm (@var{h}, @dots{})
+ ## Find the vectors normal to a meshgridded surface. The meshed gridded 
+ ## surface is defined by @var{x}, @var{y}, and @var{z}. If @var{x} and 
+ ## @var{y} are not defined, then it is assumed that they are given by
+ ##
+ ## @example
+ ## address@hidden, @var{y}] = meshgrid (1:size(@var{z}, 1), 
+ ##                      1:size(@var{z}, 2));
+ ## @end example
+ ##
+ ## If no return arguments are requested, a surface plot with the normal 
+ ## vectors to the surface is plotted. Otherwise the componets of the normal
+ ## vectors at the mesh gridded points are returned in @var{nx}, @var{ny},
+ ## and @var{nz}.
+ ##
+ ## The normal vectors are calculated by taking the cross product of the 
+ ## diagonals of eash of teh quadrilaterals in the meshgrid to find the 
+ ## normal vectors of the centers of these quadrilaterals. The four nearest
+ ## normal vectors to the meshgrid points are then averaged to obtain the 
+ ## normal to the surface at the meshgridded points.
+ ##
+ ## An example of the use of @code{surfnorm} is
+ ##
+ ## @example
+ ## surfnorm (peaks (25));
+ ## @end example
+ ## @seealso{surf, quiver3}
+ ## @end deftypefn
+ 
+ function varargout = surfnorm (varargin)
+ 
+   if (nargout > 0)
+     varargout = cell (nargout, 1);
+   else
+     varargout = cell (0, 0);
+   endif
+   if (isscalar (varargin{1}) && ishandle (varargin{1}))
+     h = varargin {1};
+     if (! strcmp (get (h, "type"), "axes"))
+       error ("surfnorm: expecting first argument to be an axes object");
+     endif
+     if (nargin != 2 && nargin != 4)
+       print_usage ();
+     endif
+     oldh = gca ();
+     unwind_protect
+       axes (h);
+       [varargout{:}] = __surfnorm__ (h, varargin{2:end});
+     unwind_protect_cleanup
+       axes (oldh);
+     end_unwind_protect
+   else
+     if (nargin != 1 && nargin != 3)
+       print_usage ();
+     endif
+     [varargout{:}] = __surfnorm__ (gca (), varargin{:});
+   endif
+ 
+ endfunction
+ 
+ function [Nx, Ny, Nz] = __surfnorm__ (h, varargin)
+ 
+   if (nargin == 2)
+     z = varargin{1};
+     [x, y] = meshgrid (1:size(z,1), 1:size(z,2));
+     ioff = 2;
+   else
+     x = varargin{1};
+     y = varargin{2};
+     z = varargin{3};
+     ioff = 4;
+   endif
+ 
+   if (nargout == 0)
+     newplot();
+     surf (x, y, z, varargin{ioff:end});
+     hold on;
+   endif
+ 
+   ## Make life easier, and avoid having to do the extrapolation later, do
+   ## a simpler linear extrapolation here. This is approximative, and works
+   ## badly for closed surfaces like spheres.
+   xx = [2 .* x(:,1) - x(:,2), x, 2 .* x(:,end) - x(:,end-1)];
+   xx = [2 .* xx(1,:) - xx(2,:); xx; 2 .* xx(end,:) - xx(end-1,:)];
+   yy = [2 .* y(:,1) - y(:,2), y, 2 .* y(:,end) - y(:,end-1)];
+   yy = [2 .* yy(1,:) - yy(2,:); yy; 2 .* yy(end,:) - yy(end-1,:)];
+   zz = [2 .* z(:,1) - z(:,2), z, 2 .* z(:,end) - z(:,end-1)];
+   zz = [2 .* zz(1,:) - zz(2,:); zz; 2 .* zz(end,:) - zz(end-1,:)];
+ 
+   u.x = xx(1:end-1,1:end-1) - xx(2:end,2:end);
+   u.y = yy(1:end-1,1:end-1) - yy(2:end,2:end);
+   u.z = zz(1:end-1,1:end-1) - zz(2:end,2:end);
+   v.x = xx(1:end-1,2:end) - xx(2:end,1:end-1);
+   v.y = yy(1:end-1,2:end) - yy(2:end,1:end-1);
+   v.z = zz(1:end-1,2:end) - zz(2:end,1:end-1);
+ 
+   c = cross ([u.x(:), u.y(:), u.z(:)], [v.x(:), v.y(:), v.z(:)]);
+   w.x = reshape (c(:,1), size(u.x));
+   w.y = reshape (c(:,2), size(u.y));
+   w.z = reshape (c(:,3), size(u.z));
+ 
+   ## Create normal vectors as mesh vectices from normals at mesh centers
+   nx = (w.x(1:end-1,1:end-1) + w.x(1:end-1,2:end) +
+       w.x(2:end,1:end-1) + w.x(2:end,2:end)) ./ 4; 
+   ny = (w.y(1:end-1,1:end-1) + w.y(1:end-1,2:end) +
+       w.y(2:end,1:end-1) + w.y(2:end,2:end)) ./ 4; 
+   nz = (w.z(1:end-1,1:end-1) + w.z(1:end-1,2:end) +
+       w.z(2:end,1:end-1) + w.z(2:end,2:end)) ./ 4; 
+ 
+   ## Normalize the normal vectors
+   len = sqrt (nx.^2 + ny.^2 + nz.^2);
+   nx = nx ./ len;
+   ny = ny ./ len;
+   nz = nz ./ len;
+ 
+   if (nargout == 0)
+     plot3 ([x(:)'; x(:).' + nx(:).' ; NaN(size(x(:).'))](:),
+          [y(:)'; y(:).' + ny(:).' ; NaN(size(y(:).'))](:),
+          [z(:)'; z(:).' + nz(:).' ; NaN(size(z(:).'))](:), 
+          varargin{ioff:end});
+   else
+     Nx = nx;
+     Ny = ny;
+     Nz = nz;
+   endif
+ endfunction
+ 
+ %!demo
+ %! [x, y, z] = peaks(10);
+ %! surfnorm (x, y, z);
+ 
+ %!demo
+ %! surfnorm (peaks(10));

2007-11-23  David Bateman  <address@hidden>

        * plot/quiver3.m, plot/surfnorm.m: New functions.
        * plot/Makefile.in (SOURCES): Add them to the sources.