## 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");