octave_packages/m/plot/isosurface.m

   1 ## Copyright (C) 2009-2012 Martin Helm
   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{fv}] =} isosurface (@var{val}, @var{iso})
  21 ## @deftypefnx {Function File} {[@var{fv}] =} isosurface (@var{x}, @var{y}, @var{z}, @var{val}, @var{iso})
  22 ## @deftypefnx {Function File} {[@var{fv}] =} isosurface (@dots{}, "noshare", "verbose")
  23 ## @deftypefnx {Function File} {[@var{fvc}] =} isosurface (@dots{}, @var{col})
  24 ## @deftypefnx {Function File} {[@var{f}, @var{v}] =} isosurface (@var{x}, @var{y}, @var{z}, @var{val}, @var{iso})
  25 ## @deftypefnx {Function File} {[@var{f}, @var{v}, @var{c}] =} isosurface (@var{x}, @var{y}, @var{z}, @var{val}, @var{iso}, @var{col})
  26 ## @deftypefnx {Function File} {} isosurface (@var{x}, @var{y}, @var{z}, @var{val}, @var{iso}, @var{col}, @var{opt})
  27 ##
  28 ## If called with one output argument and the first input argument
  29 ## @var{val} is a three-dimensional array that contains the data of an
  30 ## isosurface geometry and the second input argument @var{iso} keeps the
  31 ## isovalue as a scalar value then return a structure array @var{fv}
  32 ## that contains the fields @var{Faces} and @var{Vertices} at computed
  33 ## points @command{[x, y, z] = meshgrid (1:l, 1:m, 1:n)}.  The output
  34 ## argument @var{fv} can directly be taken as an input argument for the
  35 ## @command{patch} function.
  36 ##
  37 ## If called with further input arguments @var{x}, @var{y} and @var{z}
  38 ## which are three--dimensional arrays with the same size than @var{val}
  39 ## then the volume data is taken at those given points.
  40 ##
  41 ## The string input argument "noshare" is only for compatibility and
  42 ## has no effect.  If given the string input argument
  43 ## "verbose" then print messages to the command line interface about the
  44 ## current progress.
  45 ##
  46 ## If called with the input argument @var{col} which is a
  47 ## three-dimensional array of the same size than @var{val} then take
  48 ## those values for the interpolation of coloring the isosurface
  49 ## geometry.  Add the field @var{FaceVertexCData} to the structure
  50 ## array @var{fv}.
  51 ##
  52 ## If called with two or three output arguments then return the
  53 ## information about the faces @var{f}, vertices @var{v} and color data
  54 ## @var{c} as seperate arrays instead of a single structure array.
  55 ##
  56 ## If called with no output argument then directly process the
  57 ## isosurface geometry with the @command{patch} command.
  58 ##
  59 ## For example,
  60 ##
  61 ## @example
  62 ## @group
  63 ## [x, y, z] = meshgrid (1:5, 1:5, 1:5);
  64 ## val = rand (5, 5, 5);
  65 ## isosurface (x, y, z, val, .5);
  66 ## @end group
  67 ## @end example
  68 ##
  69 ## @noindent
  70 ## will directly draw a random isosurface geometry in a graphics window.
  71 ## Another example for an isosurface geometry with different additional
  72 ## coloring
  73 ## @c Set example in small font to prevent overfull line
  74 ##
  75 ## @smallexample
  76 ## N = 15;    # Increase number of vertices in each direction
  77 ## iso = .4;  # Change isovalue to .1 to display a sphere
  78 ## lin = linspace (0, 2, N);
  79 ## [x, y, z] = meshgrid (lin, lin, lin);
  80 ## c = abs ((x-.5).^2 + (y-.5).^2 + (z-.5).^2);
  81 ## figure (); # Open another figure window
  82 ##
  83 ## subplot (2,2,1); view (-38, 20);
  84 ## [f, v] = isosurface (x, y, z, c, iso);
  85 ## p = patch ("Faces", f, "Vertices", v, "EdgeColor", "none");
  86 ## set (gca, "PlotBoxAspectRatioMode", "manual", ...
  87 ##           "PlotBoxAspectRatio", [1 1 1]);
  88 ## # set (p, "FaceColor", "green", "FaceLighting", "phong");
  89 ## # light ("Position", [1 1 5]); # Available with the JHandles package
  90 ##
  91 ## subplot (2,2,2); view (-38, 20);
  92 ## p = patch ("Faces", f, "Vertices", v, "EdgeColor", "blue");
  93 ## set (gca, "PlotBoxAspectRatioMode", "manual", ...
  94 ##           "PlotBoxAspectRatio", [1 1 1]);
  95 ## # set (p, "FaceColor", "none", "FaceLighting", "phong");
  96 ## # light ("Position", [1 1 5]);
  97 ##
  98 ## subplot (2,2,3); view (-38, 20);
  99 ## [f, v, c] = isosurface (x, y, z, c, iso, y);
 100 ## p = patch ("Faces", f, "Vertices", v, "FaceVertexCData", c, ...
 101 ##            "FaceColor", "interp", "EdgeColor", "none");
 102 ## set (gca, "PlotBoxAspectRatioMode", "manual", ...
 103 ##           "PlotBoxAspectRatio", [1 1 1]);
 104 ## # set (p, "FaceLighting", "phong");
 105 ## # light ("Position", [1 1 5]);
 106 ##
 107 ## subplot (2,2,4); view (-38, 20);
 108 ## p = patch ("Faces", f, "Vertices", v, "FaceVertexCData", c, ...
 109 ##            "FaceColor", "interp", "EdgeColor", "blue");
 110 ## set (gca, "PlotBoxAspectRatioMode", "manual", ...
 111 ##           "PlotBoxAspectRatio", [1 1 1]);
 112 ## # set (p, "FaceLighting", "phong");
 113 ## # light ("Position", [1 1 5]);
 114 ## @end smallexample
 115 ##
 116 ## @seealso{isonormals, isocolors}
 117 ## @end deftypefn
 118
 119 ## Author: Martin Helm <martin@mhelm.de>
 120
 121 function varargout = isosurface(varargin)
 122
 123   if (nargin < 2 || nargin > 8 || nargout > 3)
 124     print_usage ();
 125   endif
 126
 127   calc_colors = false;
 128   f = v = c = [];
 129   verbose = false;
 130   noshare = false;
 131   if (nargin >= 5)
 132     x = varargin{1};
 133     y = varargin{2};
 134     z = varargin{3};
 135     val = varargin{4};
 136     iso = varargin{5};
 137     if (nargin >= 6 && ismatrix (varargin{6}))
 138       colors = varargin{6};
 139       calc_colors = true;
 140     endif
 141   else
 142     val = varargin{1};
 143     [n2, n1, n3] = size (val);
 144     [x, y, z] = meshgrid (1:n1, 1:n2, 1:n3);
 145     iso = varargin{2};
 146     if (nargin >= 3 && ismatrix (varargin{3}))
 147         colors = varargin{3};
 148         calc_colors = true;
 149     endif
 150   endif
 151   if (calc_colors)
 152     if (nargout == 2)
 153       warning ( "Colors will be calculated, but you did not specify an output argument for it!" );
 154     endif
 155     [fvc.faces, fvc.vertices, fvc.facevertexcdata] = __marching_cube__ (x, y, z, val, iso, colors);
 156   else
 157     [fvc.faces, fvc.vertices] = __marching_cube__ (x, y, z, val, iso);
 158   endif
 159
 160   if (isempty (fvc.vertices) || isempty (fvc.faces))
 161     warning ( "The resulting triangulation is empty" );
 162   endif
 163
 164   switch (nargout)
 165     case 0
 166       ## plot the calculated surface
 167       newplot ();
 168       if (calc_colors)
 169         pa = patch ("Faces", fvc.faces, "Vertices", fvc.vertices,
 170                     "FaceVertexCData", fvc.facevertexcdata,
 171                     "FaceColor", "flat", "EdgeColor", "none");
 172       else
 173         pa = patch ("Faces", fvc.faces, "Vertices", fvc.vertices,
 174                     "FaceColor", "g", "EdgeColor", "k");
 175       endif
 176       if (! ishold ())
 177         set (gca(), "view", [-37.5, 30],
 178              "xgrid", "on", "ygrid", "on", "zgrid", "on");
 179       endif
 180     case 1
 181       varargout = {fvc};
 182     case 2
 183       varargout = {fvc.faces, fvc.vertices};
 184     case 3
 185       varargout = {fvc.faces, fvc.vertices, fvc.facevertexcdata};
 186     otherwise
 187       print_usage ();
 188   endswitch
 189
 190 endfunction
 191
 192
 193 %!shared x, y, z, val
 194 %!  [x, y, z]  = meshgrid (0:1, 0:1, 0:1); ## Points for single
 195 %!  val        = [0, 0; 0, 0];             ## cube and a 3--dim
 196 %!  val(:,:,2) = [0, 0; 1, 0];             ## array of values
 197 %!test
 198 %!  fv = isosurface (x, y, z, val, 0.3);
 199 %!  assert (isfield (fv, "vertices"), true);
 200 %!  assert (isfield (fv, "faces"), true);
 201 %!  assert (size (fv.vertices), [3 3]);
 202 %!  assert (size (fv.faces), [1 3]);
 203 %!test
 204 %!  fvc = isosurface (x, y, z, val, .3, y);
 205 %!  assert (isfield (fvc, "vertices"), true);
 206 %!  assert (isfield (fvc, "faces"), true);
 207 %!  assert (isfield (fvc, "facevertexcdata"), true);
 208 %!  assert (size (fvc.vertices), [3 3]);
 209 %!  assert (size (fvc.faces), [1 3]);
 210 %!  assert (size (fvc.facevertexcdata), [3 1]);
 211 %!test
 212 %!  [f, v] = isosurface (x, y, z, val, .3);
 213 %!  assert (size (f), [1 3]);
 214 %!  assert (size (v), [3 3]);
 215 %!test
 216 %!  [f, v, c] = isosurface (x, y, z, val, .3, y);
 217 %!  assert (size (f), [1 3]);
 218 %!  assert (size (v), [3 3]);
 219 %!  assert (size (c), [3 1]);
 220
 221 %!demo
 222 %! clf
 223 %! [x,y,z] = meshgrid(-2:0.5:2, -2:0.5:2, -2:0.5:2);
 224 %! v = x.^2 + y.^2 + z.^2;
 225 %! isosurface (x, y, z, v, 1)