# HG changeset patch # User Rik # Date 1470766178 25200 # Node ID a8fd02bc895b4672f7d95133b6856a93b7495d37 # Parent 01ba6ebc52e475ee602c824ca31abc45eaaf6c7f Overhaul isosurface.m and associated functions. * isosurface.m: Rewrite docstring and examples. Don't validate nargout. Wrap long lines to < 80 characters. Use names from docstring in code. Use "EdgeColor", "none" for Matlab compatibility, but preserve exception for gnuplot. Use blank lines between case statements for readability. Update demo #2 to properly show annotation textbox. Rewrite %!error tests. * __calc_isovalue_from_data__.m: Give private function a docstring. Use '_' in numbers as thousands separator for readability. Wrap long lines to < 80 characters. * __unite_shared_vertices__.m: Reformat docstring. Use rows() rather than size (XXX, 1). Wrap long lines to < 80 characters. Use space after function name and before parenthesis. diff -r 01ba6ebc52e4 -r a8fd02bc895b scripts/plot/draw/isosurface.m --- a/scripts/plot/draw/isosurface.m Mon Jul 18 19:49:39 2016 +0200 +++ b/scripts/plot/draw/isosurface.m Tue Aug 09 11:09:38 2016 -0700 @@ -17,58 +17,65 @@ ## . ## -*- texinfo -*- -## @deftypefn {} {[@var{fv}] =} isosurface (@var{val}, @var{iso}) -## @deftypefnx {} {[@var{fv}] =} isosurface (@var{val}) -## @deftypefnx {} {[@var{fv}] =} isosurface (@var{x}, @var{y}, @var{z}, @var{val}, @var{iso}) -## @deftypefnx {} {[@var{fv}] =} isosurface (@var{x}, @var{y}, @var{z}, @var{val}) -## @deftypefnx {} {[@var{fvc}] =} isosurface (@dots{}, @var{col}) -## @deftypefnx {} {[@var{fv}] =} isosurface (@dots{}, "noshare") -## @deftypefnx {} {[@var{fv}] =} isosurface (@dots{}, "verbose") +## @deftypefn {} {@var{fv} =} isosurface (@var{v}, @var{isoval}) +## @deftypefnx {} {@var{fv} =} isosurface (@var{v}) +## @deftypefnx {} {@var{fv} =} isosurface (@var{x}, @var{y}, @var{z}, @var{v}, @var{isoval}) +## @deftypefnx {} {@var{fv} =} isosurface (@var{x}, @var{y}, @var{z}, @var{v}) +## @deftypefnx {} {@var{fvc} =} isosurface (@dots{}, @var{col}) +## @deftypefnx {} {@var{fv} =} isosurface (@dots{}, "noshare") +## @deftypefnx {} {@var{fv} =} isosurface (@dots{}, "verbose") ## @deftypefnx {} {[@var{f}, @var{v}] =} isosurface (@dots{}) ## @deftypefnx {} {[@var{f}, @var{v}, @var{c}] =} isosurface (@dots{}) ## @deftypefnx {} {} isosurface (@dots{}) ## -## Calculate isosurface of 3-D data. +## Calculate isosurface of 3-D volume data. +## +## An isosurface connects points with the same value and is analogous to a +## contour plot, but in three dimensions. ## -## If called with one output argument and the first input argument -## @var{val} is a three-dimensional array that contains the data of an -## isosurface geometry and the second input argument @var{iso} keeps the -## isovalue as a scalar value then return a structure array @var{fv} -## that contains the fields @var{faces} and @var{vertices} at computed -## points @command{[x, y, z] = meshgrid (1:l, 1:m, 1:n)}. The output -## argument @var{fv} can directly be taken as an input argument for the -## @command{patch} function. +## The input argument @var{v} is a three-dimensional array that contains data +## sampled over a volume. +## +## The input @var{isoval} is a scalar that specifies the value for the +## isosurface. If @var{isoval} is omitted or empty, a @qcode{"good"} value +## for an isosurface is determined from @var{v}. ## -## If @var{iso} is omitted or empty, a "good" value for an isosurface is -## determined from @var{val}. +## When called with a single output argument @code{isosurface} returns a +## structure array @var{fv} that contains the fields @var{faces} and +## @var{vertices} computed at the points +## @code{[@var{x}, @var{y}, @var{z}] = meshgrid (1:l, 1:m, 1:n)} where +## @code{[l, m, n] = size (@var{v})}. The output @var{fv} can be +## used directly as input to the @code{patch} function. ## -## If called with further input arguments @var{x}, @var{y} and @var{z} -## which are three--dimensional arrays with the same size as @var{val} or -## vectors with lengths corresponding to the dimensions of @var{val}, then the -## volume data is taken at those given points. If @var{x}, @var{y} or @var{z} -## are empty, the grid corresponds to the indices in the respective direction -## (see @command{meshgrid}). +## If called with additional input arguments @var{x}, @var{y}, and @var{z} +## that are three-dimensional arrays with the same size as @var{v} or +## vectors with lengths corresponding to the dimensions of @var{v}, then the +## volume data is taken at the specified points. If @var{x}, @var{y}, or +## @var{z} are empty, the grid corresponds to the indices (@code{1:n}) in +## the respective direction (@pxref{XREFmeshgrid,,meshgrid}). ## -## If called with the input argument @var{col} which is a -## three-dimensional array of the same size as @var{val}, take -## those values for the interpolation of coloring the isosurface -## geometry. In this case, the structure array @var{fv} has the additional field +## The optional input argument @var{col}, which is a three-dimensional array +## of the same size as @var{v}, specifies coloring of the isosurface. The +## color data is interpolated, as necessary, to match @var{isoval}. The +## output structure array, in this case, has the additional field ## @var{facevertexcdata}. ## -## If given the string input argument @qcode{"noshare"}, vertices might be -## returned multiple times for different faces. The default behavior is to -## search vertices shared by adjacent faces with @command{unique} which might be +## If given the string input argument @qcode{"noshare"}, vertices may be +## returned multiple times for different faces. The default behavior is to +## eliminate vertices shared by adjacent faces with @code{unique} which may be ## time consuming. -## The string input argument @qcode{"verbose"} is only for compatibility and -## has no effect. -## The string input arguments must be passed after the other arguments. +## +## The string input argument @qcode{"verbose"} is supported for @sc{matlab} +## compatibility, but has no effect. +## +## Any string arguments must be passed after the other arguments. ## ## If called with two or three output arguments, return the information about -## the faces @var{f}, vertices @var{v} and color data @var{c} as separate arrays -## instead of a single structure array. +## the faces @var{f}, vertices @var{v}, and color data @var{c} as separate +## arrays instead of a single structure array. ## ## If called with no output argument, the isosurface geometry is directly -## processed with the @command{patch} command and a light object is added to +## plotted with the @code{patch} command and a light object is added to ## the axes if not yet present. ## ## For example, @@ -76,16 +83,15 @@ ## @example ## @group ## [x, y, z] = meshgrid (1:5, 1:5, 1:5); -## val = rand (5, 5, 5); -## isosurface (x, y, z, val, .5); +## v = rand (5, 5, 5); +## isosurface (x, y, z, v, .5); ## @end group ## @end example ## ## @noindent ## will directly draw a random isosurface geometry in a graphics window. ## -## Another example for an isosurface geometry with different additional -## coloring +## An example of an isosurface geometry with different additional coloring: ## @c Set example in small font to prevent overfull line ## ## @smallexample @@ -93,42 +99,38 @@ ## iso = .4; # Change isovalue to .1 to display a sphere ## lin = linspace (0, 2, N); ## [x, y, z] = meshgrid (lin, lin, lin); -## val = abs ((x-.5).^2 + (y-.5).^2 + (z-.5).^2); -## figure (); # Open another figure window +## v = abs ((x-.5).^2 + (y-.5).^2 + (z-.5).^2); +## figure (); ## ## subplot (2,2,1); view (-38, 20); -## [f, v] = isosurface (x, y, z, val, iso); -## p = patch ("Faces", f, "Vertices", v, "EdgeColor", "none"); -## set (gca, "PlotBoxAspectRatioMode", "manual", ... -## "PlotBoxAspectRatio", [1 1 1]); -## isonormals (x, y, z, val, p) +## [f, vert] = isosurface (x, y, z, v, iso); +## p = patch ("Faces", f, "Vertices", vert, "EdgeColor", "none"); +## pbaspect ([1 1 1]); +## isonormals (x, y, z, v, p) ## set (p, "FaceColor", "green", "FaceLighting", "gouraud"); ## light ("Position", [1 1 5]); ## ## subplot (2,2,2); view (-38, 20); -## p = patch ("Faces", f, "Vertices", v, "EdgeColor", "blue"); -## set (gca, "PlotBoxAspectRatioMode", "manual", ... -## "PlotBoxAspectRatio", [1 1 1]); -## isonormals (x, y, z, val, p) +## p = patch ("Faces", f, "Vertices", vert, "EdgeColor", "blue"); +## pbaspect ([1 1 1]); +## isonormals (x, y, z, v, p) ## set (p, "FaceColor", "none", "EdgeLighting", "gouraud"); ## light ("Position", [1 1 5]); ## ## subplot (2,2,3); view (-38, 20); -## [f, v, c] = isosurface (x, y, z, val, iso, y); -## p = patch ("Faces", f, "Vertices", v, "FaceVertexCData", c, ... +## [f, vert, c] = isosurface (x, y, z, v, iso, y); +## p = patch ("Faces", f, "Vertices", vert, "FaceVertexCData", c, ... ## "FaceColor", "interp", "EdgeColor", "none"); -## set (gca, "PlotBoxAspectRatioMode", "manual", ... -## "PlotBoxAspectRatio", [1 1 1]); -## isonormals (x, y, z, val, p) +## pbaspect ([1 1 1]); +## isonormals (x, y, z, v, p) ## set (p, "FaceLighting", "gouraud"); ## light ("Position", [1 1 5]); ## ## subplot (2,2,4); view (-38, 20); -## p = patch ("Faces", f, "Vertices", v, "FaceVertexCData", c, ... +## p = patch ("Faces", f, "Vertices", vert, "FaceVertexCData", c, ... ## "FaceColor", "interp", "EdgeColor", "blue"); -## set (gca, "PlotBoxAspectRatioMode", "manual", ... -## "PlotBoxAspectRatio", [1 1 1]); -## isonormals (x, y, z, val, p) +## pbaspect ([1 1 1]); +## isonormals (x, y, z, v, p) ## set (p, "FaceLighting", "gouraud"); ## light ("Position", [1 1 5]); ## @end smallexample @@ -143,29 +145,31 @@ function varargout = isosurface (varargin) - if (nargin < 1 || nargin > 8 || nargout > 3) + if (nargin < 1 || nargin > 8) print_usage (); endif - [x, y, z, val, iso, colors, noshare, verbose] = __get_check_isosurface_args__ (nargout, varargin{:}); + [x, y, z, v, isoval, colors, noshare, verbose] = ... + __get_check_isosurface_args__ (nargout, varargin{:}); calc_colors = ! isempty (colors); if (calc_colors) [fvc.faces, fvc.vertices, fvc.facevertexcdata] = ... - __marching_cube__ (x, y, z, val, iso, colors); + __marching_cube__ (x, y, z, v, isoval, colors); else - [fvc.faces, fvc.vertices] = __marching_cube__ (x, y, z, val, iso); + [fvc.faces, fvc.vertices] = __marching_cube__ (x, y, z, v, isoval); endif if (isempty (fvc.vertices) || isempty (fvc.faces)) warning ("isosurface: triangulation is empty"); endif - if (!noshare) - [fvc.faces, fvc.vertices, J] = __unite_shared_vertices__ (fvc.faces, fvc.vertices); + if (! noshare) + [fvc.faces, fvc.vertices, J] = __unite_shared_vertices__ (fvc.faces, + fvc.vertices); if (calc_colors) - fvc.facevertexcdata(J) = []; # share very close vertices + fvc.facevertexcdata(J) = []; # share very close vertices endif endif @@ -175,42 +179,47 @@ if (calc_colors) fc = fvc.facevertexcdata; else - fc = iso; + fc = isoval; endif - ## FIXME: Matlab uses "EdgeColor", "none". But that would look odd - ## with gnuplot. + ## Matlab uses "EdgeColor", "none", but that looks odd in gnuplot. + hax = gca (); + if (strcmp (get (gcf, "__graphics_toolkit__"), "gnuplot")) + ec = "k"; + else + ec = "none"; + endif pa = patch ("Faces", fvc.faces, "Vertices", fvc.vertices, "FaceVertexCData", fc, - "FaceColor", "flat", "EdgeColor", "k", + "FaceColor", "flat", "EdgeColor", ec, "FaceLighting", "gouraud"); - hax = gca (); if (! ishold ()) set (hax, "View", [-37.5, 30], "Box", "off"); endif - isonormals (x, y, z, val, pa); + isonormals (x, y, z, v, pa); lights = findobj (hax, "Type", "light"); if (isempty (lights)) ## FIXME: Matlab seems to use camlight (patch #9014) here light (); endif + case 1 varargout = {fvc}; + case 2 varargout = {fvc.faces, fvc.vertices}; - otherwise ## 3 + + otherwise # 3 args or more varargout = {fvc.faces, fvc.vertices, fvc.facevertexcdata}; + endswitch endfunction -function [x, y, z, data, iso, colors, noshare, verbose] = __get_check_isosurface_args__ (nout, varargin) -## get arguments from input and check values - - x = []; - y = []; - z = []; - data = []; - iso = []; +function [x, y, z, v, isoval, colors, noshare, verbose] = __get_check_isosurface_args__ (nout, varargin) + ## get arguments from input and check values + x = y = z = []; + v = []; + isoval = []; colors = []; ## default values @@ -219,121 +228,132 @@ nin = length (varargin); num_string_inputs = 0; + + ## check whether last 2 input arguments are strings and assign parameters for i_arg = (nin:-1:nin-1) - ## check whether last maximum 2 input arguments are strings and assign parameters - if (!ischar (varargin{i_arg}) || i_arg < 1) - ## string arguments must be at the end - break - end + if (! ischar (varargin{i_arg}) || i_arg < 1) + break; # no string arguments at end, exit checking + endif switch (tolower (varargin{i_arg})) case {"v", "verbose"} verbose = true; num_string_inputs++; + case {"n", "noshare"} noshare = true; num_string_inputs++; + otherwise - error ("isosurface: parameter '%s' not supported", varargin{i_arg}) + error ("isosurface: parameter '%s' not supported", varargin{i_arg}); + endswitch endfor ## assign arguments switch (nin - num_string_inputs) - case 1 ## isosurface (val, ...) - data = varargin{1}; - case 2 ## isosurface (val, iso, ...) or isosurface (val, col, ...) - data = varargin{1}; + case 1 # isosurface (v, ...) + v = varargin{1}; + + case 2 # isosurface (v, isoval, ...) or isosurface (v, col, ...) + v = varargin{1}; if (isscalar (varargin{2}) || isempty (varargin{2})) - iso = varargin{2}; + isoval = varargin{2}; else colors = varargin{2}; endif - case 3 ## isosurface (val, iso, col, ...) - data = varargin{1}; - iso = varargin{2}; + + case 3 # isosurface (v, isoval, col, ...) + v = varargin{1}; + isoval = varargin{2}; colors = varargin{3}; - case 4 ## isosurface (x, y, z, val, ...) + + case 4 # isosurface (x, y, z, v, ...) x = varargin{1}; y = varargin{2}; z = varargin{3}; - data = varargin{4}; - case 5 ## isosurface (x, y, z, val, iso, ...) or isosurface (x, y, z, val, col, ...) + v = varargin{4}; + + case 5 # isosurface (x, y, z, v, isoval, ...) or + # isosurface (x, y, z, v, col, ...) x = varargin{1}; y = varargin{2}; z = varargin{3}; - data = varargin{4}; + v = varargin{4}; if (isscalar (varargin{5}) || isempty (varargin{5})) - iso = varargin{5}; + isoval = varargin{5}; else colors = varargin{5}; endif - case 6 ## isosurface (x, y, z, val, iso, col, ...) + + case 6 # isosurface (x, y, z, v, isoval, col, ...) x = varargin{1}; y = varargin{2}; z = varargin{3}; - data = varargin{4}; - iso = varargin{5}; + v = varargin{4}; + isoval = varargin{5}; colors = varargin{6}; + otherwise - error ("isosurface: wrong number of input arguments") + error ("isosurface: incorrect number of input arguments") + endswitch - ## check dimensions of data - data_size = size (data); - if (ndims (data) != 3 || any (data_size(1:3) < 2)) - error ("isosurface: VAL must be a non-singleton 3-dimensional matrix"); + ## check dimensions of v + v_sz = size (v); + if (ndims (v) != 3 || any (v_sz(1:3) < 2)) + error ("isosurface: V must be a non-singleton 3-dimensional matrix"); endif if (isempty (x)) - x = 1:size (data, 2); + x = 1:size (v, 2); endif if (isempty (y)) - y = 1:size (data, 1); + y = 1:size (v, 1); endif if (isempty (z)) - z = 1:size (data, 3); + z = 1:size (v, 3); endif ## check x - if (isvector (x) && length (x) == data_size(2)) - x = repmat (x(:)', [data_size(1) 1 data_size(3)]); - elseif (! size_equal (data, x)) - error ("isosurface: X must match the size of VAL"); + if (isvector (x) && length (x) == v_sz(2)) + x = repmat (x(:)', [v_sz(1) 1 v_sz(3)]); + elseif (! size_equal (v, x)) + error ("isosurface: X must match the size of V"); endif ## check y - if (isvector (y) && length (y) == data_size(1)) - y = repmat (y(:), [1 data_size(2) data_size(3)]); - elseif (! size_equal (data, y)) - error ("isosurface: Y must match the size of VAL"); + if (isvector (y) && length (y) == v_sz(1)) + y = repmat (y(:), [1 v_sz(2) v_sz(3)]); + elseif (! size_equal (v, y)) + error ("isosurface: Y must match the size of V"); endif ## check z - if (isvector (z) && length (z) == data_size(3)) - z = repmat (reshape (z(:), [1 1 length(z)]), [data_size(1) data_size(2) 1]); - elseif (! size_equal (data, z)) - error ("isosurface: Z must match the size of VAL"); + if (isvector (z) && length (z) == v_sz(3)) + z = repmat (reshape (z(:), [1 1 length(z)]), [v_sz(1) v_sz(2) 1]); + elseif (! size_equal (v, z)) + error ("isosurface: Z must match the size of V"); endif - ## check iso - if (isempty (iso)) - ## calculate "good" iso value from data - iso = __calc_isovalue_from_data__ (data); + ## check isoval + if (isempty (isoval)) + ## calculate "good" isoval value from v + isoval = __calc_isovalue_from_data__ (v); endif - if ~isscalar (iso) - error ("isosurface: ISO must be a scalar") + if (! isscalar (isoval)) + error ("isosurface: ISOVAL must be a scalar") endif ## check colors if (! isempty (colors)) - if (! size_equal (data, colors)) - error ("isosurface: COL must match the size of VAL") + if (! size_equal (v, colors)) + error ("isosurface: COL must match the size of V") endif if (nout == 2) warning ("isosurface: colors will be calculated, but no output argument to receive it."); endif - elseif (nout == 3) + elseif (nout >= 3) error ("isosurface: COL must be passed to return C") endif @@ -343,30 +363,30 @@ %!demo %! clf; %! [x,y,z] = meshgrid (-2:0.5:2, -2:0.5:2, -2:0.5:2); -%! val = x.^2 + y.^2 + z.^2; -%! isosurface (x, y, z, val, 3); -%! isosurface (x, y, z, val, 5); +%! v = x.^2 + y.^2 + z.^2; +%! isosurface (x, y, z, v, 3); +%! isosurface (x, y, z, v, 5); %! axis equal; -%! title ('isosurfaces of two nested spheres'); +%! title ("isosurfaces of two nested spheres"); %!demo %! x = 0:2; %! y = 0:3; %! z = 0:1; -%! [xx, yy, zz] = meshgrid (x, y, z); -%! val = [0, 0, 0; 0, 0, 0; 0, 0, 1; 0, 0, 1]; -%! val(:,:,2) = [0, 0, 0; 0, 0, 1; 0, 1, 2; 0, 1, 2]; +%! [xx, yy, zz] = meshgrid (x, y, z); +%! v = [0, 0, 0; 0, 0, 0; 0, 0, 1; 0, 0, 1]; +%! v(:,:,2) = [0, 0, 0; 0, 0, 1; 0, 1, 2; 0, 1, 2]; %! iso = 0.8; %% three arguments, no output %! figure; -%! subplot (2, 2, 1); isosurface (val, iso, yy); view(3) +%! subplot (2, 2, 1); isosurface (v, iso, yy); view(3) %% six arguments, no output (x, y, z are vectors) -%! subplot (2, 2, 2); isosurface (x, y, z, val, iso, yy); view (3) +%! subplot (2, 2, 2); isosurface (x, y, z, v, iso, yy); view (3) %% six arguments, no output (x, y, z are matrices) -%! subplot (2, 2, 3); isosurface (xx, yy, zz, val, iso, yy); view (3) +%! subplot (2, 2, 3); isosurface (xx, yy, zz, v, iso, yy); view (3) %% six arguments, no output (mixed x, y, z) and option "noshare" -%! subplot (2, 2, 4); isosurface (x, yy, z, val, iso, yy, "noshare"); view (3) -%! annotation("textbox", [0 0.95 1 0.1], ... +%! subplot (2, 2, 4); isosurface (x, yy, z, v, iso, yy, "noshare"); view (3) +%! annotation ("textbox", [0.01 0.93 1 0.1], ... %! "String", ["Apart from the first plot having a different scale, " ... %! "all four plots must look the same.\n" ... %! "The last plot might have different colors but must have " ... @@ -496,25 +516,46 @@ %! assert (size (fvc.facevertexcdata), [7 1]); ## test for each error and warning -%!test -%!error x = 1:2:24; fvc = isosurface (x, y, z, val, iso); -%!error y = -14:6:11; fvc = isosurface (x, y, z, val, iso); -%!error z = linspace (16, 18, 5); fvc = isosurface (x, y, z, val, iso); -%!error x = 1:2:24; [xx, yy, zz] = meshgrid (x, y, z); fvc = isosurface (xx, yy, zz, val, iso); -%!error y = -14:6:11; [xx, yy, zz] = meshgrid (x, y, z); fvc = isosurface (xx, yy, zz, val, iso); -%!error z = linspace (16, 18, 3); [xx, yy, zz] = meshgrid (x, y, z); fvc = isosurface (xx, yy, zz, val, iso); -%!error val = reshape(1:6*8, [6 8]); fvc = isosurface (val, iso); -%!error val = reshape(1:6*8, [6 1 8]); fvc = isosurface (val, iso); -%!error fvc = isosurface (val, [iso iso], yy); -%!error fvc = isosurface (val, [iso iso]); -%!warning [f, v] = isosurface (val, iso, yy); -%!error [f, v, c] = isosurface (val, iso); -%!error fvc = isosurface (); -%!error [f, v, c, a] = isosurface (val, iso); -%!error fvc = isosurface (xx, yy, zz, val, iso, yy, 5); -%!error fvc = isosurface (val, iso, "test_string"); +%!error isosurface () +%!error isosurface (1,2,3,4,5,6,7,8,9) +%!error +%! fvc = isosurface (val, iso, "foobar"); +%!error +%! fvc = isosurface (xx, yy, zz, val, iso, yy, 5); +%!error +%! v = reshape (1:6*8, [6 8]); +%! fvc = isosurface (v, iso); +%!error +%! v = reshape(1:6*8, [6 1 8]); fvc = isosurface (v, iso); +%!error +%! x = 1:2:24; +%! fvc = isosurface (x, y, z, val, iso); +%!error +%! y = -14:6:11; +%! fvc = isosurface (x, y, z, val, iso); +%!error +%! z = linspace (16, 18, 5); +%! fvc = isosurface (x, y, z, val, iso); +%!error +%! x = 1:2:24; +%! [xx, yy, zz] = meshgrid (x, y, z); +%! fvc = isosurface (xx, yy, zz, val, iso); +%!error +%! y = -14:6:11; +%! [xx, yy, zz] = meshgrid (x, y, z); +%! fvc = isosurface (xx, yy, zz, val, iso); +%!error +%! z = linspace (16, 18, 3); +%! [xx, yy, zz] = meshgrid (x, y, z); +%! fvc = isosurface (xx, yy, zz, val, iso); +%!error fvc = isosurface (val, [iso iso], yy) +%!error fvc = isosurface (val, [iso iso]); +%!error [f, v, c] = isosurface (val, iso) +%!warning +%! [f, v] = isosurface (val, iso, yy); ## test for __calc_isovalue_from_data__ ## FIXME: private function cannot be tested, unless bug #38776 is resolved. %!xtest %! assert (__calc_isovalue_from_data__ (1:5), 3.02); + diff -r 01ba6ebc52e4 -r a8fd02bc895b scripts/plot/draw/private/__calc_isovalue_from_data__.m --- a/scripts/plot/draw/private/__calc_isovalue_from_data__.m Mon Jul 18 19:49:39 2016 +0200 +++ b/scripts/plot/draw/private/__calc_isovalue_from_data__.m Tue Aug 09 11:09:38 2016 -0700 @@ -19,23 +19,21 @@ ## Undocumented internal function. ## -*- texinfo -*- -## @deftypefn {Function File} {@var{iso} =} __calc_isovalue_from_data__ (@var{data}) -## Undocumented internal function. +## @deftypefn {} {@var{isoval} =} __calc_isovalue_from_data__ (@var{data}) +## Calculate a "good" iso value from histogram of data. ## @end deftypefn -## calculate a "good" iso value from histogram of data ## called from isocaps, isosurface +function isoval = __calc_isovalue_from_data__ (data) -function iso = __calc_isovalue_from_data__ (data) - - ## use a maximum of 10000-20000 samples to limit run-time of hist + ## use a maximum of 10,000-20,000 samples to limit runtime of hist step = 1; - data_numel = numel (data); - if (data_numel > 20000) - step = floor (data_numel / 10000); + ndata = numel (data); + if (ndata > 20_000) + step = floor (ndata / 10_000); data = data(1:step:end); - data_numel = numel (data); + ndata = numel (data); endif num_bins = 100; @@ -43,18 +41,20 @@ ## if one of the first two bins contains more than 10 times the count as ## compared to equally distributed data, remove both (zero-padded + noise) - if (any (bin_count(1:2) > 10 * (data_numel / num_bins))) + if (any (bin_count(1:2) > 10 * (ndata / num_bins))) bin_count(1:2) = []; bin_centers(1:2) = []; endif - ## if bins have low counts, remove them (but keep them if we would loose more than 90 % of bins) + ## if bins have low counts, remove them (but keep them if we would lose + ## more than 90% of bins) bins_to_remove = find (bin_count < max (bin_count)/50); - if length (bins_to_remove) < .9 * num_bins + if (length (bins_to_remove) < .9 * num_bins) bin_centers(bins_to_remove) = []; endif ## select middle bar of histogram with previous conditions - iso = bin_centers(floor (length (bin_centers) / 2)); + isoval = bin_centers(floor (numel (bin_centers) / 2)); endfunction + diff -r 01ba6ebc52e4 -r a8fd02bc895b scripts/plot/draw/private/__unite_shared_vertices__.m --- a/scripts/plot/draw/private/__unite_shared_vertices__.m Mon Jul 18 19:49:39 2016 +0200 +++ b/scripts/plot/draw/private/__unite_shared_vertices__.m Tue Aug 09 11:09:38 2016 -0700 @@ -19,12 +19,12 @@ ## -*- texinfo -*- ## @deftypefn {} {[@var{faces}, @var{vertices}, @var{J}] =} __unite_shared_vertices__ (@var{faces}, @var{vertices}) ## -## Detect and unite shared vertices in patches +## Detect and unite shared vertices in patches. ## -## Vertices of neighboring faces are detected and united to shared vertices. For -## this, the mutual squared distances between all vertices are calculated. If -## all coordinates are closer than @command{2 * eps (max (abs (vertices(:))))}, -## the vertices are united to one. +## Vertices of neighboring faces are detected and united to shared vertices. +## For this, the mutual squared distances between all vertices are +## calculated. If all coordinates are closer than +## @code{2 * eps (max (abs (vertices(:))))}, the vertices are united to one. ## ## @var{J} holds the indices of the deleted vertices. ## @@ -34,26 +34,31 @@ ## Author: mmuetzel function [faces, vertices, J] = __unite_shared_vertices__ (faces, vertices) - ### unite shared vertices J = []; + ## Calculate the mutual differences of all vertex coordinates close_points = zeros (0, 2); - num_vertices = size (vertices, 1); + num_vertices = rows (vertices); skip_point = false (num_vertices, 1); + ## FIXME: Can this be vectorized in some way to increase performance? + ## Regardless, should probably allocate close_points to be the + ## same size as the number of vertices and then truncate the + ## array at the end of the calculation. Extending an array + ## involves a copy operation every time. for (i_point1 = 1:num_vertices - 1) if (skip_point(i_point1)) ## points already detected as duplicates can be skipped - continue + continue; endif diff = vertices(i_point1,:) - vertices(i_point1 + 1:end,:); - is_close_point = all (abs (diff) <= sqrt(3) * eps * ... + is_close_point = all (abs (diff) <= sqrt (3) * eps * ... (max (abs (vertices(i_point1,:)), abs (vertices(i_point1 + 1:end,:)))), 2); if (any (is_close_point)) close_points_idx = find (is_close_point) + i_point1; - new_close_points_num = size (close_points_idx, 1); + new_close_points_num = rows (close_points_idx); close_points(end + 1:end + new_close_points_num,1) = i_point1; close_points(end - new_close_points_num + 1:end,2) = close_points_idx; skip_point(close_points_idx) = true; @@ -76,11 +81,12 @@ faces = sort (faces, 2); faces = unique (faces, "rows"); - ## eliminate faces with zero area - is_zero_area = (faces(:,1) == faces(:,2)) | (faces(:,2) == faces(:,3)); # vertices in faces are sorted - faces = faces(!is_zero_area,:); + ## eliminate faces with zero area. Vertices in faces are sorted. + is_zero_area = (faces(:,1) == faces(:,2)) | (faces(:,2) == faces(:,3)); + faces = faces(! is_zero_area, :); J = close_points(:,2); endif endfunction +