ITK  4.9.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 #if defined( ITK_HAVE_EMMINTRIN_H ) && !defined( ITK_WRAPPING_PARSER )
43 #include <emmintrin.h> // sse 2 intrinsics
44 #endif /* ITK_HAVE_EMMINTRIN_H && ! ITK_WRAPPING_PARSER */
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__ ) && !defined( ITK_WRAPPING_PARSER )
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( ITK_WRAPPING_PARSER )
76 # undef USE_SSE2_32IMPL
77 # define USE_SSE2_32IMPL 1
78 # endif
79 # if defined(ITK_COMPILER_SUPPORTS_SSE2_64) && !defined( ITK_WRAPPING_PARSER )
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( ITK_WRAPPING_PARSER ) ) && ( 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( ITK_WRAPPING_PARSER ) ) && !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 CLANG_PRAGMA_PUSH
122 CLANG_SUPPRESS_Wfloat_equal
123 
124 template< typename TReturn, typename TInput >
125 inline TReturn RoundHalfIntegerToEven_base(TInput x)
126 {
128  {
129  x += static_cast< TInput >( 0.5 );
130  }
131  else
132  {
133  x -= static_cast< TInput >( 0.5 );
134  }
135 
136  const TReturn r = static_cast< TReturn >( x );
137  return ( x != static_cast< TInput >( r ) ) ? r : static_cast< TReturn >( 2 * ( r / 2 ) );
138 }
139 
140 template< typename TReturn, typename TInput >
141 inline TReturn RoundHalfIntegerUp_base(TInput x)
142 {
143  x += static_cast< TInput >( 0.5 );
144  const TReturn r = static_cast< TReturn >( x );
146  r :
147  ( x == static_cast< TInput >( r ) ? r : r - static_cast< TReturn >( 1 ) );
148 }
149 
150 template< typename TReturn, typename TInput >
151 inline TReturn Floor_base(TInput x)
152 {
153  const TReturn r = static_cast< TReturn >( x );
154 
156  r :
157  ( x == static_cast< TInput >( r ) ? r : r - static_cast< TReturn >( 1 ) );
158 }
159 
160 template< typename TReturn, typename TInput >
161 inline TReturn Ceil_base(TInput x)
162 {
163  const TReturn r = static_cast< TReturn >( x );
164 
166  r :
167  ( x == static_cast< TInput >( r ) ? r : r + static_cast< TReturn >( 1 ) );
168 }
169 
170 CLANG_PRAGMA_POP
171 
173 // 32 bits versions
174 
175 #if USE_SSE2_32IMPL // sse2 implementation
176 
177 inline int32_t RoundHalfIntegerToEven_32(double x)
178 {
179  #if defined( ITK_CHECK_FPU_ROUNDING_MODE ) && defined( HAVE_FENV_H )
180  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
181  #endif
182  return _mm_cvtsd_si32( _mm_set_sd(x) );
183 }
184 
185 inline int32_t RoundHalfIntegerToEven_32(float x)
186 {
187  #if defined( ITK_CHECK_FPU_ROUNDING_MODE ) && defined( HAVE_FENV_H )
188  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
189  #endif
190  return _mm_cvtss_si32( _mm_set_ss(x) );
191 }
192 
193 #elif GCC_USE_ASM_32IMPL // gcc asm implementation
194 
195 inline int32_t RoundHalfIntegerToEven_32(double x)
196 {
197  #if defined( ITK_CHECK_FPU_ROUNDING_MODE ) && defined( HAVE_FENV_H )
198  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
199  #endif
200  int32_t r;
201  __asm__ __volatile__ ( "fistpl %0" : "=m" ( r ) : "t" ( x ) : "st" );
202  return r;
203 }
204 
205 inline int32_t RoundHalfIntegerToEven_32(float x)
206 {
207  #if defined( ITK_CHECK_FPU_ROUNDING_MODE ) && defined( HAVE_FENV_H )
208  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
209  #endif
210  int32_t r;
211  __asm__ __volatile__ ( "fistpl %0" : "=m" ( r ) : "t" ( x ) : "st" );
212  return r;
213 }
214 
215 #elif VC_USE_ASM_32IMPL // msvc asm implementation
216 
217 inline int32_t RoundHalfIntegerToEven_32(double x)
218 {
219  #if defined( ITK_CHECK_FPU_ROUNDING_MODE ) && defined( HAVE_FENV_H )
220  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
221  #endif
222  int32_t r;
223  __asm
224  {
225  fld x
226  fistp r
227  }
228  return r;
229 }
230 
231 inline int32_t RoundHalfIntegerToEven_32(float x)
232 {
233  #if defined( ITK_CHECK_FPU_ROUNDING_MODE ) && defined( HAVE_FENV_H )
234  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
235  #endif
236  int32_t r;
237  __asm
238  {
239  fld x
240  fistp r
241  }
242  return r;
243 }
244 
245 #else // Base implementation
246 
247 inline int32_t RoundHalfIntegerToEven_32(double x) { return RoundHalfIntegerToEven_base< int32_t, double >(x); }
248 inline int32_t RoundHalfIntegerToEven_32(float x) { return RoundHalfIntegerToEven_base< int32_t, float >(x); }
249 
250 #endif
251 
252 #if USE_SSE2_32IMPL || GCC_USE_ASM_32IMPL || VC_USE_ASM_32IMPL
253 
254 inline int32_t RoundHalfIntegerUp_32(double x) { return RoundHalfIntegerToEven_32(2 * x + 0.5) >> 1; }
255 inline int32_t RoundHalfIntegerUp_32(float x) { return RoundHalfIntegerToEven_32(2 * x + 0.5f) >> 1; }
256 
257 inline int32_t Floor_32(double x) { return RoundHalfIntegerToEven_32(2 * x - 0.5) >> 1; }
258 inline int32_t Floor_32(float x) { return RoundHalfIntegerToEven_32(2 * x - 0.5f) >> 1; }
259 
260 inline int32_t Ceil_32(double x) { return -( RoundHalfIntegerToEven_32(-0.5 - 2 * x) >> 1 ); }
261 inline int32_t Ceil_32(float x) { return -( RoundHalfIntegerToEven_32(-0.5f - 2 * x) >> 1 ); }
262 
263 #else // Base implementation
264 
265 inline int32_t RoundHalfIntegerUp_32(double x) { return RoundHalfIntegerUp_base< int32_t, double >(x); }
266 inline int32_t RoundHalfIntegerUp_32(float x) { return RoundHalfIntegerUp_base< int32_t, float >(x); }
267 
268 inline int32_t Floor_32(double x) { return Floor_base< int32_t, double >(x); }
269 inline int32_t Floor_32(float x) { return Floor_base< int32_t, float >(x); }
270 
271 inline int32_t Ceil_32(double x) { return Ceil_base< int32_t, double >(x); }
272 inline int32_t Ceil_32(float x) { return Ceil_base< int32_t, float >(x); }
273 
274 #endif // USE_SSE2_32IMPL || GCC_USE_ASM_32IMPL || VC_USE_ASM_32IMPL
275 
277 // 64 bits versions
278 
279 #if USE_SSE2_64IMPL // sse2 implementation
280 
281 inline int64_t RoundHalfIntegerToEven_64(double x)
282 {
283  #if defined( ITK_CHECK_FPU_ROUNDING_MODE ) && defined( HAVE_FENV_H )
284  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
285  #endif
286  return _mm_cvtsd_si64( _mm_set_sd(x) );
287 }
288 
289 inline int64_t RoundHalfIntegerToEven_64(float x)
290 {
291  #if defined( ITK_CHECK_FPU_ROUNDING_MODE ) && defined( HAVE_FENV_H )
292  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
293  #endif
294  return _mm_cvtss_si64( _mm_set_ss(x) );
295 }
296 
297 #elif GCC_USE_ASM_64IMPL // gcc asm implementation
298 
299 inline int64_t RoundHalfIntegerToEven_64(double x)
300 {
301  #if defined( ITK_CHECK_FPU_ROUNDING_MODE ) && defined( HAVE_FENV_H )
302  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
303  #endif
304  int64_t r;
305  __asm__ __volatile__ ( "fistpll %0" : "=m" ( r ) : "t" ( x ) : "st" );
306  return r;
307 }
308 
309 inline int64_t RoundHalfIntegerToEven_64(float x)
310 {
311  #if defined( ITK_CHECK_FPU_ROUNDING_MODE ) && defined( HAVE_FENV_H )
312  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
313  #endif
314  int64_t r;
315  __asm__ __volatile__ ( "fistpll %0" : "=m" ( r ) : "t" ( x ) : "st" );
316  return r;
317 }
318 
319 #elif VC_USE_ASM_64IMPL // msvc asm implementation
320 
321 inline int64_t RoundHalfIntegerToEven_64(double x)
322 {
323  #if defined( ITK_CHECK_FPU_ROUNDING_MODE ) && defined( HAVE_FENV_H )
324  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
325  #endif
326  int64_t r;
327  __asm
328  {
329  fld x
330  fistp r
331  }
332  return r;
333 }
334 
335 inline int64_t RoundHalfIntegerToEven_64(float x)
336 {
337  #if defined( ITK_CHECK_FPU_ROUNDING_MODE ) && defined( HAVE_FENV_H )
338  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
339  #endif
340  int64_t r;
341  __asm
342  {
343  fld x
344  fistp r
345  }
346  return r;
347 }
348 
349 #else // Base implementation
350 
351 inline int64_t RoundHalfIntegerToEven_64(double x) { return RoundHalfIntegerToEven_base< int64_t, double >(x); }
352 inline int64_t RoundHalfIntegerToEven_64(float x) { return RoundHalfIntegerToEven_base< int64_t, float >(x); }
353 
354 #endif
355 
356 #if USE_SSE2_64IMPL || GCC_USE_ASM_64IMPL || VC_USE_ASM_64IMPL
357 
358 inline int64_t RoundHalfIntegerUp_64(double x) { return RoundHalfIntegerToEven_64(2 * x + 0.5) >> 1; }
359 inline int64_t RoundHalfIntegerUp_64(float x) { return RoundHalfIntegerToEven_64(2 * x + 0.5f) >> 1; }
360 
361 inline int64_t Floor_64(double x) { return RoundHalfIntegerToEven_64(2 * x - 0.5) >> 1; }
362 inline int64_t Floor_64(float x) { return RoundHalfIntegerToEven_64(2 * x - 0.5f) >> 1; }
363 
364 inline int64_t Ceil_64(double x) { return -( RoundHalfIntegerToEven_64(-0.5 - 2 * x) >> 1 ); }
365 inline int64_t Ceil_64(float x) { return -( RoundHalfIntegerToEven_64(-0.5f - 2 * x) >> 1 ); }
366 
367 #else // Base implementation
368 
369 inline int64_t RoundHalfIntegerUp_64(double x) { return RoundHalfIntegerUp_base< int64_t, double >(x); }
370 inline int64_t RoundHalfIntegerUp_64(float x) { return RoundHalfIntegerUp_base< int64_t, float >(x); }
371 
372 inline int64_t Floor_64(double x) { return Floor_base< int64_t, double >(x); }
373 inline int64_t Floor_64(float x) { return Floor_base< int64_t, float >(x); }
374 
375 inline int64_t Ceil_64(double x) { return Ceil_base< int64_t, double >(x); }
376 inline int64_t Ceil_64(float x) { return Ceil_base< int64_t, float >(x); }
377 
378 #endif // USE_SSE2_64IMPL || GCC_USE_ASM_64IMPL || VC_USE_ASM_64IMPL
379 
380 template <typename T>
382 
383 template <>
384 struct FloatIEEETraits<float>
385 {
386  typedef int32_t IntType;
388 };
389 
390 template <>
391 struct FloatIEEETraits<double>
392 {
393  typedef int64_t IntType;
395 };
396 
397 template <typename T>
399 {
400  typedef T FloatType;
403 
407 
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 >
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
int64_t RoundHalfIntegerUp_64(double x)
int32_t Ceil_32(double x)
CLANG_PRAGMA_PUSH CLANG_SUPPRESS_Wfloat_equal TReturn RoundHalfIntegerToEven_base(TInput x)
TReturn Ceil_base(TInput x)
CLANG_PRAGMA_POP int32_t RoundHalfIntegerToEven_32(double x)
KWIML_INT_uint32_t uint32_t
Definition: itkIntTypes.h:87
KWIML_INT_int64_t int64_t
Definition: itkIntTypes.h:88
int64_t Floor_64(double x)
TReturn RoundHalfIntegerUp_base(TInput x)
int64_t Ceil_64(double x)
FloatIEEETraits< T >::IntType IntType
KWIML_INT_uint64_t uint64_t
Definition: itkIntTypes.h:89
FloatIEEETraits< T >::UIntType UIntType
Define additional traits for native types such as int or float.
int64_t RoundHalfIntegerToEven_64(double x)
KWIML_INT_int32_t int32_t
Definition: itkIntTypes.h:86
TReturn Floor_base(TInput x)
int32_t Floor_32(double x)
int32_t RoundHalfIntegerUp_32(double x)