## 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);