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 =========================================================================*/
15 /**
16 * @class vtkInterpolatorInternals
17 * @brief internals for vtkImageInterpolator
18 */
19
20 #ifndef vtkImageInterpolatorInternals_h
21 #define vtkImageInterpolatorInternals_h
22
23 #include "vtkMath.h"
24
25 // The interpolator info struct
26 struct vtkInterpolationInfo
27 {
28 const void *Pointer;
29 int Extent[6];
30 vtkIdType Increments[3];
31 int ScalarType;
32 int NumberOfComponents;
33 int BorderMode;
34 int InterpolationMode;
35 void *ExtraInfo;
36 };
37
38 // The interpolation weights struct
39 struct vtkInterpolationWeights : public vtkInterpolationInfo
40 {
41 vtkIdType *Positions[3];
42 void *Weights[3];
43 int WeightExtent[6];
44 int KernelSize[3];
45 int WeightType; // VTK_FLOAT or VTK_DOUBLE
46 void *Workspace;
47 int LastY;
48 int LastZ;
49
50 // partial copy constructor from superclass
vtkInterpolationWeightsvtkInterpolationWeights51 vtkInterpolationWeights(const vtkInterpolationInfo &info) :
52 vtkInterpolationInfo(info), Workspace(nullptr) {}
53 };
54
55 // The internal math functions for the interpolators
56 struct vtkInterpolationMath
57 {
58 // floor with remainder (remainder can be double or float),
59 // includes a small tolerance for values just under an integer
60 template<class F>
61 static int Floor(double x, F &f);
62
63 // round function optimized for various architectures
64 static int Round(double x);
65
66 // border-handling functions for keeping index a with in bounds b, c
67 static int Clamp(int a, int b, int c);
68 static int Wrap(int a, int b, int c);
69 static int Mirror(int a, int b, int c);
70 };
71
72 //--------------------------------------------------------------------------
73 // The 'floor' function is slow, so we want to do an integer
74 // cast but keep the "floor" behavior of always rounding down,
75 // rather than truncating, i.e. we want -0.6 to become -1.
76 // The easiest way to do this is to add a large value in
77 // order to make the value "unsigned", then cast to int, and
78 // then subtract off the large value.
79
80 // On the old i386 architecture even a cast to int is very
81 // expensive because it requires changing the rounding mode
82 // on the FPU. So we use a bit-trick similar to the one
83 // described at http://www.stereopsis.com/FPU.html
84
85 #if defined ia64 || defined __ia64__ || defined _M_IA64
86 #define VTK_INTERPOLATE_64BIT_FLOOR
87 #elif defined __ppc64__ || defined __x86_64__ || defined _M_X64
88 #define VTK_INTERPOLATE_64BIT_FLOOR
89 #elif defined __ppc__ || defined sparc || defined mips
90 #define VTK_INTERPOLATE_32BIT_FLOOR
91 #elif defined i386 || defined _M_IX86
92 #define VTK_INTERPOLATE_I386_FLOOR
93 #endif
94
95 // We add a tolerance of 2^-17 (around 7.6e-6) so that float
96 // values that are just less than the closest integer are
97 // rounded up. This adds robustness against rounding errors.
98
99 #define VTK_INTERPOLATE_FLOOR_TOL 7.62939453125e-06
100
101 template<class F>
Floor(double x,F & f)102 inline int vtkInterpolationMath::Floor(double x, F &f)
103 {
104 #if defined VTK_INTERPOLATE_64BIT_FLOOR
105 x += (103079215104.0 + VTK_INTERPOLATE_FLOOR_TOL);
106 long long i = static_cast<long long>(x);
107 f = static_cast<F>(x - i);
108 return static_cast<int>(i - 103079215104LL);
109 #elif defined VTK_INTERPOLATE_32BIT_FLOOR
110 x += (2147483648.0 + VTK_INTERPOLATE_FLOOR_TOL);
111 unsigned int i = static_cast<unsigned int>(x);
112 f = x - i;
113 return static_cast<int>(i - 2147483648U);
114 #elif defined VTK_INTERPOLATE_I386_FLOOR
115 union { double d; unsigned short s[4]; unsigned int i[2]; } dual;
116 dual.d = x + 103079215104.0; // (2**(52-16))*1.5
117 f = dual.s[0]*0.0000152587890625; // 2**(-16)
118 return static_cast<int>((dual.i[1]<<16)|((dual.i[0])>>16));
119 #else
120 x += VTK_INTERPOLATE_FLOOR_TOL;
121 int i = vtkMath::Floor(x);
122 f = x - i;
123 return i;
124 #endif
125 }
126
127
Round(double x)128 inline int vtkInterpolationMath::Round(double x)
129 {
130 #if defined VTK_INTERPOLATE_64BIT_FLOOR
131 x += (103079215104.5 + VTK_INTERPOLATE_FLOOR_TOL);
132 long long i = static_cast<long long>(x);
133 return static_cast<int>(i - 103079215104LL);
134 #elif defined VTK_INTERPOLATE_32BIT_FLOOR
135 x += (2147483648.5 + VTK_INTERPOLATE_FLOOR_TOL);
136 unsigned int i = static_cast<unsigned int>(x);
137 return static_cast<int>(i - 2147483648U);
138 #elif defined VTK_INTERPOLATE_I386_FLOOR
139 union { double d; unsigned int i[2]; } dual;
140 dual.d = x + 103079215104.5; // (2**(52-16))*1.5
141 return static_cast<int>((dual.i[1]<<16)|((dual.i[0])>>16));
142 #else
143 return vtkMath::Floor(x + (0.5 + VTK_INTERPOLATE_FLOOR_TOL));
144 #endif
145 }
146
147 //----------------------------------------------------------------------------
148 // Perform a clamp to limit an index to [b, c] and subtract b.
149
Clamp(int a,int b,int c)150 inline int vtkInterpolationMath::Clamp(int a, int b, int c)
151 {
152 a = (a <= c ? a : c);
153 a -= b;
154 a = (a >= 0 ? a : 0);
155 return a;
156 }
157
158 //----------------------------------------------------------------------------
159 // Perform a wrap to limit an index to [b, c] and subtract b.
160
Wrap(int a,int b,int c)161 inline int vtkInterpolationMath::Wrap(int a, int b, int c)
162 {
163 int range = c - b + 1;
164 a -= b;
165 a %= range;
166 // required for some % implementations
167 a = (a >= 0 ? a : a + range);
168 return a;
169 }
170
171 //----------------------------------------------------------------------------
172 // Perform a mirror to limit an index to [b, c] and subtract b.
173
Mirror(int a,int b,int c)174 inline int vtkInterpolationMath::Mirror(int a, int b, int c)
175 {
176 #ifndef VTK_IMAGE_BORDER_LEGACY_MIRROR
177 int range = c - b;
178 int ifzero = (range == 0);
179 int range2 = 2*range + ifzero;
180 a -= b;
181 a = (a >= 0 ? a : -a);
182 a %= range2;
183 a = (a <= range ? a : range2 - a);
184 return a;
185 #else
186 int range = c - b + 1;
187 int range2 = 2*range;
188 a -= b;
189 a = (a >= 0 ? a : -a - 1);
190 a %= range2;
191 a = (a < range ? a : range2 - a - 1);
192 return a;
193 #endif
194 }
195
196 #endif
197 // VTK-HeaderTest-Exclude: vtkImageInterpolatorInternals.h
198