changeset 7640:3398ce778b4b

Added support for N-dimensional convolution
author sh@sh-laptop
date Tue, 25 Mar 2008 21:28:02 -0400
parents b2fbb393a072
children 115563ecbdc9
files doc/interpreter/poly.txi scripts/ChangeLog scripts/polynomial/Makefile.in scripts/polynomial/convn.m src/ChangeLog src/DLD-FUNCTIONS/__convn__.cc src/Makefile.in
diffstat 7 files changed, 241 insertions(+), 6 deletions(-) [+]
line wrap: on
line diff
--- 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)
--- 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  <hauberg@gmail.com>
+
+	* polynomial/convn.m: New function.
+	* polynomial/Makefile.in (SOURCES): Add it to the list.
+
 2008-03-25  David Bateman  <dbateman@free.fr>
 
 	* image/contrast.m: New function.
-
+	
 2008-03-24  Thomas Weber  <thomas.weber.mail@gmail.com>
 
 	* pkg/pkg.m: Allow installation of already extracted packages.
--- 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))
 
--- /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
+## <http://www.gnu.org/licenses/>.
+
+## -*- 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
--- 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  <hauberg@gmail.com>
+
+	* DLD-FUNCTIONS/__convn__.cc: New file.
+	* Makefile.in: Add __convn__.cc
+
 2008-03-25  David Bateman  <dbateman@free.fr>
 
 	* DLD-FUNCTIONS/hex2num.cc: New function
--- /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
+<http://www.gnu.org/licenses/>.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <algorithm>
+
+#include "dNDArray.h"
+#include "CNDArray.h"
+
+#include "defun-dld.h"
+
+// FIXME -- this function should maybe be available in liboctave?
+template <class MT, class ST> 
+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<octave_idx_type> a_idx (idx_dim);
+  Array<octave_idx_type> b_idx (idx_dim);
+  Array<octave_idx_type> 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<NDArray, double> (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<ComplexNDArray, Complex> (a, b);
+        }
+      else
+	error ("__convn__: first and second input should be real, or complex arrays");
+    }
+  else
+    print_usage ();
+    
+  return retval;
+}
--- 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))