1 ## Copyright (C) 2009-2012 Martin Helm
3 ## This file is part of Octave.
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.
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.
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/>.
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})
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.
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.
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
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
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.
56 ## If called with no output argument then directly process the
57 ## isosurface geometry with the @command{patch} command.
63 ## [x, y, z] = meshgrid (1:5, 1:5, 1:5);
64 ## val = rand (5, 5, 5);
65 ## isosurface (x, y, z, val, .5);
70 ## will directly draw a random isosurface geometry in a graphics window.
71 ## Another example for an isosurface geometry with different additional
73 ## @c Set example in small font to prevent overfull line
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
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
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]);
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]);
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]);
116 ## @seealso{isonormals, isocolors}
119 ## Author: Martin Helm <martin@mhelm.de>
121 function varargout = isosurface(varargin)
123 if (nargin < 2 || nargin > 8 || nargout > 3)
137 if (nargin >= 6 && ismatrix (varargin{6}))
138 colors = varargin{6};
143 [n2, n1, n3] = size (val);
144 [x, y, z] = meshgrid (1:n1, 1:n2, 1:n3);
146 if (nargin >= 3 && ismatrix (varargin{3}))
147 colors = varargin{3};
153 warning ( "Colors will be calculated, but you did not specify an output argument for it!" );
155 [fvc.faces, fvc.vertices, fvc.facevertexcdata] = __marching_cube__ (x, y, z, val, iso, colors);
157 [fvc.faces, fvc.vertices] = __marching_cube__ (x, y, z, val, iso);
160 if (isempty (fvc.vertices) || isempty (fvc.faces))
161 warning ( "The resulting triangulation is empty" );
166 ## plot the calculated surface
169 pa = patch ("Faces", fvc.faces, "Vertices", fvc.vertices,
170 "FaceVertexCData", fvc.facevertexcdata,
171 "FaceColor", "flat", "EdgeColor", "none");
173 pa = patch ("Faces", fvc.faces, "Vertices", fvc.vertices,
174 "FaceColor", "g", "EdgeColor", "k");
177 set (gca(), "view", [-37.5, 30],
178 "xgrid", "on", "ygrid", "on", "zgrid", "on");
183 varargout = {fvc.faces, fvc.vertices};
185 varargout = {fvc.faces, fvc.vertices, fvc.facevertexcdata};
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
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]);
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]);
212 %! [f, v] = isosurface (x, y, z, val, .3);
213 %! assert (size (f), [1 3]);
214 %! assert (size (v), [3 3]);
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]);
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)