OFFIS DCMTK  Version 3.6.0
displint.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: DiCubicSpline Function/Interpolation (Header/Implementation)
19  *
20  * Last Update: $Author: joergr $
21  * Update Date: $Date: 2010-10-14 13:16:27 $
22  * CVS/RCS Revision: $Revision: 1.21 $
23  * Status: $State: Exp $
24  *
25  * CVS/RCS Log at end of file
26  *
27  */
28 
29 
30 #ifndef DISPLINT_H
31 #define DISPLINT_H
32 
33 #include "dcmtk/config/osconfig.h"
34 #include "dcmtk/ofstd/ofcast.h"
35 
36 #define INCLUDE_CSTDDEF /* For NULL */
37 #include "dcmtk/ofstd/ofstdinc.h"
38 
39 
40 /*--------------------*
41  * macro definition *
42  *--------------------*/
43 
44 // SunCC 4.x does not support default values for template types :-/
45 #define T3_ double
46 
47 
48 /*------------------*
49  * template class *
50  *------------------*/
51 
54 template <class T1, class T2 /*, class T3 = double*/>
56 {
57 
58  public:
59 
74  static int Function(const T1 *x,
75  const T2 *y,
76  const unsigned int n,
77  T3_ *y2,
78  const T3_ yp1 = 1.0e30,
79  const T3_ ypn = 1.0e30)
80  {
81  if ((x != NULL) && (y != NULL) && (n > 0) && (y2 != NULL))
82  {
83  T3_ *u = new T3_[n]; // temporary vector
84  if (u != NULL)
85  {
86  register unsigned int i;
87  T3_ p, qn, sig, un;
88  if (yp1 > 0.99e30) // ignore value for first derivative at point 1
89  y2[0] = u[0] = 0.0;
90  else
91  {
92  y2[0] = -0.5;
93  u[0] = (3.0 / (OFstatic_cast(T3_, x[1]) - OFstatic_cast(T3_, x[0]))) *
94  ((OFstatic_cast(T3_, y[1]) - OFstatic_cast(T3_, y[0])) /
95  (OFstatic_cast(T3_, x[1]) - OFstatic_cast(T3_, x[0])) - yp1);
96  }
97  for (i = 1; i < n - 1; ++i)
98  {
99  sig = (OFstatic_cast(T3_, x[i]) - OFstatic_cast(T3_, x[i - 1])) /
100  (OFstatic_cast(T3_, x[i + 1]) - OFstatic_cast(T3_, x[i - 1]));
101  p = sig * y2[i - 1] + 2.0;
102  y2[i] = (sig - 1.0) / p;
103  u[i] = (OFstatic_cast(T3_, y[i + 1]) - OFstatic_cast(T3_, y[i])) /
104  (OFstatic_cast(T3_, x[i + 1]) - OFstatic_cast(T3_, x[i])) -
105  (OFstatic_cast(T3_, y[i]) - OFstatic_cast(T3_, y[i - 1])) /
106  (OFstatic_cast(T3_, x[i]) - OFstatic_cast(T3_, x[i - 1]));
107  u[i] = (6.0 * u[i] / (OFstatic_cast(T3_, x[i + 1]) -
108  OFstatic_cast(T3_, x[i - 1])) - sig * u[i - 1]) / p;
109  }
110  if (ypn > 0.99e30) // ignore value for first derivative at point 1
111  qn = un = 0.0;
112  else
113  {
114  qn = 0.5;
115  un = (3.0 / (OFstatic_cast(T3_, x[n - 1]) - OFstatic_cast(T3_, x[n - 2]))) *
116  (ypn - (OFstatic_cast(T3_, y[n - 1]) - OFstatic_cast(T3_, y[n - 2])) /
117  (OFstatic_cast(T3_, x[n - 1]) - OFstatic_cast(T3_, x[n - 2])));
118  }
119  y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.0);
120  for (i = n - 1; i > 0; --i)
121  y2[i - 1] = y2[i - 1] * y2[i] + u[i - 1];
122  delete[] u;
123  return 1;
124  }
125  }
126  return 0;
127  }
128 
129 
145  static int Interpolation(const T1 *xa,
146  const T2 *ya,
147  const T3_ *y2a,
148  const unsigned int na,
149  const T1 *x,
150  T2 *y,
151  const unsigned int n)
152  {
153  if ((xa != NULL) && (ya != NULL) && (y2a != NULL) && (na > 0) && (x != NULL) && (y != NULL) && (n > 0))
154  {
155  register unsigned int k, i;
156  register unsigned int klo = 0;
157  register unsigned int khi = na - 1;
158  T3_ h, b, a;
159  for (i = 0; i < n; ++i)
160  {
161  if ((xa[klo] > x[i]) || (xa[khi] < x[i])) // optimization
162  {
163  klo = 0;
164  khi = na - 1;
165  }
166  while (khi - klo > 1) // search right place in the table, if necessary
167  {
168  k = (khi + klo) >> 1;
169  if (xa[k] > x[i])
170  khi = k;
171  else
172  klo = k;
173  }
174  if (xa[khi] == x[i]) // optimization: use known values
175  y[i] = ya[khi];
176  else
177  {
178  h = OFstatic_cast(T3_, xa[khi]) - OFstatic_cast(T3_, xa[klo]);
179  if (h == 0.0) // bad xa input, values must be distinct !
180  return 0;
181  a = (OFstatic_cast(T3_, xa[khi]) - OFstatic_cast(T3_, x[i])) / h;
182  b = (OFstatic_cast(T3_, x[i]) - OFstatic_cast(T3_, xa[klo])) / h;
183  y[i] = OFstatic_cast(T2, a * OFstatic_cast(T3_, ya[klo]) + b * OFstatic_cast(T3_, ya[khi]) +
184  ((a * a * a - a) * y2a[klo] + (b * b * b - b) * y2a[khi]) * (h * h) / 6.0);
185  }
186  }
187  return 1;
188  }
189  return 0;
190  }
191 };
192 
193 
194 #endif
195 
196 
197 /*
198  *
199  * CVS/RCS Log:
200  * $Log: displint.h,v $
201  * Revision 1.21 2010-10-14 13:16:27 joergr
202  * Updated copyright header. Added reference to COPYRIGHT file.
203  *
204  * Revision 1.20 2010-03-01 09:08:47 uli
205  * Removed some unnecessary include directives in the headers.
206  *
207  * Revision 1.19 2005-12-08 16:48:10 meichel
208  * Changed include path schema for all DCMTK header files
209  *
210  * Revision 1.18 2003/12/23 15:53:22 joergr
211  * Replaced post-increment/decrement operators by pre-increment/decrement
212  * operators where appropriate (e.g. 'i++' by '++i').
213  *
214  * Revision 1.17 2003/12/08 19:20:47 joergr
215  * Adapted type casts to new-style typecast operators defined in ofcast.h.
216  * Removed leading underscore characters from preprocessor symbols (reserved
217  * symbols). Updated copyright header.
218  *
219  * Revision 1.16 2002/07/18 12:30:59 joergr
220  * Removed unused code.
221  *
222  * Revision 1.15 2001/06/01 15:49:51 meichel
223  * Updated copyright header
224  *
225  * Revision 1.14 2000/03/08 16:24:24 meichel
226  * Updated copyright header.
227  *
228  * Revision 1.13 2000/02/02 14:33:54 joergr
229  * Replaced 'delete' statements by 'delete[]' for objects created with 'new[]'.
230  *
231  * Revision 1.12 1999/10/21 08:29:42 joergr
232  * Renamed template type definition from 'T3' to 'T3_' to avoid naming
233  * conflicts.
234  *
235  * Revision 1.11 1999/10/20 18:38:50 joergr
236  * Eliminated default values for template types since this features is not
237  * supported by SunCC 4.x (temporarily introduced '#define' instead).
238  *
239  * Revision 1.10 1999/10/15 09:38:31 joergr
240  * Fixed typos.
241  *
242  * Revision 1.9 1999/10/14 19:05:17 joergr
243  * Fixed typo.
244  *
245  * Revision 1.8 1999/10/01 13:25:35 joergr
246  * Enhanced template class for cubic spline interpolation to support
247  * non-floating point classes/types as y-coordinates.
248  *
249  * Revision 1.7 1999/07/23 14:11:25 joergr
250  * Added preliminary support for 2D bi-cubic spline interpolation (currently
251  * not used).
252  *
253  * Revision 1.6 1999/05/03 11:09:31 joergr
254  * Minor code purifications to keep Sun CC 2.0.1 quiet.
255  *
256  * Revision 1.5 1999/04/29 13:49:08 joergr
257  * Renamed class CubicSpline to DiCubicSpline.
258  *
259  * Revision 1.4 1999/03/24 17:20:26 joergr
260  * Added/Modified comments and formatting.
261  *
262  * Revision 1.3 1999/03/22 08:52:43 joergr
263  * Added/Changed comments.
264  *
265  * Revision 1.2 1999/02/25 16:17:16 joergr
266  * Initialize local variables to avoid compiler warnings (reported by gcc
267  * 2.7.2.1 on Linux).
268  *
269  * Revision 1.1 1999/02/11 16:36:29 joergr
270  * Renamed file to indicate the use of templates.
271  *
272  * Revision 1.2 1999/02/09 14:21:54 meichel
273  * Removed default parameters from template functions, required for Sun CC 4.2
274  *
275  * Revision 1.1 1999/02/04 17:59:23 joergr
276  * Added support for calibration according to Barten transformation (incl.
277  * a DISPLAY file describing the monitor characteristic).
278  *
279  *
280  */


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