1######################################################################## 2## 3## Copyright (C) 2009-2021 The Octave Project Developers 4## 5## See the file COPYRIGHT.md in the top-level directory of this 6## distribution or <https://octave.org/copyright/>. 7## 8## This file is part of Octave. 9## 10## Octave is free software: you can redistribute it and/or modify it 11## under the terms of the GNU General Public License as published by 12## the Free Software Foundation, either version 3 of the License, or 13## (at your option) any later version. 14## 15## Octave is distributed in the hope that it will be useful, but 16## WITHOUT ANY WARRANTY; without even the implied warranty of 17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18## GNU General Public License for more details. 19## 20## You should have received a copy of the GNU General Public License 21## along with Octave; see the file COPYING. If not, see 22## <https://www.gnu.org/licenses/>. 23## 24######################################################################## 25 26## -*- texinfo -*- 27## @deftypefn {} {@var{fv} =} isosurface (@var{v}, @var{isoval}) 28## @deftypefnx {} {@var{fv} =} isosurface (@var{v}) 29## @deftypefnx {} {@var{fv} =} isosurface (@var{x}, @var{y}, @var{z}, @var{v}, @var{isoval}) 30## @deftypefnx {} {@var{fv} =} isosurface (@var{x}, @var{y}, @var{z}, @var{v}) 31## @deftypefnx {} {@var{fvc} =} isosurface (@dots{}, @var{col}) 32## @deftypefnx {} {@var{fv} =} isosurface (@dots{}, "noshare") 33## @deftypefnx {} {@var{fv} =} isosurface (@dots{}, "verbose") 34## @deftypefnx {} {[@var{f}, @var{v}] =} isosurface (@dots{}) 35## @deftypefnx {} {[@var{f}, @var{v}, @var{c}] =} isosurface (@dots{}) 36## @deftypefnx {} {} isosurface (@dots{}) 37## 38## Calculate isosurface of 3-D volume data. 39## 40## An isosurface connects points with the same value and is analogous to a 41## contour plot, but in three dimensions. 42## 43## The input argument @var{v} is a three-dimensional array that contains data 44## sampled over a volume. 45## 46## The input @var{isoval} is a scalar that specifies the value for the 47## isosurface. If @var{isoval} is omitted or empty, a @nospell{"good"} value 48## for an isosurface is determined from @var{v}. 49## 50## When called with a single output argument @code{isosurface} returns a 51## structure array @var{fv} that contains the fields @var{faces} and 52## @var{vertices} computed at the points 53## @code{[@var{x}, @var{y}, @var{z}] = meshgrid (1:l, 1:m, 1:n)} where 54## @code{[l, m, n] = size (@var{v})}. The output @var{fv} can be 55## used directly as input to the @code{patch} function. 56## 57## If called with additional input arguments @var{x}, @var{y}, and @var{z} 58## that are three-dimensional arrays with the same size as @var{v} or 59## vectors with lengths corresponding to the dimensions of @var{v}, then the 60## volume data is taken at the specified points. If @var{x}, @var{y}, or 61## @var{z} are empty, the grid corresponds to the indices (@code{1:n}) in 62## the respective direction (@pxref{XREFmeshgrid,,meshgrid}). 63## 64## The optional input argument @var{col}, which is a three-dimensional array 65## of the same size as @var{v}, specifies coloring of the isosurface. The 66## color data is interpolated, as necessary, to match @var{isoval}. The 67## output structure array, in this case, has the additional field 68## @var{facevertexcdata}. 69## 70## If given the string input argument @qcode{"noshare"}, vertices may be 71## returned multiple times for different faces. The default behavior is to 72## eliminate vertices shared by adjacent faces. 73## 74## The string input argument @qcode{"verbose"} is supported for @sc{matlab} 75## compatibility, but has no effect. 76## 77## Any string arguments must be passed after the other arguments. 78## 79## If called with two or three output arguments, return the information about 80## the faces @var{f}, vertices @var{v}, and color data @var{c} as separate 81## arrays instead of a single structure array. 82## 83## If called with no output argument, the isosurface geometry is directly 84## plotted with the @code{patch} command and a light object is added to 85## the axes if not yet present. 86## 87## For example, 88## 89## @example 90## @group 91## [x, y, z] = meshgrid (1:5, 1:5, 1:5); 92## v = rand (5, 5, 5); 93## isosurface (x, y, z, v, .5); 94## @end group 95## @end example 96## 97## @noindent 98## will directly draw a random isosurface geometry in a graphics window. 99## 100## An example of an isosurface geometry with different additional coloring: 101## @c Set example in small font to prevent overfull line 102## 103## @smallexample 104## N = 15; # Increase number of vertices in each direction 105## iso = .4; # Change isovalue to .1 to display a sphere 106## lin = linspace (0, 2, N); 107## [x, y, z] = meshgrid (lin, lin, lin); 108## v = abs ((x-.5).^2 + (y-.5).^2 + (z-.5).^2); 109## figure (); 110## 111## subplot (2,2,1); view (-38, 20); 112## [f, vert] = isosurface (x, y, z, v, iso); 113## p = patch ("Faces", f, "Vertices", vert, "EdgeColor", "none"); 114## pbaspect ([1 1 1]); 115## isonormals (x, y, z, v, p) 116## set (p, "FaceColor", "green", "FaceLighting", "gouraud"); 117## light ("Position", [1 1 5]); 118## 119## subplot (2,2,2); view (-38, 20); 120## p = patch ("Faces", f, "Vertices", vert, "EdgeColor", "blue"); 121## pbaspect ([1 1 1]); 122## isonormals (x, y, z, v, p) 123## set (p, "FaceColor", "none", "EdgeLighting", "gouraud"); 124## light ("Position", [1 1 5]); 125## 126## subplot (2,2,3); view (-38, 20); 127## [f, vert, c] = isosurface (x, y, z, v, iso, y); 128## p = patch ("Faces", f, "Vertices", vert, "FaceVertexCData", c, ... 129## "FaceColor", "interp", "EdgeColor", "none"); 130## pbaspect ([1 1 1]); 131## isonormals (x, y, z, v, p) 132## set (p, "FaceLighting", "gouraud"); 133## light ("Position", [1 1 5]); 134## 135## subplot (2,2,4); view (-38, 20); 136## p = patch ("Faces", f, "Vertices", vert, "FaceVertexCData", c, ... 137## "FaceColor", "interp", "EdgeColor", "blue"); 138## pbaspect ([1 1 1]); 139## isonormals (x, y, z, v, p) 140## set (p, "FaceLighting", "gouraud"); 141## light ("Position", [1 1 5]); 142## @end smallexample 143## 144## @seealso{isonormals, isocolors, isocaps, smooth3, reducevolume, reducepatch, patch} 145## @end deftypefn 146 147## FIXME: Add support for string input argument "verbose" 148## (needs changes to __marching_cube__.m) 149 150function varargout = isosurface (varargin) 151 152 if (nargin < 1 || nargin > 8) 153 print_usage (); 154 endif 155 156 [x, y, z, v, isoval, colors, noshare, verbose] = ... 157 __get_check_isosurface_args__ (nargout, varargin{:}); 158 159 calc_colors = ! isempty (colors); 160 if (calc_colors) 161 [fvc.faces, fvc.vertices, fvc.facevertexcdata] = ... 162 __marching_cube__ (x, y, z, v, isoval, colors); 163 else 164 [fvc.faces, fvc.vertices] = __marching_cube__ (x, y, z, v, isoval); 165 endif 166 167 if (isempty (fvc.vertices) || isempty (fvc.faces)) 168 warning ("isosurface: triangulation is empty"); 169 endif 170 171 # remove faces for which at least one of the vertices is NaN 172 vert_nan = 1:size (fvc.vertices, 1); 173 vert_nan(any (isnan (fvc.vertices), 2)) = NaN; 174 fvc.faces = vert_nan(fvc.faces); 175 fvc.faces(any (isnan (fvc.faces), 2), :) = []; 176 177 if (! noshare) 178 [fvc.faces, fvc.vertices, J] = __unite_shared_vertices__ (fvc.faces, 179 fvc.vertices); 180 181 if (calc_colors) 182 fvc.facevertexcdata = fvc.facevertexcdata(J); # share very close vertices 183 endif 184 endif 185 186 switch (nargout) 187 case 0 188 ## plot the calculated surface 189 if (calc_colors) 190 fc = fvc.facevertexcdata; 191 else 192 fc = isoval; 193 endif 194 ## Matlab uses "EdgeColor", "none", but that looks odd in gnuplot. 195 hax = gca (); 196 if (strcmp (get (gcf, "__graphics_toolkit__"), "gnuplot")) 197 ec = "k"; 198 else 199 ec = "none"; 200 endif 201 pa = patch ("Faces", fvc.faces, "Vertices", fvc.vertices, 202 "FaceVertexCData", fc, 203 "FaceColor", "flat", "EdgeColor", ec, 204 "FaceLighting", "gouraud"); 205 if (! ishold ()) 206 set (hax, "View", [-37.5, 30]); 207 endif 208 isonormals (x, y, z, v, pa); 209 lights = findobj (hax, "Type", "light"); 210 if (isempty (lights)) 211 camlight (); 212 endif 213 214 case 1 215 varargout = {fvc}; 216 217 case 2 218 varargout = {fvc.faces, fvc.vertices}; 219 220 otherwise # 3 args or more 221 varargout = {fvc.faces, fvc.vertices, fvc.facevertexcdata}; 222 223 endswitch 224 225endfunction 226 227function [x, y, z, v, isoval, colors, noshare, verbose] = __get_check_isosurface_args__ (nout, varargin) 228 ## get arguments from input and check values 229 x = y = z = []; 230 v = []; 231 isoval = []; 232 colors = []; 233 234 ## default values 235 noshare = false; 236 verbose = false; 237 238 nin = length (varargin); 239 num_string_inputs = 0; 240 241 ## check whether last 2 input arguments are strings and assign parameters 242 for i_arg = (nin:-1:nin-1) 243 if (! ischar (varargin{i_arg}) || i_arg < 1) 244 break; # no string arguments at end, exit checking 245 endif 246 switch (tolower (varargin{i_arg})) 247 case {"v", "verbose"} 248 verbose = true; 249 num_string_inputs++; 250 251 case {"n", "noshare"} 252 noshare = true; 253 num_string_inputs++; 254 255 case "" 256 num_string_inputs++; 257 ## silently ignore empty strings 258 259 otherwise 260 error ("isosurface: parameter '%s' not supported", varargin{i_arg}); 261 262 endswitch 263 endfor 264 265 ## assign arguments 266 switch (nin - num_string_inputs) 267 case 1 # isosurface (v, ...) 268 v = varargin{1}; 269 270 case 2 # isosurface (v, isoval, ...) or isosurface (v, col, ...) 271 v = varargin{1}; 272 if (isscalar (varargin{2}) || isempty (varargin{2})) 273 isoval = varargin{2}; 274 else 275 colors = varargin{2}; 276 endif 277 278 case 3 # isosurface (v, isoval, col, ...) 279 v = varargin{1}; 280 isoval = varargin{2}; 281 colors = varargin{3}; 282 283 case 4 # isosurface (x, y, z, v, ...) 284 x = varargin{1}; 285 y = varargin{2}; 286 z = varargin{3}; 287 v = varargin{4}; 288 289 case 5 # isosurface (x, y, z, v, isoval, ...) or 290 # isosurface (x, y, z, v, col, ...) 291 x = varargin{1}; 292 y = varargin{2}; 293 z = varargin{3}; 294 v = varargin{4}; 295 if (isscalar (varargin{5}) || isempty (varargin{5})) 296 isoval = varargin{5}; 297 else 298 colors = varargin{5}; 299 endif 300 301 case 6 # isosurface (x, y, z, v, isoval, col, ...) 302 x = varargin{1}; 303 y = varargin{2}; 304 z = varargin{3}; 305 v = varargin{4}; 306 isoval = varargin{5}; 307 colors = varargin{6}; 308 309 otherwise 310 error ("isosurface: incorrect number of input arguments") 311 312 endswitch 313 314 ## check dimensions of v 315 v_sz = size (v); 316 if (ndims (v) != 3 || any (v_sz(1:3) < 2)) 317 error ("isosurface: V must be a non-singleton 3-dimensional matrix"); 318 endif 319 320 if (isempty (x)) 321 x = 1:size (v, 2); 322 endif 323 if (isempty (y)) 324 y = 1:size (v, 1); 325 endif 326 if (isempty (z)) 327 z = 1:size (v, 3); 328 endif 329 330 ## check x 331 if (isvector (x) && length (x) == v_sz(2)) 332 x = repmat (x(:)', [v_sz(1) 1 v_sz(3)]); 333 elseif (! size_equal (v, x)) 334 error ("isosurface: X must match the size of V"); 335 endif 336 337 ## check y 338 if (isvector (y) && length (y) == v_sz(1)) 339 y = repmat (y(:), [1 v_sz(2) v_sz(3)]); 340 elseif (! size_equal (v, y)) 341 error ("isosurface: Y must match the size of V"); 342 endif 343 344 ## check z 345 if (isvector (z) && length (z) == v_sz(3)) 346 z = repmat (reshape (z(:), [1 1 length(z)]), [v_sz(1) v_sz(2) 1]); 347 elseif (! size_equal (v, z)) 348 error ("isosurface: Z must match the size of V"); 349 endif 350 351 ## check isoval 352 if (isempty (isoval)) 353 ## calculate "good" isoval value from v 354 isoval = __calc_isovalue_from_data__ (v); 355 endif 356 357 if (! isscalar (isoval)) 358 error ("isosurface: ISOVAL must be a scalar") 359 endif 360 361 ## check colors 362 if (! isempty (colors)) 363 if (! size_equal (v, colors)) 364 error ("isosurface: COL must match the size of V") 365 endif 366 if (nout == 2) 367 warning ("isosurface: colors will be calculated, but no output argument to receive it."); 368 endif 369 elseif (nout >= 3) 370 error ("isosurface: COL must be passed to return C") 371 endif 372 373endfunction 374 375 376%!demo 377%! clf; 378%! [x,y,z] = meshgrid (-2:0.5:2, -2:0.5:2, -2:0.5:2); 379%! v = x.^2 + y.^2 + z.^2; 380%! isosurface (x, y, z, v, 1); 381%! axis equal; 382%! title ("isosurface() of a sphere"); 383 384%!demo 385%! clf; 386%! [x,y,z] = meshgrid (-2:0.5:2, -2:0.5:2, -2:0.5:2); 387%! v = x.^2 + y.^2 + z.^2; 388%! isosurface (x, y, z, v, 3); 389%! isosurface (x, y, z, v, 5); 390%! axis equal; 391%! title ("isosurface() of two nested spheres"); 392 393%!demo 394%! clf; 395%! x = 0:2; 396%! y = 0:3; 397%! z = 0:1; 398%! [xx, yy, zz] = meshgrid (x, y, z); 399%! v = [0, 0, 0; 0, 0, 0; 0, 0, 1; 0, 0, 1]; 400%! v(:,:,2) = [0, 0, 0; 0, 0, 1; 0, 1, 2; 0, 1, 2]; 401%! iso = 0.8; 402%! 403%! ## Three arguments, no output 404%! subplot (2, 2, 1); 405%! fvc = isosurface (v, iso, yy); 406%! patch (fvc); 407%! shading faceted; 408%! view (110, 40); 409%! ## six arguments, no output (x, y, z are vectors) 410%! subplot (2, 2, 2); 411%! fvc = isosurface (x, y, z, v, iso, yy); 412%! patch (fvc); 413%! shading faceted; 414%! view (110, 40); 415%! ## six arguments, no output (x, y, z are matrices) 416%! subplot (2, 2, 3); 417%! fvc = isosurface (xx, yy, zz, v, iso, yy); 418%! patch (fvc); 419%! shading faceted; 420%! view (110, 40); 421%! ## six arguments, no output (mixed x, y, z) and option "noshare" 422%! subplot (2, 2, 4); 423%! fvc = isosurface (x, yy, z, v, iso, yy, "noshare"); 424%! patch (fvc); 425%! shading faceted; 426%! view (110, 40); 427%! annotation ("textbox", [0.41 0.9 0.9 0.1], ... 428%! "String", "isosurface() called 4 ways", ... 429%! "HorizontalAlignment", "center", ... 430%! "FontSize", 12); 431%! annotation ("textbox", [0.1 0.45 0.9 0.1], ... 432%! "String", {["Apart from the first plot having a different scale, " ... 433%! "all four plots must look the same."], 434%! ["The last plot might have different colors but must " ... 435%! "have the same shape."]}, ... 436%! "HorizontalAlignment", "left", ... 437%! "FontSize", 12); 438 439%!shared x, y, z, xx, yy, zz, val, iso 440%! x = 0:2; 441%! y = 0:3; 442%! z = 0:1; 443%! [xx, yy, zz] = meshgrid (x, y, z); 444%! val = [0, 0, 0; 0, 0, 0; 0, 0, 1; 0, 0, 1]; 445%! val(:,:,2) = [0, 0, 0; 0, 0, 1; 0, 1, 2; 0, 1, 2]; 446%! iso = 0.8; 447 448## one argument, one output 449%!test 450%! fv = isosurface (val); 451%! assert (isfield (fv, "vertices"), true); 452%! assert (isfield (fv, "faces"), true); 453%! assert (size (fv.vertices), [5 3]); 454%! assert (size (fv.faces), [3 3]); 455 456## two arguments (second is ISO), one output 457%!test 458%! fv = isosurface (val, iso); 459%! assert (isfield (fv, "vertices"), true); 460%! assert (isfield (fv, "faces"), true); 461%! assert (size (fv.vertices), [11 3]); 462%! assert (size (fv.faces), [10 3]); 463 464## two arguments (second is COL), one output 465%!test 466%! fvc = isosurface (val, yy); 467%! assert (isfield (fvc, "vertices"), true); 468%! assert (isfield (fvc, "faces"), true); 469%! assert (isfield (fvc, "facevertexcdata"), true); 470%! assert (size (fvc.vertices), [5 3]); 471%! assert (size (fvc.faces), [3 3]); 472%! assert (size (fvc.facevertexcdata), [5 1]); 473 474## three arguments, one output 475%!test 476%! fvc = isosurface (val, iso, yy); 477%! assert (isfield (fvc, "vertices"), true); 478%! assert (isfield (fvc, "faces"), true); 479%! assert (isfield (fvc, "facevertexcdata"), true); 480%! assert (size (fvc.vertices), [11 3]); 481%! assert (size (fvc.faces), [10 3]); 482%! assert (size (fvc.facevertexcdata), [11 1]); 483 484## four arguments, one output 485%!test 486%! fv = isosurface (x, [], z, val); 487%! assert (isfield (fv, "vertices"), true); 488%! assert (isfield (fv, "faces"), true); 489%! assert (size (fv.vertices), [5 3]); 490%! assert (size (fv.faces), [3 3]); 491 492## five arguments (fifth is ISO), one output 493%!test 494%! fv = isosurface (xx, y, [], val, iso); 495%! assert (isfield (fv, "vertices"), true); 496%! assert (isfield (fv, "faces"), true); 497%! assert (size (fv.vertices), [11 3]); 498%! assert (size (fv.faces), [10 3]); 499 500## five arguments (fifth is COL), one output 501%!test 502%! fvc = isosurface ([], yy, z, val, yy); 503%! assert (isfield (fvc, "vertices"), true); 504%! assert (isfield (fvc, "faces"), true); 505%! assert (isfield (fvc, "facevertexcdata"), true); 506%! assert (size (fvc.vertices), [5 3]); 507%! assert (size (fvc.faces), [3 3]); 508%! assert (size (fvc.facevertexcdata), [5 1]); 509 510## six arguments, one output 511%!test 512%! fvc = isosurface (xx, yy, zz, val, iso, yy); 513%! assert (isfield (fvc, "vertices"), true); 514%! assert (isfield (fvc, "faces"), true); 515%! assert (isfield (fvc, "facevertexcdata"), true); 516%! assert (size (fvc.vertices), [11 3]); 517%! assert (size (fvc.faces), [10 3]); 518%! assert (size (fvc.facevertexcdata), [11 1]); 519 520## five arguments (fifth is ISO), two outputs 521%!test 522%! [f, v] = isosurface (x, y, z, val, iso); 523%! assert (size (f), [10 3]); 524%! assert (size (v), [11 3]); 525 526## six arguments, three outputs 527%!test 528%! [f, v, c] = isosurface (x, y, z, val, iso, yy); 529%! assert (size (f), [10 3]); 530%! assert (size (v), [11 3]); 531%! assert (size (c), [11 1]); 532 533## two arguments (second is ISO) and one string, one output 534%!test 535%! fv = isosurface (val, iso, "verbose"); 536%! assert (isfield (fv, "vertices"), true); 537%! assert (isfield (fv, "faces"), true); 538%! assert (size (fv.vertices), [11 3]); 539%! assert (size (fv.faces), [10 3]); 540 541## six arguments and two strings, one output 542%!test 543%! fvc = isosurface (xx, yy, zz, val, iso, yy, "v", "noshare"); 544%! assert (isfield (fvc, "vertices"), true); 545%! assert (isfield (fvc, "faces"), true); 546%! assert (isfield (fvc, "facevertexcdata"), true); 547%! assert (size (fvc.vertices), [20 3]); 548%! assert (size (fvc.faces), [10 3]); 549%! assert (size (fvc.facevertexcdata), [20 1]); 550 551## five arguments (fifth is COL) and two strings (different order), one output 552%!test 553%! fvc = isosurface (xx, yy, zz, val, yy, "n", "v"); 554%! assert (isfield (fvc, "vertices"), true); 555%! assert (isfield (fvc, "faces"), true); 556%! assert (isfield (fvc, "facevertexcdata"), true); 557%! assert (size (fvc.vertices), [7 3]); 558%! assert (size (fvc.faces), [3 3]); 559%! assert (size (fvc.facevertexcdata), [7 1]); 560 561## test for each error and warning 562%!error isosurface () 563%!error isosurface (1,2,3,4,5,6,7,8,9) 564%!error <parameter 'foobar' not supported> 565%! fvc = isosurface (val, iso, "foobar"); 566%!error <incorrect number of input arguments> 567%! fvc = isosurface (xx, yy, zz, val, iso, yy, 5); 568%!error <V must be a non-singleton 3-dimensional matrix> 569%! v = reshape (1:6*8, [6 8]); 570%! fvc = isosurface (v, iso); 571%!error <V must be a non-singleton 3-dimensional matrix> 572%! v = reshape(1:6*8, [6 1 8]); fvc = isosurface (v, iso); 573%!error <X must match the size of V> 574%! x = 1:2:24; 575%! fvc = isosurface (x, y, z, val, iso); 576%!error <Y must match the size of V> 577%! y = -14:6:11; 578%! fvc = isosurface (x, y, z, val, iso); 579%!error <Z must match the size of V> 580%! z = linspace (16, 18, 5); 581%! fvc = isosurface (x, y, z, val, iso); 582%!error <X must match the size of V> 583%! x = 1:2:24; 584%! [xx, yy, zz] = meshgrid (x, y, z); 585%! fvc = isosurface (xx, yy, zz, val, iso); 586%!error <X must match the size of V> 587%! y = -14:6:11; 588%! [xx, yy, zz] = meshgrid (x, y, z); 589%! fvc = isosurface (xx, yy, zz, val, iso); 590%!error <X must match the size of V> 591%! z = linspace (16, 18, 3); 592%! [xx, yy, zz] = meshgrid (x, y, z); 593%! fvc = isosurface (xx, yy, zz, val, iso); 594%!error <ISOVAL must be a scalar> fvc = isosurface (val, [iso iso], yy) 595%!error <COL must match the size of V> fvc = isosurface (val, [iso iso]); 596%!error <COL must be passed to return C> [f, v, c] = isosurface (val, iso) 597%!warning <colors will be calculated, but no output argument to receive it.> 598%! [f, v] = isosurface (val, iso, yy); 599 600## test for __calc_isovalue_from_data__ 601## FIXME: private function cannot be tested, unless bug #38776 is resolved. 602%!test <38776> 603%! assert (__calc_isovalue_from_data__ (1:5), 3.02); 604