mtx_conv: one or two dimensional convolution for matrices. help-patch

added; object repaired.

svn path=/trunk/externals/iem/iemmatrix/; revision=5717

2 files changed, 141 insertions, 110 deletions

diff --git a/doc/mtx_conv-help.pd b/doc/mtx_conv-help.pd
new file mode 100644
index 0000000..31d2905
--- /dev/null
+++ b/doc/mtx_conv-help.pd
@@ -0,0 +1,33 @@
+#N canvas 166 141 740 545 10;
+#X text 458 17 updated for;
+#X obj 546 17 iemmatrix 0.2;
+#X obj 595 43 matrix;
+#X text 465 42 see also help for;
+#X text 90 15 [mtx_conv];
+#X text 47 34 2-dimensional convolution;
+#X text 18 63 you can calculate the convolution of a matrix with a
+filter matrix kernel (2 dimensional FIR filtering). Of course \, it
+can also be used for 1-dimensional FIR convolutions.;
+#X obj 63 196 matrix 1 10;
+#X msg 64 176 element 1 1 1;
+#X msg 42 127 bang;
+#X obj 42 151 t b b b;
+#X obj 64 215 t a a;
+#X obj 94 248 mtx_print orig;
+#X obj 64 274 mtx_conv;
+#X obj 173 197 t a a;
+#X obj 200 217 mtx_print filter_kernel;
+#X obj 64 304 mtx_print result;
+#X msg 173 174 matrix 2 2 1 1 1 1;
+#X connect 7 0 11 0;
+#X connect 8 0 7 0;
+#X connect 9 0 10 0;
+#X connect 10 0 7 0;
+#X connect 10 1 8 0;
+#X connect 10 2 17 0;
+#X connect 11 0 13 0;
+#X connect 11 1 12 0;
+#X connect 13 0 16 0;
+#X connect 14 0 13 1;
+#X connect 14 1 15 0;
+#X connect 17 0 14 0;
diff --git a/src/mtx_conv.c b/src/mtx_conv.c
index e41597a..e1b4f03 100644
--- a/src/mtx_conv.c
+++ b/src/mtx_conv.c
@@ -23,96 +23,69 @@ struct _MTXConv_
    int size;
    int rows;
    int columns;
+   
    int rows_k;
    int columns_k;
    int size_k;
+   
    int rows_y;
    int columns_y;
    int size_y;
-   t_float renorm_fac;
 
    t_float **x;
+   t_float *x_array;
    t_float **k;
+   t_float *k_array;
    t_float **y;
+   t_float *y_array;
 
    t_outlet *list_outlet;
    
    t_atom *list;
 };
 
-static t_float **getTFloatMatrix (int rows, int columns)
+static void getTFloatMatrix (int rows, int columns, t_float ***mtx, t_float **array)
 {
-   t_float **matrix = (t_float **) getbytes 
-      (sizeof (t_float *) * columns);
-   t_float **mtx = matrix;
-
-   while (columns--)
-      if (!(*matrix++ = (t_float *) getbytes
-	       (sizeof (t_float) * rows)))
-	 post("an error has occured :)");
-   return mtx;
-}
+   int size = rows*columns;
+   t_float *ptr;
+   t_float **dptr;
 
