OFFIS DCMTK: dcmimgle/include/dicrvfit.h Source File

00001 /* 00002 * 00003 * Copyright (C) 1996-2003, OFFIS 00004 * 00005 * This software and supporting documentation were developed by 00006 * 00007 * Kuratorium OFFIS e.V. 00008 * Healthcare Information and Communication Systems 00009 * Escherweg 2 00010 * D-26121 Oldenburg, Germany 00011 * 00012 * THIS SOFTWARE IS MADE AVAILABLE, AS IS, AND OFFIS MAKES NO WARRANTY 00013 * REGARDING THE SOFTWARE, ITS PERFORMANCE, ITS MERCHANTABILITY OR 00014 * FITNESS FOR ANY PARTICULAR USE, FREEDOM FROM ANY COMPUTER DISEASES OR 00015 * ITS CONFORMITY TO ANY SPECIFICATION. THE ENTIRE RISK AS TO QUALITY AND 00016 * PERFORMANCE OF THE SOFTWARE IS WITH THE USER. 00017 * 00018 * Module: dcmimgle 00019 * 00020 * Author: Joerg Riesmeier 00021 * 00022 * Purpose: DiCurveFitting (header/implementation) 00023 * 00024 * Last Update: $Author: joergr $ 00025 * Update Date: $Date: 2003/12/23 15:53:22 $ 00026 * CVS/RCS Revision: $Revision: 1.16 $ 00027 * Status: $State: Exp $ 00028 * 00029 * CVS/RCS Log at end of file 00030 * 00031 */ 00032 00033 00034 #ifndef DICRVFIT_H 00035 #define DICRVFIT_H 00036 00037 #include "osconfig.h" 00038 #include "ofcast.h" 00039 00040 #define INCLUDE_CMATH 00041 #include "ofstdinc.h" 00042 00043 00044 /*---------------------* 00045 * macro definitions * 00046 *---------------------*/ 00047 00048 // SunCC 4.x does not support default values for template types :-/ 00049 #define T3_ double 00050 00051 00052 /*------------------* 00053 * template class * 00054 *------------------*/ 00055 00058 template <class T1, class T2 /*, class T3 = double*/> 00059 class DiCurveFitting 00060 { 00061 00062 public: 00063 00077 static int calculateCoefficients(const T1 *x, 00078 const T2 *y, 00079 const unsigned int n, 00080 const unsigned int o, 00081 T3_ *c) 00082 { 00083 int result = 0; 00084 if ((x != NULL) && (y != NULL) && (c !=NULL) && (n > 0)) 00085 { 00086 const unsigned int order = o + 1; 00087 const unsigned int order2 = order * order; 00088 T3_ *basis = new T3_[order * n]; 00089 T3_ *alpha = new T3_[order2]; 00090 T3_ *beta = new T3_[order]; 00091 if ((basis != NULL) && (alpha != NULL) && (beta != NULL)) 00092 { 00093 register unsigned int i; 00094 register unsigned int j; 00095 register unsigned int k; 00096 for (i = 0; i < order; ++i) 00097 { 00098 for (j = 0; j < n; ++j) 00099 { 00100 k = i + j * order; 00101 if (i == 0) 00102 basis[k] = 1; 00103 else 00104 basis[k] = OFstatic_cast(T3_, x[j]) * basis[k - 1]; 00105 } 00106 } 00107 T3_ sum; 00108 for (i = 0; i < order; ++i) 00109 { 00110 const unsigned int i_order = i * order; 00111 for (j = 0; j <= i; ++j) 00112 { 00113 sum = 0; 00114 for (k = 0; k < n; ++k) 00115 sum += basis[i + k * order] * basis[j + k * order]; 00116 alpha[i + j * order] = sum; 00117 if (i != j) 00118 alpha[j + i_order] = sum; 00119 } 00120 } 00121 for (i = 0; i < order; ++i) 00122 { 00123 sum = 0; 00124 for (j = 0; j < n; ++j) 00125 sum += OFstatic_cast(T3_, y[j]) * basis[i + j * order]; 00126 beta[i] = sum; 00127 } 00128 if (solve(alpha, beta, order)) 00129 { 00130 for (i = 0; i < order; ++i) 00131 c[i] = beta[i]; 00132 result = 1; 00133 } 00134 } 00135 delete[] basis; 00136 delete[] alpha; 00137 delete[] beta; 00138 } 00139 return result; 00140 } 00141 00142 00158 static int calculateValues(const T1 xs, 00159 const T1 xe, 00160 T2 *y, 00161 const unsigned int n, 00162 const unsigned int o, 00163 const T3_ *c) 00164 { 00165 int result = 0; 00166 if ((y != NULL) && (c != NULL) && (n > 0) && (xe > xs)) 00167 { 00168 register unsigned int i; 00169 register unsigned int j; 00170 T3_ x; 00171 T3_ x2; 00172 T3_ w; 00173 const T3_ xo = OFstatic_cast(T3_, xs); 00174 const T3_ xi = OFstatic_cast(T3_, (OFstatic_cast(T3_, xe) - OFstatic_cast(T3_, xs)) / (n - 1)); 00175 for (i = 0; i < n; ++i) 00176 { 00177 x = xo + OFstatic_cast(T3_, i) * xi; 00178 x2 = 1; 00179 w = 0; 00180 for (j = 0; j <= o; ++j) 00181 { 00182 w += c[j] * x2; 00183 x2 *= x; 00184 } 00185 convertValue(w, y[i]); // cut value if necessary 00186 } 00187 result = 1; 00188 } 00189 return result; 00190 } 00191 00192 00193 private: 00194 00202 static void convertValue(const T3_ input, Uint8 &output) 00203 { 00204 output = (input <= 0) ? 0 : ((input >= 255) ? 255 : OFstatic_cast(Uint8, input)); 00205 } 00206 00214 static void convertValue(const T3_ input, Sint8 &output) 00215 { 00216 output = (input <= -128) ? -128 : ((input >= 127) ? 127 : OFstatic_cast(Sint8, input)); 00217 } 00218 00226 static void convertValue(const T3_ input, Uint16 &output) 00227 { 00228 output = (input <= 0) ? 0 : ((input >= 65535) ? 65535 : OFstatic_cast(Uint16, input)); 00229 } 00230 00238 static void convertValue(const T3_ input, Sint16 &output) 00239 { 00240 output = (input <= -32768) ? -32768 : ((input >= 32767) ? 32767 : OFstatic_cast(Sint16, input)); 00241 } 00242 00250 static inline void convertValue(const T3_ input, double &output) 00251 { 00252 output = OFstatic_cast(double, input); 00253 } 00254 00264 static int solve(T3_ *a, 00265 T3_ *b, 00266 const unsigned int n) 00267 { 00268 int result = 0; 00269 if ((a != NULL) && (b != NULL) && (n > 0)) 00270 { 00271 register unsigned int i; 00272 register unsigned int j; 00273 register unsigned int k; 00274 signed int pivot; 00275 T3_ mag; 00276 T3_ mag2; 00277 T3_ temp; 00278 for (i = 0; i < n; ++i) 00279 { 00280 mag = 0; 00281 pivot = -1; 00282 for (j = i; j < n; ++j) 00283 { 00284 mag2 = fabs(a[i + j * n]); 00285 if (mag2 > mag) 00286 { 00287 mag = mag2; 00288 pivot = j; 00289 } 00290 } 00291 if ((pivot == -1) || (mag == 0)) 00292 break; 00293 else 00294 { 00295 const unsigned int piv = OFstatic_cast(unsigned int, pivot); 00296 const unsigned int i_n = i * n; 00297 if (piv != i) 00298 { 00299 const unsigned int piv_n = piv * n; 00300 for (j = i; j < n; ++j) 00301 { 00302 temp = a[j + i_n]; 00303 a[j + i_n] = a[j + piv_n]; 00304 a[j + piv_n] = temp; 00305 } 00306 temp = b[i]; 00307 b[i] = b[piv]; 00308 b[piv] = temp; 00309 } 00310 mag = a[i + i_n]; 00311 for (j = i; j < n; ++j) 00312 a[j + i_n] /= mag; 00313 b[i] /= mag; 00314 for (j = 0; j < n; ++j) 00315 { 00316 if (i == j) 00317 continue; 00318 const unsigned int j_n = j * n; 00319 mag2 = a[i + j_n]; 00320 for (k = i; k < n; ++k) 00321 a[k + j_n] -= mag2 * a[k + i_n]; 00322 b[j] -= mag2 * b[i]; 00323 } 00324 result = 1; 00325 } 00326 00327 } 00328 } 00329 return result; 00330 } 00331 }; 00332 00333 00334 #endif 00335 00336 00337 /* 00338 * 00339 * CVS/RCS Log: 00340 * $Log: dicrvfit.h,v $ 00341 * Revision 1.16 2003/12/23 15:53:22 joergr 00342 * Replaced post-increment/decrement operators by pre-increment/decrement 00343 * operators where appropriate (e.g. 'i++' by '++i'). 00344 * 00345 * Revision 1.15 2003/12/08 18:54:16 joergr 00346 * Adapted type casts to new-style typecast operators defined in ofcast.h. 00347 * Removed leading underscore characters from preprocessor symbols (reserved 00348 * symbols). Updated copyright header. 00349 * 00350 * Revision 1.14 2002/11/27 14:08:03 meichel 00351 * Adapted module dcmimgle to use of new header file ofstdinc.h 00352 * 00353 * Revision 1.13 2002/11/26 18:18:35 joergr 00354 * Replaced include for "math.h" with <math.h> to avoid inclusion of math.h in 00355 * the makefile dependencies. 00356 * 00357 * Revision 1.12 2002/10/31 10:10:45 meichel 00358 * Added workaround for a bug in the Sparc optimizer in gcc 3.2 00359 * 00360 * Revision 1.11 2002/07/19 08:23:12 joergr 00361 * Added missing doc++ comments. 00362 * 00363 * Revision 1.10 2002/07/18 12:28:11 joergr 00364 * Added explicit type casts to keep Sun CC 2.0.1 quiet. 00365 * 00366 * Revision 1.9 2001/06/01 15:49:40 meichel 00367 * Updated copyright header 00368 * 00369 * Revision 1.8 2000/03/08 16:24:14 meichel 00370 * Updated copyright header. 00371 * 00372 * Revision 1.7 2000/03/06 15:58:39 meichel 00373 * Changed static template functions to methods. Required for xlC 1.0 on AIX 3.2. 00374 * 00375 * Revision 1.6 1999/10/21 08:29:41 joergr 00376 * Renamed template type definition from 'T3' to 'T3_' to avoid naming 00377 * conflicts. 00378 * 00379 * Revision 1.5 1999/10/20 18:38:49 joergr 00380 * Eliminated default values for template types since this features is not 00381 * supported by SunCC 4.x (temporarily introduced '#define' instead). 00382 * 00383 * Revision 1.4 1999/10/20 10:32:44 joergr 00384 * Added generic specification for template function convertValue to avoid 00385 * compiler warnings reported by MSVC (with additional options?). 00386 * 00387 * Revision 1.3 1999/10/18 10:15:15 joergr 00388 * Fixed typos. 00389 * 00390 * Revision 1.2 1999/10/15 09:38:31 joergr 00391 * Fixed typos. 00392 * 00393 * Revision 1.1 1999/10/14 19:04:09 joergr 00394 * Added new template class that supports polynomial curve fitting algorithm. 00395 * 00396 * 00397 */