octave_packages/m/plot/slice.m

   1 ## Copyright (C) 2007-2012 Kai Habel, David Bateman
   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} {} slice (@var{x}, @var{y}, @var{z}, @var{v}, @var{sx}, @var{sy}, @var{sz})
  21 ## @deftypefnx {Function File} {} slice (@var{x}, @var{y}, @var{z}, @var{v}, @var{xi}, @var{yi}, @var{zi})
  22 ## @deftypefnx {Function File} {} slice (@var{v}, @var{sx}, @var{sy}, @var{sz})
  23 ## @deftypefnx {Function File} {} slice (@var{v}, @var{xi}, @var{yi}, @var{zi})
  24 ## @deftypefnx {Function File} {@var{h} =} slice (@dots{})
  25 ## @deftypefnx {Function File} {@var{h} =} slice (@dots{}, @var{method})
  26 ## Plot slices of 3-D data/scalar fields.  Each element of the 3-dimensional
  27 ## array @var{v} represents a scalar value at a location given by the
  28 ## parameters @var{x}, @var{y}, and @var{z}.  The parameters @var{x},
  29 ## @var{x}, and @var{z} are either 3-dimensional arrays of the same size
  30 ## as the array @var{v} in the "meshgrid" format or vectors.  The
  31 ## parameters @var{xi}, etc. respect a similar format to @var{x}, etc.,
  32 ## and they represent the points at which the array @var{vi} is
  33 ## interpolated using interp3.  The vectors @var{sx}, @var{sy}, and
  34 ## @var{sz} contain points of orthogonal slices of the respective axes.
  35 ##
  36 ## If @var{x}, @var{y}, @var{z} are omitted, they are assumed to be
  37 ## @code{x = 1:size (@var{v}, 2)}, @code{y = 1:size (@var{v}, 1)} and
  38 ## @code{z = 1:size (@var{v}, 3)}.
  39 ##
  40 ## @var{Method} is one of:
  41 ##
  42 ## @table @asis
  43 ## @item "nearest"
  44 ## Return the nearest neighbor.
  45 ##
  46 ## @item "linear"
  47 ## Linear interpolation from nearest neighbors.
  48 ##
  49 ## @item "cubic"
  50 ## Cubic interpolation from four nearest neighbors (not implemented yet).
  51 ##
  52 ## @item "spline"
  53 ## Cubic spline interpolation---smooth first and second derivatives
  54 ## throughout the curve.
  55 ## @end table
  56 ##
  57 ## The default method is @code{"linear"}.
  58 ##
  59 ## The optional return value @var{h} is a graphics handle to the created
  60 ## surface object.
  61 ##
  62 ## Examples:
  63 ##
  64 ## @example
  65 ## @group
  66 ## [x, y, z] = meshgrid (linspace (-8, 8, 32));
  67 ## v = sin (sqrt (x.^2 + y.^2 + z.^2)) ./ (sqrt (x.^2 + y.^2 + z.^2));
  68 ## slice (x, y, z, v, [], 0, []);
  69 ## [xi, yi] = meshgrid (linspace (-7, 7));
  70 ## zi = xi + yi;
  71 ## slice (x, y, z, v, xi, yi, zi);
  72 ## @end group
  73 ## @end example
  74 ## @seealso{interp3, surface, pcolor}
  75 ## @end deftypefn
  76
  77 ## Author: Kai Habel <kai.habel@gmx.de>
  78
  79 function h = slice (varargin)
  80
  81   method = "linear";
  82   nargs = nargin;
  83
  84   if (ischar (varargin{end}))
  85     method = varargin{end};
  86     nargs -= 1;
  87   endif
  88
  89   if (nargs == 4)
  90     v = varargin{1};
  91     if (ndims (v) != 3)
  92       error ("slice: expect 3-dimensional array of values");
  93     endif
  94     [nx, ny, nz] = size (v);
  95     [x, y, z] = meshgrid (1:nx, 1:ny, 1:nz);
  96     sx = varargin{2};
  97     sy = varargin{3};
  98     sz = varargin{4};
  99   elseif (nargs == 7)
 100     v = varargin{4};
 101     if (ndims (v) != 3)
 102       error ("slice: expect 3-dimensional array of values");
 103     endif
 104     x = varargin{1};
 105     y = varargin{2};
 106     z = varargin{3};
 107     if (all ([isvector(x), isvector(y), isvector(z)]))
 108       [x, y, z] = meshgrid (x, y, z);
 109     elseif (ndims (x) == 3 && size_equal (x, y, z))
 110       ## Do nothing.
 111     else
 112       error ("slice: X, Y, Z size mismatch");
 113     endif
 114     sx = varargin{5};
 115     sy = varargin{6};
 116     sz = varargin{7};
 117   else
 118     print_usage ();
 119   endif
 120
 121   if (any ([isvector(sx), isvector(sy), isvector(sz)]))
 122     have_sval = true;
 123   elseif (ndims(sx) == 2 && size_equal (sx, sy, sz))
 124     have_sval = false;
 125   else
 126     error ("slice: dimensional mismatch for (XI, YI, ZI) or (SX, SY, SZ)");
 127   endif
 128
 129   newplot ();
 130   ax = gca ();
 131   sidx = 1;
 132   maxv = max (v(:));
 133   minv = min (v(:));
 134   set (ax, "clim", [minv, maxv]);
 135
 136   if (have_sval)
 137     ns = length (sx) + length (sy) + length (sz);
 138     hs = zeros(ns,1);
 139     [ny, nx, nz] = size (v);
 140     if (length(sz) > 0)
 141       for i = 1:length(sz)
 142         [xi, yi, zi] = meshgrid (squeeze (x(1,:,1)),
 143                                  squeeze (y(:,1,1)), sz(i));
 144         vz = squeeze (interp3 (x, y, z, v, xi, yi, zi, method));
 145         tmp(sidx++) = surface (xi, yi, sz(i) * ones (size (yi)), vz);
 146       endfor
 147     endif
 148
 149     if (length (sy) > 0)
 150       for i = length(sy):-1:1
 151         [xi, yi, zi] = meshgrid (squeeze (x(1,:,1)), sy(i), squeeze (z(1,1,:)));
 152         vy = squeeze (interp3 (x, y, z, v, xi, yi, zi, method));
 153         tmp(sidx++) = surface (squeeze (xi),
 154                                squeeze (sy(i) * ones (size (zi))),
 155                                squeeze (zi), vy);
 156       endfor
 157     endif
 158
 159     if (length (sx) > 0)
 160       for i = length(sx):-1:1
 161         [xi, yi, zi] = meshgrid (sx(i), squeeze (y(:,1,1)), squeeze (z(1,1,:)));
 162         vx = squeeze (interp3 (x, y, z, v, xi, yi, zi, method));
 163         tmp(sidx++) = surface (squeeze (sx(i) * ones (size (zi))),
 164                                squeeze (yi), squeeze(zi), vx);
 165       endfor
 166     endif
 167   else
 168     vi = interp3 (x, y, z, v, sx, sy, sz);
 169     tmp = surface (sx, sy, sz, vi);
 170   endif
 171
 172   if (! ishold ())
 173     set (ax, "view", [-37.5, 30.0], "box", "off", "xgrid", "on",
 174          "ygrid", "on", "zgrid", "on");
 175   endif
 176
 177   if (nargout > 0)
 178     h = tmp;
 179   endif
 180
 181 endfunction
 182
 183
 184 %!demo
 185 %! clf
 186 %! [x, y, z] = meshgrid (linspace (-8, 8, 32));
 187 %! v = sin (sqrt (x.^2 + y.^2 + z.^2)) ./ (sqrt (x.^2 + y.^2 + z.^2));
 188 %! slice (x, y, z, v, [], 0, []);
 189
 190 %!demo
 191 %! clf
 192 %! [x, y, z] = meshgrid (linspace (-8, 8, 32));
 193 %! v = sin (sqrt (x.^2 + y.^2 + z.^2)) ./ (sqrt (x.^2 + y.^2 + z.^2));
 194 %! [xi, yi] = meshgrid (linspace (-7, 7));
 195 %! zi = xi + yi;
 196 %! slice (x, y, z, v, xi, yi, zi);
 197