ITK  4.13.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 #ifdef ITK_HAVE_FENV_H
35 // The Sun Studio CC compiler seems to have a bug where if cstdio is
36 // included stdio.h must also be included before fenv.h
37 #include <cstdio>
38 #include <fenv.h> // should this be cfenv?
39 #endif /* ITK_HAVE_FENV_H */
40 
41 #if defined( ITK_HAVE_EMMINTRIN_H ) && !defined( ITK_WRAPPING_PARSER )
42 #include <emmintrin.h> // sse 2 intrinsics
43 #endif /* ITK_HAVE_EMMINTRIN_H && ! ITK_WRAPPING_PARSER */
44 
45 // assume no SSE2:
46 #define USE_SSE2_64IMPL 0
47 #define USE_SSE2_32IMPL 0
48 
49 // For apple assume sse2 is on for all intel builds, check for 64 and 32
50 // bit versions
51 #if defined(__APPLE__) && defined( __SSE2__ ) && !defined( ITK_WRAPPING_PARSER )
52 
53 # if defined( __i386__ )
54 # undef USE_SSE2_32IMPL
55 # define USE_SSE2_32IMPL 1
56 # endif
57 
58 # if defined( __x86_64 )
59 // Turn on the 64 bits implementation
60 # undef USE_SSE2_64IMPL
61 # define USE_SSE2_64IMPL 1
62 // Turn on also the 32 bits implementation
63 // since it is available in 64 bits versions.
64 # undef USE_SSE2_32IMPL
65 # define USE_SSE2_32IMPL 1
66 # endif
67 
68 #else
69 
70 // For non-apple (no universal binary possible) just use the
71 // try-compile set ITK_COMPILER_SUPPORTS_SSE2_32 and
72 // ITK_COMPILER_SUPPORTS_SSE2_64 to set values:
73 
74 # if defined(ITK_COMPILER_SUPPORTS_SSE2_32) && !defined( ITK_WRAPPING_PARSER )
75 # undef USE_SSE2_32IMPL
76 # define USE_SSE2_32IMPL 1
77 # endif
78 # if defined(ITK_COMPILER_SUPPORTS_SSE2_64) && !defined( ITK_WRAPPING_PARSER )
79 # undef USE_SSE2_64IMPL
80 # define USE_SSE2_64IMPL 1
81 # endif
82 
83 #endif
84 
85 
86 // Turn on 32-bit and 64-bit asm impl when using GCC on x86 platform with the
87 // following exception:
88 // GCCXML
89 #if defined( __GNUC__ ) && ( !defined( ITK_WRAPPING_PARSER ) ) && ( defined( __i386__ ) || defined( __i386 ) \
90  || defined( __x86_64__ ) || defined( __x86_64 ) )
91 #define GCC_USE_ASM_32IMPL 1
92 #define GCC_USE_ASM_64IMPL 1
93 #else
94 #define GCC_USE_ASM_32IMPL 0
95 #define GCC_USE_ASM_64IMPL 0
96 #endif
97 // Turn on 32-bit and 64-bit asm impl when using msvc on 32 bits windows
98 #if defined( VCL_VC ) && ( !defined( ITK_WRAPPING_PARSER ) ) && !defined( _WIN64 )
99 #define VC_USE_ASM_32IMPL 1
100 #define VC_USE_ASM_64IMPL 1
101 #else
102 #define VC_USE_ASM_32IMPL 0
103 #define VC_USE_ASM_64IMPL 0
104 #endif
105 
106 namespace itk
107 {
108 namespace Math
109 {
110 namespace Detail
111 {
112 // The functions defined in this namespace are not meant to be used directly
113 // and thus do not adhere to the standard backward-compatibility
114 // policy of ITK, as any Detail namespace should be considered private.
115 // Please use the functions from the itk::Math namespace instead
116 
118 // Base versions
119 
120 CLANG_PRAGMA_PUSH
121 CLANG_SUPPRESS_Wfloat_equal
122 
123 template< typename TReturn, typename TInput >
124 inline TReturn RoundHalfIntegerToEven_base(TInput x)
125 {
126  if ( NumericTraits< TInput >::IsNonnegative(x) )
127  {
128  x += static_cast< TInput >( 0.5 );
129  }
130  else
131  {
132  x -= static_cast< TInput >( 0.5 );
133  }
134 
135  const TReturn r = static_cast< TReturn >( x );
136  return ( x != static_cast< TInput >( r ) ) ? r : static_cast< TReturn >( 2 * ( r / 2 ) );
137 }
138 
139 template< typename TReturn, typename TInput >
140 inline TReturn RoundHalfIntegerUp_base(TInput x)
141 {
142  x += static_cast< TInput >( 0.5 );
143  const TReturn r = static_cast< TReturn >( x );
144  return ( NumericTraits< TInput >::IsNonnegative(x) ) ?
145  r :
146  ( x == static_cast< TInput >( r ) ? r : r - static_cast< TReturn >( 1 ) );
147 }
148 
149 template< typename TReturn, typename TInput >
150 inline TReturn Floor_base(TInput x)
151 {
152  const TReturn r = static_cast< TReturn >( x );
153 
154  return ( NumericTraits< TInput >::IsNonnegative(x) ) ?
155  r :
156  ( x == static_cast< TInput >( r ) ? r : r - static_cast< TReturn >( 1 ) );
157 }
158 
159 template< typename TReturn, typename TInput >
160 inline TReturn Ceil_base(TInput x)
161 {
162  const TReturn r = static_cast< TReturn >( x );
163 
164  return ( NumericTraits< TInput >::IsNegative(x) ) ?
165  r :
166  ( x == static_cast< TInput >( r ) ? r : r + static_cast< TReturn >( 1 ) );
167 }
168 
169 CLANG_PRAGMA_POP
170 
172 // 32 bits versions
173 
174 #if USE_SSE2_32IMPL // sse2 implementation
175 
176 inline int32_t RoundHalfIntegerToEven_32(double x)
177 {
178  #if defined( ITK_CHECK_FPU_ROUNDING_MODE ) && defined( HAVE_FENV_H )
179  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
180  #endif
181  return _mm_cvtsd_si32( _mm_set_sd(x) );
182 }
183 
184 inline int32_t RoundHalfIntegerToEven_32(float x)
185 {
186  #if defined( ITK_CHECK_FPU_ROUNDING_MODE ) && defined( HAVE_FENV_H )
187  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
188  #endif
189  return _mm_cvtss_si32( _mm_set_ss(x) );
190 }
191 
192 #elif GCC_USE_ASM_32IMPL // gcc asm implementation
193 
194 inline int32_t RoundHalfIntegerToEven_32(double x)
195 {
196  #if defined( ITK_CHECK_FPU_ROUNDING_MODE ) && defined( HAVE_FENV_H )
197  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
198  #endif
199  int32_t r;
200  __asm__ __volatile__ ( "fistpl %0" : "=m" ( r ) : "t" ( x ) : "st" );
201  return r;
202 }
203 
204 inline int32_t RoundHalfIntegerToEven_32(float x)
205 {
206  #if defined( ITK_CHECK_FPU_ROUNDING_MODE ) && defined( HAVE_FENV_H )
207  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
208  #endif
209  int32_t r;
210  __asm__ __volatile__ ( "fistpl %0" : "=m" ( r ) : "t" ( x ) : "st" );
211  return r;
212 }
213 
214 #elif VC_USE_ASM_32IMPL // msvc asm implementation
215 
216 inline int32_t RoundHalfIntegerToEven_32(double x)
217 {
218  #if defined( ITK_CHECK_FPU_ROUNDING_MODE ) && defined( HAVE_FENV_H )
219  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
220  #endif
221  int32_t r;
222  __asm
223  {
224  fld x
225  fistp r
226  }
227  return r;
228 }
229 
230 inline int32_t RoundHalfIntegerToEven_32(float x)
231 {
232  #if defined( ITK_CHECK_FPU_ROUNDING_MODE ) && defined( HAVE_FENV_H )
233  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
234  #endif
235  int32_t r;
236  __asm
237  {
238  fld x
239  fistp r
240  }
241  return r;
242 }
243 
244 #else // Base implementation
245 
246 inline int32_t RoundHalfIntegerToEven_32(double x) { return RoundHalfIntegerToEven_base< int32_t, double >(x); }
247 inline int32_t RoundHalfIntegerToEven_32(float x) { return RoundHalfIntegerToEven_base< int32_t, float >(x); }
248 
249 #endif
250 
251 #if USE_SSE2_32IMPL || GCC_USE_ASM_32IMPL || VC_USE_ASM_32IMPL
252 
253 inline int32_t RoundHalfIntegerUp_32(double x) { return RoundHalfIntegerToEven_32(2 * x + 0.5) >> 1; }
254 inline int32_t RoundHalfIntegerUp_32(float x) { return RoundHalfIntegerToEven_32(2 * x + 0.5f) >> 1; }
255 
256 inline int32_t Floor_32(double x) { return RoundHalfIntegerToEven_32(2 * x - 0.5) >> 1; }
257 inline int32_t Floor_32(float x) { return RoundHalfIntegerToEven_32(2 * x - 0.5f) >> 1; }
258 
259 inline int32_t Ceil_32(double x) { return -( RoundHalfIntegerToEven_32(-0.5 - 2 * x) >> 1 ); }
260 inline int32_t Ceil_32(float x) { return -( RoundHalfIntegerToEven_32(-0.5f - 2 * x) >> 1 ); }
261 
262 #else // Base implementation
263 
264 inline int32_t RoundHalfIntegerUp_32(double x) { return RoundHalfIntegerUp_base< int32_t, double >(x); }
265 inline int32_t RoundHalfIntegerUp_32(float x) { return RoundHalfIntegerUp_base< int32_t, float >(x); }
266 
267 inline int32_t Floor_32(double x) { return Floor_base< int32_t, double >(x); }
268 inline int32_t Floor_32(float x) { return Floor_base< int32_t, float >(x); }
269 
270 inline int32_t Ceil_32(double x) { return Ceil_base< int32_t, double >(x); }
271 inline int32_t Ceil_32(float x) { return Ceil_base< int32_t, float >(x); }
272 
273 #endif // USE_SSE2_32IMPL || GCC_USE_ASM_32IMPL || VC_USE_ASM_32IMPL
274 
276 // 64 bits versions
277 
278 #if USE_SSE2_64IMPL // sse2 implementation
279 
280 inline int64_t RoundHalfIntegerToEven_64(double x)
281 {
282  #if defined( ITK_CHECK_FPU_ROUNDING_MODE ) && defined( HAVE_FENV_H )
283  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
284  #endif
285  return _mm_cvtsd_si64( _mm_set_sd(x) );
286 }
287 
288 inline int64_t RoundHalfIntegerToEven_64(float x)
289 {
290  #if defined( ITK_CHECK_FPU_ROUNDING_MODE ) && defined( HAVE_FENV_H )
291  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
292  #endif
293  return _mm_cvtss_si64( _mm_set_ss(x) );
294 }
295 
296 #elif GCC_USE_ASM_64IMPL // gcc asm implementation
297 
298 inline int64_t RoundHalfIntegerToEven_64(double x)
299 {
300  #if defined( ITK_CHECK_FPU_ROUNDING_MODE ) && defined( HAVE_FENV_H )
301  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
302  #endif
303  int64_t r;
304  __asm__ __volatile__ ( "fistpll %0" : "=m" ( r ) : "t" ( x ) : "st" );
305  return r;
306 }
307 
308 inline int64_t RoundHalfIntegerToEven_64(float x)
309 {
310  #if defined( ITK_CHECK_FPU_ROUNDING_MODE ) && defined( HAVE_FENV_H )
311  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
312  #endif
313  int64_t r;
314  __asm__ __volatile__ ( "fistpll %0" : "=m" ( r ) : "t" ( x ) : "st" );
315  return r;
316 }
317 
318 #elif VC_USE_ASM_64IMPL // msvc asm implementation
319 
320 inline int64_t RoundHalfIntegerToEven_64(double x)
321 {
322  #if defined( ITK_CHECK_FPU_ROUNDING_MODE ) && defined( HAVE_FENV_H )
323  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
324  #endif
325  int64_t r;
326  __asm
327  {
328  fld x
329  fistp r
330  }
331  return r;
332 }
333 
334 inline int64_t RoundHalfIntegerToEven_64(float x)
335 {
336  #if defined( ITK_CHECK_FPU_ROUNDING_MODE ) && defined( HAVE_FENV_H )
337  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
338  #endif
339  int64_t r;
340  __asm
341  {
342  fld x
343  fistp r
344  }
345  return r;
346 }
347 
348 #else // Base implementation
349 
350 inline int64_t RoundHalfIntegerToEven_64(double x) { return RoundHalfIntegerToEven_base< int64_t, double >(x); }
351 inline int64_t RoundHalfIntegerToEven_64(float x) { return RoundHalfIntegerToEven_base< int64_t, float >(x); }
352 
353 #endif
354 
355 #if USE_SSE2_64IMPL || GCC_USE_ASM_64IMPL || VC_USE_ASM_64IMPL
356 
357 inline int64_t RoundHalfIntegerUp_64(double x) { return RoundHalfIntegerToEven_64(2 * x + 0.5) >> 1; }
358 inline int64_t RoundHalfIntegerUp_64(float x) { return RoundHalfIntegerToEven_64(2 * x + 0.5f) >> 1; }
359 
360 inline int64_t Floor_64(double x) { return RoundHalfIntegerToEven_64(2 * x - 0.5) >> 1; }
361 inline int64_t Floor_64(float x) { return RoundHalfIntegerToEven_64(2 * x - 0.5f) >> 1; }
362 
363 inline int64_t Ceil_64(double x) { return -( RoundHalfIntegerToEven_64(-0.5 - 2 * x) >> 1 ); }
364 inline int64_t Ceil_64(float x) { return -( RoundHalfIntegerToEven_64(-0.5f - 2 * x) >> 1 ); }
365 
366 #else // Base implementation
367 
368 inline int64_t RoundHalfIntegerUp_64(double x) { return RoundHalfIntegerUp_base< int64_t, double >(x); }
369 inline int64_t RoundHalfIntegerUp_64(float x) { return RoundHalfIntegerUp_base< int64_t, float >(x); }
370 
371 inline int64_t Floor_64(double x) { return Floor_base< int64_t, double >(x); }
372 inline int64_t Floor_64(float x) { return Floor_base< int64_t, float >(x); }
373 
374 inline int64_t Ceil_64(double x) { return Ceil_base< int64_t, double >(x); }
375 inline int64_t Ceil_64(float x) { return Ceil_base< int64_t, float >(x); }
376 
377 #endif // USE_SSE2_64IMPL || GCC_USE_ASM_64IMPL || VC_USE_ASM_64IMPL
378 
379 template <typename T>
380 struct FloatIEEETraits;
381 
382 template <>
383 struct FloatIEEETraits<float>
384 {
385  typedef int32_t IntType;
386  typedef uint32_t UIntType;
387 };
388 
389 template <>
390 struct FloatIEEETraits<double>
391 {
392  typedef int64_t IntType;
393  typedef uint64_t UIntType;
394 };
395 
396 template <typename T>
397 union FloatIEEE
398 {
399  typedef T FloatType;
400  typedef typename FloatIEEETraits<T>::IntType IntType;
401  typedef typename FloatIEEETraits<T>::UIntType UIntType;
402 
403  FloatType asFloat;
404  IntType asInt;
405  UIntType asUInt;
406 
407  FloatIEEE(FloatType f): asFloat(f) {}
408  FloatIEEE(IntType i): asInt(i) {}
409  bool Sign() const
410  {
411  return (asUInt >> (sizeof(asUInt)*8-1)) != 0;
412  }
413  IntType AsULP() const
414  {
415  return this->Sign()? IntType(~(~UIntType(0) >> 1) - asUInt) : asInt;
416  }
417 };
418 
419 } // end namespace Detail
420 } // end namespace Math
421 
422 // move to itkConceptChecking?
423 namespace Concept
424 {
425 template< typename T >
426 struct FloatOrDouble;
427 template< >
428 struct FloatOrDouble< float > {};
429 template< >
430 struct FloatOrDouble< double > {};
431 } // end namespace Concept
432 } // end namespace itk
433 
434 #undef USE_SSE2_32IMPL
435 #undef GCC_USE_ASM_32IMPL
436 #undef VC_USE_ASM_32IMPL
437 
438 #undef USE_SSE2_64IMPL
439 #undef GCC_USE_ASM_64IMPL
440 #undef VC_USE_ASM_64IMPL
441 
442 #endif // end of itkMathDetail.h
itk::Math::Detail::RoundHalfIntegerUp_64
int64_t RoundHalfIntegerUp_64(double x)
Definition: itkMathDetail.h:368
itk::Math::Detail::Ceil_32
int32_t Ceil_32(double x)
Definition: itkMathDetail.h:270
itk::Math::Detail::RoundHalfIntegerToEven_base
CLANG_PRAGMA_PUSH CLANG_SUPPRESS_Wfloat_equal TReturn RoundHalfIntegerToEven_base(TInput x)
Definition: itkMathDetail.h:124
itk::Math::Detail::Ceil_base
TReturn Ceil_base(TInput x)
Definition: itkMathDetail.h:160
itk::Math::Detail::RoundHalfIntegerToEven_32
CLANG_PRAGMA_POP int32_t RoundHalfIntegerToEven_32(double x)
Definition: itkMathDetail.h:246
itk::uint32_t
KWIML_INT_uint32_t uint32_t
Definition: itkIntTypes.h:87
itk::int64_t
KWIML_INT_int64_t int64_t
Definition: itkIntTypes.h:88
itk::Math::Detail::Floor_64
int64_t Floor_64(double x)
Definition: itkMathDetail.h:371
itk::Math::Detail::RoundHalfIntegerUp_base
TReturn RoundHalfIntegerUp_base(TInput x)
Definition: itkMathDetail.h:140
itk::Math::Detail::Ceil_64
int64_t Ceil_64(double x)
Definition: itkMathDetail.h:374
itk::Math::Detail::FloatIEEE::IntType
FloatIEEETraits< T >::IntType IntType
Definition: itkMathDetail.h:400
itk::uint64_t
KWIML_INT_uint64_t uint64_t
Definition: itkIntTypes.h:89
itk::Math::Detail::FloatIEEE::UIntType
FloatIEEETraits< T >::UIntType UIntType
Definition: itkMathDetail.h:401
itk::NumericTraits
Define additional traits for native types such as int or float.
Definition: itkNumericTraits.h:78
itk::Math::Detail::RoundHalfIntegerToEven_64
int64_t RoundHalfIntegerToEven_64(double x)
Definition: itkMathDetail.h:350
itk::int32_t
KWIML_INT_int32_t int32_t
Definition: itkIntTypes.h:86
itk::Math::Detail::Floor_base
TReturn Floor_base(TInput x)
Definition: itkMathDetail.h:150
itk::Math::Detail::Floor_32
int32_t Floor_32(double x)
Definition: itkMathDetail.h:267
itk::Math::Detail::RoundHalfIntegerUp_32
int32_t RoundHalfIntegerUp_32(double x)
Definition: itkMathDetail.h:264
Tarballs of the nightly generated Doxygen documentation are available for the html, xml, and tag file.
Generated on Tue Dec 19 2017 04:35:40 for ITK by   doxygen 1.8.5