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 DICRVFIT_H
00031 #define DICRVFIT_H
00032
00033 #include "dcmtk/config/osconfig.h"
00034 #include "dcmtk/ofstd/oftypes.h"
00035 #include "dcmtk/ofstd/ofcast.h"
00036
00037 #define INCLUDE_CMATH
00038 #define INCLUDE_CSTDDEF
00039 #include "dcmtk/ofstd/ofstdinc.h"
00040
00041
00042
00043
00044
00045
00046
00047 #define T3_ double
00048
00049
00050
00051
00052
00053
00056 template <class T1, class T2 >
00057 class DiCurveFitting
00058 {
00059
00060 public:
00061
00075 static int calculateCoefficients(const T1 *x,
00076 const T2 *y,
00077 const unsigned int n,
00078 const unsigned int o,
00079 T3_ *c)
00080 {
00081 int result = 0;
00082 if ((x != NULL) && (y != NULL) && (c !=NULL) && (n > 0))
00083 {
00084 const unsigned int order = o + 1;
00085 const unsigned int order2 = order * order;
00086 T3_ *basis = new T3_[order * n];
00087 T3_ *alpha = new T3_[order2];
00088 T3_ *beta = new T3_[order];
00089 if ((basis != NULL) && (alpha != NULL) && (beta != NULL))
00090 {
00091 register unsigned int i;
00092 register unsigned int j;
00093 register unsigned int k;
00094 for (i = 0; i < order; ++i)
00095 {
00096 for (j = 0; j < n; ++j)
00097 {
00098 k = i + j * order;
00099 if (i == 0)
00100 basis[k] = 1;
00101 else
00102 basis[k] = OFstatic_cast(T3_, x[j]) * basis[k - 1];
00103 }
00104 }
00105 T3_ sum;
00106 for (i = 0; i < order; ++i)
00107 {
00108 const unsigned int i_order = i * order;
00109 for (j = 0; j <= i; ++j)
00110 {
00111 sum = 0;
00112 for (k = 0; k < n; ++k)
00113 sum += basis[i + k * order] * basis[j + k * order];
00114 alpha[i + j * order] = sum;
00115 if (i != j)
00116 alpha[j + i_order] = sum;
00117 }
00118 }
00119 for (i = 0; i < order; ++i)
00120 {
00121 sum = 0;
00122 for (j = 0; j < n; ++j)
00123 sum += OFstatic_cast(T3_, y[j]) * basis[i + j * order];
00124 beta[i] = sum;
00125 }
00126 if (solve(alpha, beta, order))
00127 {
00128 for (i = 0; i < order; ++i)
00129 c[i] = beta[i];
00130 result = 1;
00131 }
00132 }
00133 delete[] basis;
00134 delete[] alpha;
00135 delete[] beta;
00136 }
00137 return result;
00138 }
00139
00140
00156 static int calculateValues(const T1 xs,
00157 const T1 xe,
00158 T2 *y,
00159 const unsigned int n,
00160 const unsigned int o,
00161 const T3_ *c)
00162 {
00163 int result = 0;
00164 if ((y != NULL) && (c != NULL) && (n > 0) && (xe > xs))
00165 {
00166 register unsigned int i;
00167 register unsigned int j;
00168 T3_ x;
00169 T3_ x2;
00170 T3_ w;
00171 const T3_ xo = OFstatic_cast(T3_, xs);
00172 const T3_ xi = OFstatic_cast(T3_, (OFstatic_cast(T3_, xe) - OFstatic_cast(T3_, xs)) / (n - 1));
00173 for (i = 0; i < n; ++i)
00174 {
00175 x = xo + OFstatic_cast(T3_, i) * xi;
00176 x2 = 1;
00177 w = 0;
00178 for (j = 0; j <= o; ++j)
00179 {
00180 w += c[j] * x2;
00181 x2 *= x;
00182 }
00183 convertValue(w, y[i]);
00184 }
00185 result = 1;
00186 }
00187 return result;
00188 }
00189
00190
00191 private:
00192
00200 static void convertValue(const T3_ input, Uint8 &output)
00201 {
00202 output = (input <= 0) ? 0 : ((input >= 255) ? 255 : OFstatic_cast(Uint8, input));
00203 }
00204
00212 static void convertValue(const T3_ input, Sint8 &output)
00213 {
00214 output = (input <= -128) ? -128 : ((input >= 127) ? 127 : OFstatic_cast(Sint8, input));
00215 }
00216
00224 static void convertValue(const T3_ input, Uint16 &output)
00225 {
00226 output = (input <= 0) ? 0 : ((input >= 65535) ? 65535 : OFstatic_cast(Uint16, input));
00227 }
00228
00236 static void convertValue(const T3_ input, Sint16 &output)
00237 {
00238 output = (input <= -32768) ? -32768 : ((input >= 32767) ? 32767 : OFstatic_cast(Sint16, input));
00239 }
00240
00248 static inline void convertValue(const T3_ input, double &output)
00249 {
00250 output = OFstatic_cast(double, input);
00251 }
00252
00262 static int solve(T3_ *a,
00263 T3_ *b,
00264 const unsigned int n)
00265 {
00266 int result = 0;
00267 if ((a != NULL) && (b != NULL) && (n > 0))
00268 {
00269 register unsigned int i;
00270 register unsigned int j;
00271 register unsigned int k;
00272 signed int pivot;
00273 T3_ mag;
00274 T3_ mag2;
00275 T3_ temp;
00276 for (i = 0; i < n; ++i)
00277 {
00278 mag = 0;
00279 pivot = -1;
00280 for (j = i; j < n; ++j)
00281 {
00282 mag2 = fabs(a[i + j * n]);
00283 if (mag2 > mag)
00284 {
00285 mag = mag2;
00286 pivot = j;
00287 }
00288 }
00289 if ((pivot == -1) || (mag == 0))
00290 break;
00291 else
00292 {
00293 const unsigned int piv = OFstatic_cast(unsigned int, pivot);
00294 const unsigned int i_n = i * n;
00295 if (piv != i)
00296 {
00297 const unsigned int piv_n = piv * n;
00298 for (j = i; j < n; ++j)
00299 {
00300 temp = a[j + i_n];
00301 a[j + i_n] = a[j + piv_n];
00302 a[j + piv_n] = temp;
00303 }
00304 temp = b[i];
00305 b[i] = b[piv];
00306 b[piv] = temp;
00307 }
00308 mag = a[i + i_n];
00309 for (j = i; j < n; ++j)
00310 a[j + i_n] /= mag;
00311 b[i] /= mag;
00312 for (j = 0; j < n; ++j)
00313 {
00314 if (i == j)
00315 continue;
00316 const unsigned int j_n = j * n;
00317 mag2 = a[i + j_n];
00318 for (k = i; k < n; ++k)
00319 a[k + j_n] -= mag2 * a[k + i_n];
00320 b[j] -= mag2 * b[i];
00321 }
00322 result = 1;
00323 }
00324
00325 }
00326 }
00327 return result;
00328 }
00329 };
00330
00331
00332 #endif
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404