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

Linux Cross Reference
Tina6/tina-libs/tina/math/mathSpl_spline2.c

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

  1 /**********
  2  *
  3  * Copyright (c) 2003, Division of Imaging Science and Biomedical Engineering,
  4  * University of Manchester, UK.  All rights reserved.
  5  * 
  6  * Redistribution and use in source and binary forms, with or without modification, 
  7  * are permitted provided that the following conditions are met:
  8  * 
  9  *   . Redistributions of source code must retain the above copyright notice, 
 10  *     this list of conditions and the following disclaimer.
 11  *    
 12  *   . Redistributions in binary form must reproduce the above copyright notice,
 13  *     this list of conditions and the following disclaimer in the documentation 
 14  *     and/or other materials provided with the distribution.
 15  * 
 16  *   . Neither the name of the University of Manchester nor the names of its
 17  *     contributors may be used to endorse or promote products derived from this 
 18  *     software without specific prior written permission.
 19  * 
 20  * 
 21  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
 22  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
 23  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
 24  * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
 25  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
 26  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
 27  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 28  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
 29  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
 30  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
 31  * POSSIBILITY OF SUCH DAMAGE.
 32  *
 33  **********
 34  *
 35  * Program :    TINA
 36  * File    :  $Source: /home/tina/cvs/tina-libs/tina/math/mathSpl_spline2.c,v $
 37  * Date    :  $Date: 2003/09/23 11:32:20 $
 38  * Version :  $Revision: 1.5 $
 39  * CVS Id  :  $Id: mathSpl_spline2.c,v 1.5 2003/09/23 11:32:20 matts Exp $
 40  *
 41  * Author  :  Legacy TINA
 42  *
 43  * Notes : 2D uniform cubic B-spline functions.
 44  *
 45  *********
 46 */
 47 /** 
 48  *  @file
 49  *  @brief 2D uniform cubic B-spline functions. 
 50  *
 51  *  Includes functions for;
 52  *
 53  *  - allocating/copying/freeing a Spline2.
 54  *  - Spline2 I/O.
 55  *  - evaluating/interpolating/finding values on a Spline2.
 56  *  - adding/removing etc a point from a Spline2.
 57 */
 58 
 59 #include "mathSpl_spline2.h"
 60 
 61 #if HAVE_CONFIG_H
 62 #include <config.h>
 63 #endif
 64 
 65 #include <stdio.h>
 66 #include <math.h>
 67 #include <float.h>
 68 #include <tina/sys/sysDef.h>
 69 #include <tina/sys/sysPro.h>
 70 #include <tina/math/math_SplDef.h>
 71 #include <tina/math/math_GeomDef.h>
 72 #include <tina/math/math_GeomPro.h>
 73 #include <tina/math/math_NumPro.h>
 74 #include <tina/math/mathSpl_spline.h>
 75 
 76 /*
 77 #include <tina/vision.h>
 78 #include <tina/visionfuncs.h>
 79 */
 80 
 81 
 82 /*
 83 Allocate 2D uniform cubic spline with n independent control points.
 84 Type specifies end conditions:
 85 SPLINE_PERIODIC: periodic spline, period n,
 86 SPLINE_NATURAL : zero second derivative end conditions,
 87 SPLINE_TANGENT : given tangent end conditions.
 88 */
 89 Spline2        *spline2_make(int type, int n)
 90 {
 91 
 92 
 93         Spline2        *spline = talloc(Spline2);
 94         spline->type = type;
 95         spline->n = n;
 96         spline->x = spline_make(type, n);
 97         spline->y = spline_make(type, n);
 98         return (spline);
 99 }
