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

Linux Cross Reference
Tina5/tina-libs/tina/math/mathSpl_spline.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_spline.c,v $
 37  * Date    :  $Date: 2005/01/23 19:10:21 $
 38  * Version :  $Revision: 1.6 $
 39  * CVS Id  :  $Id: mathSpl_spline.c,v 1.6 2005/01/23 19:10:21 paul Exp $
 40  *
 41  * Author  :  Legacy TINA
 42  *
 43  * Notes : Uniform cubic B-spline functions.
 44  *
 45  *********
 46 */
 47 /** 
 48  *  @file
 49  *  @brief   1D Uniform cubic B-spline functions.  
 50  *
 51  *  Includes functions for;
 52  *
 53  *  - allocating, copying, freeing a Spline structure.
 54  *  - evaluating, interpolating etc a Spline.
 55  *  - adding/removing a point on a Spline.
 56 */
 57 
 58 #include "mathSpl_spline.h"
 59 
 60 #if HAVE_CONFIG_H
 61 #include <config.h>
 62 #endif
 63 
 64 #include <stdio.h>
 65 #include <math.h>
 66 #include <float.h>
 67 #include <tina/sys/sysDef.h>
 68 #include <tina/sys/sysPro.h>
 69 #include <tina/math/math_SplDef.h>
 70 #include <tina/math/math_NumPro.h>
 71 #include <tina/math/mathSpl_bspline.h>
 72 
 73 
 74 /*
 75 Allocate 1D uniform cubic spline with n independent control points.
 76 Type specifies end conditions:
 77 SPLINE_PERIODIC: periodic spline, period n,
 78 SPLINE_NATURAL : zero second derivative end conditions,
 79 SPLINE_TANGENT : given tangent end conditions.
 80 */
 81 Spline         *spline_make(int type, int n)
 82 {
 83         Spline         *spline = talloc(Spline);
 84         spline->type = type;
 85         spline->n = n;
 86         spline->p = tvector_alloc(-1, n + 2, double);
 87         return (spline);
 88 }
 89 
 90 /*
 91 Allocate and return copy of a 1D-spline.
 92 */
 93 Spline         *spline_copy(Spline * spline)
 94 {
 95         int             n;
 96         Spline         *newspline;
 97         if (spline == NULL)
 98                 return (NULL);
 99         newspline = talloc(Spline);
100         newspline->type = spline->type;
101         n = newspline->n = spline->n;
102         newspline->p = tvector_copy(spline->p, -1, n + 2, double);
103         return (newspline);
104 }
105 
106 /*
107 Free 1D-spline and asssociated memory.
108 */
109 void            spline_free(Spline * spline)
110 {
111         if (spline == NULL)
112                 return;
113         tvector_free(spline->p, -1, double);
114         rfree(spline);
115 }
116 
117 /*
118 1D B-spline segment evaluation (0<u<1).
119 */
120 static double   bspline(double c0, double c1, double c2, double c3, double u)
121 {
122         double          u2 = u * u, u3 = u * u * u;
123         double          b0 = u3 / 6;
124         double          b1 = (1 + 3 * u + 3 * u2 - 3 * u3) / 6;
125         double          b2 = (4 - 6 * u2 + 3 * u3) / 6;
126         double          b3 = (1 - 3 * u + 3 * u2 - u3) / 6;
127         return (c0 * b3 + c1 * b2 + c2 * b1 + c3 * b0);
128 }
129 
130 /*
131 Evaluate 1D-spline at parameter t.
132 */
133 double          spline_eval(Spline * spline, double t)
134 {
135         int             i, n;
136         double          u, *p;
137         if (spline == NULL)
138                 return (0.0);
139         n = spline->n;
140         p = spline->p;
141         switch (spline->type)
142         {
143         case SPLINE_NATURAL:
144         case SPLINE_TANGENT:
145                 if (t <= 1.0)
146                         i = 0;
147                 else if (t >= n - 1.0)
148                         i = n - 2;
149                 else
150                         i = floor(t);
151                 u = t - i;
152                 break;
153         case SPLINE_PERIODIC:
154                 if (t < 0.0)
155                         t += ceil(-t / n) * n;
156                 if (t >= n)
157                         t -= floor(t / n) * n;
158                 i = floor(t);
159                 break;
160         default:
161                 return 0.0;
162         }
163         u = t - i;
164         return (bspline(p[i - 1], p[i], p[i + 1], p[i + 2], u));
165 }
166 
167 static double   p00;                    /* static data! */
168 static Spline  *spline00;               /* static data! */
169 static double   dist(double t)
170 {
171         return fabs(spline_eval(spline00, t) - p00);
172 }
173 
174 static int      tol(double a, double b)
175 {
176         return fabs(a - b) > 0.001;
177 }
178 
179 /*
180 Parameter of point on open 1D-spline closest to point p.
181 */
182 static double   spline_open_param(Spline * spline, double p)
183 {
184         int             n = spline->n;
185         double          t = 0.0, t1, t2, tmin=0.0, d, dmin = FLT_MAX;
186 
187         for (t = 0.0; t < n - 1.0; t++)
188         {
189                 d = fabs(spline_eval(spline, t) - p);
190                 if (dmin > d)
191                 {
192                         dmin = d;
193                         tmin = t;
194                 }
195         }
196         if (tmin == 0.0)
197         {
198                 t1 = 0.0;
199                 t2 = 1.0;
200         } else if (tmin == n - 1.0)
201         {
202                 t1 = n - 2.0;
203                 t2 = n - 1.0;
204         } else
205         {
206                 t1 = tmin - 1.0;
207                 t2 = tmin + 1.0;
208         }
209         spline00 = spline;
210         p00 = p;
211         t = golden_section(dist, t1, t2, tol, NULL, NULL);
212         return t;
213 }
214 
215 /*
216 Parameter of point on periodic 1D-spline closest to point p.
217 */
218 double          spline_periodic_param(Spline * spline, double p)
219 {
220         int             n = spline->n;
221         double          t = 0.0, t1, t2, tmin=0.0, d, dmin = FLT_MAX;
222 
223         for (t = 0.0; t < n; t++)
224         {
225                 d = fabs(spline_eval(spline, t) - p);
226                 if (dmin > d)
227                 {
228                         dmin = d;
229                         tmin = t;
230                 }
231         }
232 
233         t1 = tmin - 1.0;
234         t2 = tmin + 1.0;
235         spline00 = spline;
236         p00 = p;
237         t = golden_section(dist, t1, t2, tol, NULL, NULL);
238 
239         if (t < 0.0)
240         {
241                 t += n;
242         } else if (t >= n)
243         {
244                 t -= n;
245         }
246         return t;
247 }
248 
249 /*
250 Return parameter of point on 1D-spline closest to point p.
251 */
252 double          spline_param(Spline * spline, double p)
253 {
254         switch (spline->type)
255         {
256                 case SPLINE_NATURAL:
257                 case SPLINE_TANGENT:
258                 return (spline_open_param(spline, p));
259         case SPLINE_PERIODIC:
260                 return (spline_periodic_param(spline, p));
261         }
262         return (0.0);
263 }
264 
265 /*
266 Interplate an (already allocated) 1D-spline through data values
267 q[0], ... , q[n-1].
268 */
269 void            spline_interpolate(Spline * spline, double *q)
270 {
271         int             i, n = spline->n;
272         double         *a = tvector_alloc(0, n, double);
273         double         *b = tvector_alloc(0, n, double);
274         double         *c = tvector_alloc(0, n, double);
275         double         *p = spline->p;
276         for (i = 0; i < n; i++)
277         {
278                 p[i] = q[i];
279                 a[i] = 1.0 / 6.0;
280                 b[i] = 2.0 / 3.0;
281                 c[i] = 1.0 / 6.0;
282         }
283         switch (spline->type)
284         {
285         case SPLINE_NATURAL:
286                 a[0] = c[0] = 0;
287                 b[0] = 1.0;
288                 a[n - 1] = c[n - 1] = 0;
289                 b[n - 1] = 1.0;
290                 hrs_solve(n, a, b, c, p);
291                 p[-1] = 2 * p[0] - p[1];
292                 p[n] = 2 * p[n - 1] - p[n - 2];
293                 break;
294         case SPLINE_TANGENT:
295                 a[0] = 0;
296                 c[0] = 1.0 / 3.0;
297                 p[0] += q[-1] / 3.0;
298                 a[n - 1] = 1.0 / 3.0;
299                 c[n - 1] = 0;
300                 p[n - 1] -= q[n] / 3.0;
301                 hrs_solve(n, a, b, c, p);
302                 p[-1] = p[1] - 2 * q[-1];
303                 p[n] = p[n - 2] + 2 * q[n];
304                 break;
305         case SPLINE_PERIODIC:
306                 hrs_solve(n, a, b, c, p);
307                 p[-1] = p[n - 1];
308                 p[n] = p[0];
309                 p[n + 1] = p[1];
310                 break;
311         }
312         tvector_free(a, 0, double);
313         tvector_free(b, 0, double);
314 }
315 
316 /*
317 Makes 1D-spline interpolate a different value p at knot ip.
318 */
319 void            spline_replace_point(Spline * spline, int ip, double p)
320 {
321         int             i, n = spline->n;
322         double         *pnew = tvector_alloc(0, n, double);
323         for (i = 0; i < n; i++)
324                 pnew[i] = spline_eval(spline, i);
325         pnew[ip] = p;
326         spline_interpolate(spline, pnew);
327         tvector_free(pnew, 0, double);
328 }
329 
330 /*
331 Makes 1D-spline interpolate a new value p after knot ibelow.
332 */
333 void            spline_add_point(Spline * spline, int ibelow, double p)
334 {
335         int             i, n = spline->n;
336         double         *pnew = tvector_alloc(0, n + 1, double);
337 
338         for (i = 0; i <= ibelow; i++)
339                 pnew[i] = spline_eval(spline, i);
340         pnew[ibelow + 1] = p;
341         for (i = ibelow + 2; i < n + 1; i++)
342                 pnew[i] = spline_eval(spline, i - 1);
343 
344         tvector_free(spline->p, -1, double);
345         spline->n++;
346         spline->p = tvector_alloc(-1, n + 3, double);
347         spline_interpolate(spline, pnew);
348         tvector_free(pnew, 0, double);
349 }
350 
351 /*
352 Deletes interpolation point at knot ip from 1D-spline.
353 */
354 void            spline_delete_point(Spline * spline, int ip)
355 {
356         int             i, n = spline->n;
357         double         *pnew;
358 
359         if (spline->n < 3)
360                 return;
361         pnew = tvector_alloc(0, n - 1, double);
362         for (i = 0; i < ip; i++)
363                 pnew[i] = spline_eval(spline, i);
364         for (i = ip; i < n - 1; i++)
365                 pnew[i] = spline_eval(spline, i + 1);
366         tvector_free(spline->p, -1, double);
367         spline->n--;
368         spline->p = tvector_alloc(-1, n + 1, double);
369         spline_interpolate(spline, pnew);
370         tvector_free(pnew, 0, double);
371 }
372 
373 
374 /*
375 Given an array of n-1D-splines in stack, interpolate a 1D-spline
376 between them at parameter 0<=t<=n-1;
377 */
378 Spline         *spline_stack_interp(int n, Spline ** stack, double t)
379 {
380         Spline         *horiz, *vert;
381         double         *p, *q;
382         int             i, j, m = stack[0]->n;
383         int             type = stack[0]->type;
384 
385         p = tvector_alloc(0, n, double);
386         q = tvector_alloc(0, m, double);
387         vert = spline_make(SPLINE_NATURAL, n);
388         for (j = 0; j < m; j++)
389         {
390                 for (i = 0; i < n; i++)
391                         p[i] = stack[i]->p[j];
392                 spline_interpolate(vert, p);
393                 q[j] = spline_eval(vert, t);
394         }
395         horiz = spline_make(type, m);
396         spline_interpolate(horiz, q);
397         tvector_free(p, 0, double);
398         tvector_free(q, 0, double);
399         spline_free(vert);
400         return (horiz);
401 }
402 

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