-static void deleteTFloatMatrix (t_float **matrix, int rows, int columns)
-{
-   int n = columns;
-   t_float **matr = matrix;
-   if (matr) {
-      while (n--)
-	 if (*matr)
-	    freebytes (*matr++, sizeof(t_float) * rows);
-      freebytes (matrix, sizeof(t_float*) * columns);
+   if (!size)
+      return;
+   
+   if (*array=ptr=(t_float *)calloc(sizeof(t_float),size)) {
+      if (*mtx=dptr=(t_float **)calloc(sizeof(t_float *),rows)) {
+	 for(;rows-- ; ptr+=columns) {
+	    *dptr++ = ptr;
+	 }
+      } else {
+	 freebytes (*array,sizeof(t_float)*size);
+	 array=0;
+      }
    }
 }
 
-
-static t_float **resizeTFloatMatrix (t_float **old, int rows_old, int columns_old,
-      int rows_new, int columns_new)
+static void deleteTFloatMatrix (int rows, int columns, t_float ***mtx, t_float **array)
 {
-   t_float **mtx = old;
-   int count1;
-   post("resizing from %dx%d to %dx%d", rows_old, columns_old, rows_new, columns_new);
-   
-   if ((rows_new == 0)||(columns_new == 0)) {
-      deleteTFloatMatrix (old, rows_old, columns_old);
-      old = 0;
-      return old;
-   }
-   // 1. if rows_old>rows_new: old row disposal
-   if (rows_old>rows_new)
-      for (count1 = (rows_old - rows_new), mtx += rows_new;
-	    count1--; mtx++) 
-	 freebytes (*mtx, sizeof(t_float) * columns_old);
-   // 2. resize row (double) pointer
-   mtx = old = (t_float **) resizebytes(old, sizeof(t_float*) * rows_old, 
-	 sizeof(t_float*) * rows_new);
-   // 3. resizing rows if new size is smaller
-   if (rows_old>rows_new)
-      for (count1 = rows_new; count1--; mtx++) 
-	 *mtx = (t_float *) resizebytes (*mtx, sizeof(t_float) * columns_old,
-	       sizeof(t_float) * columns_new);
-   else { // 4. resizing old rows if new size is bigger, allocating new rows
-      for (count1 = rows_old; count1--; mtx++)
-	 *mtx = (t_float *) resizebytes (*mtx, sizeof(t_float) * columns_old,
-	       sizeof(t_float) * columns_new);
-      for (count1 = (rows_new - rows_old); count1--; mtx++)
-	 *mtx = (t_float *) getbytes (sizeof(t_float) * columns_new);
-   }
-   /* post("return resize"); */
-   return old;
+   int size = rows*columns;
+
+   if (*mtx) 
+      freebytes (*mtx, sizeof(t_float*) * columns);
+   if (*array)
+      freebytes (*array, sizeof(t_float) * size);
+   *mtx=0;
+   *array=0;
 }
 
+
 static void deleteMTXConv (MTXConv *mtx_conv_obj) 
 {
-   deleteTFloatMatrix (mtx_conv_obj->k, mtx_conv_obj->rows_k, mtx_conv_obj->columns_k);
-   deleteTFloatMatrix (mtx_conv_obj->x, mtx_conv_obj->rows, mtx_conv_obj->columns);
-   deleteTFloatMatrix (mtx_conv_obj->y, mtx_conv_obj->rows_y, mtx_conv_obj->columns_y);
+   deleteTFloatMatrix (mtx_conv_obj->rows_k, mtx_conv_obj->columns_k, &mtx_conv_obj->k, &mtx_conv_obj->k_array);
+   deleteTFloatMatrix (mtx_conv_obj->rows, mtx_conv_obj->columns, &mtx_conv_obj->x, &mtx_conv_obj->x_array);
+   deleteTFloatMatrix (mtx_conv_obj->rows_y, mtx_conv_obj->columns_y, &mtx_conv_obj->y, &mtx_conv_obj->y_array);
    if (mtx_conv_obj->list)
       freebytes (mtx_conv_obj->list, sizeof(t_float) * (mtx_conv_obj->size_y + 2));
          
-   mtx_conv_obj->k = 0;
-   mtx_conv_obj->x = 0;
-   mtx_conv_obj->y = 0;
    mtx_conv_obj->list = 0;
 }
 
@@ -162,13 +135,17 @@ static void readFloatFromList (int n, t_atom *l, t_float *f)
 }
 static void readMatrixFromList (int rows, int columns, t_atom *l, t_float **mtx) 
 {
-   for (;rows--; l+=columns)
-      readFloatFromList (columns, l, *mtx++);
+   int n,m;
+   for (n=0;n<rows; n++)
+      for (m=0;m<columns; m++)
+	 mtx[n][m]=atom_getfloat (l++);
 }
 static void writeMatrixIntoList (int rows, int columns, t_atom *l, t_float **mtx)
 {
-   for (;rows--; l+=columns)
-      writeFloatIntoList (columns, l, *mtx++);
+   int n,m;
+   for (n=0;n<rows; n++)
+      for (m=0;m<columns; m++, l++)
+	 SETFLOAT(l,mtx[n][m]);
 }
 
 static void mTXConvKernelMatrix (MTXConv *mtx_conv_obj, t_symbol *s, int argc,
@@ -178,30 +155,35 @@ static void mTXConvKernelMatrix (MTXConv *mtx_conv_obj, t_symbol *s, int argc,
    int columns_k = atom_getint (argv++);
    int in_size = argc-2;
    int size_k = rows_k * columns_k;
-   t_float **k = mtx_conv_obj->k;
 
-   if (!size_k) 
+   if (!size_k) {
       post ("mtx_conv: invalid matrix dimensions!");
-   else if (in_size < size_k)
+      return;
+   }
+   if (in_size < size_k) {
       post("mtx_conv: sparse matrix not yet supported: use \"mtx_check\"");
-   else if ((rows_k != mtx_conv_obj->rows_k) || (columns_k != mtx_conv_obj->columns_k)) {
-      if (k)
-	 k = resizeTFloatMatrix (k, mtx_conv_obj->rows_k, mtx_conv_obj->columns_k,
-	       rows_k, columns_k);
-      else
-	 k = getTFloatMatrix (rows_k, columns_k);
+      return;
+   }
 
+   if ((rows_k != mtx_conv_obj->rows_k) || (columns_k != mtx_conv_obj->columns_k)) {
+      if (mtx_conv_obj->k)
+	 deleteTFloatMatrix (mtx_conv_obj->rows_k, mtx_conv_obj->columns_k,
+	       &mtx_conv_obj->k, &mtx_conv_obj->k_array);
+      post("getting %d rows and %d columns for k",rows_k,columns_k);
+      getTFloatMatrix (rows_k, columns_k, &mtx_conv_obj->k, &mtx_conv_obj->k_array);
+      if ((!mtx_conv_obj->k)||(!mtx_conv_obj->k_array)) {
+	 post("mtx_conv: memory allocation failed!");
+	 return;
+      }
       mtx_conv_obj->rows_k = rows_k;
       mtx_conv_obj->columns_k = columns_k;
-
-      readMatrixFromList (rows_k, columns_k, argv, k);
-
-      mtx_conv_obj->k = k;
       mtx_conv_obj->size_k = size_k;
+
    }
-   else
-      readMatrixFromList (rows_k, columns_k, argv, k);
+
+   readMatrixFromList (rows_k, columns_k, argv, mtx_conv_obj->k);
 }
+
 static void convolveRow (int columns, int columns_c, t_float *x, t_float *c, t_float *y)
 {
    int n,k,count;
@@ -210,14 +192,20 @@ static void convolveRow (int columns, int columns_c, t_float *x, t_float *c, t_f
 	 y[n] += x[n-k] * c[k];
 }
 
-static void convolveMtx (int rows, int columns, int rows_c, int columns_c, 
-      t_float **x, t_float **c, t_float **y)
+static void convolveMtx (int rows, int columns, int rows_h, int columns_h, 
+      t_float **x, t_float **h, t_float **y)
 {
-   int n,k,count;
-   zeroTFloatMatrix (y, rows+rows_c-1, columns+columns_c-1);
-   for (k = 0; k < rows_c; k++)
-      for (n = k, count = rows; count--; n++) 
-	 convolveRow (columns, columns_c, x[n-k], c[k], y[n]);
+   int n,m,k,l;
+   int rows_y=rows+rows_h-1;
+   int cols_y=columns+columns_h-1;
+   zeroTFloatMatrix (y, rows_y, cols_y);
+
+
+   for (n=0; n<rows_y; n++)
+      for (m=0; m<cols_y; m++) 
+	 for (k=n; (k<rows_h)&&(k-n<rows); k++)
+	    for (l=m; (l<columns_h)&&(l-m<columns); l++)
+	       y[n][m]+=x[k-n][l-m]*h[k][l];
 }
 
 
@@ -235,9 +223,6 @@ static void mTXConvMatrix (MTXConv *mtx_conv_obj, t_symbol *s,
    int columns_y;
    int size_y = mtx_conv_obj->size_y;
    t_atom *list_ptr = mtx_conv_obj->list;
-   t_float **x = mtx_conv_obj->x;
-   t_float **y = mtx_conv_obj->y;
-   t_float **k = mtx_conv_obj->k;
 
    // fftsize check
    if (!size){
@@ -254,12 +239,15 @@ static void mTXConvMatrix (MTXConv *mtx_conv_obj, t_symbol *s,
    //   post("1");
 
    if ((mtx_conv_obj->rows != rows)||(mtx_conv_obj->columns != columns)) { 
-     if (x)
-       x = resizeTFloatMatrix (x, mtx_conv_obj->rows, mtx_conv_obj->columns,
-                               rows, columns);
-     else
-       x = getTFloatMatrix (rows, columns);
-     mtx_conv_obj->x = x;
+     if (mtx_conv_obj->x)
+       deleteTFloatMatrix (mtx_conv_obj->rows, mtx_conv_obj->columns, 
+	     &mtx_conv_obj->x, &mtx_conv_obj->x_array);
+     post("getting %d rows and %d columns for x",rows,columns);
+     getTFloatMatrix (rows, columns, &mtx_conv_obj->x, &mtx_conv_obj->x_array);
+     if ((!mtx_conv_obj->x)||(!mtx_conv_obj->x_array)) {
+	post("mtx_conv: memory allocation failed!");
+	return;
+     }
      mtx_conv_obj->size = size;
      mtx_conv_obj->rows = rows;
      mtx_conv_obj->columns = columns;
@@ -269,29 +257,39 @@ static void mTXConvMatrix (MTXConv *mtx_conv_obj, t_symbol *s,
    columns_y = columns+columns_k-1;
    if ((mtx_conv_obj->rows_y != rows_y)||(mtx_conv_obj->columns_y != columns_y)) { 
      size_y = rows_y * columns_y;      
-     if (y)
-       y = resizeTFloatMatrix (y, mtx_conv_obj->rows_y, mtx_conv_obj->columns_y,
-                               rows_y, columns_y);
-     else
-       y = getTFloatMatrix (rows_y, columns_y);
-     if (list_ptr)
-       list_ptr = (t_atom *) resizebytes (list_ptr, sizeof(t_atom) * (mtx_conv_obj->size_y+2),
-                                          sizeof (t_atom) * (size_y+2));
-     else
-       list_ptr = (t_atom *) getbytes (sizeof (t_atom) * (size_y+2));
+     if (mtx_conv_obj->y)
+        deleteTFloatMatrix (mtx_conv_obj->rows_y, mtx_conv_obj->columns_y,
+                               &mtx_conv_obj->y, &mtx_conv_obj->y_array);
+     post("getting %d rows and %d columns for y",rows_y,columns_y);
+     getTFloatMatrix (rows_y, columns_y, &mtx_conv_obj->y, &mtx_conv_obj->y_array);
+     if ((!mtx_conv_obj->y)||(!mtx_conv_obj->y_array)) {
+	post("mtx_conv: memory allocation failed!");
+	return;
+     }
      mtx_conv_obj->size_y = size_y;
      mtx_conv_obj->rows_y = rows_y;
      mtx_conv_obj->columns_y = columns_y;
-     mtx_conv_obj->y = y;
+
+     if (list_ptr)
+	list_ptr = (t_atom *) resizebytes (list_ptr, sizeof(t_atom) * (mtx_conv_obj->size_y+2),
+                                          sizeof (t_atom) * (size_y+2));
+     else
+	list_ptr = (t_atom *) getbytes (sizeof (t_atom) * (size_y+2));
      mtx_conv_obj->list = list_ptr;
+     if (!list_ptr) {
+	post("mtx_conv: memory allocation failed!");
+	return;
+     }
+
    }
    //post("3");
    // main part
-   readMatrixFromList (rows, columns, argv, x); 
+   readMatrixFromList (rows, columns, argv, mtx_conv_obj->x); 
    //post("4");
-   convolveMtx (rows, columns, rows_k, columns_k, x, k, y);
+   convolveMtx (rows, columns, rows_k, columns_k, 
+	 mtx_conv_obj->x, mtx_conv_obj->k, mtx_conv_obj->y);
    //post("5");
-   writeMatrixIntoList (rows_y, columns_y, list_ptr+2, y);
+   writeMatrixIntoList (rows_y, columns_y, list_ptr+2, mtx_conv_obj->y);
    //post("6");
    SETSYMBOL(list_ptr, gensym("matrix"));
    SETFLOAT(list_ptr, rows_y);