00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030 #ifndef DISPLINT_H
00031 #define DISPLINT_H
00032
00033 #include "dcmtk/config/osconfig.h"
00034 #include "dcmtk/ofstd/ofcast.h"
00035
00036 #define INCLUDE_CSTDDEF
00037 #include "dcmtk/ofstd/ofstdinc.h"
00038
00039
00040
00041
00042
00043
00044
00045 #define T3_ double
00046
00047
00048
00049
00050
00051
00054 template <class T1, class T2 >
00055 class DiCubicSpline
00056 {
00057
00058 public:
00059
00074 static int Function(const T1 *x,
00075 const T2 *y,
00076 const unsigned int n,
00077 T3_ *y2,
00078 const T3_ yp1 = 1.0e30,
00079 const T3_ ypn = 1.0e30)
00080 {
00081 if ((x != NULL) && (y != NULL) && (n > 0) && (y2 != NULL))
00082 {
00083 T3_ *u = new T3_[n];
00084 if (u != NULL)
00085 {
00086 register unsigned int i;
00087 T3_ p, qn, sig, un;
00088 if (yp1 > 0.99e30)
00089 y2[0] = u[0] = 0.0;
00090 else
00091 {
00092 y2[0] = -0.5;
00093 u[0] = (3.0 / (OFstatic_cast(T3_, x[1]) - OFstatic_cast(T3_, x[0]))) *
00094 ((OFstatic_cast(T3_, y[1]) - OFstatic_cast(T3_, y[0])) /
00095 (OFstatic_cast(T3_, x[1]) - OFstatic_cast(T3_, x[0])) - yp1);
00096 }
00097 for (i = 1; i < n - 1; ++i)
00098 {
00099 sig = (OFstatic_cast(T3_, x[i]) - OFstatic_cast(T3_, x[i - 1])) /
00100 (OFstatic_cast(T3_, x[i + 1]) - OFstatic_cast(T3_, x[i - 1]));
00101 p = sig * y2[i - 1] + 2.0;
00102 y2[i] = (sig - 1.0) / p;
00103 u[i] = (OFstatic_cast(T3_, y[i + 1]) - OFstatic_cast(T3_, y[i])) /
00104 (OFstatic_cast(T3_, x[i + 1]) - OFstatic_cast(T3_, x[i])) -
00105 (OFstatic_cast(T3_, y[i]) - OFstatic_cast(T3_, y[i - 1])) /
00106 (OFstatic_cast(T3_, x[i]) - OFstatic_cast(T3_, x[i - 1]));
00107 u[i] = (6.0 * u[i] / (OFstatic_cast(T3_, x[i + 1]) -
00108 OFstatic_cast(T3_, x[i - 1])) - sig * u[i - 1]) / p;
00109 }
00110 if (ypn > 0.99e30)
00111 qn = un = 0.0;
00112 else
00113 {
00114 qn = 0.5;
00115 un = (3.0 / (OFstatic_cast(T3_, x[n - 1]) - OFstatic_cast(T3_, x[n - 2]))) *
00116 (ypn - (OFstatic_cast(T3_, y[n - 1]) - OFstatic_cast(T3_, y[n - 2])) /
00117 (OFstatic_cast(T3_, x[n - 1]) - OFstatic_cast(T3_, x[n - 2])));
00118 }
00119 y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.0);
00120 for (i = n - 1; i > 0; --i)
00121 y2[i - 1] = y2[i - 1] * y2[i] + u[i - 1];
00122 delete[] u;
00123 return 1;
00124 }
00125 }
00126 return 0;
00127 }
00128
00129
00145 static int Interpolation(const T1 *xa,
00146 const T2 *ya,
00147 const T3_ *y2a,
00148 const unsigned int na,
00149 const T1 *x,
00150 T2 *y,
00151 const unsigned int n)
00152 {
00153 if ((xa != NULL) && (ya != NULL) && (y2a != NULL) && (na > 0) && (x != NULL) && (y != NULL) && (n > 0))
00154 {
00155 register unsigned int k, i;
00156 register unsigned int klo = 0;
00157 register unsigned int khi = na - 1;
00158 T3_ h, b, a;
00159 for (i = 0; i < n; ++i)
00160 {
00161 if ((xa[klo] > x[i]) || (xa[khi] < x[i]))
00162 {
00163 klo = 0;
00164 khi = na - 1;
00165 }
00166 while (khi - klo > 1)
00167 {
00168 k = (khi + klo) >> 1;
00169 if (xa[k] > x[i])
00170 khi = k;
00171 else
00172 klo = k;
00173 }
00174 if (xa[khi] == x[i])
00175 y[i] = ya[khi];
00176 else
00177 {
00178 h = OFstatic_cast(T3_, xa[khi]) - OFstatic_cast(T3_, xa[klo]);
00179 if (h == 0.0)
00180 return 0;
00181 a = (OFstatic_cast(T3_, xa[khi]) - OFstatic_cast(T3_, x[i])) / h;
00182 b = (OFstatic_cast(T3_, x[i]) - OFstatic_cast(T3_, xa[klo])) / h;
00183 y[i] = OFstatic_cast(T2, a * OFstatic_cast(T3_, ya[klo]) + b * OFstatic_cast(T3_, ya[khi]) +
00184 ((a * a * a - a) * y2a[klo] + (b * b * b - b) * y2a[khi]) * (h * h) / 6.0);
00185 }
00186 }
00187 return 1;
00188 }
00189 return 0;
00190 }
00191 };
00192
00193
00194 #endif
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280