dcmimgle/include/dcmtk/dcmimgle/dicrvfit.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: 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  */


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