Mercurial > hg > octave-nkf

comparison scripts/plot/isosurface.m @ 9110:22ae6b3411a7

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add isocolor, isonormals and isosurface functions (For Martin Helm). Add 3D filled triangular patches and the trisurf function
author David Bateman <dbateman@free.fr>
date Sat, 11 Apr 2009 16:26:01 +0200
parents
children 09da0bd91412
comparison
equal deleted inserted replaced
9109:978c863bc8e5 9110:22ae6b3411a7
1 ## Copyright (C) 2009 Martin Helm
2 ##
3 ## This program is free software; you can redistribute it and/or modify
4 ## it under the terms of the GNU General Public License as published by
5 ## the Free Software Foundation; either version 3 of the License, or
6 ## (at your option) any later version.
7 ##
8 ## This program is distributed in the hope that it will be useful,
9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 ## GNU General Public License for more details.
12 ##
13 ## You should have received a copy of the GNU General Public License
14 ## along with this program; if not, see http://www.gnu.org/licenses/gpl.html.
15
16 ## -*- texinfo -*-
17 ## @deftypefn {Function File} {[@var{fv}] =} isosurface (@var{val}, @var{iso})
18 ## @deftypefnx {Function File} {[@var{fv}] =} isosurface (@var{x}, @var{y}, @var{z}, @var{val}, @var{iso})
19 ## @deftypefnx {Function File} {[@var{fv}] =} isosurface (@dots{}, "noshare", "verbose")
20 ## @deftypefnx {Function File} {[@var{fvc}] =} isosurface (@dots{}, @var{col})
21 ## @deftypefnx {Function File} {[@var{f}, @var{v}] =} isosurface (@var{x}, @var{y}, @var{z}, @var{val}, @var{iso})
22 ## @deftypefnx {Function File} {[@var{f}, @var{v}, @var{c}] =} isosurface (@var{x}, @var{y}, @var{z}, @var{val}, @var{iso}, @var{col})
23 ## @deftypefnx {Function File} {} isosurface (@var{x}, @var{y}, @var{z}, @var{val}, @var{iso}, @var{col}, @var{opt})
24 ##
25 ## If called with one output argument and the first input argument
26 ## @var{val} is a three--dimensional array that contains the data of an
27 ## isosurface geometry and the second input argument @var{iso} keeps the
28 ## isovalue as a scalar value then return a structure array @var{fv}
29 ## that contains the fields @var{Faces} and @var{Vertices} at computed
30 ## points @command{[x, y, z] = meshgrid (1:l, 1:m, 1:n)}. The output
31 ## argument @var{fv} can directly be taken as an input argument for the
32 ## @command{patch} function.
33 ##
34 ## If called with further input arguments @var{x}, @var{y} and @var{z}
35 ## which are three--dimensional arrays with the same size than @var{val}
36 ## then the volume data is taken at those given points.
37 ##
38 ## The string input argument "noshare" is only for compatibility and
39 ## has no effect. If given the string input argument
40 ## "verbose" then print messages to the command line interface about the
41 ## current progress.
42 ##
43 ## If called with the input argument @var{col} which is a
44 ## three-dimensional array of the same size than @var{val} then take
45 ## those values for the interpolation of coloring the isosurface
46 ## geometry. Add the field @var{FaceVertexCData} to the structure
47 ## array @var{fv}.
48 ##
49 ## If called with two or three output arguments then return the
50 ## information about the faces @var{f}, vertices @var{v} and color data
51 ## @var{c} as seperate arrays instead of a single structure array.
52 ##
53 ## If called with no output argument then directly process the
54 ## isosurface geometry with the @command{patch} command.
55 ##
56 ## For example
57 ##
58 ## @example
59 ## [x, y, z] = meshgrid (1:5, 1:5, 1:5);
60 ## val = rand (5, 5, 5);
61 ## isosurface (x, y, z, val, .5);
62 ## @end example
63 ##
64 ## will directly draw a random isosurface geometry in a graphics window.
65 ## Another example for an isosurface geometry with different additional
66 ## coloring
67 ##
68 ## @example
69 ## N = 15; ## Increase number of vertices in each direction
70 ## iso = .4; ## Change isovalue to .1 to display a sphere
71 ## lin = linspace (0, 2, N);
72 ## [x, y, z] = meshgrid (lin, lin, lin);
73 ## c = abs ((x-.5).^2 + (y-.5).^2 + (z-.5).^2);
74 ## figure (); ## Open another figure window
75 ##
76 ## subplot (2, 2, 1); view (-38, 20);
77 ## [f, v] = isosurface (x, y, z, c, iso);
78 ## p = patch ("Faces", f, "Vertices", v, "EdgeColor", "none");
79 ## set (gca, "DataAspectRatioMode","manual", "DataAspectRatio", [1 1 1]);
80 ## # set (p, "FaceColor", "green", "FaceLighting", "phong");
81 ## # light ("Position", [1 1 5]); ## Available with the JHandles package
82 ##
83 ## subplot (2, 2, 2); view (-38, 20);
84 ## p = patch ("Faces", f, "Vertices", v, "EdgeColor", "blue");
85 ## set (gca, "DataAspectRatioMode","manual", "DataAspectRatio", [1 1 1]);
86 ## # set (p, "FaceColor", "none", "FaceLighting", "phong");
87 ## # light ("Position", [1 1 5]);
88 ##
89 ## subplot (2, 2, 3); view (-38, 20);
90 ## [f, v, c] = isosurface (x, y, z, c, iso, y);
91 ## p = patch ("Faces", f, "Vertices", v, "FaceVertexCData", c, \
92 ## "FaceColor", "interp", "EdgeColor", "none");
93 ## set (gca, "DataAspectRatioMode","manual", "DataAspectRatio", [1 1 1]);
94 ## # set (p, "FaceLighting", "phong");
95 ## # light ("Position", [1 1 5]);
96 ##
97 ## subplot (2, 2, 4); view (-38, 20);
98 ## p = patch ("Faces", f, "Vertices", v, "FaceVertexCData", c, \
99 ## "FaceColor", "interp", "EdgeColor", "blue");
100 ## set (gca, "DataAspectRatioMode","manual", "DataAspectRatio", [1 1 1]);
101 ## # set (p, "FaceLighting", "phong");
102 ## # light ("Position", [1 1 5]);
103 ## @end example
104 ##
105 ## @seealso{isocolors, isonormals, isocaps}
106 ##
107 ## @end deftypefn
108
109 ## Author: Martin Helm <martin@mhelm.de>
110
111 function varargout = isosurface(varargin)
112
113 if (nargin < 2 || nargin > 8 || nargout > 3)
114 print_usage ();
115 endif
116
117 calc_colors = false;
118 f = v = c = [];
119 verbose = false;
120 noshare = false;
121 if (nargin >= 5)
122 x = varargin{1};
123 y = varargin{2};
124 z = varargin{3};
125 val = varargin{4};
126 iso = varargin{5};
127 if (nargin >= 6 && ismatrix (varargin{6}))
128 colors = varargin{6};
129 calc_colors = true;
130 endif
131 else
132 val = varargin{1};
133 [n1, n2, n3] = size (val);
134 [x, y, z] = meshgrid (1:n1, 1:n2, 1:n3);
135 iso = varargin{2};
136 if (nargin >= 3 && ismatrix (varargin{3}))
137 colors = varargin{3};
138 calc_colors = true;
139 endif
140 endif
141 if (calc_colors)
142 if (nargout == 2)
143 warning ( "Colors will be calculated, but you did not specify an output argument for it!" );
144 endif
145 [fvc.faces, fvc.vertices, fvc.facevertexcdata] = __marching_cube__ (x, y, z, val, iso, colors);
146 else
147 [fvc.faces, fvc.vertices] = __marching_cube__ (x, y, z, val, iso);
148 endif
149
150 if (isempty (fvc.vertices) || isempty (fvc.faces))
151 warning ( "The resulting triangulation is empty" );
152 endif
153
154 switch (nargout)
155 case 0
156 ## plot the calculated surface
157 newplot ();
158 if (calc_colors)
159 pa = patch ("Faces", fvc.faces, "Vertices", fvc.vertices,
160 "FaceVertexCData", fvc.facevertexcdata,
161 "FaceColor", "flat", "EdgeColor", "none");
162 else
163 pa = patch ("Faces", fvc.faces, "Vertices", fvc.vertices,
164 "FaceColor", "g", "EdgeColor", "k");
165 endif
166 if (! ishold ())
167 set (gca(), "view", [-37.5, 30],
168 "xgrid", "on", "ygrid", "on", "zgrid", "on");
169 endif
170 case 1
171 varargout = {fvc};
172 case 2
173 varargout = {fvc.faces, fvc.vertices};
174 case 3
175 varargout = {fvc.faces, fvc.vertices, fvc.facevertexcdata};
176 otherwise
177 print_usage ();
178 endswitch
179
180 endfunction
181
182
183 %!shared x, y, z, val
184 %! [x, y, z] = meshgrid (0:1, 0:1, 0:1); ## Points for single
185 %! val = [0, 0; 0, 0]; ## cube and a 3--dim
186 %! val(:,:,2) = [0, 0; 1, 0]; ## array of values
187 %!test
188 %! fv = isosurface (x, y, z, val, 0.3);
189 %! assert (isfield (fv, "vertices"), true);
190 %! assert (isfield (fv, "faces"), true);
191 %! assert (size (fv.vertices), [3 3]);
192 %! assert (size (fv.faces), [1 3]);
193 %!test
194 %! fvc = isosurface (x, y, z, val, .3, y);
195 %! assert (isfield (fvc, "vertices"), true);
196 %! assert (isfield (fvc, "faces"), true);
197 %! assert (isfield (fvc, "facevertexcdata"), true);
198 %! assert (size (fvc.vertices), [3 3]);
199 %! assert (size (fvc.faces), [1 3]);
200 %! assert (size (fvc.facevertexcdata), [3 1]);
201 %!test
202 %! [f, v] = isosurface (x, y, z, val, .3);
203 %! assert (size (f), [1 3]);
204 %! assert (size (v), [3 3]);
205 %!test
206 %! [f, v, c] = isosurface (x, y, z, val, .3, y);
207 %! assert (size (f), [1 3]);
208 %! assert (size (v), [3 3]);
209 %! assert (size (c), [3 1]);
210
211 %!demo
212 %! clf
213 %! [x,y,z] = meshgrid(-2:0.5:2, -2:0.5:2, -2:0.5:2);
214 %! v = x.^2 + y.^2 + z.^2;
215 %! isosurface (x, y, z, v, 1)