## 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{nx}, @var{ny}, @var{nz}, @var{nv}] =} reducevolume (@var{v}, @var{r})
## @deftypefnx {Function File} {@var{nv} =} reducevolume (@dots{})
## @deftypefnx {Function File} {[@var{nx}, @var{ny}, @var{nz}, @var{nv}] =} reducevolume (@var{x}, @var{y}, @var{z}, @var{v}, @var{r})
## @deftypefnx {Function File} {@var{nv} =} reducevolume (@dots{})
##
## Reduce the volume of the dataset in @var{v} according to the values in
## @var{r}.
##
## @var{v} is a matrix that is non-singleton in the first 3 dimensions.
##
## @var{r} can be either a vector of 3 elements representing the reduction
## factors in x-, y- and z-direction or a scalar in which case the same
## reduction factor is used in all three dimensions.
##
## @code{reducevolume} reduces the number of elements of @var{v} by taking only
## every @var{r}-th element in the respective dimension.
##
## Optionally, @var{x}, @var{y} and @var{z} can be supplied to represent the
## set of coordinates of @var{v}. They can either be matrices of the same size
## as @var{v} or vectors with sizes according to the dimensions of @var{v} in
## which case they are expanded to matrices (see: @command{meshgrid}).
##
## The reduced matrix is returned in @var{nv}.
##
## Optionally, the reduced set of coordinates are returned in @var{nx}, @var{ny}
## and @var{nz}, respectively.
##
## If @command{reducevolume} is called with two arguments, @var{x}, @var{y} and
## @var{z} are assumed to match the respective indices of @var{v}.
##
## Examples:
## @example
## @group
## @var{v} = reshape (1:6*8*4, [6 8 4]);
## [@var{nv}] = reducevolume (@var{v}, [4 3 2]);
## @end group
## @end example
##
## @example
## @group
## @var{v} = reshape (1:6*8*4, [6 8 4]);
## @var{x} = 1:3:24; @var{y} = -14:5:11; @var{z} = linspace (16, 18, 4);
## [@var{nx}, @var{ny}, @var{nz}, @var{nv}] = reducevolume (@var{x}, @var{y}, @var{z}, @var{v}, [4 3 2]);
## @end group
## @end example
##
## @seealso{isosurface, isonormals, isocaps, smooth3}
## @end deftypefn
## Author: mmuetzel
function [reduced_x, reduced_y, reduced_z, reduced_data] = reducevolume (varargin)
if (nargin < 2 || nargin > 5 || (nargout != 1 && nargout != 4))
print_usage ();
endif
[x, y, z, data, r] = __get_check_reducevolume_args__ (nargout, varargin{:});
[reduced_x, reduced_y, reduced_z, reduced_data] = __reducevolume__ (x, y, z, data, r);
if (nargout == 1)
reduced_x = reduced_data;
endif
endfunction
function [x, y, z, data, r] = __get_check_reducevolume_args__ (naout, varargin)
x = [];
y = [];
z = [];
switch (nargin)
case 3
data = varargin{1};
r = varargin{2};
case 6
if (naout == 4)
x = varargin{1};
y = varargin{2};
z = varargin{3};
endif
data = varargin{4};
r = varargin{5};
otherwise
error ("reducevolume: wrong number of arguments");
endswitch
## check reduction values r
if (isscalar (r))
r = [r r r];
endif
if (numel (r) != 3)
error ("reducevolume: reduction value R must either be a scalar or a vector of length 3");
endif
if (any (r(:) < 1 | r(:) != fix (r(:))))
error ("reducevolume: reduction values R must be positive integers");
endif
## check dimensions of data
if (ndims (data) < 3)
error ("reducevolume: data must have at least 3 dimensions");
endif
data_size = size (data);
if (any (data_size(1:3) < 2))
error ("reducevolume: data must be a non-singleton 3-dimensional matrix");
endif
if (naout == 4)
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 ("reducevolume: 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 ("reducevolume: 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 ("reducevolume: Z must match the size of data");
endif
endif
endfunction
function [reduced_x, reduced_y, reduced_z, reduced_data] = __reducevolume__ (x, y, z, data, r)
data_size = size (data);
reduced_data = data(1:r(2):end, 1:r(1):end, 1:r(3):end, :);
reduced_data_size = size (reduced_data);
if (length (reduced_data_size) < 3 || min (reduced_data_size) < 2)
error ("reducevolume: reduction value R is too high");
endif
if (length (data_size) > 3)
reduced_data = reshape (reduced_data, [reduced_data_size(1:3) data_size(4:end)]);
endif
if (isempty (x))
[reduced_x, reduced_y, reduced_z] = deal ([]);
else
reduced_x = x(1:r(2):end, 1:r(1):end, 1:r(3):end);
reduced_y = y(1:r(2):end, 1:r(1):end, 1:r(3):end);
reduced_z = z(1:r(2):end, 1:r(1):end, 1:r(3):end);
endif
endfunction
%!shared v, x, y, z, xx, yy, zz
%! v = reshape (1:6*8*4, [6 8 4]);
%! x = 1:3:22; y = -14:5:11; z = linspace (16, 18, 4);
%! [xx, yy, zz] = meshgrid (x, y, z);
## two inputs, one output
%!test
%! [nv] = reducevolume (v, [4 3 2]);
%! nv_expected = [1 25; 4 28]; nv_expected(:,:,2) = [97 121; 100 124];
%! assert (nv, nv_expected);
## two inputs, four outputs
%!test
%! [nx, ny, nz, nv] = reducevolume (v, [4 3 2]);
%! nx_expected(1:2,1,1:2) = 1; nx_expected(:,2,:) = 5;
%! ny_expected(1,1:2,1:2) = 1; ny_expected(2,:,:) = 4;
%! nz_expected(1:2,1:2,1) = 1; nz_expected(:,:,2) = 3;
%! nv_expected = [1 25; 4 28]; nv_expected(:,:,2) = [97 121; 100 124];
%! assert (nx, nx_expected);
%! assert (ny, ny_expected);
%! assert (nz, nz_expected);
%! assert (nv, nv_expected);
## five inputs, one output
%!test
%! [nv] = reducevolume (x, y, z, v, [4 3 2]);
%! nv_expected = [1 25; 4 28]; nv_expected(:,:,2) = [97 121; 100 124];
%! assert (nv, nv_expected);
## five inputs, four outputs (coordinates are vectors)
%!test
%! [nx, ny, nz, nv] = reducevolume (x, y, z, v, [4 3 2]);
%! nx_expected(1:2,1,1:2) = x(1); nx_expected(:,2,:) = x(5);
%! ny_expected(1,1:2,1:2) = y(1); ny_expected(2,:,:) = y(4);
%! nz_expected(1:2,1:2,1) = z(1); nz_expected(:,:,2) = z(3);
%! nv_expected = [1 25; 4 28]; nv_expected(:,:,2) = [97 121; 100 124];
%! assert (nx, nx_expected);
%! assert (ny, ny_expected);
%! assert (nz, nz_expected);
%! assert (nv, nv_expected);
## five inputs, four outputs (coordinates are matrices)
%!test
%! [nx, ny, nz, nv] = reducevolume (xx, yy, zz, v, [4 3 2]);
%! nx_expected(1:2,1,1:2) = x(1); nx_expected(:,2,:) = x(5);
%! ny_expected(1,1:2,1:2) = y(1); ny_expected(2,:,:) = y(4);
%! nz_expected(1:2,1:2,1) = z(1); nz_expected(:,:,2) = z(3);
%! nv_expected = [1 25; 4 28]; nv_expected(:,:,2) = [97 121; 100 124];
%! assert (nx, nx_expected);
%! assert (ny, ny_expected);
%! assert (nz, nz_expected);
%! assert (nv, nv_expected);
## five inputs, four outputs (coordinates are matrices, R is scalar)
%!test
%! [nx, ny, nz, nv] = reducevolume (xx, yy, zz, v, 3);
%! nx_expected(1:2,1,1:2) = x(1); nx_expected(:,2,:) = x(4); nx_expected(:,3,:) = x(7);
%! ny_expected(1,1:3,1:2) = y(1); ny_expected(2,:,:) = y(4);
%! nz_expected(1:2,1:3,1) = z(1); nz_expected(:,:,2) = z(4);
%! nv_expected = [1 19 37; 4 22 40]; nv_expected(:,:,2) = [145 163 181; 148 166 184];
%! assert (nx, nx_expected);
%! assert (ny, ny_expected);
%! assert (nz, nz_expected);
%! assert (nv, nv_expected);
## test for each error
%!test
%!error x = 1:2:24; [nx, ny, nz, nv] = reducevolume (x, y, z, v, [4 3 2]);
%!error y = -14:6:11; [nx, ny, nz, nv] = reducevolume (x, y, z, v, [4 3 2]);
%!error z = linspace (16, 18, 5); [nx, ny, nz, nv] = reducevolume (x, y, z, v, [4 3 2]);
%!error x = 1:2:24; [xx, yy, zz] = meshgrid (x, y, z); [nx, ny, nz, nv] = reducevolume (xx, yy, zz, v, [4 3 2]);
%!error y = -14:6:11; [xx, yy, zz] = meshgrid (x, y, z); [nx, ny, nz, nv] = reducevolume (xx, yy, zz, v, [4 3 2]);
%!error z = linspace (16, 18, 3); [xx, yy, zz] = meshgrid (x, y, z); [nx, ny, nz, nv] = reducevolume (xx, yy, zz, v, [4 3 2]);
%!error v = reshape(1:6*8, [6 8]); [nv] = reducevolume (v, [4 3 2]);
%!error v = reshape(1:6*8, [6 1 8]); [nv] = reducevolume (v, [4 3 2]);
%!error [nv] = reducevolume (v, []);
%!error [nv] = reducevolume (v, [1 2]);
%!error [nv] = reducevolume (v, 1.5);
%!error [nv] = reducevolume (v, 0);
%!error [nv] = reducevolume (v);
%!error [nx, ny] = reducevolume (v, 3);
%!error [nv] = reducevolume (x, v, 3);
%!error [nv] = reducevolume (v, 5);
%!error [nv] = reducevolume (v, [4 7 2]);