ITK  5.3.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_COMPILER_SUPPORTS_SSE2_32) || defined(ITK_COMPILER_SUPPORTS_SSE2_64)) && !defined(ITK_WRAPPING_PARSER)
37 # include <emmintrin.h> // SSE2 intrinsics
38 #endif
39 
40 // Initially assume no SSE2:
41 #define USE_SSE2_64IMPL 0
42 #define USE_SSE2_32IMPL 0
43 
44 // Turn on 32-bit and/or 64-bit SSE2 impl when using any compiler on x86 platform.
45 #if defined(ITK_COMPILER_SUPPORTS_SSE2_32) && !defined(ITK_WRAPPING_PARSER)
46 # undef USE_SSE2_32IMPL
47 # define USE_SSE2_32IMPL 1
48 #endif
49 #if defined(ITK_COMPILER_SUPPORTS_SSE2_64) && !defined(ITK_WRAPPING_PARSER)
50 # undef USE_SSE2_64IMPL
51 # define USE_SSE2_64IMPL 1
52 #endif
53 
54 // Turn on 32-bit and 64-bit asm impl when using GCC/clang on x86 platform with the
55 // following exception:
56 // GCCXML
57 #if defined(__GNUC__) && !defined(ITK_WRAPPING_PARSER) && \
58  (defined(__i386__) || defined(__i386) || defined(__x86_64__) || defined(__x86_64))
59 # define GCC_USE_ASM_32IMPL 1
60 # define GCC_USE_ASM_64IMPL 1
61 #else
62 # define GCC_USE_ASM_32IMPL 0
63 # define GCC_USE_ASM_64IMPL 0
64 #endif
65 
66 // Turn on 32-bit and 64-bit asm impl when using MSVC on 32-bit x86 Windows
67 #if defined(_MSC_VER) && !defined(ITK_WRAPPING_PARSER) && !defined(_WIN64)
68 # define VC_USE_ASM_32IMPL 1
69 # define VC_USE_ASM_64IMPL 1
70 #else
71 # define VC_USE_ASM_32IMPL 0
72 # define VC_USE_ASM_64IMPL 0
73 #endif
74 
75 namespace itk
76 {
77 namespace Math
78 {
79 namespace Detail
80 {
81 // The functions defined in this namespace are not meant to be used directly
82 // and thus do not adhere to the standard backward-compatibility
83 // policy of ITK, as any Detail namespace should be considered private.
84 // Please use the functions from the itk::Math namespace instead
85 
87 // Base versions
88 
89 CLANG_PRAGMA_PUSH
90 CLANG_SUPPRESS_Wfloat_equal
91 
92  template <typename TReturn, typename TInput>
93  inline TReturn
95 {
97  {
98  x += static_cast<TInput>(0.5);
99  }
100  else
101  {
102  x -= static_cast<TInput>(0.5);
103  }
104 
105  const auto r = static_cast<TReturn>(x);
106  return (x != static_cast<TInput>(r)) ? r : static_cast<TReturn>(2 * (r / 2));
107 }
108 
109 template <typename TReturn, typename TInput>
110 inline TReturn
112 {
113  x += static_cast<TInput>(0.5);
114  const auto r = static_cast<TReturn>(x);
116  : (x == static_cast<TInput>(r) ? r : r - static_cast<TReturn>(1));
117 }
118 
119 template <typename TReturn, typename TInput>
120 inline TReturn
121 Floor_base(TInput x)
122 {
123  const auto r = static_cast<TReturn>(x);
124 
126  : (x == static_cast<TInput>(r) ? r : r - static_cast<TReturn>(1));
127 }
128 
129 template <typename TReturn, typename TInput>
130 inline TReturn
131 Ceil_base(TInput x)
132 {
133  const auto r = static_cast<TReturn>(x);
134 
135  return (NumericTraits<TInput>::IsNegative(x)) ? r : (x == static_cast<TInput>(r) ? r : r + static_cast<TReturn>(1));
136 }
137 
138 CLANG_PRAGMA_POP
139 
141 // 32 bits versions
142 
143 #if USE_SSE2_32IMPL // SSE2 implementation
144 
145 inline int32_t
147 {
148 # if defined(ITK_CHECK_FPU_ROUNDING_MODE)
149  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
150 # endif
151  return _mm_cvtsd_si32(_mm_set_sd(x));
152 }
153 
154 inline int32_t
156 {
157 # if defined(ITK_CHECK_FPU_ROUNDING_MODE)
158  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
159 # endif
160  return _mm_cvtss_si32(_mm_set_ss(x));
161 }
162 
163 #elif GCC_USE_ASM_32IMPL // GCC/clang x86 asm implementation
164 
165 inline int32_t
167 {
168 # if defined(ITK_CHECK_FPU_ROUNDING_MODE)
169  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
170 # endif
171  int32_t r;
172  __asm__ __volatile__("fistpl %0" : "=m"(r) : "t"(x) : "st");
173  return r;
174 }
175 
176 inline int32_t
178 {
179 # if defined(ITK_CHECK_FPU_ROUNDING_MODE)
180  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
181 # endif
182  int32_t r;
183  __asm__ __volatile__("fistpl %0" : "=m"(r) : "t"(x) : "st");
184  return r;
185 }
186 
187 #elif VC_USE_ASM_32IMPL // msvc asm implementation
188 
189 inline int32_t
191 {
192 # if defined(ITK_CHECK_FPU_ROUNDING_MODE)
193  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
194 # endif
195  int32_t r;
196  __asm
197  {
198  fld x
199  fistp r
200  }
201  return r;
202 }
203 
204 inline int32_t
206 {
207 # if defined(ITK_CHECK_FPU_ROUNDING_MODE)
208  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
209 # endif
210  int32_t r;
211  __asm
212  {
213  fld x
214  fistp r
215  }
216  return r;
217 }
218 
219 #else // Base implementation
220 
221 inline int32_t
223 {
224  return RoundHalfIntegerToEven_base<int32_t, double>(x);
225 }
226 inline int32_t
228 {
229  return RoundHalfIntegerToEven_base<int32_t, float>(x);
230 }
231 
232 #endif
233 
234 #if USE_SSE2_32IMPL || GCC_USE_ASM_32IMPL || VC_USE_ASM_32IMPL
235 
236 inline int32_t
237 RoundHalfIntegerUp_32(double x)
238 {
239  return RoundHalfIntegerToEven_32(2 * x + 0.5) >> 1;
240 }
241 inline int32_t
242 RoundHalfIntegerUp_32(float x)
243 {
244  return RoundHalfIntegerToEven_32(2 * x + 0.5f) >> 1;
245 }
246 
247 inline int32_t
248 Floor_32(double x)
249 {
250  return RoundHalfIntegerToEven_32(2 * x - 0.5) >> 1;
251 }
252 inline int32_t
253 Floor_32(float x)
254 {
255  return RoundHalfIntegerToEven_32(2 * x - 0.5f) >> 1;
256 }
257 
258 inline int32_t
259 Ceil_32(double x)
260 {
261  return -(RoundHalfIntegerToEven_32(-0.5 - 2 * x) >> 1);
262 }
263 inline int32_t
264 Ceil_32(float x)
265 {
266  return -(RoundHalfIntegerToEven_32(-0.5f - 2 * x) >> 1);
267 }
268 
269 #else // Base implementation
270 
271 inline int32_t
273 {
274  return RoundHalfIntegerUp_base<int32_t, double>(x);
275 }
276 inline int32_t
278 {
279  return RoundHalfIntegerUp_base<int32_t, float>(x);
280 }
281 
282 inline int32_t
283 Floor_32(double x)
284 {
285  return Floor_base<int32_t, double>(x);
286 }
287 inline int32_t
288 Floor_32(float x)
289 {
290  return Floor_base<int32_t, float>(x);
291 }
292 
293 inline int32_t
294 Ceil_32(double x)
295 {
296  return Ceil_base<int32_t, double>(x);
297 }
298 inline int32_t
299 Ceil_32(float x)
300 {
301  return Ceil_base<int32_t, float>(x);
302 }
303 
304 #endif // USE_SSE2_32IMPL || GCC_USE_ASM_32IMPL || VC_USE_ASM_32IMPL
305 
307 // 64 bits versions
308 
309 #if USE_SSE2_64IMPL // SSE2 implementation
310 
311 inline int64_t
313 {
314 # if defined(ITK_CHECK_FPU_ROUNDING_MODE)
315  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
316 # endif
317  return _mm_cvtsd_si64(_mm_set_sd(x));
318 }
319 
320 inline int64_t
322 {
323 # if defined(ITK_CHECK_FPU_ROUNDING_MODE)
324  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
325 # endif
326  return _mm_cvtss_si64(_mm_set_ss(x));
327 }
328 
329 #elif GCC_USE_ASM_64IMPL // GCC/clang x86 asm implementation
330 
331 inline int64_t
333 {
334 # if defined(ITK_CHECK_FPU_ROUNDING_MODE)
335  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
336 # endif
337  int64_t r;
338  __asm__ __volatile__("fistpll %0" : "=m"(r) : "t"(x) : "st");
339  return r;
340 }
341 
342 inline int64_t
344 {
345 # if defined(ITK_CHECK_FPU_ROUNDING_MODE)
346  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
347 # endif
348  int64_t r;
349  __asm__ __volatile__("fistpll %0" : "=m"(r) : "t"(x) : "st");
350  return r;
351 }
352 
353 #elif VC_USE_ASM_64IMPL // msvc asm implementation
354 
355 inline int64_t
357 {
358 # if defined(ITK_CHECK_FPU_ROUNDING_MODE)
359  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
360 # endif
361  int64_t r;
362  __asm
363  {
364  fld x
365  fistp r
366  }
367  return r;
368 }
369 
370 inline int64_t
372 {
373 # if defined(ITK_CHECK_FPU_ROUNDING_MODE)
374  itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
375 # endif
376  int64_t r;
377  __asm
378  {
379  fld x
380  fistp r
381  }
382  return r;
383 }
384 
385 #else // Base implementation
386 
387 inline int64_t
389 {
390  return RoundHalfIntegerToEven_base<int64_t, double>(x);
391 }
392 inline int64_t
394 {
395  return RoundHalfIntegerToEven_base<int64_t, float>(x);
396 }
397 
398 #endif
399 
400 #if USE_SSE2_64IMPL || GCC_USE_ASM_64IMPL || VC_USE_ASM_64IMPL
401 
402 inline int64_t
403 RoundHalfIntegerUp_64(double x)
404 {
405  return RoundHalfIntegerToEven_64(2 * x + 0.5) >> 1;
406 }
407 inline int64_t
408 RoundHalfIntegerUp_64(float x)
409 {
410  return RoundHalfIntegerToEven_64(2 * x + 0.5f) >> 1;
411 }
412 
413 inline int64_t
414 Floor_64(double x)
415 {
416  return RoundHalfIntegerToEven_64(2 * x - 0.5) >> 1;
417 }
418 inline int64_t
419 Floor_64(float x)
420 {
421  return RoundHalfIntegerToEven_64(2 * x - 0.5f) >> 1;
422 }
423 
424 inline int64_t
425 Ceil_64(double x)
426 {
427  return -(RoundHalfIntegerToEven_64(-0.5 - 2 * x) >> 1);
428 }
429 inline int64_t
430 Ceil_64(float x)
431 {
432  return -(RoundHalfIntegerToEven_64(-0.5f - 2 * x) >> 1);
433 }
434 
435 #else // Base implementation
436 
437 inline int64_t
439 {
440  return RoundHalfIntegerUp_base<int64_t, double>(x);
441 }
442 inline int64_t
444 {
445  return RoundHalfIntegerUp_base<int64_t, float>(x);
446 }
447 
448 inline int64_t
449 Floor_64(double x)
450 {
451  return Floor_base<int64_t, double>(x);
452 }
453 inline int64_t
454 Floor_64(float x)
455 {
456  return Floor_base<int64_t, float>(x);
457 }
458 
459 inline int64_t
460 Ceil_64(double x)
461 {
462  return Ceil_base<int64_t, double>(x);
463 }
464 inline int64_t
465 Ceil_64(float x)
466 {
467  return Ceil_base<int64_t, float>(x);
468 }
469 
470 #endif // USE_SSE2_64IMPL || GCC_USE_ASM_64IMPL || VC_USE_ASM_64IMPL
471 
472 template <typename T>
474 
475 template <>
476 struct FloatIEEETraits<float>
477 {
478  using IntType = int32_t;
480 };
481 
482 template <>
483 struct FloatIEEETraits<double>
484 {
485  using IntType = int64_t;
487 };
488 
489 template <typename T>
491 {
492  using FloatType = T;
495 
499 
501  : asFloat(f)
502  {}
504  : asInt(i)
505  {}
506  bool
507  Sign() const
508  {
509  return (asUInt >> (sizeof(asUInt) * 8 - 1)) != 0;
510  }
511  IntType
512  AsULP() const
513  {
514  return this->Sign() ? IntType(~(~UIntType(0) >> 1) - asUInt) : asInt;
515  }
516 };
517 
518 } // end namespace Detail
519 } // end namespace Math
520 
521 // move to itkConceptChecking?
522 namespace Concept
523 {
524 template <typename T>
526 template <>
527 struct FloatOrDouble<float>
528 {};
529 template <>
530 struct FloatOrDouble<double>
531 {};
532 } // end namespace Concept
533 } // end namespace itk
534 
535 #undef USE_SSE2_32IMPL
536 #undef GCC_USE_ASM_32IMPL
537 #undef VC_USE_ASM_32IMPL
538 
539 #undef USE_SSE2_64IMPL
540 #undef GCC_USE_ASM_64IMPL
541 #undef VC_USE_ASM_64IMPL
542 
543 #endif // end of itkMathDetail.h
itk::Math::Detail::FloatIEEETraits< float >::IntType
int32_t IntType
Definition: itkMathDetail.h:478
itk::Math::Detail::RoundHalfIntegerUp_base
TReturn RoundHalfIntegerUp_base(TInput x)
Definition: itkMathDetail.h:111
itk::Math::Detail::FloatIEEE::FloatType
T FloatType
Definition: itkMathDetail.h:492
itk::Math::Detail::Floor_base
TReturn Floor_base(TInput x)
Definition: itkMathDetail.h:121
itk::Math::Detail::FloatIEEE::IntType
typename FloatIEEETraits< T >::IntType IntType
Definition: itkMathDetail.h:493
itk::Math::Detail::RoundHalfIntegerUp_32
int32_t RoundHalfIntegerUp_32(double x)
Definition: itkMathDetail.h:272
itk::Math::Detail::FloatIEEE
Definition: itkMathDetail.h:490
itk::Math::Detail::FloatIEEE::AsULP
IntType AsULP() const
Definition: itkMathDetail.h:512
itk::Math::Detail::FloatIEEE::Sign
bool Sign() const
Definition: itkMathDetail.h:507
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:438
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:460
itk::Math::Detail::RoundHalfIntegerToEven_base
CLANG_PRAGMA_PUSH CLANG_SUPPRESS_Wfloat_equal TReturn RoundHalfIntegerToEven_base(TInput x)
Definition: itkMathDetail.h:94
itk::Math::Detail::FloatIEEE::asInt
IntType asInt
Definition: itkMathDetail.h:497
itk::Math::Detail::FloatIEEETraits
Definition: itkMathDetail.h:473
itk::Math::Detail::FloatIEEE::asUInt
UIntType asUInt
Definition: itkMathDetail.h:498
itk::Math::Detail::FloatIEEETraits< double >::IntType
int64_t IntType
Definition: itkMathDetail.h:485
itk::Math::Detail::Floor_32
int32_t Floor_32(double x)
Definition: itkMathDetail.h:283
itk::Math::Detail::FloatIEEE::FloatIEEE
FloatIEEE(IntType i)
Definition: itkMathDetail.h:503
itk::Math::Detail::FloatIEEETraits< float >::UIntType
uint32_t UIntType
Definition: itkMathDetail.h:479
itk::Concept::FloatOrDouble
Definition: itkMathDetail.h:525
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:486
itk::Math::Detail::Floor_64
int64_t Floor_64(double x)
Definition: itkMathDetail.h:449
itk::Math::Detail::FloatIEEE::UIntType
typename FloatIEEETraits< T >::UIntType UIntType
Definition: itkMathDetail.h:494
itk::Math::Detail::RoundHalfIntegerToEven_32
CLANG_PRAGMA_POP int32_t RoundHalfIntegerToEven_32(double x)
Definition: itkMathDetail.h:222
itk::Math::Detail::FloatIEEE::FloatIEEE
FloatIEEE(FloatType f)
Definition: itkMathDetail.h:500
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:294
itk::Math::Detail::RoundHalfIntegerToEven_64
int64_t RoundHalfIntegerToEven_64(double x)
Definition: itkMathDetail.h:388
itk::Math::Detail::Ceil_base
TReturn Ceil_base(TInput x)
Definition: itkMathDetail.h:131
itk::Math::Detail::FloatIEEE::asFloat
FloatType asFloat
Definition: itkMathDetail.h:496