100 
101 /*
102 As spline2_make but does not allocate control points.
103 */
104 Spline2        *spline2_alloc(int type, int n)
105 {
106         Spline2        *spline = talloc(Spline2);
107 
108         spline->type = type;
109         spline->n = n;
110         spline->x = NULL;
111         spline->y = NULL;
112         return (spline);
113 }
114 
115 void            spline2_print(FILE * pf, Spline2 * spline)
116 {
117         int             i;
118         if (spline == NULL || pf == NULL)
119                 return;
120         fprintf(pf, "type %d points %d \n", spline->type, spline->n);
121         for (i = -1; i < spline->n + 2; i++)
122                 fprintf(pf, "x %3.2f y %3.2f \n", spline->x->p[i], spline->y->p[i]);
123 }
124 
125 Spline2        *spline2_read(FILE * pf)
126 {
127         int             i, n, type;
128         Spline2        *newspline;
129         if (pf == NULL)
130                 return (NULL);
131 
132         fscanf(pf, "%*s %d %*s %d", &type, &n);
133         newspline = spline2_make(type, n);
134         for (i = -1; i < n + 2; i++)
135                 fscanf(pf, "%*s %lf %*s %lf", &(newspline->x->p[i]), &(newspline->y->p[i]));
136 
137         return (newspline);
138 }
139 
140 /*
141 Allocate and return copy of a 2D-spline.
142 */
143 Spline2        *spline2_copy(Spline2 * spline)
144 {
145         Spline2        *newspline;
146         if (spline == NULL)
147                 return (NULL);
148         newspline = talloc(Spline2);
149         newspline->type = spline->type;
150         newspline->n = spline->n;
151         newspline->x = spline_copy(spline->x);
152         newspline->y = spline_copy(spline->y);
153         return (newspline);
154 }
155 
156 /*
157 Free 2D-spline and asssociated memory.
158 */
159 void            spline2_free(Spline2 * spline)
160 {
161         if (spline == NULL)
162                 return;
163         spline_free(spline->x);
164         spline_free(spline->y);
165         rfree(spline);
166 }
167 
168 /*
169 Evaluate 2D-spline at parameter t.
170 */
171 Vec2            spline2_eval(Spline2 * spline, double t)
172 {
173         double          x = spline_eval(spline->x, t);
174         double          y = spline_eval(spline->y, t);
175         return (vec2(x, y));
176 }
177 
178 static Vec2     p00;                    /* static data! */
179 static Spline2 *spline00;               /* static data! */
180 static double   dist(double t)  /* static data! */
181 {
182         return (vec2_dist(spline2_eval(spline00, t), p00));
183 }
184 
185 static int      tol(double a, double b)
186 {
187         return (fabs(a - b) > 0.01);
188 }
189 
190 double          spline2_natural_param(Spline2 * spline, Vec2 p)
191 {
192         int             n = spline->n;
193         double          t, t1, t2, tmin = 0.0, d, dmin = FLT_MAX;
194 
195         for (t = 0.0; t < n - 1.0; t++)
196         {
197                 d = vec2_dist(spline2_eval(spline, t), p);
198                 if (dmin > d)
199                 {
200                         dmin = d;
201                         tmin = t;
202                 }
203         }
204         if (tmin == 0.0)
205         {
206                 t1 = 0.0;
207                 t2 = 1.0;
208         } else if (tmin == n - 1.0)
209         {
210                 t1 = n - 2.0;
211                 t2 = n - 1.0;
212         } else
213         {
214                 t1 = tmin - 1.0;
215                 t2 = tmin + 1.0;
216         }
217         spline00 = spline;
218         p00 = p;
219         t = golden_section(dist, t1, t2, tol, NULL, NULL);
220         return (t);
221 }
222 
223 static double   spline2_periodic_param(Spline2 * spline, Vec2 p)
224 {
225         int             n = spline->n;
226         double          t, t1, t2, tmin = 0.0, d, dmin = FLT_MAX;
227 
228         for (t = 0.0; t < n; t++)
229         {
230                 d = vec2_dist(spline2_eval(spline, t), p);
231                 if (dmin > d)
232                 {
233                         dmin = d;
234                         tmin = t;
235                 }
236         }
237         t1 = tmin - 1.0;
238         t2 = tmin + 1.0;
239         spline00 = spline;
240         p00 = p;
241         t = golden_section(dist, t1, t2, tol, NULL, NULL);
242         if (t < 0.0)
243                 t += n;
244         if (t >= n)
245                 t -= n;
246         return (t);
247 }
248 
249 /*
250 Return parameter of point on 2D-spline closest to point p.
251 */
252 double          spline2_param(Spline2 * spline, Vec2 p)
253 {
254         switch (spline->type)
255         {
256                 case SPLINE_NATURAL:
257                 case SPLINE_TANGENT:
258                 return (spline2_natural_param(spline, p));
259         case SPLINE_PERIODIC:
260                 return (spline2_periodic_param(spline, p));
261         }
262         return (0.0);
263 }
264 
265 /*
266 Return closest point on spline to p.
267 */
268 Vec2            spline2_closest(Spline2 * spline, Vec2 p)
269 {
270         double          t = spline2_param(spline, p);
271         return (spline2_eval(spline, t));
272 }
273 
274 
275 /*
276 Return distance from point p of closest point on 2D-spline.
277 */
278 double          spline2_dist(Spline2 * spline, Vec2 p)
279 {
280         double          t = spline2_param(spline, p);
281         return (vec2_dist(p, spline2_eval(spline, t)));
282 }
283 
284 /*
285 Interplate an (already allocated) 2D-spline through data values
286 q[0], ... , q[n-1].
287 */
288 void            spline2_interpolate(Spline2 * spline, Vec2 * p)
289 {
290 
291         int             i, l, u;
292         double         *x, *y;
293         switch (spline->type)
294         {
295         case SPLINE_TANGENT:
296                 l = -1;
297                 u = spline->n + 1;
298                 break;
299         case SPLINE_NATURAL:
300         case SPLINE_PERIODIC:
301                 l = 0;
302                 u = spline->n;
303                 break;
304         default:
305                 return;
306         }
307         x = tvector_alloc(l, u, double);
308         y = tvector_alloc(l, u, double);
309         for (i = l; i < u; i++)
310         {
311                 x[i] = vec2_x(p[i]);
312                 y[i] = vec2_y(p[i]);
313         }
314         spline_interpolate(spline->x, x);
315         spline_interpolate(spline->y, y);
316         tvector_free(x, l, double);
317         tvector_free(y, l, double);
318 }
319 
320 
321 /*
322 Interplate an (already allocated) 2D-spline through data points
323 stored in list points.
324 */
325 Spline2        *spline2_interpolate_list(int type, List * points)
326 {
327 
328         List           *ptr;
329         Spline2        *spline;
330         int             i, n, l, u;
331         double         *x, *y;
332         switch (type)
333         {
334         case SPLINE_TANGENT:
335                 n = list_length(points) - 2;
336                 l = -1;
337                 u = n + 1;
338                 break;
339         case SPLINE_NATURAL:
340         case SPLINE_PERIODIC:
341                 n = list_length(points);
342                 l = 0;
343                 u = n;
344                 break;
345         default:
346                 return (NULL);
347         }
348         x = tvector_alloc(l, u, double);
349         y = tvector_alloc(l, u, double);
350         spline = spline2_make(type, n);
351         for (i = l, ptr = points; i < u; i++, ptr = ptr->next)
352         {
353                 Vec2           *p = (Vec2 *) ptr->to;
354                 x[i] = vec2_x(*p);
355                 y[i] = vec2_y(*p);
356         }
357         spline_interpolate(spline->x, x);
358         spline_interpolate(spline->y, y);
359         tvector_free(x, l, double);
360         tvector_free(y, l, double);
361         return (spline);
362 }
363 
364 Spline2        *spline2_interpolate_ddlist(int type, List * points)
365 {
366 
367         List           *ptr;
368         Spline2        *spline;
369         int             i, n, l, u;
370         double         *x, *y;
371         switch (type)
372         {
373         case SPLINE_TANGENT:
374                 n = dd_list_length(points) - 2;
375                 l = -1;
376                 u = n + 1;
377                 break;
378         case SPLINE_NATURAL:
379         case SPLINE_PERIODIC:
380                 n = dd_list_length(points);
381                 l = 0;
382                 u = n;
383                 break;
384         default:
385                 return (NULL);
386         }
387         x = tvector_alloc(l, u, double);
388         y = tvector_alloc(l, u, double);
389         spline = spline2_make(type, n);
390         for (i = l, ptr = points; i < u; i++, ptr = ptr->next)
391         {
392                 Vec2           *p = (Vec2 *) ptr->to;
393                 x[i] = vec2_x(*p);
394                 y[i] = vec2_y(*p);
395         }
396         spline_interpolate(spline->x, x);
397         spline_interpolate(spline->y, y);
398         tvector_free(x, l, double);
399         tvector_free(y, l, double);
400         return (spline);
401 }
402 
403 /*
404 Makes 2D-spline interpolate a different point p at knot ip.
405 */
406 void            spline2_replace_point(Spline2 * spline, int ip, Vec2 p)
407 {
408         spline_replace_point(spline->x, ip, vec2_x(p));
409         spline_replace_point(spline->y, ip, vec2_y(p));
410 }
411 
412 /*
413 Makes 2D-spline interpolate a new point p after knot ibelow.
414 */
415 void            spline2_add_point(Spline2 * spline, int ibelow, Vec2 p)
416 {
417         spline_add_point(spline->x, ibelow, vec2_x(p));
418         spline_add_point(spline->y, ibelow, vec2_y(p));
419         spline->n++;
420 }
421 
422 /*
423 Deletes interpolation point at knot ip from 2D-spline.
424 */
425 void            spline2_delete_point(Spline2 * spline, int i)
426 {
427         if (spline->n < 3)
428                 return;
429         spline_delete_point(spline->x, i);
430         spline_delete_point(spline->y, i);
431         spline->n--;
432 }
433 
434 /*
435 Given an array of n-2D-splines in stack, interpolate a 2D-spline
436 between them at parameter 0<=t<=n-1;
437 */
438 Spline2        *spline2_stack_interp(int n, Spline2 ** spline, double t)
439 {
440         Spline        **x, **y;
441         Spline2        *s;
442         int             i;
443 
444         x = tvector_alloc(0, n, Spline *);
445         y = tvector_alloc(0, n, Spline *);
446         for (i = 0; i < n; i++)
447         {
448                 x[i] = spline[i]->x;
449                 y[i] = spline[i]->y;
450         }
451         s = spline2_alloc(spline[0]->type, spline[0]->n);
452         s->x = spline_stack_interp(n, x, t);
453         s->y = spline_stack_interp(n, y, t);
454         tvector_free(x, 0, Spline *);
455         tvector_free(y, 0, Spline *);
456 
457         return (s);
458 }
459 
460 
461 /*
462 Apply rot translate and scale to 2D spline.
463 */
464 void            spline2_rts(Spline2 * spline, Mat2 r, Vec2 t, double s)
465 {
466         int             i;
467         if (spline == NULL)
468                 return;
469         for (i = 0; i < spline->n; i++)
470         {
471                 Vec2            p = spline2_eval(spline, i);
472                 p = vec2_sum(t, vec2_times(s, mat2_vprod(r, p)));
473                 spline2_replace_point(spline, i, p);
474         }
475 }
476 

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