VTK  9.2.6
vtkImageInterpolatorInternals.h
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: Visualization Toolkit
4  Module: vtkInterpolatorInternals.h
5 
6  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7  All rights reserved.
8  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9 
10  This software is distributed WITHOUT ANY WARRANTY; without even
11  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12  PURPOSE. See the above copyright notice for more information.
13 
14 =========================================================================*/
20 #ifndef vtkImageInterpolatorInternals_h
21 #define vtkImageInterpolatorInternals_h
22 
24 #include "vtkEndian.h"
25 #include "vtkMath.h"
26 
27 // The interpolator info struct
29 {
30  const void* Pointer;
31  int Extent[6];
37  void* ExtraInfo;
38 
41 };
42 
43 // The interpolation weights struct
45 {
47  void* Weights[3];
48  int WeightExtent[6];
49  int KernelSize[3];
50  int WeightType; // VTK_FLOAT or VTK_DOUBLE
51  void* Workspace;
52  int LastY;
53  int LastZ;
54 
55  // partial copy constructor from superclass
57  : vtkInterpolationInfo(info)
58  , Workspace(nullptr)
59  {
60  }
61 };
62 
63 // The internal math functions for the interpolators
65 {
66  // floor with remainder (remainder can be double or float),
67  // includes a small tolerance for values just under an integer
68  template <class F>
69  static int Floor(double x, F& f);
70 
71  // round function optimized for various architectures
72  static int Round(double x);
73 
74  // border-handling functions for keeping index a with in bounds b, c
75  static int Clamp(int a, int b, int c);
76  static int Wrap(int a, int b, int c);
77  static int Mirror(int a, int b, int c);
78 };
79 
80 //--------------------------------------------------------------------------
81 // The 'floor' function is slow, so we want a faster replacement.
82 // The goal is to cast double to integer, but round down instead of
83 // rounding towards zero. In other words, we want -0.1 to become -1.
84 
85 // The easiest way to do this is to add a large value (a bias)
86 // to the input of our 'fast floor' in order to ensure that the
87 // double that we cast to integer is positive. This ensures the
88 // cast will round the value down. After the cast, we can subtract
89 // the bias from the integer result.
90 
91 // We choose a bias of 103079215104 because it has a special property
92 // with respect to ieee-754 double-precision floats. It uses 37 bits
93 // of the 53 significant bits available, leaving 16 bits of precision
94 // after the radix. And the same is true for any number in the range
95 // [-34359738368,34359738367] when added to this bias. This is a
96 // very large range, 16 times the range of a 32-bit int. Essentially,
97 // this bias allows us to use the mantissa of a 'double' as a 52-bit
98 // (36.16) fixed-point value. Hence, we can use our floating-point
99 // hardware for fixed-point math, with float-to-fixed and vice-versa
100 // conversions achieved by simply by adding or subtracting the bias.
101 // See http://www.stereopsis.com/FPU.html for further explanation.
102 
103 // The advantage of fixed (absolute) precision over float (relative)
104 // precision is that when we do math on a coordinate (x,y,z) in the
105 // image, the available precision will be the same regardless of
106 // whether x, y, z are close to (0,0,0) or whether they are far away.
107 // This protects against a common problem in computer graphics where
108 // there is lots of available precision near the origin, but less
109 // precision far from the origin. Instead of relying on relative
110 // precision, we have enforced the use of fixed precision. As a
111 // trade-off, we are limited to the range [-34359738368,34359738367].
112 
113 // The value 2^-17 (around 7.6e-6) is exactly half the value of the
114 // 16th bit past the decimal, so it is a useful tolerance to apply in
115 // our calculations. For our 'fast floor', floating-point values
116 // that are within this tolerance from the closest integer will always
117 // be rounded to the integer, even when the value is less than the
118 // integer. Values further than this tolerance from an integer will
119 // always be rounded down.
120 
121 #define VTK_INTERPOLATE_FLOOR_TOL 7.62939453125e-06
122 
123 // A fast replacement for 'floor' that provides fixed precision:
124 template <class F>
125 inline int vtkInterpolationMath::Floor(double x, F& f)
126 {
127 #if VTK_SIZEOF_VOID_P >= 8
128  // add the bias and then subtract it to achieve the desired result
129  x += (103079215104.0 + VTK_INTERPOLATE_FLOOR_TOL);
130  long long i = static_cast<long long>(x);
131  f = static_cast<F>(x - i);
132  return static_cast<int>(i - 103079215104LL);
133 #elif !defined VTK_WORDS_BIGENDIAN
134  // same as above, but avoid doing any 64-bit integer arithmetic
135  union {
136  double d;
137  unsigned short s[4];
138  unsigned int i[2];
139  } dual;
140  dual.d = x + 103079215104.0; // (2**(52-16))*1.5
141  f = dual.s[0] * 0.0000152587890625; // 2**(-16)
142  return static_cast<int>((dual.i[1] << 16) | ((dual.i[0]) >> 16));
143 #else
144  // and again for big-endian architectures
145  union {
146  double d;
147  unsigned short s[4];
148  unsigned int i[2];
149  } dual;
150  dual.d = x + 103079215104.0; // (2**(52-16))*1.5
151  f = dual.s[3] * 0.0000152587890625; // 2**(-16)
152  return static_cast<int>((dual.i[0] << 16) | ((dual.i[1]) >> 16));
153 #endif
154 }
155 
156 inline int vtkInterpolationMath::Round(double x)
157 {
158 #if VTK_SIZEOF_VOID_P >= 8
159  // add the bias and then subtract it to achieve the desired result
160  x += (103079215104.5 + VTK_INTERPOLATE_FLOOR_TOL);
161  long long i = static_cast<long long>(x);
162  return static_cast<int>(i - 103079215104LL);
163 #elif !defined VTK_WORDS_BIGENDIAN
164  // same as above, but avoid doing any 64-bit integer arithmetic
165  union {
166  double d;
167  unsigned int i[2];
168  } dual;
169  dual.d = x + 103079215104.5; // (2**(52-16))*1.5
170  return static_cast<int>((dual.i[1] << 16) | ((dual.i[0]) >> 16));
171 #else
172  // and again for big-endian architectures
173  union {
174  double d;
175  unsigned int i[2];
176  } dual;
177  dual.d = x + 103079215104.5; // (2**(52-16))*1.5
178  return static_cast<int>((dual.i[0] << 16) | ((dual.i[1]) >> 16));
179 #endif
180 }
181 
182 //----------------------------------------------------------------------------
183 // Perform a clamp to limit an index to [b, c] and subtract b.
184 
185 inline int vtkInterpolationMath::Clamp(int a, int b, int c)
186 {
187  a = (a <= c ? a : c);
188  a -= b;
189  a = (a >= 0 ? a : 0);
190  return a;
191 }
192 
193 //----------------------------------------------------------------------------
194 // Perform a wrap to limit an index to [b, c] and subtract b.
195 
196 inline int vtkInterpolationMath::Wrap(int a, int b, int c)
197 {
198  int range = c - b + 1;
199  a -= b;
200  a %= range;
201  // required for some % implementations
202  a = (a >= 0 ? a : a + range);
203  return a;
204 }
205 
206 //----------------------------------------------------------------------------
207 // Perform a mirror to limit an index to [b, c] and subtract b.
208 
209 inline int vtkInterpolationMath::Mirror(int a, int b, int c)
210 {
211 #ifndef VTK_IMAGE_BORDER_LEGACY_MIRROR
212  int range = c - b;
213  int ifzero = (range == 0);
214  int range2 = 2 * range + ifzero;
215  a -= b;
216  a = (a >= 0 ? a : -a);
217  a %= range2;
218  a = (a <= range ? a : range2 - a);
219  return a;
220 #else
221  int range = c - b + 1;
222  int range2 = 2 * range;
223  a -= b;
224  a = (a >= 0 ? a : -a - 1);
225  a %= range2;
226  a = (a < range ? a : range2 - a - 1);
227  return a;
228 #endif
229 }
230 
231 #endif
232 // VTK-HeaderTest-Exclude: vtkImageInterpolatorInternals.h
static int Mirror(int a, int b, int c)
static int Clamp(int a, int b, int c)
static int Wrap(int a, int b, int c)
int vtkIdType
Definition: vtkType.h:332
abstract superclass for arrays of numeric data
Definition: vtkDataArray.h:55
static int Floor(double x, F &f)
#define VTK_INTERPOLATE_FLOOR_TOL
vtkInterpolationWeights(const vtkInterpolationInfo &info)