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