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];
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:
123template <class F>
124inline 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
155inline 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
184inline 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
195inline 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
208inline 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:159
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