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

Linux Cross Reference
Tina6/tina-libs/tina/math/mathSpl_bspline.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_bspline.c,v $
 37  * Date    :  $Date: 2003/09/23 11:32:20 $
 38  * Version :  $Revision: 1.5 $
 39  * CVS Id  :  $Id: mathSpl_bspline.c,v 1.5 2003/09/23 11:32:20 matts Exp $
 40  *
 41  * Author  :  Legacy TINA
 42  *
 43  * Notes : Code to generate and interpolate 1D uniform cubic B-splines.
 44  *
 45  *********
 46 */
 47 /** 
 48  *  @file
 49  *  @brief   Code to generate and interpolate 1D uniform cubic B-splines.       
 50  *
 51  *
 52  * 
 53 */
 54 
 55 #include "mathSpl_bspline.h"
 56 
 57 #if HAVE_CONFIG_H
 58 #include <config.h>
 59 #endif
 60 
 61 
 62 #include <math.h>
 63 #include <tina/sys/sysDef.h>
 64 #include <tina/sys/sysPro.h>
 65 
 66 /*
 67 Decomposition of a cyclic-tri-diagonal n by n matrix (see hrs_solve
 68 below).
 69 Decomposition stored in previously allocated arrays l, m, k, e, f,
 70 of n doubles.
 71 */
 72 static void     hrs_tridiag(int n, double *a, double *b, double *c,
 73                       double *l, double *m, double *k, double *e, double *f)
 74 {
 75         int             i;
 76         double          s;
 77 
 78         l[0] = b[0];
 79         m[1] = a[1];
 80         k[0] = c[0] / l[0];
 81         e[0] = c[n - 1];
 82         f[0] = a[0] / l[0];
 83         s = e[0] * f[0];
 84 
 85         for (i = 1; i < n - 2; i++)
 86         {
 87                 l[i] = b[i] - m[i] * k[i - 1];
 88                 m[i + 1] = a[i + 1];
 89                 k[i] = c[i] / l[i];
 90                 e[i] = -e[i - 1] * k[i - 1];
 91                 f[i] = -f[i - 1] * m[i] / l[i];
 92                 s += e[i] * f[i];
 93         }
 94 
 95         l[n - 2] = b[n - 2] - m[n - 2] * k[n - 3];
 96         e[n - 2] = a[n - 1] - e[n - 3] * k[n - 3];
 97         f[n - 2] = (c[n - 2] - f[n - 3] * m[n - 2]) / l[n - 2];
 98         s += e[n - 2] * f[n - 2];
 99 
100         l[n - 1] = b[n - 1] - s;
101 }
102 
103 /*
104 Forward substitution stage of solution.
105 */
106 static void     hrs_fsubst(int n, double *l, double *m, double *e, double *y)
107 {
108         int             i;
109         double          s;
110         y[0] /= l[0];
111         s = e[0] * y[0];
112         for (i = 1; i < n - 1; i++)
113         {
114                 y[i] = (y[i] - m[i] * y[i - 1]) / l[i];
115                 s += e[i] * y[i];
116         }
117         y[n - 1] = (y[n - 1] - s) / l[n - 1];
118 }
119 
120 /*
121 Backward substitution stage of solution.
122 */
123 static void     hrs_bsubst(int n, double *k, double *f, double *y)
124 {
125         int             i;
126         y[n - 2] -= f[n - 2] * y[n - 1];
127         for (i = n - 3; i >= 0; i--)
128                 y[i] -= k[i] * y[i + 1] + f[i] * y[n - 1];
129 }
130 
131 /*
132 Solution of cyclic-tr-diagonal system (a b c) x = y.
133 Entries on i-th row of matrix are a[i] = lower band, b[i] = diagonal,
134 c[i] = upper band. a[0] wraps on 0-th row, c[n-1] on (n-1)-th row.
135 Solution x overwrites y.
136 */
137 void            hrs_solve(int n, double *a, double *b, double *c, double *y)
138 {
139         static int      nlast = -1;                                                                                             /* static data! */
140         static double  *l = NULL, *m = NULL, *k = NULL, *e = NULL, *f = NULL;   /* static data! */
141         if (n != nlast)
142         {
143                 /* unstable freeing fixed NAT 27/2/2000 */
144                 if (l != NULL)
145                         tvector_free(l, 0, double);
146                 if (m != NULL)
147                         tvector_free(m, 0, double);
148                 if (k != NULL)
149                         tvector_free(k, 0, double);
150                 if (e != NULL)
151                         tvector_free(e, 0, double);
152                 if (f != NULL)
153                         tvector_free(f, 0, double);
154                 l = tvector_alloc(0, n, double);
155                 m = tvector_alloc(0, n, double);
156                 k = tvector_alloc(0, n, double);
157                 e = tvector_alloc(0, n, double);
158                 f = tvector_alloc(0, n, double);
159         }
160         hrs_tridiag(n, a, b, c, l, m, k, e, f);
161         hrs_fsubst(n, l, m, e, y);
162         hrs_bsubst(n, k, f, y);
163         nlast = n;
164 }
165 

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