## Copyright (C) 2016 Markus Mützel
##
## 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
## <http://www.gnu.org/licenses/>.

## -*- 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: @code{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 @code{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
## @end deftypefn

## Author: mm

function [reduced_x, reduced_y, reduced_z, reduced_data] = reducevolume (varargin)

  [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

## two inputs, one output
%!test
%! v = reshape (1:6*8*4, [6 8 4]);
%! [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
%! v = reshape (1:6*8*4, [6 8 4]);
%! [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
%! v = reshape (1:6*8*4, [6 8 4]);
%! x = 1:3:24; y = -14:5:11; z = linspace (16, 18, 4);
%! [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
%! v = reshape (1:6*8*4, [6 8 4]);
%! x = 1:3:24; y = -14:5:11; z = linspace (16, 18, 4);
%! [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
%! v = reshape (1:6*8*4, [6 8 4]);
%! x = 1:3:24; y = -14:5:11; z = linspace (16, 18, 4);
%! [xx, yy, zz] = meshgrid (x, y, z);
%! [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
%! v = reshape(1:6*8*4, [6 8 4]);
%! x = 1:3:24; y = -14:5:11; z = linspace (16, 18, 4);
%! [xx, yy, zz] = meshgrid (x, y, z);
%! [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
%!error <X must match the size of data> v = reshape (1:6*8*4, [6 8 4]); x = 1:2:24; y = -14:5:11; z = linspace (16, 18, 4); [nx, ny, nz, nv] = reducevolume (x, y, z, v, [4 3 2]);
%!error <Y must match the size of data> v = reshape (1:6*8*4, [6 8 4]); x = 1:3:24; y = -14:6:11; z = linspace (16, 18, 4); [nx, ny, nz, nv] = reducevolume (x, y, z, v, [4 3 2]);
%!error <Z must match the size of data> v = reshape (1:6*8*4, [6 8 4]); x = 1:3:24; y = -14:5:11; z = linspace (16, 18, 5); [nx, ny, nz, nv] = reducevolume (x, y, z, v, [4 3 2]);
%!error <X must match the size of data> v = reshape (1:6*8*4, [6 8 4]); x = 1:2:24; y = -14:5:11; z = linspace (16, 18, 4); [xx, yy, zz] = meshgrid (x, y, z); [nx, ny, nz, nv] = reducevolume (xx, yy, zz, v, [4 3 2]);
%!error <X must match the size of data> v = reshape (1:6*8*4, [6 8 4]); x = 1:3:24; y = -14:6:11; z = linspace (16, 18, 4); [xx, yy, zz] = meshgrid (x, y, z); [nx, ny, nz, nv] = reducevolume (xx, yy, zz, v, [4 3 2]);
%!error <X must match the size of data> v = reshape (1:6*8*4, [6 8 4]); x = 1:3:24; y = -14:5:11; 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 <data must have at least 3 dimensions> v = reshape(1:6*8, [6 8]); [nv] = reducevolume (v, [4 3 2]);
%!error <data must be a non-singleton 3-dimensional matrix> v = reshape(1:6*8, [6 1 8]); [nv] = reducevolume (v, [4 3 2]);
%!error <reduction value R must either be a scalar or a vector of length 3> v = reshape(1:6*8*4, [6 8 4]); [nv] = reducevolume (v, []);
%!error <reduction value R must either be a scalar or a vector of length 3> v = reshape(1:6*8*4, [6 8 4]); [nv] = reducevolume (v, [1 2]);
%!error <reduction values R must be positive integers> v = reshape(1:6*8*4, [6 8 4]); [nv] = reducevolume (v, 1.5);
%!error <reduction values R must be positive integers> v = reshape(1:6*8*4, [6 8 4]); [nv] = reducevolume (v, 0);
%!error <wrong number of arguments> v = reshape(1:6*8*4, [6 8 4]); [nv] = reducevolume (v);
%!error <wrong number of arguments> x = 1:2:24; v = reshape(1:6*8*4, [6 8 4]); [nv] = reducevolume (x, v, 3);
%!error <reduction value R is too high> v = reshape(1:6*8*4, [6 8 4]); [nv] = reducevolume (v, 5);
%!error <reduction value R is too high> v = reshape(1:6*8*4, [6 8 4]); [nv] = reducevolume (v, [4 7 2]);
