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