# HG changeset patch # User address@hidden # Date 1206397313 -3600 # Node ID 87b52a332769a83a519ba4f1460b81d445cb8df0 # Parent 205e43d9d9da89a3841c9e669969aa4f5b7d2a10 Added support for N-dimensional convolution diff -r 205e43d9d9da -r 87b52a332769 ChangeLog --- a/ChangeLog Sat Mar 22 14:59:41 2008 +0100 +++ b/ChangeLog Mon Mar 24 23:21:53 2008 +0100 @@ -1,3 +1,10 @@ 2008-03-21 David Bateman + + * src/DLD-FUNCTIONS/__convn__.cc: new file that implements the 'heavy + lifting' of N-dimensional convolution. + * scripts/polynomial/convn.m: new function for N-dimensional convolution. + * src/Makefile.in: add __convn__.cc to DLD_XSRC. + 2008-03-21 David Bateman * configure.in (HAVE_AMD): Complete test for presence of amd. diff -r 205e43d9d9da -r 87b52a332769 scripts/polynomial/convn.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/polynomial/convn.m Mon Mar 24 23:21:53 2008 +0100 @@ -0,0 +1,83 @@ +## Copyright (C) 2008 Soren Hauberg +## +## This program 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. +## +## This program 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 this program; if not, write to the Free Software +## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +## -*- texinfo -*- +## @deftypefn {Function File} @var{c} = convn (@var{a}, @var{b}, @var{shape}) +## @math{N}-dimensional convolution of matrices @var{a} and @var{b}. +## +## The size of the output is determined by the @var{shape} argument. This can be +## any of the following strings +## +## @table @asis +## @item "full" +## The full convolution result is returned. The size out of the output is +## @code{size (@var{a}) + size (@var{b})-1}. This is the default behaviour. +## @item "same" +## The central part of the convolution result is returned. The size out of the +## output is the same as @var{a}. +## @item "valid" +## The valid part of the convolution is returned. The size of the result is +## @code{max (size (@var{a}) - size (@var{b})+1, 0)}. +## @end table +## +## @seealso{conv, conv2} +## @end deftypefn + +function c = convn (a, b, shape = "full") + ## Check input + if (nargin < 2) + error ("convn: not enough input arguments"); + endif + if (!ismatrix (a) || !ismatrix (b) || ndims (a) != ndims (b)) + error ("convn: first and second arguments must be matrices of the same dimensionality"); + endif + if (!ischar (shape)) + error ("convn: third input argument must be a string"); + endif + if (!any (strcmpi (shape, {"full", "same", "valid"}))) + error ("convn: invalid shape argument: '%s'", shape); + endif + + ## Should we swap 'a' and 'b'? + ## XXX: Should we also swap in any of the non-full cases? + if (numel (b) > numel (a) && strcmpi (shape, "full")) + tmp = a; + a = b; + b = tmp; + endif + + ## Pad A + switch (lower(shape)) + case "full" + a = pad (a, size (b)-1, size (b)-1); + case "same" + a = pad (a, floor ((size (b)-1)/2), ceil ((size (b)-1)/2)); + endswitch + + ## Perform convolution + c = __convn__ (a, b); +endfunction + +## Helper function that performs the padding +function a = pad (a, left, right) + cl = class (a); + for dim = 1:ndims (a) + l = r = size (a); + l (dim) = left (dim); + r (dim) = right (dim); + a = cat (dim, zeros (l, cl), a, zeros (r, cl)); + endfor +endfunction diff -r 205e43d9d9da -r 87b52a332769 src/DLD-FUNCTIONS/__convn__.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/DLD-FUNCTIONS/__convn__.cc Mon Mar 24 23:21:53 2008 +0100 @@ -0,0 +1,115 @@ +// Copyright (C) 2008 Soren Hauberg +// +// This program 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. +// +// This program 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 this program; if not, see . + +#include + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include "defun-dld.h" + +template +octave_value +convn (const MT &A, const MT &B) +{ + // Get sizes + const octave_idx_type ndims = A.ndims (); + const octave_idx_type B_numel = B.numel (); + const dim_vector A_size = A.dims (); + const dim_vector B_size = B.dims (); + + // Check input + if (ndims != B.ndims ()) + { + error ("__convn__: first and second argument must have same dimensionality"); + return octave_value (); + } + + // Allocate output + dim_vector out_size (A_size); + for (octave_idx_type n = 0; n < ndims; n++) + out_size(n) = std::max (A_size (n) - B_size (n) + 1, 0); + MT out = MT (out_size); + const octave_idx_type out_numel = out.numel (); + + // Iterate over every element of 'out'. + dim_vector idx_dim (ndims); + Array A_idx (idx_dim); + Array B_idx (idx_dim); + Array out_idx (idx_dim, 0); + for (octave_idx_type i = 0; i < out_numel; i++) + { + OCTAVE_QUIT; + + // For each neighbour + ST sum = 0; + for (octave_idx_type n = 0; n < ndims; n++) + B_idx(n) = 0; + for (octave_idx_type j = 0; j < B_numel; j++) + { + for (octave_idx_type n = 0; n < ndims; n++) + A_idx (n) = out_idx (n) + (B_size (n) - 1 - B_idx (n)); + sum += A (A_idx) * B (B_idx); + B.increment_index (B_idx, B_size); + } + + // Compute filter result + out (out_idx) = sum; + + // Prepare for next iteration + out.increment_index (out_idx, out_size); + } + + return octave_value (out); +} + +DEFUN_DLD (__convn__, args, , + "-*- texinfo -*-\n\ address@hidden {Loadable Function} {} __convn__(@var{a}, @var{b})\n\ +Undocumented internal function.\n\ address@hidden deftypefn\n\ +") +{ + octave_value_list retval; + if (args.length () == 2) + { + // Take action depending on input type + if (args (0).is_real_matrix () && args (1).is_real_matrix ()) + { + const NDArray A = args (0).array_value(); + const NDArray B = args (1).array_value(); + retval (0) = convn (A, B); + } + else if (args (0).is_complex_matrix () && args (1).is_complex_matrix ()) + { + const ComplexNDArray A = args (0).complex_matrix_value (); + const ComplexNDArray B = args (1).complex_matrix_value (); + retval (0) = convn (A, B); + } + else + { + error ("__convn__: first and second input should be real, or complex arrays"); + } + } + else // nargin != 2 + { + print_usage (); + } + + return retval; +} diff -r 205e43d9d9da -r 87b52a332769 src/Makefile.in --- a/src/Makefile.in Sat Mar 22 14:59:41 2008 +0100 +++ b/src/Makefile.in Mon Mar 24 23:21:53 2008 +0100 @@ -74,7 +74,7 @@ DLD_XSRC := amd.cc balance.cc besselj.cc time.cc tsearch.cc typecast.cc \ urlwrite.cc __contourc__.cc __delaunayn__.cc __dsearchn__.cc \ __glpk__.cc __lin_interpn__.cc __pchip_deriv__.cc \ - __qp__.cc __voronoi__.cc + __qp__.cc __voronoi__.cc __convn__.cc DLD_SRC := $(addprefix DLD-FUNCTIONS/, $(DLD_XSRC))