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

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

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

  1 /*
  2  *
  3  * spline2.c
  4  *
  5  * 2D uniform cubic B-spline functions.
  6  *
  7  */
  8 
  9 #include <values.h>
 10 #include <tina/all_tina.h>
 11 #include <tina/brain.h>
 12 #include <tina/brainfuncs.h>
 13 
 14 /*
 15 Allocate 2D uniform cubic spline with n independent control points.
 16 Type specifies end conditions:
 17 SPLINE_PERIODIC: periodic spline, period n,
 18 SPLINE_NATURAL : zero second derivative end conditions,
 19 SPLINE_TANGENT : given tangent end conditions.
 20 */
 21 Spline2 *spline2_make(int type, int n)
 22 {
 23     Spline2 *spline = talloc(Spline2);
 24     spline->type = type;
 25     spline->n = n;
 26     spline->x = spline_make(type, n);
 27     spline->y = spline_make(type, n);
 28     return (spline);
 29 }
 30 
 31 /*
 32 As spline2_make but does not allocate control points.
 33 */
 34 Spline2 *spline2_alloc(int type, int n)
 35 {
 36     Spline2 *spline = talloc(Spline2);
 37 
 38     spline->type = type;
 39     spline->n = n;
 40     spline->x = NULL;
 41     spline->y = NULL;
 42     return (spline);
 43 }
 44 
 45 void spline2_print(FILE *pf, Spline2 *spline)
 46 {
 47     int i;
 48     if (spline==NULL||pf==NULL) return;
 49     fprintf(pf,"type %d points %d \n",spline->type,spline->n);
 50     for (i=-1;i<spline->n+2;i++)
 51        fprintf(pf,"x %3.2f y %3.2f \n",spline->x->p[i],spline->y->p[i]);
 52 }
 53 
 54 Spline2 *spline2_read(FILE *pf)
 55 {
 56     int i, n, type;
 57     Spline2 *newspline;
 58     if (pf==NULL) return(NULL);
 59 
 60     fscanf(pf,"%*s %d %*s %d",&type,&n);  
 61     newspline = spline2_make(type, n);
 62     for (i=-1;i<n+2;i++)
 63        fscanf(pf,"%*s %lf %*s %lf",&(newspline->x->p[i]),&(newspline->y->p[i]));
 64 
 65     return(newspline);
 66 }
 67 
 68 /*
 69 Allocate and return copy of a 2D-spline.
 70 */
 71 Spline2 *spline2_copy(Spline2 * spline)
 72 {
 73     Spline2 *newspline;
 74     if (spline == NULL)
 75         return (NULL);
 76     newspline = talloc(Spline2);
 77     newspline->type = spline->type;
 78     newspline->n = spline->n;
 79     newspline->x = spline_copy(spline->x);
 80     newspline->y = spline_copy(spline->y);
 81     return (newspline);
 82 }
 83 
 84 /*
 85 Free 2D-spline and asssociated memory.
 86 */
 87 void spline2_free(Spline2 * spline)
 88 {
 89     if (spline == NULL)
 90         return;
 91     spline_free(spline->x);
 92     spline_free(spline->y);
 93     rfree(spline);
 94 }
 95 
 96 /*
 97 Evaluate 2D-spline at parameter t.
 98 */
 99 Vec2 spline2_eval(Spline2 * spline, double t)
