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