## 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
## <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn  {Function File} {@var{smoothed_data} =} smooth3 (@var{data})
## @deftypefnx {Function File} {@var{smoothed_data} =} smooth3 (@var{data}, @var{method})
## @deftypefnx {Function File} {@var{smoothed_data} =} smooth3 (@var{data}, @var{method}, @var{sz})
## @deftypefnx {Function File} {@var{smoothed_data} =} smooth3 (@var{data}, @var{method}, @var{sz}, @var{st_dev})
## Smooth values of 3-dimensional matrix @var{data}.
##
## @var{data} must be a non-singleton 3-dimensional matrix. The smoothed data of
## this matrix is returned in @var{smoothed_data} which is of the same size as
## @var{data}.
##
## Optionally, @var{method} determines which convolution kernel is used for the
## smoothing process. Choose one of:
## @table @asis
## @item @code{box}, @code{b} (default)
## to use a convolution kernel with sharp edges.
## @item @code{gaussian}, @code{g}
## to use a convolution kernel that is represented by a non-correlated
## trivariate normal distribution function.
## @end table
##
## Optionally, @var{sz} can be either a vector of 3 elements representing the
## size of the convolution kernel in x-, y- and z-direction or a scalar in which
## case the same size is used in all three dimensions. The default value is 3.
##
## When @var{method} is @code{gaussian}, @var{st_dev} can define the standard
## deviation of the trivariate normal distribution function. @var{st_dev} can be
## either a vector of 3 elements representing the standard deviation of the
## gaussian convolution kernel in x-, y- and z-direction or a scalar in which
## case the same value is used in all three dimensions. The default value is
## 0.65.
##
## Example:
## @example
## @group
## data = rand (10, 10, 10);
## smoothed_data = smooth3 (data);
## patch (isosurface (data, .5), "FaceColor", "blue", "EdgeColor", "none");
## @end group
## @end example
##
## @seealso{isosurface, reducevolume}
## @end deftypefn

## Author: mmuetzel

function smoothed_data = smooth3 (data, method, sz, st_dev)

  narginchk (1, 4);

  if (nargin < 2), method = "box"; endif
  if (nargin < 3), sz = 3; endif
  if (nargin < 4), st_dev = .65; endif

  [data, conv_kernel, sz, st_dev] = __get_check_smooth3_args (data, method, sz, st_dev);

  ## manually pad data by replicating the values at the edges (convn would pad with zeros)
  idx = cell (3, 1);
  for i_dim = 1:3
    sz_dim = size (data, i_dim);
    pad_vec = ones (1, (sz(i_dim)-1)/2);
    idx{i_dim} = [pad_vec 1:sz_dim sz_dim*pad_vec];
  endfor
  data_padded = data(idx{:});
  
  ## actual smoothing
  smoothed_data = convn (data_padded, conv_kernel, "valid");
  
endfunction

function [data, conv_kernel, sz, st_dev] = __get_check_smooth3_args (data, method, sz, st_dev);

  ## check data
  if (ndims (data) != 3)
    error("smooth3: DATA must have 3 dimensions");
  endif

%  data_size = size (data);
%  if (any (data_size < 2))
%    error ("smooth3: DATA must be a non-singleton 3-dimensional matrix");
%  endif

  ## check sz
  if (isscalar (sz))
    sz(1:3) = sz;
  endif
  
  if (numel (sz) != 3)
    error ("smooth3: the size SZ of the convolution kernel must either be a scalar or a vector of length 3");
  endif
  
  if (any (sz < 1) || any (rem (sz, 2) != 1))
    error ("smooth3: the size SZ of the convolution kernel must consist of positive odd integers");
  endif
  
  ## check method
  switch lower(method)
    case {"g", "gaussian"}
      ## check st_dev
      if (isscalar (st_dev))
        st_dev(1:3) = st_dev;
      endif
      
      if (numel (st_dev) != 3)
        error ("smooth3: the standard deviation of the gaussian convolution kernel must either be a scalar or a vector of length 3");
      endif
      
      conv_kernel = __smooth3_gaussian3__ (sz, st_dev);
    case {"b", "box"}
      conv_kernel = ones (sz) / prod (sz);
    otherwise
      error ("smooth3: METHOD '%s' unknown", method)
  endswitch
  
