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

Linux Cross Reference
Tina6/tina-libs/tina/math/mathSpl_ics.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_ics.c,v $
 37  * Date    :  $Date: 2005/01/23 19:10:21 $
 38  * Version :  $Revision: 1.5 $
 39  * CVS Id  :  $Id: mathSpl_ics.c,v 1.5 2005/01/23 19:10:21 paul Exp $
 40  *
 41  * Author  :  Legacy TINA
 42  *
 43  * Notes :Interpolation Cubic Spline. Functions converted from Numeric
 44  *        Recipies to perform 1d cubic spline interpolation
 45  *
 46  *********
 47 */
 48 /** 
 49  *  @file
 50  *  @brief  Interpolation Cubic Spline. 
 51  *
 52  *  Functions converted from Numeric Recipies to perform 1d cubic spline interpolation.
 53  *  See 3rd Ed ch 3.3 pp 113-116 for more details.
 54 */
 55 
 56 #include "mathSpl_ics.h"
 57 
 58 #if HAVE_CONFIG_H
 59 #include <config.h>
 60 #endif
 61 
 62 
 63 #include <tina/sys/sysDef.h>
 64 #include <tina/sys/sysPro.h>
 65 
 66 void            spl_ci_gen(float *x, float *y, int n, double *ypf, double *ypl, float *y2)
 67 {
 68         int             i, k;
 69         float           p, qn, sig, un, *u;
 70 
 71         u = fvector_alloc(0, n - 1);    /* intermediate values */
 72 
 73         if (ypf == NULL)        /* natural splines: 0 second derivative at
 74                                  * ends */
 75                 y2[0] = u[0] = 0.0;
 76         else
 77         {
 78                 y2[0] = -0.5;
 79                 u[0] = (3.0 / (x[1] - x[0])) * ((y[1] - y[0]) / (x[1] - x[0]) - *ypf);
 80         }
 81 
 82         for (i = 1; i <= n - 2; i++)
 83         {
 84                 sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
 85                 p = sig * y2[i - 1] + 2.0;
 86                 y2[i] = (sig - 1.0) / p;
 87                 u[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]) - (y[i] - y[i - 1]) / (x[i] - x[i - 1]);
 88                 u[i] = (6.0 * u[i] / (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p;
 89         }
 90 
 91         if (ypl == NULL)
 92                 qn = un = 0.0;
 93         else
 94         {
 95                 qn = 0.5;
 96                 un = (3.0 / (x[n - 1] - x[n - 2])) * (*ypl - (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2]));
 97         }
 98         y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.0);
 99 
100         for (k = n - 2; k >= 0; k--)
101                 y2[k] = y2[k] * y2[k + 1] + u[k];
102 
103         fvector_free(u, 0);
104 }
105 
106 /* ARGSUSED Quieten Lint */
107 void            spl_ci_sbp(float *x, float *y, int n, double *ypf, double *ypl, float *y2)
108 {
109         int             i;
110         float          *u, *g;
111 
112         u = fvector_alloc(0, n);/* intermediate values */
113         g = fvector_alloc(0, n);/* intermediate values */
114 
115         u[0] = g[0] = 0.0;      /* start condition */
116 
117         for (i = 1; i <= n - 2; i++)
118         {
119                 double          a, b, c, d, temp;
120 
121                 a = x[i] - x[i - 1];
122                 b = (x[i + 1] - x[i - 1]) / 3;
123                 c = x[i + 1] - x[i];
124                 d = (y[i + 1] - y[i]) / c - (y[i] - y[i - 1]) / a;
125                 a /= 6;
126                 c /= 6;
127 
128                 temp = b - a * g[i - 1];
129                 u[i] = (d - a * u[i - 1]) / temp;
130                 g[i] = c / temp;
131         }
132 
133         u[n - 1] = g[n - 1] = 0.0;      /* end condition */
134 
135         y2[n - 1] = u[n - 1];
136         for (i = n - 2; i >= 0; i--)
137                 y2[i] = u[i] + y2[i + 1] * g[i];
138 
139         fvector_free(u, 0);
140         fvector_free(g, 0);
141 }
142 
143 double          spl_ci_val(float *xa, float *ya, float *y2a, int n, double x)
144 {
145         int             klo, khi, k;
146         float           h, b, a, y;
147 
148         klo = 0;
149         khi = n - 1;
150         while (khi - klo > 1)
151         {                       /* binary search */
152                 k = (khi + klo) >> 1;   /* mid point */
153                 if (xa[k] > x)
154                         khi = k;
155                 else
156                         klo = k;
157         }
158         h = xa[khi] - xa[klo];  /* should be strictly increasing */
159         if (h == 0.0)
160         {
161                 error("spl_ci_val: non monotonic absisa", warning);
162                 return (0.0);
163         }
164         a = (xa[khi] - x) / h;
165         b = (x - xa[klo]) / h;
166         y = a * ya[klo] + b * ya[khi] + ((a * a * a - a) * y2a[klo] + (b * b * b - b) * y2a[khi]) * (h * h) / 6.0;
167         return (y);
168 }
169 

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