~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina4/src/covira/splinesurf.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 /*
  2  *
  3  * splinesurf.c
  4  *
  5  * Uniform cubic B-spline surface z=z(s,t) functions.
  6  *
  7  */
  8 
  9 #include <tina/all_tina.h>
 10 #include <tina/brain.h>
 11 #include <tina/brainfuncs.h>
 12 
 13 /*
 14 Allocate and return a surface spline with nu by nv independent knots.
 15 utype and vtype specify spline types (SPLINE_NATURAL, SPLINE_PERIODIC, ...)
 16 in u and v parameters, see spline_make.
 17 */
 18 Splinesurf *splinesurf_make(int utype, int vtype, int nu, int nv)
 19 {
 20     Splinesurf *ss = talloc(Splinesurf);
 21     ss->utype = utype;
 22     ss->vtype = vtype;
 23     ss->nu = nu;
 24     ss->nv = nv;
 25     ss->p = tarray_alloc(-1, -1, nu+2, nv+2, double);
 26     return (ss);
 27 }
 28 
 29 /*
 30 Allocate and return a copy of a surface spline.
 31 */
 32 Splinesurf *splinesurf_copy(Splinesurf * ss)
 33 {
 34     int nu, nv;
 35     Splinesurf *newss;
 36     if (ss == NULL)
 37         return(NULL);
 38     newss = talloc(Splinesurf);
 39     newss->utype = ss->utype;
 40     newss->vtype = ss->vtype;
 41     nu = newss->nu = ss->nu;
 42     nv = newss->nv = ss->nv;
 43     newss->p = tarray_copy(ss->p, -1, -1, nu+2, nv+2, double);
 44     return (newss);
 45 }
 46 
 47 /*
 48 Free a spline surface and associated memory.
 49 */
 50 void splinesurf_free(Splinesurf * ss)
 51 {
 52     if (ss == NULL)
 53         return;
 54     tarray_free(ss->p, -1, -1, ss->nu+2, ss->nv+2, double);
 55     rfree(ss);
 56 }
 57 
 58 /*
 59 Basic B-spline patch evaluation function, 0<=u,v<=1.
 60 */
 61 static double bsplinesurf(double **c, int i, int j, double u, double v)
 62 {
 63     double u2 = u*u, u3 = u*u*u;
 64     double v2 = v*v, v3 = v*v*v;
 65     double bu0 = u3/6;
 66     double bu1 = (1+3*u+3*u2-3*u3)/6;
 67     double bu2 = (4-6*u2+3*u3)/6;
 68     double bu3 = (1-3*u+3*u2-u3)/6;
 69     double bv0 = v3/6;
 70     double bv1 = (1+3*v+3*v2-3*v3)/6;
 71     double bv2 = (4-6*v2+3*v3)/6;
 72     double bv3 = (1-3*v+3*v2-v3)/6;
 73     return((c[i-1][j-1]*bv3+c[i-1][j]*bv2+c[i-1][j+1]*bv1+c[i-1][j+2]*bv0)*bu3+
 74            (c[i  ][j-1]*bv3+c[i  ][j]*bv2+c[i  ][j+1]*bv1+c[i  ][j+2]*bv0)*bu2+
 75            (c[i+1][j-1]*bv3+c[i+1][j]*bv2+c[i+1][j+1]*bv1+c[i+1][j+2]*bv0)*bu1+
 76            (c[i+2][j-1]*bv3+c[i+2][j]*bv2+c[i+2][j+1]*bv1+c[i+2][j+2]*bv0)*bu0);
 77 }
 78 
 79 /*
 80 Basic B-spline patch u-derivative evaluation function, 0<=u,v<=1.
 81 */
 82 static double bsplinesurf_du(double **c, int i, int j, double u, double v)
 83 {
 84     double u2 = u*u;
 85     double v2 = v*v, v3 = v*v*v;
 86     double bu0 = u2/2;
 87     double bu1 = (1+2*u-3*u2)/2;
 88     double bu2 = (-4*u+3*u2)/2;
 89     double bu3 = (-1+2*u-u2)/2;
 90     double bv0 = v3/6;
 91     double bv1 = (1+3*v+3*v2-3*v3)/6;
 92     double bv2 = (4-6*v2+3*v3)/6;
 93     double bv3 = (1-3*v+3*v2-v3)/6;
 94     return((c[i-1][j-1]*bv3+c[i-1][j]*bv2+c[i-1][j+1]*bv1+c[i-1][j+2]*bv0)*bu3+
 95            (c[i  ][j-1]*bv3+c[i  ][j]*bv2+c[i  ][j+1]*bv1+c[i  ][j+2]*bv0)*bu2+
 96            (c[i+1][j-1]*bv3+c[i+1][j]*bv2+c[i+1][j+1]*bv1+c[i+1][j+2]*bv0)*bu1+
 97            (c[i+2][j-1]*bv3+c[i+2][j]*bv2+c[i+2][j+1]*bv1+c[i+2][j+2]*bv0)*bu0);
 98 }
 99 
