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: DiCubicSpline Function/Interpolation (Header/Implementation) 00019 * 00020 * Last Update: $Author: joergr $ 00021 * Update Date: $Date: 2010-10-14 13:16:27 $ 00022 * CVS/RCS Revision: $Revision: 1.21 $ 00023 * Status: $State: Exp $ 00024 * 00025 * CVS/RCS Log at end of file 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 /* For NULL */ 00037 #include "dcmtk/ofstd/ofstdinc.h" 00038 00039 00040 /*--------------------* 00041 * macro definition * 00042 *--------------------*/ 00043 00044 // SunCC 4.x does not support default values for template types :-/ 00045 #define T3_ double 00046 00047 00048 /*------------------* 00049 * template class * 00050 *------------------*/ 00051 00054 template <class T1, class T2 /*, class T3 = double*/> 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]; // temporary vector 00084 if (u != NULL) 00085 { 00086 register unsigned int i; 00087 T3_ p, qn, sig, un; 00088 if (yp1 > 0.99e30) // ignore value for first derivative at point 1 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) // ignore value for first derivative at point 1 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])) // optimization 00162 { 00163 klo = 0; 00164 khi = na - 1; 00165 } 00166 while (khi - klo > 1) // search right place in the table, if necessary 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]) // optimization: use known values 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) // bad xa input, values must be distinct ! 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 * CVS/RCS Log: 00200 * $Log: displint.h,v $ 00201 * Revision 1.21 2010-10-14 13:16:27 joergr 00202 * Updated copyright header. Added reference to COPYRIGHT file. 00203 * 00204 * Revision 1.20 2010-03-01 09:08:47 uli 00205 * Removed some unnecessary include directives in the headers. 00206 * 00207 * Revision 1.19 2005-12-08 16:48:10 meichel 00208 * Changed include path schema for all DCMTK header files 00209 * 00210 * Revision 1.18 2003/12/23 15:53:22 joergr 00211 * Replaced post-increment/decrement operators by pre-increment/decrement 00212 * operators where appropriate (e.g. 'i++' by '++i'). 00213 * 00214 * Revision 1.17 2003/12/08 19:20:47 joergr 00215 * Adapted type casts to new-style typecast operators defined in ofcast.h. 00216 * Removed leading underscore characters from preprocessor symbols (reserved 00217 * symbols). Updated copyright header. 00218 * 00219 * Revision 1.16 2002/07/18 12:30:59 joergr 00220 * Removed unused code. 00221 * 00222 * Revision 1.15 2001/06/01 15:49:51 meichel 00223 * Updated copyright header 00224 * 00225 * Revision 1.14 2000/03/08 16:24:24 meichel 00226 * Updated copyright header. 00227 * 00228 * Revision 1.13 2000/02/02 14:33:54 joergr 00229 * Replaced 'delete' statements by 'delete[]' for objects created with 'new[]'. 00230 * 00231 * Revision 1.12 1999/10/21 08:29:42 joergr 00232 * Renamed template type definition from 'T3' to 'T3_' to avoid naming 00233 * conflicts. 00234 * 00235 * Revision 1.11 1999/10/20 18:38:50 joergr 00236 * Eliminated default values for template types since this features is not 00237 * supported by SunCC 4.x (temporarily introduced '#define' instead). 00238 * 00239 * Revision 1.10 1999/10/15 09:38:31 joergr 00240 * Fixed typos. 00241 * 00242 * Revision 1.9 1999/10/14 19:05:17 joergr 00243 * Fixed typo. 00244 * 00245 * Revision 1.8 1999/10/01 13:25:35 joergr 00246 * Enhanced template class for cubic spline interpolation to support 00247 * non-floating point classes/types as y-coordinates. 00248 * 00249 * Revision 1.7 1999/07/23 14:11:25 joergr 00250 * Added preliminary support for 2D bi-cubic spline interpolation (currently 00251 * not used). 00252 * 00253 * Revision 1.6 1999/05/03 11:09:31 joergr 00254 * Minor code purifications to keep Sun CC 2.0.1 quiet. 00255 * 00256 * Revision 1.5 1999/04/29 13:49:08 joergr 00257 * Renamed class CubicSpline to DiCubicSpline. 00258 * 00259 * Revision 1.4 1999/03/24 17:20:26 joergr 00260 * Added/Modified comments and formatting. 00261 * 00262 * Revision 1.3 1999/03/22 08:52:43 joergr 00263 * Added/Changed comments. 00264 * 00265 * Revision 1.2 1999/02/25 16:17:16 joergr 00266 * Initialize local variables to avoid compiler warnings (reported by gcc 00267 * 2.7.2.1 on Linux). 00268 * 00269 * Revision 1.1 1999/02/11 16:36:29 joergr 00270 * Renamed file to indicate the use of templates. 00271 * 00272 * Revision 1.2 1999/02/09 14:21:54 meichel 00273 * Removed default parameters from template functions, required for Sun CC 4.2 00274 * 00275 * Revision 1.1 1999/02/04 17:59:23 joergr 00276 * Added support for calibration according to Barten transformation (incl. 00277 * a DISPLAY file describing the monitor characteristic). 00278 * 00279 * 00280 */