endfunction

function gaussian3 = __smooth3_gaussian3__ (sz, st_dev)
  ## trivariate non-correlated gaussian distribution function
  
  x = (-(sz(2)-1)/2:(sz(2)-1)/2) / st_dev(2);
  y = (-(sz(1)-1)/2:(sz(1)-1)/2) / st_dev(1);
  z = (-(sz(3)-1)/2:(sz(3)-1)/2) / st_dev(3);
  
  [xx, yy, zz] = meshgrid (x, y, z);
  
  gaussian3 = exp (-(xx.*xx + yy.*yy + zz.*zz)/2);
  
  gaussian3 = gaussian3 / sum (gaussian3(:)); ## normalize

endfunction

% one input argument (method: "box")
%!test
%! a = rand(10, 8, 7);
%! b = smooth3(a);
%! assert (size_equal (a, b), true);

% two input argument (method: "gaussian")
%!test
%! a = rand(5, 8, 7);
%! b = smooth3(a, "gaussian");
%! assert (size_equal (a, b), true);

% three input argument (method: "box")
%!test
%! a = rand(3, 8, 7);
%! b = smooth3(a, "box", 5);
%! assert (size_equal (a, b), true);

% three input argument (method: "gaussian")
%!test
%! a = rand(3, 8, 7);
%! b = smooth3(a, "gaussian", 7);
%! assert (size_equal (a, b), true);

% size of convolution kernel = 1: no smoothing (method: "box")
%!test
%! a = rand(9, 8, 7);
%! b = smooth3(a, "box", 1);
%! assert (a, b);

% size of convolution kernel = 1: no smoothing (method: "gaussian")
%!test
%! a = rand(9, 8, 7);
%! b = smooth3(a, "gaussian", 1);
%! assert (a, b);

% four input arguments (method: "gaussian")
%!test
%! a = rand(3, 8, 7);
%! b = smooth3(a, "gaussian", 7, .5);
%! assert (size_equal (a, b), true);

% size of convolution kernel is different in x, y and z (method: "box")
%!test
%! a = rand(3, 8, 7);
%! b = smooth3(a, "box", [5 3 7]);
%! assert (size_equal (a, b), true);

% size of convolution kernel is different in x, y and z (method: "gaussian")
%!test
%! a = rand(3, 8, 7);
%! b = smooth3(a, "gaussian", [5 3 7]);
%! assert (size_equal (a, b), true);

% size and width of gaussian convolution kernel is different in x, y and z (method: "gaussian")
%!test
%! a = rand(3, 8, 7);
%! b = smooth3(a, "gaussian", [7 3 5], [.3 .5 .4]);
%! assert (size_equal (a, b), true);

% test for each error
%!test
%!error <DATA must have 3 dimensions> a = rand(3, 8); b = smooth3(a);
%!error <the size SZ of the convolution kernel must either be a scalar or a vector of length 3> a = rand(3, 8, 3); b = smooth3(a, "box", [3 5]);
%!error <the size SZ of the convolution kernel must consist of positive odd integers> a = rand(3, 8, 3); b = smooth3(a, "box", [3 2 5]);
%!error <the size SZ of the convolution kernel must consist of positive odd integers> a = rand(3, 8, 3); b = smooth3(a, "box", [3 0 5]);
%!error <the size SZ of the convolution kernel must consist of positive odd integers> a = rand(3, 8, 3); b = smooth3(a, "box", [3 2.5 5]);
%!error <the standard deviation of the gaussian convolution kernel must either be a scalar or a vector of length 3> a = rand(3, 8, 3); b = smooth3(a, "gaussian", 3, [.3 .4]);
%!error <METHOD 'other' unknown> a = rand(3, 8, 3); b = smooth3(a, "other");
