33 #include "dcmtk/config/osconfig.h"
34 #include "dcmtk/ofstd/ofcast.h"
36 #define INCLUDE_CSTDDEF
37 #include "dcmtk/ofstd/ofstdinc.h"
54 template <
class T1,
class T2 >
78 const T3_ yp1 = 1.0e30,
79 const T3_ ypn = 1.0e30)
81 if ((x != NULL) && (y != NULL) && (n > 0) && (y2 != NULL))
86 register unsigned int i;
93 u[0] = (3.0 / (OFstatic_cast(T3_, x[1]) - OFstatic_cast(T3_, x[0]))) *
94 ((OFstatic_cast(T3_, y[1]) - OFstatic_cast(T3_, y[0])) /
95 (OFstatic_cast(T3_, x[1]) - OFstatic_cast(T3_, x[0])) - yp1);
97 for (i = 1; i < n - 1; ++i)
99 sig = (OFstatic_cast(T3_, x[i]) - OFstatic_cast(T3_, x[i - 1])) /
100 (OFstatic_cast(T3_, x[i + 1]) - OFstatic_cast(T3_, x[i - 1]));
101 p = sig * y2[i - 1] + 2.0;
102 y2[i] = (sig - 1.0) / p;
103 u[i] = (OFstatic_cast(T3_, y[i + 1]) - OFstatic_cast(T3_, y[i])) /
104 (OFstatic_cast(T3_, x[i + 1]) - OFstatic_cast(T3_, x[i])) -
105 (OFstatic_cast(T3_, y[i]) - OFstatic_cast(T3_, y[i - 1])) /
106 (OFstatic_cast(T3_, x[i]) - OFstatic_cast(T3_, x[i - 1]));
107 u[i] = (6.0 * u[i] / (OFstatic_cast(T3_, x[i + 1]) -
108 OFstatic_cast(T3_, x[i - 1])) - sig * u[i - 1]) / p;
115 un = (3.0 / (OFstatic_cast(T3_, x[n - 1]) - OFstatic_cast(T3_, x[n - 2]))) *
116 (ypn - (OFstatic_cast(T3_, y[n - 1]) - OFstatic_cast(T3_, y[n - 2])) /
117 (OFstatic_cast(T3_, x[n - 1]) - OFstatic_cast(T3_, x[n - 2])));
119 y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.0);
120 for (i = n - 1; i > 0; --i)
121 y2[i - 1] = y2[i - 1] * y2[i] + u[i - 1];
148 const unsigned int na,
151 const unsigned int n)
153 if ((xa != NULL) && (ya != NULL) && (y2a != NULL) && (na > 0) && (x != NULL) && (y != NULL) && (n > 0))
155 register unsigned int k, i;
156 register unsigned int klo = 0;
157 register unsigned int khi = na - 1;
159 for (i = 0; i < n; ++i)
161 if ((xa[klo] > x[i]) || (xa[khi] < x[i]))
166 while (khi - klo > 1)
168 k = (khi + klo) >> 1;
178 h = OFstatic_cast(T3_, xa[khi]) - OFstatic_cast(T3_, xa[klo]);
181 a = (OFstatic_cast(T3_, xa[khi]) - OFstatic_cast(T3_, x[i])) / h;
182 b = (OFstatic_cast(T3_, x[i]) - OFstatic_cast(T3_, xa[klo])) / h;
183 y[i] = OFstatic_cast(T2, a * OFstatic_cast(T3_, ya[klo]) + b * OFstatic_cast(T3_, ya[khi]) +
184 ((a * a * a - a) * y2a[klo] + (b * b * b - b) * y2a[khi]) * (h * h) / 6.0);