OFFIS DCMTK  Version 3.6.0
dicrvfit.h
1 /*
2  *
3  * Copyright (C) 1996-2010, OFFIS e.V.
4  * All rights reserved. See COPYRIGHT file for details.
5  *
6  * This software and supporting documentation were developed by
7  *
8  * OFFIS e.V.
9  * R&D Division Health
10  * Escherweg 2
11  * D-26121 Oldenburg, Germany
12  *
13  *
14  * Module: dcmimgle
15  *
16  * Author: Joerg Riesmeier
17  *
18  * Purpose: DiCurveFitting (header/implementation)
19  *
20  * Last Update: $Author: joergr $
21  * Update Date: $Date: 2010-10-14 13:16:25 $
22  * CVS/RCS Revision: $Revision: 1.19 $
23  * Status: $State: Exp $
24  *
25  * CVS/RCS Log at end of file
26  *
27  */
28 
29 
30 #ifndef DICRVFIT_H
31 #define DICRVFIT_H
32 
33 #include "dcmtk/config/osconfig.h"
34 #include "dcmtk/ofstd/oftypes.h"
35 #include "dcmtk/ofstd/ofcast.h"
36 
37 #define INCLUDE_CMATH
38 #define INCLUDE_CSTDDEF /* For NULL */
39 #include "dcmtk/ofstd/ofstdinc.h"
40 
41 
42 /*---------------------*
43  * macro definitions *
44  *---------------------*/
45 
46 // SunCC 4.x does not support default values for template types :-/
47 #define T3_ double
48 
49 
50 /*------------------*
51  * template class *
52  *------------------*/
53 
56 template <class T1, class T2 /*, class T3 = double*/>
58 {
59 
60  public:
61 
75  static int calculateCoefficients(const T1 *x,
76  const T2 *y,
77  const unsigned int n,
78  const unsigned int o,
79  T3_ *c)
80  {
81  int result = 0;
82  if ((x != NULL) && (y != NULL) && (c !=NULL) && (n > 0))
83  {
84  const unsigned int order = o + 1;
85  const unsigned int order2 = order * order;
86  T3_ *basis = new T3_[order * n];
87  T3_ *alpha = new T3_[order2];
88  T3_ *beta = new T3_[order];
89  if ((basis != NULL) && (alpha != NULL) && (beta != NULL))
90  {
91  register unsigned int i;
92  register unsigned int j;
93  register unsigned int k;
94  for (i = 0; i < order; ++i)
95  {
96  for (j = 0; j < n; ++j)
97  {
98  k = i + j * order;
99  if (i == 0)
100  basis[k] = 1;
101  else
102  basis[k] = OFstatic_cast(T3_, x[j]) * basis[k - 1];
103  }
104  }
105  T3_ sum;
106  for (i = 0; i < order; ++i)
107  {
108  const unsigned int i_order = i * order;
109  for (j = 0; j <= i; ++j)
110  {
111  sum = 0;
112  for (k = 0; k < n; ++k)
113  sum += basis[i + k * order] * basis[j + k * order];
114  alpha[i + j * order] = sum;
115  if (i != j)
116  alpha[j + i_order] = sum;
117  }
118  }
119  for (i = 0; i < order; ++i)
120  {
121  sum = 0;
122  for (j = 0; j < n; ++j)
123  sum += OFstatic_cast(T3_, y[j]) * basis[i + j * order];
124  beta[i] = sum;
125  }
126  if (solve(alpha, beta, order))
127  {
128  for (i = 0; i < order; ++i)
129  c[i] = beta[i];
130  result = 1;
131  }
132  }
133  delete[] basis;
134  delete[] alpha;
135  delete[] beta;
136  }
137  return result;
138  }
139 
140 
156  static int calculateValues(const T1 xs,
157  const T1 xe,
158  T2 *y,
159  const unsigned int n,
160  const unsigned int o,
161  const T3_ *c)
162  {
163  int result = 0;
164  if ((y != NULL) && (c != NULL) && (n > 0) && (xe > xs))
165  {
166  register unsigned int i;
167  register unsigned int j;
168  T3_ x;
169  T3_ x2;
170  T3_ w;
171  const T3_ xo = OFstatic_cast(T3_, xs);
172  const T3_ xi = OFstatic_cast(T3_, (OFstatic_cast(T3_, xe) - OFstatic_cast(T3_, xs)) / (n - 1));
173  for (i = 0; i < n; ++i)
174  {
175  x = xo + OFstatic_cast(T3_, i) * xi;
176  x2 = 1;
177  w = 0;
178  for (j = 0; j <= o; ++j)
179  {
180  w += c[j] * x2;
181  x2 *= x;
182  }
183  convertValue(w, y[i]); // cut value if necessary
184  }
185  result = 1;
186  }
187  return result;
188  }
189 
190 
191  private:
192 
200  static void convertValue(const T3_ input, Uint8 &output)
201  {
202  output = (input <= 0) ? 0 : ((input >= 255) ? 255 : OFstatic_cast(Uint8, input));
203  }
204 
212  static void convertValue(const T3_ input, Sint8 &output)
213  {
214  output = (input <= -128) ? -128 : ((input >= 127) ? 127 : OFstatic_cast(Sint8, input));
215  }
216 
224  static void convertValue(const T3_ input, Uint16 &output)
225  {
226  output = (input <= 0) ? 0 : ((input >= 65535) ? 65535 : OFstatic_cast(Uint16, input));
227  }
228 
236  static void convertValue(const T3_ input, Sint16 &output)
237  {
238  output = (input <= -32768) ? -32768 : ((input >= 32767) ? 32767 : OFstatic_cast(Sint16, input));
239  }
240 
248  static inline void convertValue(const T3_ input, double &output)
249  {
250  output = OFstatic_cast(double, input);
251  }
252 
262  static int solve(T3_ *a,
263  T3_ *b,
264  const unsigned int n)
265  {
266  int result = 0;
267  if ((a != NULL) && (b != NULL) && (n > 0))
268  {
269  register unsigned int i;
270  register unsigned int j;
271  register unsigned int k;
272  signed int pivot;
273  T3_ mag;
274  T3_ mag2;
275  T3_ temp;
276  for (i = 0; i < n; ++i)
277  {
278  mag = 0;
279  pivot = -1;
280  for (j = i; j < n; ++j)
281  {
282  mag2 = fabs(a[i + j * n]);
283  if (mag2 > mag)
284  {
285  mag = mag2;
286  pivot = j;
287  }
288  }
289  if ((pivot == -1) || (mag == 0))
290  break;
291  else
292  {
293  const unsigned int piv = OFstatic_cast(unsigned int, pivot);
294  const unsigned int i_n = i * n;
295  if (piv != i)
296  {
297  const unsigned int piv_n = piv * n;
298  for (j = i; j < n; ++j)
299  {
300  temp = a[j + i_n];
301  a[j + i_n] = a[j + piv_n];
302  a[j + piv_n] = temp;
303  }
304  temp = b[i];
305  b[i] = b[piv];
306  b[piv] = temp;
307  }
308  mag = a[i + i_n];
309  for (j = i; j < n; ++j)
310  a[j + i_n] /= mag;
311  b[i] /= mag;
312  for (j = 0; j < n; ++j)
313  {
314  if (i == j)
315  continue;
316  const unsigned int j_n = j * n;
317  mag2 = a[i + j_n];
318  for (k = i; k < n; ++k)
319  a[k + j_n] -= mag2 * a[k + i_n];
320  b[j] -= mag2 * b[i];
321  }
322  result = 1;
323  }
324 
325  }
326  }
327  return result;
328  }
329 };
330 
331 
332 #endif
333 
334 
335 /*
336  *
337  * CVS/RCS Log:
338  * $Log: dicrvfit.h,v $
339  * Revision 1.19 2010-10-14 13:16:25 joergr
340  * Updated copyright header. Added reference to COPYRIGHT file.
341  *
342  * Revision 1.18 2010-03-01 09:08:46 uli
343  * Removed some unnecessary include directives in the headers.
344  *
345  * Revision 1.17 2005-12-08 16:47:35 meichel
346  * Changed include path schema for all DCMTK header files
347  *
348  * Revision 1.16 2003/12/23 15:53:22 joergr
349  * Replaced post-increment/decrement operators by pre-increment/decrement
350  * operators where appropriate (e.g. 'i++' by '++i').
351  *
352  * Revision 1.15 2003/12/08 18:54:16 joergr
353  * Adapted type casts to new-style typecast operators defined in ofcast.h.
354  * Removed leading underscore characters from preprocessor symbols (reserved
355  * symbols). Updated copyright header.
356  *
357  * Revision 1.14 2002/11/27 14:08:03 meichel
358  * Adapted module dcmimgle to use of new header file ofstdinc.h
359  *
360  * Revision 1.13 2002/11/26 18:18:35 joergr
361  * Replaced include for "math.h" with <math.h> to avoid inclusion of math.h in
362  * the makefile dependencies.
363  *
364  * Revision 1.12 2002/10/31 10:10:45 meichel
365  * Added workaround for a bug in the Sparc optimizer in gcc 3.2
366  *
367  * Revision 1.11 2002/07/19 08:23:12 joergr
368  * Added missing doc++ comments.
369  *
370  * Revision 1.10 2002/07/18 12:28:11 joergr
371  * Added explicit type casts to keep Sun CC 2.0.1 quiet.
372  *
373  * Revision 1.9 2001/06/01 15:49:40 meichel
374  * Updated copyright header
375  *
376  * Revision 1.8 2000/03/08 16:24:14 meichel
377  * Updated copyright header.
378  *
379  * Revision 1.7 2000/03/06 15:58:39 meichel
380  * Changed static template functions to methods. Required for xlC 1.0 on AIX 3.2.
381  *
382  * Revision 1.6 1999/10/21 08:29:41 joergr
383  * Renamed template type definition from 'T3' to 'T3_' to avoid naming
384  * conflicts.
385  *
386  * Revision 1.5 1999/10/20 18:38:49 joergr
387  * Eliminated default values for template types since this features is not
388  * supported by SunCC 4.x (temporarily introduced '#define' instead).
389  *
390  * Revision 1.4 1999/10/20 10:32:44 joergr
391  * Added generic specification for template function convertValue to avoid
392  * compiler warnings reported by MSVC (with additional options?).
393  *
394  * Revision 1.3 1999/10/18 10:15:15 joergr
395  * Fixed typos.
396  *
397  * Revision 1.2 1999/10/15 09:38:31 joergr
398  * Fixed typos.
399  *
400  * Revision 1.1 1999/10/14 19:04:09 joergr
401  * Added new template class that supports polynomial curve fitting algorithm.
402  *
403  *
404  */


Generated on Thu Dec 20 2012 for OFFIS DCMTK Version 3.6.0 by Doxygen 1.8.2