dcmimgle/include/dcmtk/dcmimgle/displint.h

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  */


Generated on 6 Jan 2011 for OFFIS DCMTK Version 3.6.0 by Doxygen 1.5.1