aboutsummaryrefslogtreecommitdiff
path: root/src/mtx_bspline.c
diff options
context:
space:
mode:
authorIOhannes m zmölnig <zmoelnig@users.sourceforge.net>2006-04-09 20:50:20 +0000
committerIOhannes m zmölnig <zmoelnig@users.sourceforge.net>2006-04-09 20:50:20 +0000
commitc57c7a2923d66eab7d09b6a94153d37b3e5f6d4e (patch)
tree825ce2647bf6ac74b8fc672c698c5b3367d49b7a /src/mtx_bspline.c
parentef1e7c91dba43cd215ea1aa98191d8326997a461 (diff)
an n-dimensional b-spline interpolator
svn path=/trunk/externals/iem/iemmatrix/; revision=4859
Diffstat (limited to 'src/mtx_bspline.c')
-rw-r--r--src/mtx_bspline.c227
1 files changed, 227 insertions, 0 deletions
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 zmlnig, forum::fr::umlute
+ * 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();
+}