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