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=========================================================================*/
19
20#ifndef vtkImageInterpolatorInternals_h
21#define vtkImageInterpolatorInternals_h
22
24#include "vtkEndian.h"
25#include "vtkMath.h"
26
27// The interpolator info struct
42
43// The interpolation weights struct
45{
47 void* Weights[3];
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
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:
124template <class F>
125inline 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
156inline 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
185inline 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
196inline 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
209inline 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 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 vtkDataArray
#define VTK_INTERPOLATE_FLOOR_TOL
int vtkIdType
Definition vtkType.h:332