ITK  5.2.0
Insight Toolkit
itkMathDetail.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright NumFOCUS
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)) && \
85  (defined(__i386__) || defined(__i386) || 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
121 {
123  {
124  x += static_cast<TInput>(0.5);
125  }
126  else
127  {
128  x -= static_cast<TInput>(0.5);
129  }
130 
131  const auto r = static_cast<TReturn>(x);
132  return (x != static_cast<TInput>(r)) ? r : static_cast<TReturn>(2 * (r / 2));
133 }
134 
135 template <typename TReturn, typename TInput>
136 inline TReturn
138 {
139  x += static_cast<TInput>(0.5);
140  const auto r = static_cast<TReturn>(x);
142  : (x == static_cast<TInput>(r) ? r : r - static_cast<TReturn>(1));
143 }
144 
145 template <typename TReturn, typename TInput>
146 inline TReturn
147 Floor_base(TInput x)
148 {
149  const auto r = static_cast<TReturn>(x);
150 
152  : (x == static_cast<TInput>(r) ? r : r - static_cast<TReturn>(1));
153 }
154 
155 template <typename TReturn, typename TInput>
156 inline TReturn
157 Ceil_base(TInput x)
158 {
159  const auto r = static_cast<TReturn>(x);
160 
161  return (NumericTraits<TInput>::IsNegative(x)) ? r : (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
173 {
174 # if defined(ITK_CHECK_FPU_ROUNDING_MODE)
175  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
176 # endif
177  return _mm_cvtsd_si32(_mm_set_sd(x));
178 }
179 
180 inline int32_t
182 {
183 # if defined(ITK_CHECK_FPU_ROUNDING_MODE)
184  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
185 # endif
186  return _mm_cvtss_si32(_mm_set_ss(x));
187 }
188 
189 #elif GCC_USE_ASM_32IMPL // gcc asm implementation
190 
191 inline int32_t
193 {
194 # if defined(ITK_CHECK_FPU_ROUNDING_MODE)
195  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
196 # endif
197  int32_t r;
198  __asm__ __volatile__("fistpl %0" : "=m"(r) : "t"(x) : "st");
199  return r;
200 }
201 
202 inline int32_t
204 {
205 # if defined(ITK_CHECK_FPU_ROUNDING_MODE)
206  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
207 # endif
208  int32_t r;
209  __asm__ __volatile__("fistpl %0" : "=m"(r) : "t"(x) : "st");
210  return r;
211 }
212 
213 #elif VC_USE_ASM_32IMPL // msvc asm implementation
214 
215 inline int32_t
217 {
218 # if defined(ITK_CHECK_FPU_ROUNDING_MODE)
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
232 {
233 # if defined(ITK_CHECK_FPU_ROUNDING_MODE)
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
249 {
250  return RoundHalfIntegerToEven_base<int32_t, double>(x);
251 }
252 inline int32_t
254 {
255  return RoundHalfIntegerToEven_base<int32_t, float>(x);
256 }
257 
258 #endif
259 
260 #if USE_SSE2_32IMPL || GCC_USE_ASM_32IMPL || VC_USE_ASM_32IMPL
261 
262 inline int32_t
263 RoundHalfIntegerUp_32(double x)
264 {
265  return RoundHalfIntegerToEven_32(2 * x + 0.5) >> 1;
266 }
267 inline int32_t
268 RoundHalfIntegerUp_32(float x)
269 {
270  return RoundHalfIntegerToEven_32(2 * x + 0.5f) >> 1;
271 }
272 
273 inline int32_t
274 Floor_32(double x)
275 {
276  return RoundHalfIntegerToEven_32(2 * x - 0.5) >> 1;
277 }
278 inline int32_t
279 Floor_32(float x)
280 {
281  return RoundHalfIntegerToEven_32(2 * x - 0.5f) >> 1;
282 }
283 
284 inline int32_t
285 Ceil_32(double x)
286 {
287  return -(RoundHalfIntegerToEven_32(-0.5 - 2 * x) >> 1);
288 }
289 inline int32_t
290 Ceil_32(float x)
291 {
292  return -(RoundHalfIntegerToEven_32(-0.5f - 2 * x) >> 1);
293 }
294 
295 #else // Base implementation
296 
297 inline int32_t
299 {
300  return RoundHalfIntegerUp_base<int32_t, double>(x);
301 }
302 inline int32_t
304 {
305  return RoundHalfIntegerUp_base<int32_t, float>(x);
306 }
307 
308 inline int32_t
309 Floor_32(double x)
310 {
311  return Floor_base<int32_t, double>(x);
312 }
313 inline int32_t
314 Floor_32(float x)
315 {
316  return Floor_base<int32_t, float>(x);
317 }
318 
319 inline int32_t
320 Ceil_32(double x)
321 {
322  return Ceil_base<int32_t, double>(x);
323 }
324 inline int32_t
325 Ceil_32(float x)
326 {
327  return Ceil_base<int32_t, float>(x);
328 }
329 
330 #endif // USE_SSE2_32IMPL || GCC_USE_ASM_32IMPL || VC_USE_ASM_32IMPL
331 
333 // 64 bits versions
334 
335 #if USE_SSE2_64IMPL // sse2 implementation
336 
337 inline int64_t
339 {
340 # if defined(ITK_CHECK_FPU_ROUNDING_MODE)
341  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
342 # endif
343  return _mm_cvtsd_si64(_mm_set_sd(x));
344 }
345 
346 inline int64_t
348 {
349 # if defined(ITK_CHECK_FPU_ROUNDING_MODE)
350  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
351 # endif
352  return _mm_cvtss_si64(_mm_set_ss(x));
353 }
354 
355 #elif GCC_USE_ASM_64IMPL // gcc asm implementation
356 
357 inline int64_t
359 {
360 # if defined(ITK_CHECK_FPU_ROUNDING_MODE)
361  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
362 # endif
363  int64_t r;
364  __asm__ __volatile__("fistpll %0" : "=m"(r) : "t"(x) : "st");
365  return r;
366 }
367 
368 inline int64_t
370 {
371 # if defined(ITK_CHECK_FPU_ROUNDING_MODE)
372  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
373 # endif
374  int64_t r;
375  __asm__ __volatile__("fistpll %0" : "=m"(r) : "t"(x) : "st");
376  return r;
377 }
378 
379 #elif VC_USE_ASM_64IMPL // msvc asm implementation
380 
381 inline int64_t
383 {
384 # if defined(ITK_CHECK_FPU_ROUNDING_MODE)
385  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
386 # endif
387  int64_t r;
388  __asm
389  {
390  fld x
391  fistp r
392  }
393  return r;
394 }
395 
396 inline int64_t
398 {
399 # if defined(ITK_CHECK_FPU_ROUNDING_MODE)
400  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
401 # endif
402  int64_t r;
403  __asm
404  {
405  fld x
406  fistp r
407  }
408  return r;
409 }
410 
411 #else // Base implementation
412 
413 inline int64_t
415 {
416  return RoundHalfIntegerToEven_base<int64_t, double>(x);
417 }
418 inline int64_t
420 {
421  return RoundHalfIntegerToEven_base<int64_t, float>(x);
422 }
423 
424 #endif
425 
426 #if USE_SSE2_64IMPL || GCC_USE_ASM_64IMPL || VC_USE_ASM_64IMPL
427 
428 inline int64_t
429 RoundHalfIntegerUp_64(double x)
430 {
431  return RoundHalfIntegerToEven_64(2 * x + 0.5) >> 1;
432 }
433 inline int64_t
434 RoundHalfIntegerUp_64(float x)
435 {
436  return RoundHalfIntegerToEven_64(2 * x + 0.5f) >> 1;
437 }
438 
439 inline int64_t
440 Floor_64(double x)
441 {
442  return RoundHalfIntegerToEven_64(2 * x - 0.5) >> 1;
443 }
444 inline int64_t
445 Floor_64(float x)
446 {
447  return RoundHalfIntegerToEven_64(2 * x - 0.5f) >> 1;
448 }
449 
450 inline int64_t
451 Ceil_64(double x)
452 {
453  return -(RoundHalfIntegerToEven_64(-0.5 - 2 * x) >> 1);
454 }
455 inline int64_t
456 Ceil_64(float x)
457 {
458  return -(RoundHalfIntegerToEven_64(-0.5f - 2 * x) >> 1);
459 }
460 
461 #else // Base implementation
462 
463 inline int64_t
465 {
466  return RoundHalfIntegerUp_base<int64_t, double>(x);
467 }
468 inline int64_t
470 {
471  return RoundHalfIntegerUp_base<int64_t, float>(x);
472 }
473 
474 inline int64_t
475 Floor_64(double x)
476 {
477  return Floor_base<int64_t, double>(x);
478 }
479 inline int64_t
480 Floor_64(float x)
481 {
482  return Floor_base<int64_t, float>(x);
483 }
484 
485 inline int64_t
486 Ceil_64(double x)
487 {
488  return Ceil_base<int64_t, double>(x);
489 }
490 inline int64_t
491 Ceil_64(float x)
492 {
493  return Ceil_base<int64_t, float>(x);
494 }
495 
496 #endif // USE_SSE2_64IMPL || GCC_USE_ASM_64IMPL || VC_USE_ASM_64IMPL
497 
498 template <typename T>
500 
501 template <>
502 struct FloatIEEETraits<float>
503 {
504  using IntType = int32_t;
506 };
507 
508 template <>
509 struct FloatIEEETraits<double>
510 {
511  using IntType = int64_t;
513 };
514 
515 template <typename T>
517 {
518  using FloatType = T;
521 
525 
527  : asFloat(f)
528  {}
530  : asInt(i)
531  {}
532  bool
533  Sign() const
534  {
535  return (asUInt >> (sizeof(asUInt) * 8 - 1)) != 0;
536  }
537  IntType
538  AsULP() const
539  {
540  return this->Sign() ? IntType(~(~UIntType(0) >> 1) - asUInt) : asInt;
541  }
542 };
543 
544 } // end namespace Detail
545 } // end namespace Math
546 
547 // move to itkConceptChecking?
548 namespace Concept
549 {
550 template <typename T>
552 template <>
553 struct FloatOrDouble<float>
554 {};
555 template <>
556 struct FloatOrDouble<double>
557 {};
558 } // end namespace Concept
559 } // end namespace itk
560 
561 #undef USE_SSE2_32IMPL
562 #undef GCC_USE_ASM_32IMPL
563 #undef VC_USE_ASM_32IMPL
564 
565 #undef USE_SSE2_64IMPL
566 #undef GCC_USE_ASM_64IMPL
567 #undef VC_USE_ASM_64IMPL
568 
569 #endif // end of itkMathDetail.h
itk::Math::Detail::FloatIEEETraits< float >::IntType
int32_t IntType
Definition: itkMathDetail.h:504
itk::Math::Detail::RoundHalfIntegerUp_base
TReturn RoundHalfIntegerUp_base(TInput x)
Definition: itkMathDetail.h:137
itk::Math::Detail::FloatIEEE::FloatType
T FloatType
Definition: itkMathDetail.h:518
itk::Math::Detail::Floor_base
TReturn Floor_base(TInput x)
Definition: itkMathDetail.h:147
itk::Math::Detail::FloatIEEE::IntType
typename FloatIEEETraits< T >::IntType IntType
Definition: itkMathDetail.h:519
itk::Math::Detail::RoundHalfIntegerUp_32
int32_t RoundHalfIntegerUp_32(double x)
Definition: itkMathDetail.h:298
itk::Math::Detail::FloatIEEE
Definition: itkMathDetail.h:516
itk::Math::Detail::FloatIEEE::AsULP
IntType AsULP() const
Definition: itkMathDetail.h:538
itk::Math::Detail::FloatIEEE::Sign
bool Sign() const
Definition: itkMathDetail.h:533
itk::int64_t
::int64_t int64_t
Definition: itkIntTypes.h:34
itk::Math::Detail::RoundHalfIntegerUp_64
int64_t RoundHalfIntegerUp_64(double x)
Definition: itkMathDetail.h:464
itk::uint64_t
::uint64_t uint64_t
Definition: itkIntTypes.h:35
itk::Math::Detail::Ceil_64
int64_t Ceil_64(double x)
Definition: itkMathDetail.h:486
itk::Math::Detail::RoundHalfIntegerToEven_base
CLANG_PRAGMA_PUSH CLANG_SUPPRESS_Wfloat_equal TReturn RoundHalfIntegerToEven_base(TInput x)
Definition: itkMathDetail.h:120
itk::Math::Detail::FloatIEEE::asInt
IntType asInt
Definition: itkMathDetail.h:523
itk::Math::Detail::FloatIEEETraits
Definition: itkMathDetail.h:499
itk::Math::Detail::FloatIEEE::asUInt
UIntType asUInt
Definition: itkMathDetail.h:524
itk::Math::Detail::FloatIEEETraits< double >::IntType
int64_t IntType
Definition: itkMathDetail.h:511
itk::Math::Detail::Floor_32
int32_t Floor_32(double x)
Definition: itkMathDetail.h:309
itk::Math::Detail::FloatIEEE::FloatIEEE
FloatIEEE(IntType i)
Definition: itkMathDetail.h:529
itk::Math::Detail::FloatIEEETraits< float >::UIntType
uint32_t UIntType
Definition: itkMathDetail.h:505
itk::Concept::FloatOrDouble
Definition: itkMathDetail.h:551
itkIntTypes.h
itk::NumericTraits
Define additional traits for native types such as int or float.
Definition: itkNumericTraits.h:58
itk
The "itk" namespace contains all Insight Segmentation and Registration Toolkit (ITK) classes....
Definition: itkAnnulusOperator.h:24
itk::Math::Detail::FloatIEEETraits< double >::UIntType
uint64_t UIntType
Definition: itkMathDetail.h:512
itk::Math::Detail::Floor_64
int64_t Floor_64(double x)
Definition: itkMathDetail.h:475
itk::Math::Detail::FloatIEEE::UIntType
typename FloatIEEETraits< T >::UIntType UIntType
Definition: itkMathDetail.h:520
itk::Math::Detail::RoundHalfIntegerToEven_32
CLANG_PRAGMA_POP int32_t RoundHalfIntegerToEven_32(double x)
Definition: itkMathDetail.h:248
itk::Math::Detail::FloatIEEE::FloatIEEE
FloatIEEE(FloatType f)
Definition: itkMathDetail.h:526
itkNumericTraits.h
itk::int32_t
::int32_t int32_t
Definition: itkIntTypes.h:32
itk::uint32_t
::uint32_t uint32_t
Definition: itkIntTypes.h:33
itk::Math::Detail::Ceil_32
int32_t Ceil_32(double x)
Definition: itkMathDetail.h:320
itk::Math::Detail::RoundHalfIntegerToEven_64
int64_t RoundHalfIntegerToEven_64(double x)
Definition: itkMathDetail.h:414
itk::Math::Detail::Ceil_base
TReturn Ceil_base(TInput x)
Definition: itkMathDetail.h:157
itk::Math::Detail::FloatIEEE::asFloat
FloatType asFloat
Definition: itkMathDetail.h:522