100 {
101     double x = spline_eval(spline->x, t);
102     double y = spline_eval(spline->y, t);
103     return (vec2(x, y));
104 }
105 
106 static Vec2 p00;
107 static Spline2 *spline00;
108 static double dist(double t)
109 {
110     return (vec2_dist(spline2_eval(spline00, t), p00));
111 }
112 
113 static int tol(double a, double b)
114 {
115     return (fabs(a - b) > 0.01);
116 }
117 
118 double spline2_natural_param(Spline2 * spline, Vec2 p)
119 {
120     int n = spline->n;
121     double t, t1, t2, tmin = 0.0, d, dmin = MAXFLOAT;
122 
123     for (t = 0.0; t < n - 1.0; t++)
124     {
125         d = vec2_dist(spline2_eval(spline, t), p);
126         if (dmin > d)
127         {
128             dmin = d;
129             tmin = t;
130         }
131     }
132     if (tmin == 0.0)
133     {
134         t1 = 0.0;
135         t2 = 1.0;
136     } else if (tmin == n - 1.0)
137     {
138         t1 = n - 2.0;
139         t2 = n - 1.0;
140     } else
141     {
142         t1 = tmin - 1.0;
143         t2 = tmin + 1.0;
144     }
145     spline00 = spline;
146     p00 = p;
147     t = golden_section(dist, t1, t2, tol, NULL, NULL);
148     return (t);
149 }
150 
151 static double spline2_periodic_param(Spline2 * spline, Vec2 p)
152 {
153     int n = spline->n;
154     double t, t1, t2, tmin = 0.0, d, dmin = MAXFLOAT;
155 
156     for (t = 0.0; t < n; t++)
157     {
158         d = vec2_dist(spline2_eval(spline, t), p);
159         if (dmin > d)
160         {
161             dmin = d;
162             tmin = t;
163         }
164     }
165     t1 = tmin - 1.0;
166     t2 = tmin + 1.0;
167     spline00 = spline;
168     p00 = p;
169     t = golden_section(dist, t1, t2, tol, NULL, NULL);
170     if (t < 0.0)
171         t += n;
172     if (t >= n)
173         t -= n;
174     return (t);
175 }
176 
177 /*
178 Return parameter of point on 2D-spline closest to point p.
179 */
180 double spline2_param(Spline2 * spline, Vec2 p)
181 {
182     switch (spline->type)
183     {
184         case SPLINE_NATURAL:
185         case SPLINE_TANGENT:
186             return (spline2_natural_param(spline, p));
187         case SPLINE_PERIODIC:
188             return (spline2_periodic_param(spline, p));
189     }
190     return (0.0);
191 }
192 
193 /*
194 Return closest point on spline to p.
195 */
196 Vec2 spline2_closest(Spline2 *spline, Vec2 p)
197 {
198     double t = spline2_param(spline, p);
199     return(spline2_eval(spline, t));
200 }
201  
202 
203 /*
204 Return distance from point p of closest point on 2D-spline.
205 */
206 double spline2_dist(Spline2 * spline, Vec2 p)
207 {
208     double t = spline2_param(spline, p);
209     return(vec2_dist(p, spline2_eval(spline, t)));
210 }
211 
212 /*
213 Interplate an (already allocated) 2D-spline through data values
214 q[0], ... , q[n-1].
215 */
216 void spline2_interpolate(Spline2 * spline, Vec2 * p)
217 {
218 
219     int i, l, u;
220     double *x, *y;
221     switch(spline->type)
222     {
223         case SPLINE_TANGENT:
224             l = -1;
225             u = spline->n+1;
226             break;
227         case SPLINE_NATURAL:
228         case SPLINE_PERIODIC:
229             l = 0;
230             u = spline->n;
231             break;
232         default:
233             return;
234     }
235     x = tvector_alloc(l, u, double);
236     y = tvector_alloc(l, u, double);
237     for (i = l; i < u; i++)
238     {
239         x[i] = vec2_x(p[i]);
240         y[i] = vec2_y(p[i]);
241     }
242     spline_interpolate(spline->x, x);
243     spline_interpolate(spline->y, y);
244     tvector_free(x, l, double);
245     tvector_free(y, l, double);
246 }
247 
248 
249 /*
250 Interplate an (already allocated) 2D-spline through data points
251 stored in list points.
252 */
253 Spline2 *spline2_interpolate_list(int type, List * points)
254 {
255 
256     List *ptr;
257     Spline2 *spline;
258     int i, n, l, u;
259     double *x, *y;
260     switch(type)
261     {
262         case SPLINE_TANGENT:
263             n = list_length(points)-2;
264             l = -1;
265             u = n+1;
266             break;
267         case SPLINE_NATURAL:
268         case SPLINE_PERIODIC:
269             n = list_length(points);
270             l = 0;
271             u = n;
272             break;
273         default:
274             return(NULL);
275     }
276     x = tvector_alloc(l, u, double);
277     y = tvector_alloc(l, u, double);
278     spline = spline2_make(type, n);
279     for (i = l, ptr = points; i < u; i++, ptr = ptr->next)
280     {
281         Vec2 *p = (Vec2 *) ptr->to;
282         x[i] = vec2_x(*p);
283         y[i] = vec2_y(*p);
284     }
285     spline_interpolate(spline->x, x);
286     spline_interpolate(spline->y, y);
287     tvector_free(x, l, double);
288     tvector_free(y, l, double);
289     return (spline);
290 }
291 
292 Spline2 *spline2_interpolate_ddlist(int type, List * points)
293 {
294 
295     List *ptr;
296     Spline2 *spline;
297     int i, n, l, u;
298     double *x, *y;
299     switch(type)
300     {
301         case SPLINE_TANGENT:
302             n = dd_list_length(points)-2;
303             l = -1;
304             u = n+1;
305             break;
306         case SPLINE_NATURAL:
307         case SPLINE_PERIODIC:
308             n = dd_list_length(points);
309             l = 0;
310             u = n;
311             break;
312         default:
313             return(NULL);
314     }
315     x = tvector_alloc(l, u, double);
316     y = tvector_alloc(l, u, double);
317     spline = spline2_make(type, n);
318     for (i = l, ptr = points; i < u; i++, ptr = ptr->next)
319     {
320         Vec2 *p = (Vec2 *) ptr->to;
321         x[i] = vec2_x(*p);
322         y[i] = vec2_y(*p);
323     }
324     spline_interpolate(spline->x, x);
325     spline_interpolate(spline->y, y);
326     tvector_free(x, l, double);
327     tvector_free(y, l, double);
328     return (spline);
329 }
330 
331 /*
332 Makes 2D-spline interpolate a different point p at knot ip.
333 */
334 void spline2_replace_point(Spline2 * spline, int ip, Vec2 p)
335 {
336     spline_replace_point(spline->x, ip, vec2_x(p));
337     spline_replace_point(spline->y, ip, vec2_y(p));
338 }
339 
340 /*
341 Makes 2D-spline interpolate a new point p after knot ibelow.
342 */
343 void spline2_add_point(Spline2 * spline, int ibelow, Vec2 p)
344 {
345     spline_add_point(spline->x, ibelow, vec2_x(p));
346     spline_add_point(spline->y, ibelow, vec2_y(p));
347     spline->n++;
348 }
349 
350 /*
351 Deletes interpolation point at knot ip from 2D-spline.
352 */
353 void spline2_delete_point(Spline2 * spline, int i)
354 {
355     if (spline->n < 3)
356         return;
357     spline_delete_point(spline->x, i);
358     spline_delete_point(spline->y, i);
359     spline->n--;
360 }
361 
362 /*
363 Given an array of n-2D-splines in stack, interpolate a 2D-spline
364 between them at parameter 0<=t<=n-1;
365 */
366 Spline2 *spline2_stack_interp(int n, Spline2 ** spline, double t)
367 {
368     Spline **x, **y;
369     Spline2 *s;
370     int i;
371 
372     x = tvector_alloc(0, n, Spline *);
373     y = tvector_alloc(0, n, Spline *);
374     for (i = 0; i < n; i++)
375     {
376         x[i] = spline[i]->x;
377         y[i] = spline[i]->y;
378     }
379     s = spline2_alloc(spline[0]->type, spline[0]->n);
380     s->x = spline_stack_interp(n, x, t);
381     s->y = spline_stack_interp(n, y, t);
382     tvector_free(x, 0, Spline *);
383     tvector_free(y, 0, Spline *);
384 
385     return (s);
386 }
387 
388 /*
389 Return 8-neighbour connected edge string occupied by spline.
390 Points returned are sub-pixel accuracy.
391 */
392 Tstring *str2_of_spline2(Spline2 * spline)
393 {
394     double t1 = 0.0, t2;
395 
396     if (spline == NULL)
397         return (NULL);
398 
399     switch (spline->type)
400     {
401         case SPLINE_NATURAL:
402         case SPLINE_TANGENT:
403             t2 = spline->n - 1;
404             break;
405         case SPLINE_PERIODIC:
406             t2 = spline->n;
407             break;
408         default:
409             return(NULL);
410     }
411     return (str2_of_curve2(spline2_eval, spline, t1, t2, (void *) NULL));
412 }
413 
414 /*
415 Returns a spline interpolating an edge string at n evenly spaced points.
416 */
417 Spline2 *spline2_of_str2(Tstring * str, int n)
418 {
419     Spline2 *spline;
420     Vec2 *p;
421 
422     if (str == NULL)
423         return (NULL);
424 
425     switch (str->type)
426     {
427         case STRING:
428             spline = spline2_make(SPLINE_NATURAL, n);
429             break;
430         case LOOP:
431             spline = spline2_make(SPLINE_PERIODIC, n);
432             break;
433     }
434     p = tvector_alloc(0, n, Vec2);
435     str2_get_interp_vec2_knots(str, n, p);
436     spline2_interpolate(spline, p);
437     tvector_free((void *) p, 0, Vec2);
438     return (spline);
439 }
440 
441 /*
442 Returns a 2D-spline approximating an edge string to an
443 accuracy dmax. Can be fooled by shapes with symmetry,
444 since accuracy is checked only at knot-mid-points.
445 Recursive call is fairly expensive.
446 */
447 Spline2 *spline2_approx_str2(Tstring * str, int *pn, double dmax)
448 {
449     Spline2 *spline;
450     Vec2 *p;
451     double dist;
452     int i, n = *pn;
453 
454     if (str == NULL)
455         return (NULL);
456     switch (str->type)
457     {
458         case STRING:
459             spline = spline2_make(SPLINE_NATURAL, n);
460             break;
461         case LOOP:
462             spline = spline2_make(SPLINE_PERIODIC, n);
463             break;
464     }
465     p = tvector_alloc(0, 2*n, Vec2);
466     str2_get_interp_vec2_knots(str, n, p);
467     spline2_interpolate(spline, p);
468 
469     str2_get_interp_vec2_knots(str, 2*n, p);
470     dist = 0;
471     for(i = 0; i < n; i++)
472     {
473         double d = spline2_dist(spline, p[2*i+1]);
474         dist = MAX(dist, d);
475     }
476     tvector_free((void *) p, 0, Vec2);
477     if(dist > dmax && *pn < 500)
478     {
479         spline2_free(spline);
480         *pn = 1.4142*n;
481         return(spline2_approx_str2(str, pn, dmax));
482     }
483     return (spline);
484 }
485 
486 /*
487 Apply rot translate and scale to 2D spline.
488 */
489 void spline2_rts(Spline2 *spline, Mat2 r, Vec2 t, double s)
490 {
491     int i;
492     if(spline == NULL)
493         return;
494     for(i = 0; i < spline->n; i++)
495     {
496         Vec2 p = spline2_eval(spline, i);
497         p = vec2_sum(t, vec2_times(s, mat2_vprod(r, p)));
498         spline2_replace_point(spline, i, p);
499     }
500 }
501 

~ [ 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.