00001 /* 00002 * 00003 * Copyright (C) 1996-2010, OFFIS e.V. 00004 * All rights reserved. See COPYRIGHT file for details. 00005 * 00006 * This software and supporting documentation were developed by 00007 * 00008 * OFFIS e.V. 00009 * R&D Division Health 00010 * Escherweg 2 00011 * D-26121 Oldenburg, Germany 00012 * 00013 * 00014 * Module: dcmimgle 00015 * 00016 * Author: Joerg Riesmeier 00017 * 00018 * Purpose: DiCurveFitting (header/implementation) 00019 * 00020 * Last Update: $Author: joergr $ 00021 * Update Date: $Date: 2010-10-14 13:16:25 $ 00022 * CVS/RCS Revision: $Revision: 1.19 $ 00023 * Status: $State: Exp $ 00024 * 00025 * CVS/RCS Log at end of file 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 /* For NULL */ 00039 #include "dcmtk/ofstd/ofstdinc.h" 00040 00041 00042 /*---------------------* 00043 * macro definitions * 00044 *---------------------*/ 00045 00046 // SunCC 4.x does not support default values for template types :-/ 00047 #define T3_ double 00048 00049 00050 /*------------------* 00051 * template class * 00052 *------------------*/ 00053 00056 template <class T1, class T2 /*, class T3 = double*/> 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]); // cut value if necessary 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 * CVS/RCS Log: 00338 * $Log: dicrvfit.h,v $ 00339 * Revision 1.19 2010-10-14 13:16:25 joergr 00340 * Updated copyright header. Added reference to COPYRIGHT file. 00341 * 00342 * Revision 1.18 2010-03-01 09:08:46 uli 00343 * Removed some unnecessary include directives in the headers. 00344 * 00345 * Revision 1.17 2005-12-08 16:47:35 meichel 00346 * Changed include path schema for all DCMTK header files 00347 * 00348 * Revision 1.16 2003/12/23 15:53:22 joergr 00349 * Replaced post-increment/decrement operators by pre-increment/decrement 00350 * operators where appropriate (e.g. 'i++' by '++i'). 00351 * 00352 * Revision 1.15 2003/12/08 18:54:16 joergr 00353 * Adapted type casts to new-style typecast operators defined in ofcast.h. 00354 * Removed leading underscore characters from preprocessor symbols (reserved 00355 * symbols). Updated copyright header. 00356 * 00357 * Revision 1.14 2002/11/27 14:08:03 meichel 00358 * Adapted module dcmimgle to use of new header file ofstdinc.h 00359 * 00360 * Revision 1.13 2002/11/26 18:18:35 joergr 00361 * Replaced include for "math.h" with <math.h> to avoid inclusion of math.h in 00362 * the makefile dependencies. 00363 * 00364 * Revision 1.12 2002/10/31 10:10:45 meichel 00365 * Added workaround for a bug in the Sparc optimizer in gcc 3.2 00366 * 00367 * Revision 1.11 2002/07/19 08:23:12 joergr 00368 * Added missing doc++ comments. 00369 * 00370 * Revision 1.10 2002/07/18 12:28:11 joergr 00371 * Added explicit type casts to keep Sun CC 2.0.1 quiet. 00372 * 00373 * Revision 1.9 2001/06/01 15:49:40 meichel 00374 * Updated copyright header 00375 * 00376 * Revision 1.8 2000/03/08 16:24:14 meichel 00377 * Updated copyright header. 00378 * 00379 * Revision 1.7 2000/03/06 15:58:39 meichel 00380 * Changed static template functions to methods. Required for xlC 1.0 on AIX 3.2. 00381 * 00382 * Revision 1.6 1999/10/21 08:29:41 joergr 00383 * Renamed template type definition from 'T3' to 'T3_' to avoid naming 00384 * conflicts. 00385 * 00386 * Revision 1.5 1999/10/20 18:38:49 joergr 00387 * Eliminated default values for template types since this features is not 00388 * supported by SunCC 4.x (temporarily introduced '#define' instead). 00389 * 00390 * Revision 1.4 1999/10/20 10:32:44 joergr 00391 * Added generic specification for template function convertValue to avoid 00392 * compiler warnings reported by MSVC (with additional options?). 00393 * 00394 * Revision 1.3 1999/10/18 10:15:15 joergr 00395 * Fixed typos. 00396 * 00397 * Revision 1.2 1999/10/15 09:38:31 joergr 00398 * Fixed typos. 00399 * 00400 * Revision 1.1 1999/10/14 19:04:09 joergr 00401 * Added new template class that supports polynomial curve fitting algorithm. 00402 * 00403 * 00404 */