## 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{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);
## figure;
## subplot (1, 2, 1)
## patch (isosurface (data, .5), "FaceColor", "blue", "EdgeColor", "k");
## title ("Original data")
## view(3)
## smoothed_data = smooth3 (data);
## subplot (1, 2, 2)
## patch (isosurface (smoothed_data, .5), "FaceColor", "blue", "EdgeColor", "k");
## title ("Smoothed data")
## view(3)
## @end group
## @end example
##
## @seealso{isosurface, isocaps, isonormals, patch, reducevolume}
## @end deftypefn
## Author: mmuetzel
function smoothed_data = smooth3 (data, method, sz, st_dev)
if (nargin < 1 || nargin > 4 || nargout > 1)
print_usage ();
endif
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
## 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 a = rand(3, 8); b = smooth3(a);
%!error a = rand(3, 8, 3); b = smooth3(a, "box", [3 5]);
%!error a = rand(3, 8, 3); b = smooth3(a, "box", [3 2 5]);
%!error a = rand(3, 8, 3); b = smooth3(a, "box", [3 0 5]);
%!error a = rand(3, 8, 3); b = smooth3(a, "box", [3 2.5 5]);
%!error a = rand(3, 8, 3); b = smooth3(a, "gaussian", 3, [.3 .4]);
%!error a = rand(3, 8, 3); b = smooth3(a, "other");
%!error a = rand(3, 8, 3); b = smooth3(a, "box", [7 3 5], .4, 1);