Update of /cvsroot/pure-data/externals/iem/iemmatrix/src In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv12857

Added Files: mtx_bspline.c Log Message: an n-dimensional b-spline interpolator

--- NEW FILE: mtx_bspline.c --- /* * 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(); }