From octave-sources-request at bevo dot che dot wisc dot edu Thu Jun 3 19:34:01 1999 Subject: conv2 2-D convolution From: A+A Adler To: octave-sources at bevo dot che dot wisc dot edu Date: Thu, 3 Jun 1999 20:32:08 -0400 (EDT) The following *.oct file replicates the conv2 functionality in Matlab (including, of course, the Matlab bugs) (ie conv2(small matrix, bigmatrix, 'valid') outputs [] , which is incorrect ) It runs about 2 x faster (on an i686 Linux, for large random matrices) While writing this code I noticed some interesting things about optimizing *.oct files. See my accompanying msg on help-octave. ______________________________________________________________ Andy Adler, adler at mondenet dot com /* * conv2: 2D convolution for octave * * Copyright (C) 1999 Andy Adler * This code has no warrany whatsoever. * Do what you like with this code as long as you * leave this copyright in place. * * $Id: conv2.cc,v 1.6 1999/06/03 17:25:58 aadler Exp aadler $ */ #include #define MAX(a,b) ((a) > (b) ? (a) : (b)) DEFUN_DLD (conv2, args, , "[...] = conv2 (...) CONV2: do 2 dimentional convolution c= conv2(a,b) -> same as c= conv2(a,b,'full') c= conv2(a,b,shape) returns 2-D convolution of a and b where the size of c is given by shape= 'full' -> returns full 2-D convolution shape= 'same' -> same size as a. 'central' part of convolution shape= 'valid' -> only parts which do not include zero-padded edges c= conv2(a,b,shape) returns 2-D convolution of a and b c= conv2(v1,v2,a) -> same as c= conv2(v1,v2,a,'full') c= conv2(v1,v2,a,shape) returns convolution of a by vector v1 in the column direction and vector v2 in the row direction ") { octave_value_list retval; octave_value tmp; int nargin = args.length (); string shape= "full"; bool separable= false; int outM, outN, edgM, edgN; if (nargin < 2 ) { print_usage ("conv2"); return retval; } else if (nargin == 3) { if ( args(2).is_string() ) shape= args(2).string_value(); else separable= true; } else if (nargin >= 4) { separable= true; shape= args(3).string_value(); } /* * Here we calculate the size of the output matrix, * in order to stay Matlab compatible, it is based * on the third parameter if its separable, and the * first if it's not */ if (separable) { /* * Check that the first two parameters are vectors * if we're doing separable */ if ( !( 1== args(0).rows() || 1== args(0).columns() ) || !( 1== args(1).rows() || 1== args(1).columns() ) ) { print_usage ("conv2"); return retval; } ColumnVector C= args(0).vector_value(); ColumnVector R= args(1).vector_value(); Matrix A= args(2).matrix_value(); int Rn= R.length(); int Cm= C.length(); int Am = A.rows(); int An = A.columns(); if ( shape == "full" ) { outM= Am + Cm - 1; outN= An + Rn - 1; edgM= Cm - 1; edgN= Rn - 1; } else if ( shape == "same" ) { outM= Am; outN= An; // Matlab seems to arbitrarily choose this convention for // 'same' with even length R, C edgM= ( Cm - 1) /2; edgN= ( Rn - 1) /2; } else if ( shape == "valid" ) { outM= Am - Cm + 1; outN= An - Rn + 1; edgM= edgN= 0; } else { // if ( shape error("Shape type not valid"); print_usage ("conv2"); return retval; } // printf("A(%d,%d) C(%d) R(%d) O(%d,%d) E(%d,%d)\n", // Am,An, Cm,Rn, outM, outN, edgM, edgN); Matrix O(outM,outN); /* * T accumulated the 1-D conv for each row, before calculating * the convolution in the other direction * There is no efficiency advantage to doing it in either direction * first */ RowVector T( An ); for( int oi=0; oi < outM; oi++ ) { for( int oj=0; oj < An; oj++ ) { double sum=0; int ci= Cm - 1 - MAX(0, edgM-oi); int ai= MAX(0, oi-edgM) ; const double* Ad= A.data() + ai + Am*oj; const double* Cd= C.data() + ci; for( ; ci >= 0 && ai < Am; ci--, Cd--, ai++, Ad++) { sum+= (*Ad) * (*Cd); } // for( int ci= T(oj)= sum; } // for( int oj=0 for( int oj=0; oj < outN; oj++ ) { double sum=0; int rj= Rn - 1 - MAX(0, edgN-oj); int aj= MAX(0, oj-edgN) ; const double* Td= T.data() + aj; const double* Rd= R.data() + rj; for( ; rj >= 0 && aj < An; rj--, Rd--, aj++, Td++) { sum+= (*Td) * (*Rd); } //for( int rj= O(oi,oj)= sum; } // for( int oj=0 } // for( int oi=0 retval(0)= O; } else { // if (separable) /* Convolution works fastest if we choose the A matrix to be * the largest. * * NOTE in order to be Matlab compatible, we give * wrong sizes for 'valid' if the smallest matrix is first */ Matrix A = args( 0 ).matrix_value(); Matrix B = args( 1 ).matrix_value(); int Am = A.rows(); int An = A.columns(); int Bm = B.rows(); int Bn = B.columns(); if ( shape == "full" ) { outM= Am + Bm - 1; outN= An + Bn - 1; edgM= Bm - 1; edgN= Bn - 1; } else if ( shape == "same" ) { outM= Am; outN= An; // Matlab seems to arbitrarily choose this convention for // 'same' with even length R, C edgM= ( Bm - 1) /2; edgN= ( Bn - 1) /2; } else if ( shape == "valid" ) { outM= Am - Bm + 1; outN= An - Bn + 1; edgM= edgN= 0; } else { // if ( shape error("Shape type not valid"); print_usage ("conv2"); return retval; } // printf("A(%d,%d) B(%d,%d) O(%d,%d) E(%d,%d)\n", // Am,An, Bm,Bn, outM, outN, edgM, edgN); Matrix O(outM,outN); for( int oi=0; oi < outM; oi++ ) { for( int oj=0; oj < outN; oj++ ) { double sum=0; for( int bj= Bn - 1 - MAX(0, edgN-oj), aj= MAX(0, oj-edgN); bj >= 0 && aj < An; bj--, aj++) { int bi= Bm - 1 - MAX(0, edgM-oi); int ai= MAX(0, oi-edgM); const double* Ad= A.data() + ai + Am*aj; const double* Bd= B.data() + bi + Bm*bj; for( ; bi >= 0 && ai < Am; bi--, Bd--, ai++, Ad++) { sum+= (*Ad) * (*Bd); /* * It seems to be about 2.5 times faster to use pointers than * to do this * sum+= A(ai,aj) * B(bi,bj); */ } // for( int bi= } //for( int bj= O(oi,oj)= sum; } // for( int oj= } // for( int oi= retval(0)= O; } // if (separable) return retval; }