ITK 6.0.0
Insight Toolkit
 
Loading...
Searching...
No Matches
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{
77
78namespace Math::Detail
79{
80// The functions defined in this namespace are not meant to be used directly
81// and thus do not adhere to the standard backward-compatibility
82// policy of ITK, as any Detail namespace should be considered private.
83// Please use the functions from the itk::Math namespace instead
84
86// Base versions
87
88ITK_GCC_PRAGMA_PUSH
89ITK_GCC_SUPPRESS_Wfloat_equal
90
91template <typename TReturn, typename TInput>
92inline TReturn
94{
96 {
97 x += static_cast<TInput>(0.5);
98 }
99 else
100 {
101 x -= static_cast<TInput>(0.5);
102 }
103
104 const auto r = static_cast<TReturn>(x);
105 return (x != static_cast<TInput>(r)) ? r : static_cast<TReturn>(2 * (r / 2));
106}
107
108template <typename TReturn, typename TInput>
109inline TReturn
111{
112 x += static_cast<TInput>(0.5);
113 const auto r = static_cast<TReturn>(x);
115 : (x == static_cast<TInput>(r) ? r : r - static_cast<TReturn>(1));
116}
117
118template <typename TReturn, typename TInput>
119inline TReturn
120Floor_base(TInput x)
121{
122 const auto r = static_cast<TReturn>(x);
123
125 : (x == static_cast<TInput>(r) ? r : r - static_cast<TReturn>(1));
126}
127
128template <typename TReturn, typename TInput>
129inline TReturn
130Ceil_base(TInput x)
131{
132 const auto r = static_cast<TReturn>(x);
133
134 return (NumericTraits<TInput>::IsNegative(x)) ? r : (x == static_cast<TInput>(r) ? r : r + static_cast<TReturn>(1));
135}
136
137ITK_GCC_PRAGMA_POP
138
140// 32 bits versions
141
142#if USE_SSE2_32IMPL // SSE2 implementation
143
144inline int32_t
146{
147# if defined(ITK_CHECK_FPU_ROUNDING_MODE)
148 itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
149# endif
150 return _mm_cvtsd_si32(_mm_set_sd(x));
151}
152
153inline int32_t
155{
156# if defined(ITK_CHECK_FPU_ROUNDING_MODE)
157 itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
158# endif
159 return _mm_cvtss_si32(_mm_set_ss(x));
160}
161
162#elif GCC_USE_ASM_32IMPL // GCC/clang x86 asm implementation
163
164inline int32_t
166{
167# if defined(ITK_CHECK_FPU_ROUNDING_MODE)
168 itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
169# endif
170 int32_t r;
171 __asm__ __volatile__("fistpl %0" : "=m"(r) : "t"(x) : "st");
172 return r;
173}
174
175inline int32_t
177{
178# if defined(ITK_CHECK_FPU_ROUNDING_MODE)
179 itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
180# endif
181 int32_t r;
182 __asm__ __volatile__("fistpl %0" : "=m"(r) : "t"(x) : "st");
183 return r;
184}
185
186#elif VC_USE_ASM_32IMPL // msvc asm implementation
187
188inline int32_t
190{
191# if defined(ITK_CHECK_FPU_ROUNDING_MODE)
192 itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
193# endif
194 int32_t r;
195 __asm
196 {
197 fld x
198 fistp r
199 }
200 return r;
201}
202
203inline int32_t
205{
206# if defined(ITK_CHECK_FPU_ROUNDING_MODE)
207 itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
208# endif
209 int32_t r;
210 __asm
211 {
212 fld x
213 fistp r
214 }
215 return r;
216}
217
218#else // Base implementation
219
220inline int32_t
225inline int32_t
230
231#endif
232
233#if USE_SSE2_32IMPL || GCC_USE_ASM_32IMPL || VC_USE_ASM_32IMPL
234
235inline int32_t
236RoundHalfIntegerUp_32(double x)
237{
238 return RoundHalfIntegerToEven_32(2 * x + 0.5) >> 1;
239}
240inline int32_t
242{
243 return RoundHalfIntegerToEven_32(2 * x + 0.5f) >> 1;
244}
245
246inline int32_t
247Floor_32(double x)
248{
249 return RoundHalfIntegerToEven_32(2 * x - 0.5) >> 1;
250}
251inline int32_t
252Floor_32(float x)
253{
254 return RoundHalfIntegerToEven_32(2 * x - 0.5f) >> 1;
255}
256
257inline int32_t
258Ceil_32(double x)
259{
260 return -(RoundHalfIntegerToEven_32(-0.5 - 2 * x) >> 1);
261}
262inline int32_t
263Ceil_32(float x)
264{
265 return -(RoundHalfIntegerToEven_32(-0.5f - 2 * x) >> 1);
266}
267
268#else // Base implementation
269
270inline int32_t
275inline int32_t
280
281inline int32_t
282Floor_32(double x)
283{
285}
286inline int32_t
287Floor_32(float x)
288{
290}
291
292inline int32_t
293Ceil_32(double x)
294{
296}
297inline int32_t
298Ceil_32(float x)
299{
301}
302
303#endif // USE_SSE2_32IMPL || GCC_USE_ASM_32IMPL || VC_USE_ASM_32IMPL
304
306// 64 bits versions
307
308#if USE_SSE2_64IMPL // SSE2 implementation
309
310inline int64_t
312{
313# if defined(ITK_CHECK_FPU_ROUNDING_MODE)
314 itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
315# endif
316 return _mm_cvtsd_si64(_mm_set_sd(x));
317}
318
319inline int64_t
321{
322# if defined(ITK_CHECK_FPU_ROUNDING_MODE)
323 itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
324# endif
325 return _mm_cvtss_si64(_mm_set_ss(x));
326}
327
328#elif GCC_USE_ASM_64IMPL // GCC/clang x86 asm implementation
329
330inline int64_t
332{
333# if defined(ITK_CHECK_FPU_ROUNDING_MODE)
334 itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
335# endif
336 int64_t r;
337 __asm__ __volatile__("fistpll %0" : "=m"(r) : "t"(x) : "st");
338 return r;
339}
340
341inline int64_t
343{
344# if defined(ITK_CHECK_FPU_ROUNDING_MODE)
345 itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
346# endif
347 int64_t r;
348 __asm__ __volatile__("fistpll %0" : "=m"(r) : "t"(x) : "st");
349 return r;
350}
351
352#elif VC_USE_ASM_64IMPL // msvc asm implementation
353
354inline int64_t
356{
357# if defined(ITK_CHECK_FPU_ROUNDING_MODE)
358 itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
359# endif
360 int64_t r;
361 __asm
362 {
363 fld x
364 fistp r
365 }
366 return r;
367}
368
369inline int64_t
371{
372# if defined(ITK_CHECK_FPU_ROUNDING_MODE)
373 itkAssertInDebugAndIgnoreInReleaseMacro(fegetround() == FE_TONEAREST);
374# endif
375 int64_t r;
376 __asm
377 {
378 fld x
379 fistp r
380 }
381 return r;
382}
383
384#else // Base implementation
385
386inline int64_t
391inline int64_t
396
397#endif
398
399#if USE_SSE2_64IMPL || GCC_USE_ASM_64IMPL || VC_USE_ASM_64IMPL
400
401inline int64_t
402RoundHalfIntegerUp_64(double x)
403{
404 return RoundHalfIntegerToEven_64(2 * x + 0.5) >> 1;
405}
406inline int64_t
408{
409 return RoundHalfIntegerToEven_64(2 * x + 0.5f) >> 1;
410}
411
412inline int64_t
413Floor_64(double x)
414{
415 return RoundHalfIntegerToEven_64(2 * x - 0.5) >> 1;
416}
417inline int64_t
418Floor_64(float x)
419{
420 return RoundHalfIntegerToEven_64(2 * x - 0.5f) >> 1;
421}
422
423inline int64_t
424Ceil_64(double x)
425{
426 return -(RoundHalfIntegerToEven_64(-0.5 - 2 * x) >> 1);
427}
428inline int64_t
429Ceil_64(float x)
430{
431 return -(RoundHalfIntegerToEven_64(-0.5f - 2 * x) >> 1);
432}
433
434#else // Base implementation
435
436inline int64_t
441inline int64_t
446
447inline int64_t
448Floor_64(double x)
449{
451}
452inline int64_t
453Floor_64(float x)
454{
456}
457
458inline int64_t
459Ceil_64(double x)
460{
462}
463inline int64_t
464Ceil_64(float x)
465{
467}
468
469#endif // USE_SSE2_64IMPL || GCC_USE_ASM_64IMPL || VC_USE_ASM_64IMPL
470
471template <typename T>
473
474template <>
475struct FloatIEEETraits<float>
476{
477 using IntType = int32_t;
478 using UIntType = uint32_t;
479};
480
481template <>
482struct FloatIEEETraits<double>
483{
484 using IntType = int64_t;
485 using UIntType = uint64_t;
486};
487
488template <typename T>
490{
491 using FloatType = T;
494
498
500 : asFloat(f)
501 {}
503 : asInt(i)
504 {}
505 [[nodiscard]] bool
506 Sign() const
507 {
508 return (asUInt >> (sizeof(asUInt) * 8 - 1)) != 0;
509 }
510 [[nodiscard]] IntType
511 AsULP() const
512 {
513 return this->Sign() ? IntType(~(~UIntType(0) >> 1) - asUInt) : asInt;
514 }
515};
516
517} // namespace Math::Detail
518// end namespace Math
519
520// move to itkConceptChecking?
521namespace Concept
522{
523template <typename T>
525template <>
526struct FloatOrDouble<float>
527{};
528template <>
529struct FloatOrDouble<double>
530{};
531} // end namespace Concept
532} // end namespace itk
533
534#undef USE_SSE2_32IMPL
535#undef GCC_USE_ASM_32IMPL
536#undef VC_USE_ASM_32IMPL
537
538#undef USE_SSE2_64IMPL
539#undef GCC_USE_ASM_64IMPL
540#undef VC_USE_ASM_64IMPL
541
542#endif // end of itkMathDetail.h
static bool IsNegative(T val)
static bool IsNonnegative(T val)
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)
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