## Copyright (C) 2016 Markus Muetzel ## ## This file is part of Octave. ## ## Octave is free software; you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by ## the Free Software Foundation; either version 3 of the License, or (at ## your option) any later version. ## ## Octave is distributed in the hope that it will be useful, but ## WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ## General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with Octave; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {@var{fvc_struct} =} isocaps (@var{val}, @var{iso}) ## @deftypefnx {Function File} {@var{fvc_struct} =} isocaps (@var{val}) ## @deftypefnx {Function File} {@var{fvc_struct} =} isocaps (@var{x}, @var{y}, @var{z}, @var{val}, @var{iso}) ## @deftypefnx {Function File} {@var{fvc_struct} =} isocaps (@var{x}, @var{y}, @var{z}, @var{val}) ## @deftypefnx {Function File} {@var{fvc_struct} =} isocaps (@dots{}, @var{which_caps}) ## @deftypefnx {Function File} {@var{fvc_struct} =} isocaps (@dots{}, @var{which_plane}) ## @deftypefnx {Function File} {@var{fvc_struct} =} isocaps (@dots{}, @qcode{"verbose"}) ## @deftypefnx {Function File} {[@var{faces}, @var{vertices}, @var{cdata}] =} isocaps (@dots{}) ## @deftypefnx {Function File} {} isocaps (@dots{}) ## ## Create end-caps for isosurfaces of 3d data ## ## @var{val} must be a 3-dimensional matrix. The end-caps for this data are ## calculated for the isosurface at value @var{iso}. ## ## If @var{iso} is omitted or empty, a "good" value for an isosurface is ## determined from @var{val}. ## ## If called with one output argument the results are returned in a structure ## with the fields: @code{vertices}, @code{faces} and @code{facevertexcdata}. ## ## Optionally, @var{x}, @var{y} and @var{z} can be supplied to represent the ## set of coordinates of @var{val}. They can either be matrices of the same size ## as @var{val} or vectors with sizes according to the dimensions of @var{val} ## in which case they are expanded to matrices (see: @command{meshgrid}). ## ## If @command{isocaps} is called without passing @var{x}, @var{y} and @var{z} ## or they are empty, they are assumed to match the respective indices of ## @var{val}. ## ## The optional parameter @var{which_caps} can have one of the following string ## values to define how the data should be enclosed: ## @table @asis ## @item @code{above}, @code{a} (default) ## for end-caps that enclose the data above @var{iso}. ## @item @code{below}, @code{b} ## for end-caps that enclose the data below @var{iso}. ## @end table ## ## The optional parameter @var{which_plane} can have one of the following string ## values to define in which plane the cap should be drawn: ## @table @asis ## @item @code{xmax} ## for the end-caps at the upper x-plain of the data. ## @item @code{xmin} ## for the end-caps at the lower x-plain of the data. ## @item @code{ymax} ## for the end-caps at the upper y-plain of the data. ## @item @code{ymin} ## for the end-caps at the lower y-plain of the data. ## @item @code{zmax} ## for the end-caps at the upper z-plain of the data. ## @item @code{zmin} ## for the end-caps at the lower z-plain of the data. ## @item @code{all} (default) ## for all of these end-caps. ## @end table ## ## The optional parameter @qcode{"verbose"} is currently ignored. ## ## If called with two or three output arguments, the data for faces @var{faces}, ## vertices @var{vertices} and the color data @var{cdata} are returned in ## separate arrays instead of a single structure. ## ## If called with no output arguments, the results are drawn in the current ## figure with the @command{patch} command. ## ## Example: ## @example ## @group ## data = rand (10, 10, 10); ## smoothed_data = smooth3 (data); ## patch (isosurface (smoothed_data, .5), "FaceColor", "blue", "EdgeColor", "k"); ## isocaps (smoothed_data, .5); ## view(3) ## @end group ## @end example ## ## @seealso{isosurface, isonormals, patch, reducevolume, smooth3} ## @end deftypefn ## Author: mmuetzel function varargout = isocaps (varargin) if (nargin < 1 || nargin > 8 || nargout > 3) print_usage (); endif faces = vertices = caps_cdata = []; [x, y, z, data, iso, which_caps, which_plane, verbose] = __get_check_isocaps_args__ (varargin{:}); ## select type of cap (above or below iso value) data_min = min([data(:); iso]); data_max = max([data(:); iso]); switch (tolower (which_caps)) case {"a", "above"} pad_val = data_min - 1; otherwise ## "b", "below" pad_val = data_max + 1; endswitch ## create patches for caps if (strcmpi (which_plane, "all")) ## get patches for all planes [faces, vertices] = __get_isocaps_patches__ (data, iso, pad_val, "xmin", verbose); [f, v] = __get_isocaps_patches__ (data, iso, pad_val, "xmax", verbose); faces = [faces; f + size(vertices, 1)]; vertices = [vertices; v]; [f, v] = __get_isocaps_patches__ (data, iso, pad_val, "ymin", verbose); faces = [faces; f + size(vertices, 1)]; vertices = [vertices; v]; [f, v] = __get_isocaps_patches__ (data, iso, pad_val, "ymax", verbose); faces = [faces; f + size(vertices, 1)]; vertices = [vertices; v]; [f, v] = __get_isocaps_patches__ (data, iso, pad_val, "zmin", verbose); faces = [faces; f + size(vertices, 1)]; vertices = [vertices; v]; [f, v] = __get_isocaps_patches__ (data, iso, pad_val, "zmax", verbose); faces = [faces; f + size(vertices, 1)]; vertices = [vertices; v]; else ## only one plain specified [faces, vertices] = __get_isocaps_patches__ (data, iso, pad_val, which_plane, verbose); endif if (!isempty (vertices)) ## interpolate data at the vertices (on a "meshgrid") caps_cdata = interp3 (data, vertices(:,1), vertices(:,2), vertices(:,3))'; ## check if vertices must be re-calculated on given grid data_size = size(data); if !(isequal (x, 1:data_size(2)) && isequal (y, 1:data_size(1)) && isequal (z, 1:data_size(3))) vertices_grid(:,1) = interp3 (x, vertices(:,1), vertices(:,2), vertices(:,3)); vertices_grid(:,2) = interp3 (y, vertices(:,1), vertices(:,2), vertices(:,3)); vertices_grid(:,3) = interp3 (z, vertices(:,1), vertices(:,2), vertices(:,3)); vertices = vertices_grid; endif endif switch nargout case 0 hp = patch ("Faces", faces, "Vertices", vertices, ... "FaceVertexCData", caps_cdata, ... "FaceColor", "interp", "EdgeColor", "none"); case 1 vfc.vertices = vertices; vfc.faces = faces; vfc.facevertexcdata = caps_cdata; varargout = {vfc}; otherwise varargout{1} = faces; varargout{2} = vertices; varargout{3} = caps_cdata; endswitch endfunction function [x, y, z, data, iso, which_caps, which_plane, verbose] = __get_check_isocaps_args__ (varargin) ## get arguments from input and check values x = []; y = []; z = []; data = []; iso = []; ## default values which_caps = "above"; which_plane = "all"; verbose = ""; num_string_inputs = 0; for i_arg = (nargin:-1:nargin-2) ## check whether last maximum 3 input arguments are strings and assign parameters if (!ischar (varargin{i_arg}) || i_arg < 1) ## string arguments must be at the end break end switch (tolower (varargin{i_arg})) case {"v", "verbose"} verbose = "verbose"; num_string_inputs++; case {"all", "xmin", "xmax", "ymin", "ymax", "zmin", "zmax"} which_plane = varargin{i_arg}; num_string_inputs++; case {"above", "a", "below", "b"} which_caps = varargin{i_arg}; num_string_inputs++; otherwise error ("isocaps: parameter '%s' not supported", varargin{i_arg}) endswitch endfor switch (nargin - num_string_inputs) case 1 ## isocaps (val, ...) data = varargin{1}; case 2 ## isocaps (val, iso, ...) data = varargin{1}; iso = varargin{2}; case 4 ## isocaps (x, y, z, val, ...) x = varargin{1}; y = varargin{2}; z = varargin{3}; data = varargin{4}; case 5 ## isocaps (x, y, z, val, iso, ...) x = varargin{1}; y = varargin{2}; z = varargin{3}; data = varargin{4}; iso = varargin{5}; otherwise error ("isocaps: wrong number of input arguments") endswitch ## check dimensions of data if (ndims (data) < 3) error ("isocaps: data must have at least 3 dimensions"); endif data_size = size (data); if (any (data_size(1:3) < 2)) error ("isocaps: data must be a non-singleton 3-dimensional matrix"); endif if (isempty (x)) x = 1:size (data, 2); endif if (isempty (y)) y = 1:size (data, 1); endif if (isempty (z)) z = 1:size (data, 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 ("isocaps: X must match the size of data"); 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 ("isocaps: Y must match the size of data"); 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 ("isocaps: Z must match the size of data"); endif ## check iso if (isempty (iso)) ## calculate "good" iso value from data iso = __calc_isovalue_from_data__ (data); endif if ~isscalar (iso) error ("isocaps: ISO must be a scalar") endif endfunction function [faces, vertices] = __get_isocaps_patches__ (data, iso, pad_val, which_plane, verbose) ## calculate patches for end-caps data_size = size (data); switch (lower(which_plane)) case "ymin" cap_data = zeros(2, data_size(2), data_size(3)); cap_data(1,:,:) = pad_val; cap_data(2,:,:) = data(1,:,:); [faces, vertices] = isosurface (cap_data, iso);##, verbose); # bug #46946 if (!isempty (vertices)) vertices(:,2) = 1; end case "ymax" cap_data = zeros(2, data_size(2), data_size(3)); cap_data(1,:,:) = data(end,:,:); cap_data(2,:,:) = pad_val; [faces, vertices] = isosurface (cap_data, iso);##, verbose); # bug #46946 if (!isempty (vertices)) vertices(:,2) = data_size(1); end case "xmin" cap_data = zeros(data_size(1), 2, data_size(3)); cap_data(:,1,:) = pad_val; cap_data(:,2,:) = data(:,1,:); [faces, vertices] = isosurface (cap_data, iso);##, verbose); # bug #46946 if (!isempty (vertices)) vertices(:,1) = 1; end case "xmax" cap_data = zeros(data_size(1), 2, data_size(3)); cap_data(:,1,:) = data(:,end,:); cap_data(:,2,:) = pad_val; [faces, vertices] = isosurface (cap_data, iso);##, verbose); # bug #46946 if (!isempty (vertices)) vertices(:,1) = data_size(2); end case "zmin" cap_data = zeros(data_size(1), data_size(2), 2); cap_data(:,:,1) = pad_val; cap_data(:,:,2) = data(:,:,1); [faces, vertices] = isosurface (cap_data, iso);##, verbose); # bug #46946 if (!isempty (vertices)) vertices(:,3) = 1; end otherwise ## "zmax" cap_data = zeros(data_size(1), data_size(2), 2); cap_data(:,:,1) = data(:,:,end); cap_data(:,:,2) = pad_val; [faces, vertices] = isosurface (cap_data, iso);##, verbose); # bug #46946 if (!isempty (vertices)) vertices(:,3) = data_size(3); end endswitch endfunction %!shared x, y, z, xx, yy, zz, val, iso %! val = rand (6, 8, 4); %! iso = .5; %! x = 1:3:22; y = -14:5:11; z = linspace (16, 18, 4); %! [xx, yy, zz] = meshgrid (x, y, z); ## two arguments, no output %!test %! figure; %! subplot (2, 2, 1); isocaps (val, iso); view (3) ## five arguments, no output (x, y, z are vectors) %!test %! subplot (2, 2, 2); isocaps (x, y, z, val, iso); view (3) ## five arguments, no output (x, y, z are matrices) %!test %! subplot (2, 2, 3); isocaps (xx, yy, zz, val, iso); view (3) ## five arguments, no output (mixed x, y, z) %!test %! subplot (2, 2, 4); isocaps (x, yy, z, val, iso); view (3) %! annotation("textbox", [0 0.95 1 0.1], ... %! "String", "Besides the first plot having a different scale, all four plots must look the same.", ... %! "HorizontalAlignment", "center"); ## check results for differently shaped input coordinates %!test %! fvc_vectors = isocaps (x, y, z, val, iso); %! fvc_matrices = isocaps (xx, yy, zz, val, iso); %! fvc_mixed = isocaps (xx, y, zz, val, iso); %! assert (fvc_vectors, fvc_matrices); %! assert (fvc_vectors, fvc_mixed); ## two arguments, one output %!test %! fvc = isocaps (val, iso); %! assert (isfield (fvc, "vertices"), true); %! assert (isfield (fvc, "faces"), true); %! assert (isfield (fvc, "facevertexcdata"), true); ## one argument(un-documented Matlab) %!test %! fvc = isocaps (val); %! fvc2 = isocaps (val, []); %! assert (fvc, fvc2); %! assert (isfield (fvc, "vertices"), true); %! assert (isfield (fvc, "faces"), true); %! assert (isfield (fvc, "facevertexcdata"), true); ## four arguments (un-documented Matlab) %!test %! fvc = isocaps (x, [], z, val); %! assert (isfield (fvc, "vertices"), true); %! assert (isfield (fvc, "faces"), true); %! assert (isfield (fvc, "facevertexcdata"), true); ## five arguments, two outputs %!test %! [faces, vertices] = isocaps ([], y, z, val, iso); %! assert (size (faces, 2), 3); %! assert (size (vertices, 2), 3); ## five arguments, three outputs %!test %! [faces, vertices, caps_cdata] = isocaps (x, y, [], val, iso); %! assert (size (faces, 2), 3); %! assert (size (vertices, 2), 3); %! assert (size (caps_cdata, 2), 1); %! assert (size (vertices, 1), size (caps_cdata, 1)); ## two arguments + one string, one output %!test %! fvc = isocaps (val, iso, "below"); %! assert (isfield (fvc, "vertices"), true); %! assert (isfield (fvc, "faces"), true); %! assert (isfield (fvc, "facevertexcdata"), true); ## two arguments + two strings, one output %!test %! fvc = isocaps (val, iso, "b", "ymax"); %! assert (isfield (fvc, "vertices"), true); %! assert (isfield (fvc, "faces"), true); %! assert (isfield (fvc, "facevertexcdata"), true); ## two arguments + three strings, one output %!test %! fvc = isocaps (val, iso, "a", "ymin", "verbose"); %! assert (isfield (fvc, "vertices"), true); %! assert (isfield (fvc, "faces"), true); %! assert (isfield (fvc, "facevertexcdata"), true); ## five arguments + one string, three outputs %!test %! [faces, vertices, caps_cdata] = isocaps (x, y, z, val, iso, "xmin"); %! assert (size (faces, 2), 3); %! assert (size (vertices, 2), 3); %! assert (size (caps_cdata, 2), 1); %! assert (size (vertices, 1), size (caps_cdata, 1)); ## five arguments + one string, three outputs %!test %! [faces, vertices, caps_cdata] = isocaps (x, y, z, val, iso, "verbose"); ## un-documented Matlab %! assert (size (faces, 2), 3); %! assert (size (vertices, 2), 3); %! assert (size (caps_cdata, 2), 1); %! assert (size (vertices, 1), size (caps_cdata, 1)); ## five arguments + two strings, three outputs %!test %! [faces, vertices, caps_cdata] = isocaps (x, y, z, val, iso, "xmax", "above"); %! assert (size (faces, 2), 3); %! assert (size (vertices, 2), 3); %! assert (size (caps_cdata, 2), 1); %! assert (size (vertices, 1), size (caps_cdata, 1)); ## five arguments + three strings, three outputs %!test %! [faces, vertices, caps_cdata] = isocaps (x, y, z, val, iso, "zmin", "b", "verbose"); %! assert (size (faces, 2), 3); %! assert (size (vertices, 2), 3); %! assert (size (caps_cdata, 2), 1); %! assert (size (vertices, 1), size (caps_cdata, 1)); ## five arguments + three strings (different order), three outputs %!test %! [faces, vertices, caps_cdata] = isocaps (x, y, z, val, iso, "below", "v", "zmax"); %! assert (size (faces, 2), 3); %! assert (size (vertices, 2), 3); %! assert (size (caps_cdata, 2), 1); %! assert (size (vertices, 1), size (caps_cdata, 1)); ## test for each error %!test %!error x = 1:2:24; fvc = isocaps (x, y, z, val, iso); %!error y = -14:6:11; fvc = isocaps (x, y, z, val, iso); %!error z = linspace (16, 18, 5); fvc = isocaps (x, y, z, val, iso); %!error x = 1:2:24; [xx, yy, zz] = meshgrid (x, y, z); fvc = isocaps (xx, yy, zz, val, iso); %!error y = -14:6:11; [xx, yy, zz] = meshgrid (x, y, z); fvc = isocaps (xx, yy, zz, val, iso); %!error z = linspace (16, 18, 3); [xx, yy, zz] = meshgrid (x, y, z); fvc = isocaps (xx, yy, zz, val, iso); %!error v2 = reshape(1:6*8, [6 8]); fvc = isocaps (v2, iso); %!error v3 = reshape(1:6*8, [6 1 8]); fvc = isocaps (v3, iso); %!error fvc = isocaps (val, [iso iso]); %!error fvc = isocaps (); %!error fvc = isocaps (xx, yy, zz, val, iso, 5); %!error fvc = isocaps (x, val, iso); %!error fvc = isocaps (val, iso, "test_string");