From c57c7a2923d66eab7d09b6a94153d37b3e5f6d4e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?IOhannes=20m=20zm=C3=B6lnig?=
 <zmoelnig@users.sourceforge.net>
Date: Sun, 9 Apr 2006 20:50:20 +0000
Subject: an n-dimensional b-spline interpolator

svn path=/trunk/externals/iem/iemmatrix/; revision=4859
---
 src/mtx_bspline.c | 227 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 227 insertions(+)
 create mode 100644 src/mtx_bspline.c

(limited to 'src/mtx_bspline.c')

diff --git a/src/mtx_bspline.c b/src/mtx_bspline.c
new file mode 100644
index 0000000..c6d45be
--- /dev/null
+++ b/src/mtx_bspline.c
@@ -0,0 +1,227 @@
+/*
+ *  iemmatrix
+ *
+ *  objects for manipulating simple matrices
+ *  mostly refering to matlab/octave matrix functions
+ *
+ * Copyright (c) IOhannes m zmölnig, forum::für::umläute
+ * IEM, Graz, Austria
+ *
+ * For information on usage and redistribution, and for a DISCLAIMER OF ALL
+ * WARRANTIES, see the file, "LICENSE.txt," in this distribution.
+ *
+ */
+#include "iemmatrix.h"
+
+/*
+  mtx_bspline: 
+  this is only in the iemmatrix library since i have to make sure that there is an x-value
+  for each y-value; this however enforces that for each point we have to define a point in all dimensions;
+*/
+
+/* mtx_bspline */
+static t_class *mtx_bspline_class;
+
+typedef struct _mtx_spline
+{
+  t_object x_obj;
+  t_outlet *x_outlet;
+
+  int x_numpoints;
+  int x_dimension;
+  t_matrixfloat x_min, x_max;
+  t_matrixfloat*x_x;
+  t_matrixfloat**x_y, **x_u, **x_p;
+  t_atom*x_result;
+} t_mtx_spline;
+
+static void mtx_bspline_resize(t_mtx_spline *x, int cols, int dim){
+  int size=x->x_numpoints*sizeof(t_matrixfloat);
+  int i=0;
+
+  if(x->x_x)freebytes(x->x_x, size); x->x_x=0;
+
+  for(i=0; i<x->x_dimension; i++){
+    if(x->x_y&&x->x_y[i])freebytes(x->x_y[i], size); x->x_y[i]=0;
+    if(x->x_u&&x->x_u[i])freebytes(x->x_u[i], size); x->x_u[i]=0;
+    if(x->x_p&&x->x_p[i])freebytes(x->x_p[i], size); x->x_p[i]=0;
+  }
+  if(x->x_y)freebytes(x->x_y, x->x_dimension*sizeof(t_matrixfloat*)); x->x_y=0;
+  if(x->x_u)freebytes(x->x_u, x->x_dimension*sizeof(t_matrixfloat*)); x->x_u=0;
+  if(x->x_p)freebytes(x->x_p, x->x_dimension*sizeof(t_matrixfloat*)); x->x_p=0;
+
+  if(x->x_result)freebytes(x->x_result, x->x_dimension*sizeof(t_atom)); x->x_p=0;
+
+  if(dim<1)dim=1;
+  if(cols<0)cols=0;
+
+  x->x_numpoints = cols;
+  x->x_dimension = dim;
+
+  if(cols>0){
+    size=x->x_numpoints*sizeof(t_matrixfloat);
+    x->x_x = (t_matrixfloat*)getbytes(size);
+    x->x_result = (t_atom*)getbytes(x->x_dimension*sizeof(t_atom));
+
+    x->x_y = (t_matrixfloat**)getbytes(dim*sizeof(t_matrixfloat*));
+    x->x_u = (t_matrixfloat**)getbytes(dim*sizeof(t_matrixfloat*));
+    x->x_p = (t_matrixfloat**)getbytes(dim*sizeof(t_matrixfloat*));
+
+    for(i=0; i<x->x_dimension; i++){
+      x->x_y[i] = (t_matrixfloat*)getbytes(size);
+      x->x_u[i] = (t_matrixfloat*)getbytes(size);
+      x->x_p[i] = (t_matrixfloat*)getbytes(size);
+    }
+  }
+}
+
+
+static void mtx_bspline_matrix2(t_mtx_spline *X, t_symbol *s, int argc, t_atom *argv)
+{
+  int row=0;
+  int col=0;
+  t_atom *m = argv+2;
+
+  t_matrixfloat *x, **y, **u, **p, *w, *d, *fp;
+  t_matrixfloat*dummy;
+  int i,j,n;
+  int N;
+
+  if (argc<2){    error("mtx_bspline: crippled matrix");    return;  }
+
+  row=atom_getfloat(argv);
+  col=atom_getfloat(argv+1);
+  
+  if ((col<2)||(row<3)) {    error("mtx_bspline: invalid dimensions");    return;  }
+  if (col*row>argc-2){    error("sparse matrix not yet supported : use \"mtx_check\"");    return;  }
+
+  col--;
+
+  mtx_bspline_resize(X, row, col);
+
+  //for(i=0; i<col*row; i++)post("mtx[%d]=%f", i, atom_getfloat(argv+2+i));
+
+  /* 1st fill the matrix into the arrays */
+  fp=matrix2float(argv);
+
+  dummy=fp;
+  x=X->x_x;
+  y=X->x_y;
+  u=X->x_u;
+  p=X->x_p;
+
+  for(i=0; i<row; i++){
+    x[i]=*dummy++;
+    for(j=0; j<col; j++){
+      y[j][i]=*dummy++;
+    }
+  }
+  X->x_min=x[0];
+  X->x_max=x[row-1];
+
+
+  w=(t_matrixfloat*)getbytes(X->x_numpoints*sizeof(t_matrixfloat));
+  d=(t_matrixfloat*)getbytes(X->x_numpoints*sizeof(t_matrixfloat));
+
+  N=row-1;
+
+  for(j=0; j<col;j++){
+    d[0]=0.0; d[1]=0.0;
+
+    for(i=1; i<N; i++)
+      d[i] = 2*(x[i+1]-x[i-1]);
+
+    for(i=0; i<N; i++)
+      u[j][i] = x[i+1]-x[i];
+
+    for(i=1; i<N; i++)
+      w[i]=6.0*((y[j][i+1]-y[j][i])/u[j][i]
+		-(y[j][i]-y[j][i-1])/u[j][i-1]);
+
+    // now solve this tridiagonal matrix
+
+    for(i=1; i<N-1; i++)
+      {
+	w[i+1] -= w[i]*u[j][i]/d[i];
+	d[i+1] -= u[j][i]*u[j][i]/d[i];
+      }
+
+    for(i=0; i<N; i++)p[j][i]=0.0;
+
+    for(i=N-1; i>0; i--)
+      p[j][i] = (w[i]-u[j][i]*p[j][i+1])/d[i];  
+  }
+
+}
+
+
+
+static void mtx_bspline_list(t_mtx_spline *x, t_symbol *s, int argc, t_atom *argv)
+{
+  // this should output a matrix, one row for each element of this list
+}
+static void mtx_bspline_float(t_mtx_spline *X, t_float f)
+{
+  t_matrixfloat t, t3, t3i;
+  int i=0, j=0;
+  int dim=X->x_dimension;
+  t_matrixfloat *x=X->x_x, **y=X->x_y, **u=X->x_u, **p=X->x_p;
+  t_atom*result=X->x_result;
+
+  if(dim<1){
+    outlet_float(X->x_outlet, 0.f);
+    return;
+  }
+
+  if(f<X->x_min)f=X->x_min;
+  if(f>X->x_max)f=X->x_max;
+
+  while(f>x[i+1])i++;
+
+  for(j=0; j<dim; j++){
+    t_matrixfloat uji=u[j][i];
+    t_matrixfloat t=(f-x[i])/uji;
+    t_matrixfloat t3=t*t*t-t;
+    t_matrixfloat t3i=(1-t)*(1-t)*(1-t)-(1-t);
+    t_float ret=(t_float)(t*y[j][i+1]+(1-t)*y[j][i]+uji*uji*(p[j][i+1]*t3+p[j][i]*t3i)/6.0);
+    SETFLOAT(result+j, ret);
+  }
+  outlet_list(X->x_outlet, 0, dim, result);
+
+}
+static void mtx_bspline_free(t_mtx_spline *x)
+{
+  mtx_bspline_resize(x, 0, 0);
+}
+static void *mtx_bspline_new(void)
+{
+  t_mtx_spline *x = (t_mtx_spline *)pd_new(mtx_bspline_class);
+  inlet_new(&x->x_obj, &x->x_obj.ob_pd, gensym("matrix"), gensym(""));
+
+  x->x_numpoints=0;
+  x->x_dimension=0;
+
+  x->x_min=0.0; 
+  x->x_max=0.0;
+
+  x->x_x=0;
+  x->x_y=x->x_u=x->x_p=0;
+  x->x_result=0;
+
+  x->x_outlet=outlet_new(&x->x_obj, 0);
+  return(x);
+}
+
+void mtx_bspline_setup(void)
+{
+  mtx_bspline_class = class_new(gensym("mtx_bspline"), (t_newmethod)mtx_bspline_new, (t_method)mtx_bspline_free,
+                            sizeof(t_mtx_spline), 0, A_NULL);
+  //  class_addmethod(mtx_bspline_class, (t_method)mtx_bspline_matrix, gensym("list"), A_GIMME, 0);
+  class_addmethod(mtx_bspline_class, (t_method)mtx_bspline_matrix2, gensym(""), A_GIMME, 0);
+  class_addfloat (mtx_bspline_class, mtx_bspline_float);
+}
+
+void iemtx_bspline_setup(void)
+{
+  mtx_bspline_setup();
+}
-- 
cgit v1.2.1