ITK  5.0.0
Insight Segmentation and Registration Toolkit
itkMathDetail.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright Insight Software Consortium
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
18 /*=========================================================================
19  *
20  * Portions of this file are subject to the VTK Toolkit Version 3 copyright.
21  *
22  * Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
23  *
24  * For complete copyright, license and disclaimer of warranty information
25  * please refer to the NOTICE file at the top of the ITK source tree.
26  *
27  *=========================================================================*/
28 #ifndef itkMathDetail_h
29 #define itkMathDetail_h
30 
31 #include "itkIntTypes.h"
32 #include "itkNumericTraits.h"
33 
34 #include <cfenv>
35 
36 #if defined( ITK_HAVE_EMMINTRIN_H ) && !defined( ITK_WRAPPING_PARSER )
37 #include <emmintrin.h> // sse 2 intrinsics
38 #endif /* ITK_HAVE_EMMINTRIN_H && ! ITK_WRAPPING_PARSER */
39 
40 // assume no SSE2:
41 #define USE_SSE2_64IMPL 0
42 #define USE_SSE2_32IMPL 0
43 
44 // For apple assume sse2 is on for all intel builds, check for 64 and 32
45 // bit versions
46 #if defined(__APPLE__) && defined( __SSE2__ ) && !defined( ITK_WRAPPING_PARSER )
47 
48 # if defined( __i386__ )
49 # undef USE_SSE2_32IMPL
50 # define USE_SSE2_32IMPL 1
51 # endif
52 
53 # if defined( __x86_64 )
54 // Turn on the 64 bits implementation
55 # undef USE_SSE2_64IMPL
56 # define USE_SSE2_64IMPL 1
57 // Turn on also the 32 bits implementation
58 // since it is available in 64 bits versions.
59 # undef USE_SSE2_32IMPL
60 # define USE_SSE2_32IMPL 1
61 # endif
62 
63 #else
64 
65 // For non-apple (no universal binary possible) just use the
66 // try-compile set ITK_COMPILER_SUPPORTS_SSE2_32 and
67 // ITK_COMPILER_SUPPORTS_SSE2_64 to set values:
68 
69 # if defined(ITK_COMPILER_SUPPORTS_SSE2_32) && !defined( ITK_WRAPPING_PARSER )
70 # undef USE_SSE2_32IMPL
71 # define USE_SSE2_32IMPL 1
72 # endif
73 # if defined(ITK_COMPILER_SUPPORTS_SSE2_64) && !defined( ITK_WRAPPING_PARSER )
74 # undef USE_SSE2_64IMPL
75 # define USE_SSE2_64IMPL 1
76 # endif
77 
78 #endif
79 
80 
81 // Turn on 32-bit and 64-bit asm impl when using GCC on x86 platform with the
82 // following exception:
83 // GCCXML
84 #if defined( __GNUC__ ) && ( !defined( ITK_WRAPPING_PARSER ) ) && ( defined( __i386__ ) || defined( __i386 ) \
85  || defined( __x86_64__ ) || defined( __x86_64 ) )
86 #define GCC_USE_ASM_32IMPL 1
87 #define GCC_USE_ASM_64IMPL 1
88 #else
89 #define GCC_USE_ASM_32IMPL 0
90 #define GCC_USE_ASM_64IMPL 0
91 #endif
92 // Turn on 32-bit and 64-bit asm impl when using msvc on 32 bits windows
93 #if defined( VCL_VC ) && ( !defined( ITK_WRAPPING_PARSER ) ) && !defined( _WIN64 )
94 #define VC_USE_ASM_32IMPL 1
95 #define VC_USE_ASM_64IMPL 1
96 #else
97 #define VC_USE_ASM_32IMPL 0
98 #define VC_USE_ASM_64IMPL 0
99 #endif
100 
101 namespace itk
102 {
103 namespace Math
104 {
105 namespace Detail
106 {
107 // The functions defined in this namespace are not meant to be used directly
108 // and thus do not adhere to the standard backward-compatibility
109 // policy of ITK, as any Detail namespace should be considered private.
110 // Please use the functions from the itk::Math namespace instead
111 
113 // Base versions
114 
115 CLANG_PRAGMA_PUSH
116 CLANG_SUPPRESS_Wfloat_equal
117 
118 template< typename TReturn, typename TInput >
119 inline TReturn RoundHalfIntegerToEven_base(TInput x)
120 {
122  {
123  x += static_cast< TInput >( 0.5 );
124  }
125  else
126  {
127  x -= static_cast< TInput >( 0.5 );
128  }
129 
130  const auto r = static_cast< TReturn >( x );
131  return ( x != static_cast< TInput >( r ) ) ? r : static_cast< TReturn >( 2 * ( r / 2 ) );
132 }
133 
134 template< typename TReturn, typename TInput >
135 inline TReturn RoundHalfIntegerUp_base(TInput x)
136 {
137  x += static_cast< TInput >( 0.5 );
138  const auto r = static_cast< TReturn >( x );
140  r :
141  ( x == static_cast< TInput >( r ) ? r : r - static_cast< TReturn >( 1 ) );
142 }
143 
144 template< typename TReturn, typename TInput >
145 inline TReturn Floor_base(TInput x)
146 {
147  const auto r = static_cast< TReturn >( x );
148 
150  r :
151  ( x == static_cast< TInput >( r ) ? r : r - static_cast< TReturn >( 1 ) );
152 }
153 
154 template< typename TReturn, typename TInput >
155 inline TReturn Ceil_base(TInput x)
156 {
157  const auto r = static_cast< TReturn >( x );
158 
160  r :
161  ( x == static_cast< TInput >( r ) ? r : r + static_cast< TReturn >( 1 ) );
162 }
163 
164 CLANG_PRAGMA_POP
165 
167 // 32 bits versions
168 
169 #if USE_SSE2_32IMPL // sse2 implementation
170 
171 inline int32_t RoundHalfIntegerToEven_32(double x)
172 {
173  #if defined( ITK_CHECK_FPU_ROUNDING_MODE )
174  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
175  #endif
176  return _mm_cvtsd_si32( _mm_set_sd(x) );
177 }
178 
179 inline int32_t RoundHalfIntegerToEven_32(float x)
180 {
181  #if defined( ITK_CHECK_FPU_ROUNDING_MODE )
182  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
183  #endif
184  return _mm_cvtss_si32( _mm_set_ss(x) );
185 }
186 
187 #elif GCC_USE_ASM_32IMPL // gcc asm implementation
188 
189 inline int32_t RoundHalfIntegerToEven_32(double x)
190 {
191  #if defined( ITK_CHECK_FPU_ROUNDING_MODE )
192  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
193  #endif
194  int32_t r;
195  __asm__ __volatile__ ( "fistpl %0" : "=m" ( r ) : "t" ( x ) : "st" );
196  return r;
197 }
198 
199 inline int32_t RoundHalfIntegerToEven_32(float x)
200 {
201  #if defined( ITK_CHECK_FPU_ROUNDING_MODE )
202  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
203  #endif
204  int32_t r;
205  __asm__ __volatile__ ( "fistpl %0" : "=m" ( r ) : "t" ( x ) : "st" );
206  return r;
207 }
208 
209 #elif VC_USE_ASM_32IMPL // msvc asm implementation
210 
211 inline int32_t RoundHalfIntegerToEven_32(double x)
212 {
213  #if defined( ITK_CHECK_FPU_ROUNDING_MODE )
214  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
215  #endif
216  int32_t r;
217  __asm
218  {
219  fld x
220  fistp r
221  }
222  return r;
223 }
224 
225 inline int32_t RoundHalfIntegerToEven_32(float x)
226 {
227  #if defined( ITK_CHECK_FPU_ROUNDING_MODE )
228  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
229  #endif
230  int32_t r;
231  __asm
232  {
233  fld x
234  fistp r
235  }
236  return r;
237 }
238 
239 #else // Base implementation
240 
241 inline int32_t RoundHalfIntegerToEven_32(double x) { return RoundHalfIntegerToEven_base< int32_t, double >(x); }
242 inline int32_t RoundHalfIntegerToEven_32(float x) { return RoundHalfIntegerToEven_base< int32_t, float >(x); }
243 
244 #endif
245 
246 #if USE_SSE2_32IMPL || GCC_USE_ASM_32IMPL || VC_USE_ASM_32IMPL
247 
248 inline int32_t RoundHalfIntegerUp_32(double x) { return RoundHalfIntegerToEven_32(2 * x + 0.5) >> 1; }
249 inline int32_t RoundHalfIntegerUp_32(float x) { return RoundHalfIntegerToEven_32(2 * x + 0.5f) >> 1; }
250 
251 inline int32_t Floor_32(double x) { return RoundHalfIntegerToEven_32(2 * x - 0.5) >> 1; }
252 inline int32_t Floor_32(float x) { return RoundHalfIntegerToEven_32(2 * x - 0.5f) >> 1; }
253 
254 inline int32_t Ceil_32(double x) { return -( RoundHalfIntegerToEven_32(-0.5 - 2 * x) >> 1 ); }
255 inline int32_t Ceil_32(float x) { return -( RoundHalfIntegerToEven_32(-0.5f - 2 * x) >> 1 ); }
256 
257 #else // Base implementation
258 
259 inline int32_t RoundHalfIntegerUp_32(double x) { return RoundHalfIntegerUp_base< int32_t, double >(x); }
260 inline int32_t RoundHalfIntegerUp_32(float x) { return RoundHalfIntegerUp_base< int32_t, float >(x); }
261 
262 inline int32_t Floor_32(double x) { return Floor_base< int32_t, double >(x); }
263 inline int32_t Floor_32(float x) { return Floor_base< int32_t, float >(x); }
264 
265 inline int32_t Ceil_32(double x) { return Ceil_base< int32_t, double >(x); }
266 inline int32_t Ceil_32(float x) { return Ceil_base< int32_t, float >(x); }
267 
268 #endif // USE_SSE2_32IMPL || GCC_USE_ASM_32IMPL || VC_USE_ASM_32IMPL
269 
271 // 64 bits versions
272 
273 #if USE_SSE2_64IMPL // sse2 implementation
274 
275 inline int64_t RoundHalfIntegerToEven_64(double x)
276 {
277  #if defined( ITK_CHECK_FPU_ROUNDING_MODE )
278  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
279  #endif
280  return _mm_cvtsd_si64( _mm_set_sd(x) );
281 }
282 
283 inline int64_t RoundHalfIntegerToEven_64(float x)
284 {
285  #if defined( ITK_CHECK_FPU_ROUNDING_MODE )
286  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
287  #endif
288  return _mm_cvtss_si64( _mm_set_ss(x) );
289 }
290 
291 #elif GCC_USE_ASM_64IMPL // gcc asm implementation
292 
293 inline int64_t RoundHalfIntegerToEven_64(double x)
294 {
295  #if defined( ITK_CHECK_FPU_ROUNDING_MODE )
296  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
297  #endif
298  int64_t r;
299  __asm__ __volatile__ ( "fistpll %0" : "=m" ( r ) : "t" ( x ) : "st" );
300  return r;
301 }
302 
303 inline int64_t RoundHalfIntegerToEven_64(float x)
304 {
305  #if defined( ITK_CHECK_FPU_ROUNDING_MODE )
306  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
307  #endif
308  int64_t r;
309  __asm__ __volatile__ ( "fistpll %0" : "=m" ( r ) : "t" ( x ) : "st" );
310  return r;
311 }
312 
313 #elif VC_USE_ASM_64IMPL // msvc asm implementation
314 
315 inline int64_t RoundHalfIntegerToEven_64(double x)
316 {
317  #if defined( ITK_CHECK_FPU_ROUNDING_MODE )
318  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
319  #endif
320  int64_t r;
321  __asm
322  {
323  fld x
324  fistp r
325  }
326  return r;
327 }
328 
329 inline int64_t RoundHalfIntegerToEven_64(float x)
330 {
331  #if defined( ITK_CHECK_FPU_ROUNDING_MODE )
332  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
333  #endif
334  int64_t r;
335  __asm
336  {
337  fld x
338  fistp r
339  }
340  return r;
341 }
342 
343 #else // Base implementation
344 
345 inline int64_t RoundHalfIntegerToEven_64(double x) { return RoundHalfIntegerToEven_base< int64_t, double >(x); }
346 inline int64_t RoundHalfIntegerToEven_64(float x) { return RoundHalfIntegerToEven_base< int64_t, float >(x); }
347 
348 #endif
349 
350 #if USE_SSE2_64IMPL || GCC_USE_ASM_64IMPL || VC_USE_ASM_64IMPL
351 
352 inline int64_t RoundHalfIntegerUp_64(double x) { return RoundHalfIntegerToEven_64(2 * x + 0.5) >> 1; }
353 inline int64_t RoundHalfIntegerUp_64(float x) { return RoundHalfIntegerToEven_64(2 * x + 0.5f) >> 1; }
354 
355 inline int64_t Floor_64(double x) { return RoundHalfIntegerToEven_64(2 * x - 0.5) >> 1; }
356 inline int64_t Floor_64(float x) { return RoundHalfIntegerToEven_64(2 * x - 0.5f) >> 1; }
357 
358 inline int64_t Ceil_64(double x) { return -( RoundHalfIntegerToEven_64(-0.5 - 2 * x) >> 1 ); }
359 inline int64_t Ceil_64(float x) { return -( RoundHalfIntegerToEven_64(-0.5f - 2 * x) >> 1 ); }
360 
361 #else // Base implementation
362 
363 inline int64_t RoundHalfIntegerUp_64(double x) { return RoundHalfIntegerUp_base< int64_t, double >(x); }
364 inline int64_t RoundHalfIntegerUp_64(float x) { return RoundHalfIntegerUp_base< int64_t, float >(x); }
365 
366 inline int64_t Floor_64(double x) { return Floor_base< int64_t, double >(x); }
367 inline int64_t Floor_64(float x) { return Floor_base< int64_t, float >(x); }
368 
369 inline int64_t Ceil_64(double x) { return Ceil_base< int64_t, double >(x); }
370 inline int64_t Ceil_64(float x) { return Ceil_base< int64_t, float >(x); }
371 
372 #endif // USE_SSE2_64IMPL || GCC_USE_ASM_64IMPL || VC_USE_ASM_64IMPL
373 
374 template <typename T>
376 
377 template <>
378 struct FloatIEEETraits<float>
379 {
380  using IntType = int32_t;
382 };
383 
384 template <>
385 struct FloatIEEETraits<double>
386 {
387  using IntType = int64_t;
389 };
390 
391 template <typename T>
393 {
394  using FloatType = T;
397 
401 
404  bool Sign() const
405  {
406  return (asUInt >> (sizeof(asUInt)*8-1)) != 0;
407  }
408  IntType AsULP() const
409  {
410  return this->Sign()? IntType(~(~UIntType(0) >> 1) - asUInt) : asInt;
411  }
412 };
413 
414 } // end namespace Detail
415 } // end namespace Math
416 
417 // move to itkConceptChecking?
418 namespace Concept
419 {
420 template< typename T >
422 template< >
423 struct FloatOrDouble< float > {};
424 template< >
425 struct FloatOrDouble< double > {};
426 } // end namespace Concept
427 } // end namespace itk
428 
429 #undef USE_SSE2_32IMPL
430 #undef GCC_USE_ASM_32IMPL
431 #undef VC_USE_ASM_32IMPL
432 
433 #undef USE_SSE2_64IMPL
434 #undef GCC_USE_ASM_64IMPL
435 #undef VC_USE_ASM_64IMPL
436 
437 #endif // end of itkMathDetail.h
int64_t RoundHalfIntegerUp_64(double x)
int32_t Ceil_32(double x)
CLANG_PRAGMA_PUSH CLANG_SUPPRESS_Wfloat_equal TReturn RoundHalfIntegerToEven_base(TInput x)
Define numeric traits for std::vector.
typename FloatIEEETraits< T >::IntType IntType
TReturn Ceil_base(TInput x)
CLANG_PRAGMA_POP int32_t RoundHalfIntegerToEven_32(double x)
int64_t Floor_64(double x)
TReturn RoundHalfIntegerUp_base(TInput x)
int64_t Ceil_64(double x)
typename FloatIEEETraits< T >::UIntType UIntType
::uint32_t uint32_t
Definition: itkIntTypes.h:33
::int64_t int64_t
Definition: itkIntTypes.h:34
int64_t RoundHalfIntegerToEven_64(double x)
::int32_t int32_t
Definition: itkIntTypes.h:32
TReturn Floor_base(TInput x)
::uint64_t uint64_t
Definition: itkIntTypes.h:35
int32_t Floor_32(double x)
int32_t RoundHalfIntegerUp_32(double x)