changelog shortlog tags changeset files revisions annotate raw

src/DLD-FUNCTIONS/conv2.cc

changeset 10289: 4b124317dc38
parent:40dfc0c99116
author: John W. Eaton <jwe@octave.org>
date: Tue Feb 09 20:58:55 2010 -0500 (71 minutes ago)
permissions: -rw-r--r--
description: base_properties::set_children: account for hidden children 
     1/*

     2

     3Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008

     4              Andy Adler

     5

     6This file is part of Octave.

     7

     8Octave is free software; you can redistribute it and/or modify it

     9under the terms of the GNU General Public License as published by the

    10Free Software Foundation; either version 3 of the License, or (at your

    11option) any later version.

    12

    13Octave is distributed in the hope that it will be useful, but WITHOUT

    14ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or

    15FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License

    16for more details.

    17

    18You should have received a copy of the GNU General Public License

    19along with Octave; see the file COPYING.  If not, see

    20<http://www.gnu.org/licenses/>.

    21

    22*/

    23

    24#ifdef HAVE_CONFIG_H

    25#include <config.h>

    26#endif

    27

    28#include "defun-dld.h"

    29#include "error.h"

    30#include "oct-obj.h"

    31#include "utils.h"

    32

    33#define MAX(a,b) ((a) > (b) ? (a) : (b))

    34

    35enum Shape { SHAPE_FULL, SHAPE_SAME, SHAPE_VALID };

    36

    37#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)

    38extern MArray2<double>

    39conv2 (MArray<double>&, MArray<double>&, MArray2<double>&, Shape);

    40

    41extern MArray2<Complex>

    42conv2 (MArray<Complex>&, MArray<Complex>&, MArray2<Complex>&, Shape);

    43

    44extern MArray2<float>

    45conv2 (MArray<float>&, MArray<float>&, MArray2<float>&, Shape);

    46

    47extern MArray2<FloatComplex>

    48conv2 (MArray<FloatComplex>&, MArray<FloatComplex>&, MArray2<FloatComplex>&, Shape);

    49#endif

    50

    51template <class T>

    52MArray2<T>

    53conv2 (MArray<T>& R, MArray<T>& C, MArray2<T>& A, Shape ishape)

    54{

    55  octave_idx_type  Rn =  R.length ();

    56  octave_idx_type  Cm =  C.length ();

    57  octave_idx_type  Am = A.rows ();

    58  octave_idx_type  An = A.columns ();

    59

    60  // Calculate the size of the output matrix:

    61  // in order to stay Matlab compatible, it is based

    62  // on the third parameter if it's separable, and the

    63  // first if it's not

    64

    65  octave_idx_type outM = 0;

    66  octave_idx_type outN = 0;

    67  octave_idx_type edgM = 0;

    68  octave_idx_type edgN = 0;

    69

    70  switch (ishape)

    71    {

    72      case SHAPE_FULL:

    73        outM = Am + Cm - 1;

    74        outN = An + Rn - 1;

    75        edgM = Cm - 1;

    76        edgN = Rn - 1;

    77        break;

    78

    79      case SHAPE_SAME:

    80        outM = Am;

    81        outN = An;

    82        // Follow the Matlab convention (ie + instead of -)

    83        edgM = (Cm - 1) /2;

    84        edgN = (Rn - 1) /2;

    85        break;

    86

    87      case SHAPE_VALID:

    88        outM = Am - Cm + 1;

    89        outN = An - Rn + 1;

    90        if (outM < 0)

    91          outM = 0;

    92        if (outN < 0)

    93          outN = 0;

    94        edgM = edgN = 0;

    95        break;

    96

    97      default:

    98        error ("conv2: invalid value of parameter ishape");

    99    }

   100

   101  MArray2<T> O (outM, outN);

   102

   103  // X accumulates the 1-D conv for each row, before calculating

   104  //    the convolution in the other direction

   105  // There is no efficiency advantage to doing it in either direction

   106  //     first

   107

   108  MArray<T> X (An);

   109

   110  for (octave_idx_type oi = 0; oi < outM; oi++)

   111    {

   112      for (octave_idx_type oj = 0; oj < An; oj++)

   113        {

   114          T sum = 0;

   115

   116          octave_idx_type ci = Cm - 1 - MAX(0, edgM-oi);

   117          octave_idx_type ai = MAX(0, oi-edgM);

   118          const T* Ad = A.data() + ai + Am*oj;

   119          const T* Cd = C.data() + ci;

   120          for ( ; ci >= 0 && ai < Am; ci--, Cd--, ai++, Ad++)

   121            sum += (*Ad) * (*Cd);

   122

   123          X(oj) = sum;

   124        }

   125

   126      for (octave_idx_type oj = 0; oj < outN; oj++)

   127        {

   128          T sum = 0;

   129

   130          octave_idx_type rj = Rn - 1 - MAX(0, edgN-oj);

   131          octave_idx_type aj = MAX(0, oj-edgN) ;

   132          const T* Xd = X.data() + aj;

   133          const T* Rd = R.data() + rj;

   134

   135          for ( ; rj >= 0 && aj < An; rj--, Rd--, aj++, Xd++)

   136            sum += (*Xd) * (*Rd);

   137

   138          O(oi,oj)=  sum;

   139        }

   140    }

   141

   142  return O;

   143}

   144

   145#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)

   146extern MArray2<double>

   147conv2 (MArray2<double>&, MArray2<double>&, Shape);

   148

   149extern MArray2<Complex>

   150conv2 (MArray2<Complex>&, MArray2<Complex>&, Shape);

   151

   152extern MArray2<float>

   153conv2 (MArray2<float>&, MArray2<float>&, Shape);

   154

   155extern MArray2<FloatComplex>

   156conv2 (MArray2<FloatComplex>&, MArray2<FloatComplex>&, Shape);

   157#endif

   158

   159template <class T>

   160MArray2<T>

   161conv2 (MArray2<T>&A, MArray2<T>&B, Shape ishape)

   162{

   163  // Convolution works fastest if we choose the A matrix to be

   164  // the largest.

   165

   166  // Here we calculate the size of the output matrix,

   167  // in order to stay Matlab compatible, it is based

   168  // on the third parameter if it's separable, and the

   169  // first if it's not

   170

   171  // NOTE in order to be Matlab compatible, we give argueably

   172  // wrong sizes for 'valid' if the smallest matrix is first

   173

   174  octave_idx_type Am = A.rows ();

   175  octave_idx_type An = A.columns ();

   176  octave_idx_type Bm = B.rows ();

   177  octave_idx_type Bn = B.columns ();

   178

   179  octave_idx_type outM = 0;

   180  octave_idx_type outN = 0;

   181  octave_idx_type edgM = 0;

   182  octave_idx_type edgN = 0;

   183

   184  switch (ishape)

   185    {

   186      case SHAPE_FULL:

   187        outM = Am + Bm - 1;

   188        outN = An + Bn - 1;

   189        edgM = Bm - 1;

   190        edgN = Bn - 1;

   191        break;

   192

   193      case SHAPE_SAME:

   194        outM = Am;

   195        outN = An;

   196        edgM = (Bm - 1) /2;

   197        edgN = (Bn - 1) /2;

   198        break;

   199

   200      case SHAPE_VALID:

   201        outM = Am - Bm + 1;

   202        outN = An - Bn + 1;

   203        if (outM < 0)

   204          outM = 0;

   205        if (outN < 0)

   206          outN = 0;

   207        edgM = edgN = 0;

   208        break;

   209    }

   210

   211  MArray2<T> O (outM, outN);

   212

   213  for (octave_idx_type oi = 0; oi < outM; oi++)

   214    {

   215      for (octave_idx_type oj = 0; oj < outN; oj++)

   216        {

   217          T sum = 0;

   218

   219          for (octave_idx_type bj = Bn - 1 - MAX (0, edgN-oj), aj= MAX (0, oj-edgN);

   220               bj >= 0 && aj < An; bj--, aj++)

   221            {

   222              octave_idx_type bi = Bm - 1 - MAX (0, edgM-oi);

   223              octave_idx_type ai = MAX (0, oi-edgM);

   224              const T* Ad = A.data () + ai + Am*aj;

   225              const T* Bd = B.data () + bi + Bm*bj;

   226

   227              for ( ; bi >= 0 && ai < Am; bi--, Bd--, ai++, Ad++)

   228                {

   229                  sum += (*Ad) * (*Bd);

   230                  // Comment: it seems to be 2.5 x faster than this:

   231                  //        sum+= A(ai,aj) * B(bi,bj);

   232                }

   233            }

   234

   235          O(oi,oj) = sum;

   236        }

   237    }

   238

   239  return O;

   240}

   241

   242/*

   243%!test

   244%! b = [0,1,2,3;1,8,12,12;4,20,24,21;7,22,25,18];

   245%! assert(conv2([0,1;1,2],[1,2,3;4,5,6;7,8,9]),b);

   246

   247%!test

   248%! b = single([0,1,2,3;1,8,12,12;4,20,24,21;7,22,25,18]);

   249%! assert(conv2(single([0,1;1,2]),single([1,2,3;4,5,6;7,8,9])),b);

   250*/

   251

   252DEFUN_DLD (conv2, args, ,

   253  "-*- texinfo -*-\n\

   254@deftypefn {Loadable Function} {y =} conv2 (@var{a}, @var{b}, @var{shape})\n\

   255@deftypefnx {Loadable Function} {y =} conv2 (@var{v1}, @var{v2}, @var{M}, @var{shape})\n\

   256\n\

   257Returns 2D convolution of @var{a} and @var{b} where the size\n\

   258of @var{c} is given by\n\

   259\n\

   260@table @asis\n\

   261@item @var{shape}= 'full'\n\

   262returns full 2-D convolution\n\

   263@item @var{shape}= 'same'\n\

   264same size as a. 'central' part of convolution\n\

   265@item @var{shape}= 'valid'\n\

   266only parts which do not include zero-padded edges\n\

   267@end table\n\

   268\n\

   269By default @var{shape} is 'full'.  When the third argument is a matrix\n\

   270returns the convolution of the matrix @var{M} by the vector @var{v1}\n\

   271in the column direction and by vector @var{v2} in the row direction\n\

   272@end deftypefn")

   273{

   274  octave_value retval;

   275  octave_value tmp;

   276  int nargin = args.length ();

   277  std::string shape= "full"; //default

   278  bool separable= false;

   279  Shape ishape;

   280

   281  if (nargin < 2)

   282    {

   283     print_usage ();

   284     return retval;

   285    }

   286  else if (nargin == 3)

   287    {

   288      if (args(2).is_string ())

   289        shape = args(2).string_value ();

   290      else

   291        separable = true;

   292    }

   293  else if (nargin >= 4)

   294    {

   295      separable = true;

   296      shape = args(3).string_value ();

   297    }

   298

   299  if (shape == "full")

   300    ishape = SHAPE_FULL;

   301  else if (shape == "same")

   302    ishape = SHAPE_SAME;

   303  else if (shape == "valid")

   304    ishape = SHAPE_VALID;

   305  else

   306    {

   307      error ("conv2: shape type not valid");

   308      print_usage ();

   309      return retval;

   310    }

   311

   312   if (separable)

   313     {

   314      // If user requests separable, check first two params are vectors

   315

   316       if (! (1 == args(0).rows () || 1 == args(0).columns ())

   317           || ! (1 == args(1).rows () || 1 == args(1).columns ()))

   318         {

   319           print_usage ();

   320           return retval;

   321         }

   322

   323       if (args(0).is_single_type () || 

   324           args(1).is_single_type () || 

   325           args(2).is_single_type ())

   326         {

   327           if (args(0).is_complex_type ()

   328               || args(1).is_complex_type ()

   329               || args(2).is_complex_type ())

   330             {

   331               FloatComplexColumnVector v1 (args(0).float_complex_vector_value ());

   332               FloatComplexColumnVector v2 (args(1).float_complex_vector_value ());

   333               FloatComplexMatrix a (args(2).float_complex_matrix_value ());

   334               FloatComplexMatrix c (conv2 (v1, v2, a, ishape));

   335               if (! error_state)

   336                 retval = c;

   337             }

   338           else

   339             {

   340               FloatColumnVector v1 (args(0).float_vector_value ());

   341               FloatColumnVector v2 (args(1).float_vector_value ());

   342               FloatMatrix a (args(2).float_matrix_value ());

   343               FloatMatrix c (conv2 (v1, v2, a, ishape));

   344               if (! error_state)

   345                 retval = c;

   346             }

   347         }

   348       else

   349         {

   350           if (args(0).is_complex_type ()

   351               || args(1).is_complex_type ()

   352               || args(2).is_complex_type ())

   353             {

   354               ComplexColumnVector v1 (args(0).complex_vector_value ());

   355               ComplexColumnVector v2 (args(1).complex_vector_value ());

   356               ComplexMatrix a (args(2).complex_matrix_value ());

   357               ComplexMatrix c (conv2 (v1, v2, a, ishape));

   358               if (! error_state)

   359                 retval = c;

   360             }

   361           else

   362             {

   363               ColumnVector v1 (args(0).vector_value ());

   364               ColumnVector v2 (args(1).vector_value ());

   365               Matrix a (args(2).matrix_value ());

   366               Matrix c (conv2 (v1, v2, a, ishape));

   367               if (! error_state)

   368                 retval = c;

   369             }

   370         }

   371     } // if (separable)

   372   else

   373     {

   374       if (args(0).is_single_type () || 

   375           args(1).is_single_type ())

   376         {

   377           if (args(0).is_complex_type ()

   378               || args(1).is_complex_type ())

   379             {

   380               FloatComplexMatrix a (args(0).float_complex_matrix_value ());

   381               FloatComplexMatrix b (args(1).float_complex_matrix_value ());

   382               FloatComplexMatrix c (conv2 (a, b, ishape));

   383               if (! error_state)

   384                 retval = c;

   385             }

   386           else

   387             {

   388               FloatMatrix a (args(0).float_matrix_value ());

   389               FloatMatrix b (args(1).float_matrix_value ());

   390               FloatMatrix c (conv2 (a, b, ishape));

   391               if (! error_state)

   392                 retval = c;

   393             }

   394         }

   395       else

   396         {

   397           if (args(0).is_complex_type ()

   398               || args(1).is_complex_type ())

   399             {

   400               ComplexMatrix a (args(0).complex_matrix_value ());

   401               ComplexMatrix b (args(1).complex_matrix_value ());

   402               ComplexMatrix c (conv2 (a, b, ishape));

   403               if (! error_state)

   404                 retval = c;

   405             }

   406           else

   407             {

   408               Matrix a (args(0).matrix_value ());

   409               Matrix b (args(1).matrix_value ());

   410               Matrix c (conv2 (a, b, ishape));

   411               if (! error_state)

   412                 retval = c;

   413             }

   414         }

   415

   416     } // if (separable)

   417

   418   return retval;

   419}

   420

   421template MArray2<double>

   422conv2 (MArray<double>&, MArray<double>&, MArray2<double>&, Shape);

   423

   424template MArray2<double>

   425conv2 (MArray2<double>&, MArray2<double>&, Shape);

   426

   427template MArray2<Complex>

   428conv2 (MArray<Complex>&, MArray<Complex>&, MArray2<Complex>&, Shape);

   429

   430template MArray2<Complex>

   431conv2 (MArray2<Complex>&, MArray2<Complex>&, Shape);