# HG changeset patch # User sh@sh-laptop # Date 1206494882 14400 # Node ID 3398ce778b4bfec105dffd515b0f162dfd45ae21 # Parent b2fbb393a072bd7feb019fa3c9756c465f8fba9c Added support for N-dimensional convolution diff -r b2fbb393a072 -r 3398ce778b4b doc/interpreter/poly.txi --- a/doc/interpreter/poly.txi Tue Mar 25 20:22:46 2008 -0400 +++ b/doc/interpreter/poly.txi Tue Mar 25 21:28:02 2008 -0400 @@ -90,6 +90,8 @@ @DOCSTRING(conv) +@DOCSTRING(convn) + @DOCSTRING(deconv) @DOCSTRING(conv2) diff -r b2fbb393a072 -r 3398ce778b4b scripts/ChangeLog --- a/scripts/ChangeLog Tue Mar 25 20:22:46 2008 -0400 +++ b/scripts/ChangeLog Tue Mar 25 21:28:02 2008 -0400 @@ -1,7 +1,12 @@ +2008-03-25 Soren Hauberg + + * polynomial/convn.m: New function. + * polynomial/Makefile.in (SOURCES): Add it to the list. + 2008-03-25 David Bateman * image/contrast.m: New function. - + 2008-03-24 Thomas Weber * pkg/pkg.m: Allow installation of already extracted packages. diff -r b2fbb393a072 -r 3398ce778b4b scripts/polynomial/Makefile.in --- a/scripts/polynomial/Makefile.in Tue Mar 25 20:22:46 2008 -0400 +++ b/scripts/polynomial/Makefile.in Tue Mar 25 21:28:02 2008 -0400 @@ -33,10 +33,10 @@ INSTALL_PROGRAM = @INSTALL_PROGRAM@ INSTALL_DATA = @INSTALL_DATA@ -SOURCES = compan.m conv.m deconv.m mkpp.m mpoles.m pchip.m poly.m \ - polyder.m polyderiv.m polyfit.m polygcd.m polyint.m \ - polyout.m polyreduce.m polyval.m polyvalm.m ppval.m residue.m \ - roots.m spline.m unmkpp.m +SOURCES = compan.m conv.m convn.m deconv.m mkpp.m mpoles.m \ + pchip.m poly.m polyder.m polyderiv.m polyfit.m polygcd.m \ + polyint.m polyout.m polyreduce.m polyval.m polyvalm.m \ + ppval.m residue.m roots.m spline.m unmkpp.m DISTFILES = $(addprefix $(srcdir)/, Makefile.in $(SOURCES)) diff -r b2fbb393a072 -r 3398ce778b4b scripts/polynomial/convn.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/polynomial/convn.m Tue Mar 25 21:28:02 2008 -0400 @@ -0,0 +1,89 @@ +## Copyright (C) 2008 SÃ?ren Hauberg +## +## 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{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 character 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") + + 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'? + ## FIXME -- 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 b2fbb393a072 -r 3398ce778b4b src/ChangeLog --- a/src/ChangeLog Tue Mar 25 20:22:46 2008 -0400 +++ b/src/ChangeLog Tue Mar 25 21:28:02 2008 -0400 @@ -1,3 +1,8 @@ +2008-03-25 Soren Hauberg + + * DLD-FUNCTIONS/__convn__.cc: New file. + * Makefile.in: Add __convn__.cc + 2008-03-25 David Bateman * DLD-FUNCTIONS/hex2num.cc: New function diff -r b2fbb393a072 -r 3398ce778b4b src/DLD-FUNCTIONS/__convn__.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/DLD-FUNCTIONS/__convn__.cc Tue Mar 25 21:28:02 2008 -0400 @@ -0,0 +1,134 @@ +/* + +Copyright (C) 2008 Soren Hauberg + +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 +. + +*/ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +#include "dNDArray.h" +#include "CNDArray.h" + +#include "defun-dld.h" + +// FIXME -- this function should maybe be available in liboctave? +template +octave_value +convn (const MT& a, const MT& b) +{ + octave_value retval; + + // 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 (); + + if (ndims != b.ndims ()) + { + error ("__convn__: first and second argument must have same dimensionality"); + return retval; + } + + // 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 out; +} + +DEFUN_DLD (__convn__, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {} __convn__(@var{a}, @var{b})\n\ +Undocumented internal function.\n\ +@end deftypefn\n\ +") +{ + octave_value retval; + + if (args.length () == 2) + { + if (args(0).is_real_type () && args(1).is_real_type ()) + { + const NDArray a = args (0).array_value (); + const NDArray b = args (1).array_value (); + + if (! error_state) + retval = convn (a, b); + } + else if (args(0).is_complex_type () && args(1).is_complex_type ()) + { + const ComplexNDArray a = args (0).complex_matrix_value (); + const ComplexNDArray b = args (1).complex_matrix_value (); + + if (! error_state) + retval = convn (a, b); + } + else + error ("__convn__: first and second input should be real, or complex arrays"); + } + else + print_usage (); + + return retval; +} diff -r b2fbb393a072 -r 3398ce778b4b src/Makefile.in --- a/src/Makefile.in Tue Mar 25 20:22:46 2008 -0400 +++ b/src/Makefile.in Tue Mar 25 21:28:02 2008 -0400 @@ -74,7 +74,7 @@ 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))