## 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 -*-
## @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} =} surface (@dots{})
## @deftypefnx {Function File} {[@var{nx}, @var{ny}, @var{nz}, @var{nv}] =} reducevolume (@var{v}, @var{r})
## @deftypefnx {Function File} {@var{nv} =} surface (@dots{})
## Reduce the volume of the dataset in @var{v} according to the values in
## @var{r}. 

## 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 ("reducvolume: 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 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 = 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);
  else
    [reduced_x, reduced_y, reduced_z] = deal([]);
  endif

endfunction
