ITK  6.0.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 * https://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
75namespace itk
76{
77namespace Math
78{
79namespace 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
89ITK_GCC_PRAGMA_PUSH
90ITK_GCC_SUPPRESS_Wfloat_equal
91
92template <typename TReturn, typename TInput>
93inline 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
109template <typename TReturn, typename TInput>
110inline 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
119template <typename TReturn, typename TInput>
120inline TReturn
121Floor_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
129template <typename TReturn, typename TInput>
130inline TReturn
131Ceil_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
138ITK_GCC_PRAGMA_POP
139
141// 32 bits versions
142
143#if USE_SSE2_32IMPL // SSE2 implementation
144
145inline 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
154inline 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
165inline 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
176inline 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
189inline 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
204inline 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
221inline int32_t
223{
224 return RoundHalfIntegerToEven_base<int32_t, double>(x);
225}
226inline 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
236inline int32_t
237RoundHalfIntegerUp_32(double x)
238{
239 return RoundHalfIntegerToEven_32(2 * x + 0.5) >> 1;
240}
241inline int32_t
243{
244 return RoundHalfIntegerToEven_32(2 * x + 0.5f) >> 1;
245}
246
247inline int32_t
248Floor_32(double x)
249{
250 return RoundHalfIntegerToEven_32(2 * x - 0.5) >> 1;
251}
252inline int32_t
253Floor_32(float x)
254{
255 return RoundHalfIntegerToEven_32(2 * x - 0.5f) >> 1;
256}
257
258inline int32_t
259Ceil_32(double x)
260{
261 return -(RoundHalfIntegerToEven_32(-0.5 - 2 * x) >> 1);
262}
263inline int32_t
264Ceil_32(float x)
265{
266 return -(RoundHalfIntegerToEven_32(-0.5f - 2 * x) >> 1);
267}
268
269#else // Base implementation
270
271inline int32_t
273{
274 return RoundHalfIntegerUp_base<int32_t, double>(x);
275}
276inline int32_t
278{
279 return RoundHalfIntegerUp_base<int32_t, float>(x);
280}
281
282inline int32_t
283Floor_32(double x)
284{
285 return Floor_base<int32_t, double>(x);
286}
287inline int32_t
288Floor_32(float x)
289{
290 return Floor_base<int32_t, float>(x);
291}
292
293inline int32_t
294Ceil_32(double x)
295{
296 return Ceil_base<int32_t, double>(x);
297}
298inline int32_t
299Ceil_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
311inline 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
320inline 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
331inline 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
342inline 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
355inline 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
370inline 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
387inline int64_t
389{
390 return RoundHalfIntegerToEven_base<int64_t, double>(x);
391}
392inline 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
402inline int64_t
403RoundHalfIntegerUp_64(double x)
404{
405 return RoundHalfIntegerToEven_64(2 * x + 0.5) >> 1;
406}
407inline int64_t
409{
410 return RoundHalfIntegerToEven_64(2 * x + 0.5f) >> 1;
411}
412
413inline int64_t
414Floor_64(double x)
415{
416 return RoundHalfIntegerToEven_64(2 * x - 0.5) >> 1;
417}
418inline int64_t
419Floor_64(float x)
420{
421 return RoundHalfIntegerToEven_64(2 * x - 0.5f) >> 1;
422}
423
424inline int64_t
425Ceil_64(double x)
426{
427 return -(RoundHalfIntegerToEven_64(-0.5 - 2 * x) >> 1);
428}
429inline int64_t
430Ceil_64(float x)
431{
432 return -(RoundHalfIntegerToEven_64(-0.5f - 2 * x) >> 1);
433}
434
435#else // Base implementation
436
437inline int64_t
439{
440 return RoundHalfIntegerUp_base<int64_t, double>(x);
441}
442inline int64_t
444{
445 return RoundHalfIntegerUp_base<int64_t, float>(x);
446}
447
448inline int64_t
449Floor_64(double x)
450{
451 return Floor_base<int64_t, double>(x);
452}
453inline int64_t
454Floor_64(float x)
455{
456 return Floor_base<int64_t, float>(x);
457}
458
459inline int64_t
460Ceil_64(double x)
461{
462 return Ceil_base<int64_t, double>(x);
463}
464inline int64_t
465Ceil_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
472template <typename T>
474
475template <>
476struct FloatIEEETraits<float>
477{
478 using IntType = int32_t;
479 using UIntType = uint32_t;
480};
481
482template <>
483struct FloatIEEETraits<double>
484{
485 using IntType = int64_t;
486 using UIntType = uint64_t;
487};
488
489template <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?
522namespace Concept
523{
524template <typename T>
526template <>
527struct FloatOrDouble<float>
528{};
529template <>
530struct 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
Define additional traits for native types such as int or float.
ITK_GCC_PRAGMA_POP int32_t RoundHalfIntegerToEven_32(double x)
int32_t Floor_32(double x)
TReturn RoundHalfIntegerUp_base(TInput x)
int32_t RoundHalfIntegerUp_32(double x)
TReturn Ceil_base(TInput x)
int64_t Floor_64(double x)
int64_t Ceil_64(double x)
int32_t Ceil_32(double x)
int64_t RoundHalfIntegerUp_64(double x)
int64_t RoundHalfIntegerToEven_64(double x)
ITK_GCC_PRAGMA_PUSH ITK_GCC_SUPPRESS_Wfloat_equal TReturn RoundHalfIntegerToEven_base(TInput x)
Definition: itkMathDetail.h:94
TReturn Floor_base(TInput x)
The "itk" namespace contains all Insight Segmentation and Registration Toolkit (ITK) classes....
typename FloatIEEETraits< T >::IntType IntType
typename FloatIEEETraits< T >::UIntType UIntType