100 /*
101 Basic B-spline patch v-derivative evaluation function, 0<=u,v<=1.
102 */
103 static double bsplinesurf_dv(double **c, int i, int j, double u, double v)
104 {
105     double u2 = u*u, u3 = u*u*u;
106     double v2 = v*u;
107     double bu0 = u3/6;
108     double bu1 = (1+3*u+3*u2-3*u3)/6;
109     double bu2 = (4-6*u2+3*u3)/6;
110     double bu3 = (1-3*u+3*u2-u3)/6;
111     double bv0 = v2/2;
112     double bv1 = (1+2*v-3*v2)/2;
113     double bv2 = (-4*v+3*v2)/2;
114     double bv3 = (-1+2*v-v2)/2;
115     return((c[i-1][j-1]*bv3+c[i-1][j]*bv2+c[i-1][j+1]*bv1+c[i-1][j+2]*bv0)*bu3+
116            (c[i  ][j-1]*bv3+c[i  ][j]*bv2+c[i  ][j+1]*bv1+c[i  ][j+2]*bv0)*bu2+
117            (c[i+1][j-1]*bv3+c[i+1][j]*bv2+c[i+1][j+1]*bv1+c[i+1][j+2]*bv0)*bu1+
118            (c[i+2][j-1]*bv3+c[i+2][j]*bv2+c[i+2][j+1]*bv1+c[i+2][j+2]*bv0)*bu0);
119 }
120 
121 /*
122 Split a parameter into its integer index i and fractional part u,
123 details depend on spline type.
124 */
125 static void split(int type, int n, double s, int *i, double *u)
126 {
127     switch (type)
128     {
129         case SPLINE_NATURAL:
130             if (s <= 1.0)
131                 *i = 0;
132             else if (s >= n - 1.0)
133                 *i = n - 2;
134             else
135                 *i = floor(s);
136             break;
137         case SPLINE_PERIODIC:
138             if(s < 0.0)
139                 s += ceil(-s/n)*n;
140             if(s >= n)
141                 s -= floor(s/n)*n;
142             *i = floor(s);
143             break;
144         default:
145             return;
146     }
147     *u = s-*i;
148 }
149 
150 /*
151 Return value of a spline surface at point (s, t).
152 */
153 double splinesurf_eval(Splinesurf * ss, double s, double t)
154 {
155     int i, j;
156     double u, v;
157     if (ss == NULL)
158         return(0.0);
159     split(ss->utype, ss->nu, s, &i, &u);
160     split(ss->vtype, ss->nv, t, &j, &v);
161     return(bsplinesurf(ss->p, i, j, u, v));
162 }
163 
164 /*
165 Return value of u-derivative of a spline surface at point (s, t).
166 */
167 double splinesurf_eval_du(Splinesurf * ss, double s, double t)
168 {
169     int i, j;
170     double u, v;
171     if (ss == NULL)
172         return(0.0);
173     split(ss->utype, ss->nu, s, &i, &u);
174     split(ss->vtype, ss->nv, t, &j, &v);
175     return(bsplinesurf_du(ss->p, i, j, u, v));
176 }
177 
178 /*
179 Return value of v-derivative of a spline surface at point (s, t).
180 */
181 double splinesurf_eval_dv(Splinesurf * ss, double s, double t)
182 {
183     int i, j;
184     double u, v;
185     if (ss == NULL)
186         return(0.0);
187     split(ss->utype, ss->nu, s, &i, &u);
188     split(ss->vtype, ss->nv, t, &j, &v);
189     return(bsplinesurf_dv(ss->p, i, j, u, v));
190 }
191 
192 /*
193 Interpolate an aready allocated spline through data values
194 data[0][0], ... , data[ss->nu-1][ss->nv-1].
195 */
196 void splinesurf_interpolate(Splinesurf *ss, double **data)
197 {
198     int nu = ss->nu, nv = ss->nv;
199     Spline *uspline = spline_make(ss->utype, nu);
200     Spline *vspline = spline_make(ss->vtype, nv);
201     double *udata = tvector_alloc(0, nu, double);
202     double *vdata = tvector_alloc(0, nv, double);
203     int u, v;
204 
205     for(v = 0; v < nv; v++)
206     {
207         for(u = 0; u < nu; u++)
208             udata[u] = data[u][v];
209         spline_interpolate(uspline, udata);
210         for(u = -1; u < nu+2; u++)
211             ss->p[u][v] = uspline->p[u];
212     }
213     for(u = -1; u < nu+2; u++)
214     {
215         for(v = 0; v < nv; v++)
216             vdata[v] = ss->p[u][v];
217         spline_interpolate(vspline, vdata);
218         for(v = -1; v < nv+2; v++)
219             ss->p[u][v] = vspline->p[v];
220     }
221     spline_free(uspline);
222     spline_free(vspline);
223     tvector_free(udata, 0, double);
224     tvector_free(vdata, 0, double);
225 }
226 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.