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 itkMath_h
29 #define itkMath_h
30
31 #include "itkMathDetail.h"
32 #include "itkConceptChecking.h"
33 #include <vnl/vnl_math.h>
34
35 /* Only maintain backwards compatibility with old versions
36 * of VXL back to the point where vnl_math:: was introduced
37 * versions of VXL where only vnl_math_ was available are not
38 * supported.
39 */
40 #include <vxl_version.h>
41
42 namespace itk
43 {
44 namespace Math
45 {
46 // These constants originate from VXL's vnl_math.h. They have been
47 // moved here to improve visibility, and to ensure that the constants
48 // are available during compile time ( as opposed to static ITK_CONSTEXPR
49 // member vaiables ).
50
51
52 /** \brief \f[e\f] The base of the natural logarithm or Euler's number */
53 static constexpr double e = vnl_math::e;
54 /** \brief \f[ \log_2 e \f] */
55 static constexpr double log2e = vnl_math::log2e;
56 /** \brief \f[ \log_{10} e \f] */
57 static constexpr double log10e = vnl_math::log10e;
58 /** \brief \f[ \log_e 2 \f] */
59 static constexpr double ln2 = vnl_math::ln2;
60 /** \brief \f[ \log_e 10 \f] */
61 static constexpr double ln10 = vnl_math::ln10;
62 /** \brief \f[ \pi \f] */
63 static constexpr double pi = vnl_math::pi;
64 /** \brief \f[ 2\pi \f] */
65 static constexpr double twopi = vnl_math::twopi;
66 /** \brief \f[ \frac{\pi}{2} \f] */
67 static constexpr double pi_over_2 = vnl_math::pi_over_2;
68 /** \brief \f[ \frac{\pi}{4} \f] */
69 static constexpr double pi_over_4 = vnl_math::pi_over_4;
70 /** \brief \f[ \frac{\pi}{180} \f] */
71 static constexpr double pi_over_180 = vnl_math::pi_over_180;
72 /** \brief \f[ \frac{1}{\pi} \f] */
73 static constexpr double one_over_pi = vnl_math::one_over_pi;
74 /** \brief \f[ \frac{2}{\pi} \f] */
75 static constexpr double two_over_pi = vnl_math::two_over_pi;
76 /** \brief \f[ \frac{180}{\pi} \f] */
77 static constexpr double deg_per_rad = vnl_math::deg_per_rad;
78 /** \brief \f[ \sqrt{2\pi} \f] */
79 static constexpr double sqrt2pi = vnl_math::sqrt2pi;
80 /** \brief \f[ \frac{2}{\sqrt{\pi}} \f] */
81 static constexpr double two_over_sqrtpi = vnl_math::two_over_sqrtpi;
82 /** \brief \f[ \frac{2}{\sqrt{2\pi}} \f] */
83 static constexpr double one_over_sqrt2pi = vnl_math::one_over_sqrt2pi;
84 /** \brief \f[ \sqrt{2} \f] */
85 static constexpr double sqrt2 = vnl_math::sqrt2;
86 /** \brief \f[ \sqrt{ \frac{1}{2}} \f] */
87 static constexpr double sqrt1_2 = vnl_math::sqrt1_2;
88 /** \brief \f[ \sqrt{ \frac{1}{3}} \f] */
89 static constexpr double sqrt1_3 = vnl_math::sqrt1_3;
90 /** \brief euler constant */
91 static constexpr double euler = vnl_math::euler;
92
93 //: IEEE double machine precision
94 static constexpr double eps = vnl_math::eps;
95 static constexpr double sqrteps = vnl_math::sqrteps;
96 //: IEEE single machine precision
97 static constexpr float float_eps = vnl_math::float_eps;
98 static constexpr float float_sqrteps = vnl_math::float_sqrteps;
99
100 /** A useful macro to generate a template floating point to integer
101 * conversion templated on the return type and using either the 32
102 * bit, the 64 bit or the vanilla version */
103 #define itkTemplateFloatingToIntegerMacro(name) \
104 template< typename TReturn, typename TInput > \
105 inline TReturn name(TInput x) \
106 { \
107 \
108 if ( sizeof( TReturn ) <= 4 ) \
109 { \
110 return static_cast< TReturn >( Detail::name##_32(x) ); \
111 } \
112 else if ( sizeof( TReturn ) <= 8 ) \
113 { \
114 return static_cast< TReturn >( Detail::name##_64(x) ); \
115 } \
116 else \
117 { \
118 return static_cast< TReturn >( Detail::name##_base< TReturn, TInput >(x) ); \
119 } \
120 }
121
122 /** \brief Round towards nearest integer
123 *
124 * \tparam TReturn must be an integer type
125 * \tparam TInput must be float or double
126 *
127 * halfway cases are rounded towards the nearest even
128 * integer, e.g.
129 \code
130 RoundHalfIntegerToEven( 1.5) == 2
131 RoundHalfIntegerToEven(-1.5) == -2
132 RoundHalfIntegerToEven( 2.5) == 2
133 RoundHalfIntegerToEven( 3.5) == 4
134 \endcode
135 *
136 * The behavior of overflow is undefined due to numerous implementations.
137 *
138 * \warning We assume that the rounding mode is not changed from the default
139 * one (or at least that it is always restored to the default one).
140 */
141 itkTemplateFloatingToIntegerMacro(RoundHalfIntegerToEven);
142
143 /** \brief Round towards nearest integer
144 *
145 * \tparam TReturn must be an integer type
146 * \tparam TInput must be float or double
147 *
148 * halfway cases are rounded upward, e.g.
149 \code
150 RoundHalfIntegerUp( 1.5) == 2
151 RoundHalfIntegerUp(-1.5) == -1
152 RoundHalfIntegerUp( 2.5) == 3
153 \endcode
154 *
155 * The behavior of overflow is undefined due to numerous implementations.
156 *
157 * \warning The argument absolute value must be less than
158 * NumbericTraits<TReturn>::max()/2 for RoundHalfIntegerUp to be
159 * guaranteed to work.
160 *
161 * \warning We also assume that the rounding mode is not changed from
162 * the default one (or at least that it is always restored to the
163 * default one).
164 */
165 itkTemplateFloatingToIntegerMacro(RoundHalfIntegerUp);
166
167 /** \brief Round towards nearest integer (This is a synonym for RoundHalfIntegerUp)
168 *
169 * \tparam TReturn must be an integer type
170 * \tparam TInput must be float or double
171 *
172 * \sa RoundHalfIntegerUp<TReturn, TInput>()
173 */
174 template< typename TReturn, typename TInput >
Round(TInput x)175 inline TReturn Round(TInput x) { return RoundHalfIntegerUp< TReturn, TInput >(x); }
176
177 /** \brief Round towards minus infinity
178 *
179 * The behavior of overflow is undefined due to numerous implementations.
180 *
181 * \warning argument absolute value must be less than
182 * NumbericTraits<TReturn>::max()/2 for vnl_math_floor to be
183 * guaranteed to work.
184 *
185 * \warning We also assume that the rounding mode is not changed from
186 * the default one (or at least that it is always restored to the
187 * default one).
188 */
189 itkTemplateFloatingToIntegerMacro(Floor);
190
191 /** \brief Round towards plus infinity
192 *
193 * The behavior of overflow is undefined due to numerous implementations.
194 *
195 * \warning argument absolute value must be less than INT_MAX/2
196 * for vnl_math_ceil to be guaranteed to work.
197 * \warning We also assume that the rounding mode is not changed from
198 * the default one (or at least that it is always restored to the
199 * default one).
200 */
201 itkTemplateFloatingToIntegerMacro(Ceil);
202
203 #undef itkTemplateFloatingToIntegerMacro
204
205 template< typename TReturn, typename TInput >
CastWithRangeCheck(TInput x)206 inline TReturn CastWithRangeCheck(TInput x)
207 {
208 #ifdef ITK_USE_CONCEPT_CHECKING
209 itkConceptMacro( OnlyDefinedForIntegerTypes1, ( itk::Concept::IsInteger< TReturn > ) );
210 itkConceptMacro( OnlyDefinedForIntegerTypes2, ( itk::Concept::IsInteger< TInput > ) );
211 #endif // ITK_USE_CONCEPT_CHECKING
212
213 auto ret = static_cast< TReturn >( x );
214 if ( sizeof( TReturn ) > sizeof( TInput )
215 && !( !itk::NumericTraits< TReturn >::is_signed && itk::NumericTraits< TInput >::is_signed ) )
216 {
217 // if the output type is bigger and we are not converting a signed
218 // integer to an unsigned integer then we have no problems
219 return ret;
220 }
221 else if ( sizeof( TReturn ) >= sizeof( TInput ) )
222 {
223 if ( itk::NumericTraits< TInput >::IsPositive(x) != itk::NumericTraits< TReturn >::IsPositive(ret) )
224 {
225 itk::RangeError _e(__FILE__, __LINE__);
226 throw _e;
227 }
228 }
229 else if ( static_cast< TInput >( ret ) != x
230 || ( itk::NumericTraits< TInput >::IsPositive(x) != itk::NumericTraits< TReturn >::IsPositive(ret) ) )
231 {
232 itk::RangeError _e(__FILE__, __LINE__);
233 throw _e;
234 }
235 return ret;
236 }
237
238 /** \brief Return the signed distance in ULPs (units in the last place) between two floats.
239 *
240 * This is the signed distance, i.e., if x1 > x2, then the result is positive.
241 *
242 * \sa FloatAlmostEqual
243 * \sa FloatAddULP
244 */
245 template <typename T>
246 inline typename Detail::FloatIEEE<T>::IntType
FloatDifferenceULP(T x1,T x2)247 FloatDifferenceULP( T x1, T x2 )
248 {
249 Detail::FloatIEEE<T> x1f(x1);
250 Detail::FloatIEEE<T> x2f(x2);
251 return x1f.AsULP() - x2f.AsULP();
252 }
253
254 /** \brief Add the given ULPs (units in the last place) to a float.
255 *
256 * If the given ULPs can are negative, they are subtracted.
257 *
258 * \sa FloatAlmostEqual
259 * \sa FloatDifferenceULP
260 */
261 template <typename T>
262 inline T
FloatAddULP(T x,typename Detail::FloatIEEE<T>::IntType ulps)263 FloatAddULP( T x, typename Detail::FloatIEEE<T>::IntType ulps )
264 {
265 Detail::FloatIEEE<T> representInput( x );
266 Detail::FloatIEEE<T> representOutput( representInput.asInt + ulps );
267 return representOutput.asFloat;
268 }
269
270 /** \brief Compare two floats and return if they are effectively equal.
271 *
272 * Determining when floats are almost equal is difficult because of their
273 * IEEE bit representation. This function uses the integer representation of
274 * the float to determine if they are almost equal.
275 *
276 * The implementation is based off the explanation in the white papers:
277 *
278 * - http://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
279 * - http://randomascii.wordpress.com/category/floating-point/
280 *
281 * This function is not a cure-all, and reading those articles is important
282 * to understand its appropriate use in the context of ULPs, zeros, subnormals,
283 * infinities, and NANs. For example, it is preferable to use this function on
284 * two floats directly instead of subtracting them and comparing them to zero.
285 *
286 * The tolerance is specified in ULPs (units in the last place), i.e. how many
287 * floats there are in between the numbers. Therefore, the tolerance depends on
288 * the magnitude of the values that are being compared. A second tolerance is
289 * a maximum difference allowed, which is important when comparing numbers close to
290 * zero.
291 *
292 * A NAN compares as not equal to a number, but two NAN's may compare as equal
293 * to each other.
294 *
295 * \param x1 first floating value to compare
296 * \param x2 second floating values to compare
297 * \param maxUlps maximum units in the last place to be considered equal
298 * \param maxAbsoluteDifference maximum absolute difference to be considered equal
299 */
300 template <typename T>
301 inline bool
302 FloatAlmostEqual( T x1, T x2,
303 typename Detail::FloatIEEE<T>::IntType maxUlps = 4,
304 typename Detail::FloatIEEE<T>::FloatType maxAbsoluteDifference = 0.1*itk::NumericTraits<T>::epsilon() )
305 {
306 // Check if the numbers are really close -- needed
307 // when comparing numbers near zero.
308 const T absDifference = std::abs(x1 - x2);
309 if ( absDifference <= maxAbsoluteDifference )
310 {
311 return true;
312 }
313
314 #if defined(__APPLE__) && (__clang_major__ == 3) && (__clang_minor__ == 0) && defined(NDEBUG) && defined(__x86_64__)
315 Detail::FloatIEEE<T> x1f(x1);
316 Detail::FloatIEEE<T> x2f(x2);
317 double x1fAsULP = static_cast<double>(x1f.AsULP());
318 double x2fAsULP = static_cast<double>(x2f.AsULP());
319 double ulps = x1fAsULP - x2fAsULP;
320 if(ulps < 0)
321 {
322 ulps = -ulps;
323 }
324 return ulps <= static_cast<double>(maxUlps);
325 #else
326 typename Detail::FloatIEEE<T>::IntType
327 ulps = FloatDifferenceULP(x1, x2);
328 if(ulps < 0)
329 {
330 ulps = -ulps;
331 }
332 return ulps <= maxUlps;
333 #endif
334 }
335
336 // The following code cannot be moved to the itkMathDetail.h file without introducing circular dependencies
337 namespace Detail // The Detail namespace holds the templates used by AlmostEquals
338 {
339 // The following structs and templates are used to choose
340 // which version of the AlmostEquals function
341 // should be implemented base on input parameter types
342
343 // Structs for choosing AlmostEquals function
344
345 struct AlmostEqualsFloatVsFloat
346 {
347 template <typename TFloatType1, typename TFloatType2>
AlmostEqualsFunctionAlmostEqualsFloatVsFloat348 static bool AlmostEqualsFunction(TFloatType1 x1, TFloatType2 x2)
349 {
350 return FloatAlmostEqual<double>(x1, x2);
351 }
352
353 template <typename TFloatType1, typename TFloatType2>
354 static bool
AlmostEqualsFunctionAlmostEqualsFloatVsFloat355 AlmostEqualsFunction(double x1, double x2)
356 {
357 return FloatAlmostEqual<double>(x1, x2);
358 }
359
360 template <typename TFloatType1, typename TFloatType2>
361 static bool
AlmostEqualsFunctionAlmostEqualsFloatVsFloat362 AlmostEqualsFunction(double x1, float x2)
363 {
364 return FloatAlmostEqual<float>(x1, x2);
365 }
366
367 template <typename TFloatType1, typename TFloatType2>
368 static bool
AlmostEqualsFunctionAlmostEqualsFloatVsFloat369 AlmostEqualsFunction(float x1, double x2)
370 {
371 return FloatAlmostEqual<float>(x1, x2);
372 }
373
374 template <typename TFloatType1, typename TFloatType2>
375 static bool
AlmostEqualsFunctionAlmostEqualsFloatVsFloat376 AlmostEqualsFunction(float x1, float x2)
377 {
378 return FloatAlmostEqual<float>(x1, x2);
379 }
380 };
381
382 struct AlmostEqualsFloatVsInteger
383 {
384 template <typename TFloatType, typename TIntType>
AlmostEqualsFunctionAlmostEqualsFloatVsInteger385 static bool AlmostEqualsFunction(TFloatType floatingVariable, TIntType integerVariable)
386 {
387 return FloatAlmostEqual<TFloatType> (floatingVariable, integerVariable);
388 }
389 };
390
391 struct AlmostEqualsIntegerVsFloat
392 {
393 template <typename TIntType, typename TFloatType>
AlmostEqualsFunctionAlmostEqualsIntegerVsFloat394 static bool AlmostEqualsFunction(TIntType integerVariable, TFloatType floatingVariable)
395 {
396 return AlmostEqualsFloatVsInteger::AlmostEqualsFunction(floatingVariable, integerVariable);
397 }
398 };
399
400 struct AlmostEqualsSignedVsUnsigned
401 {
402 template <typename TSignedInt, typename TUnsignedInt>
AlmostEqualsFunctionAlmostEqualsSignedVsUnsigned403 static bool AlmostEqualsFunction(TSignedInt signedVariable, TUnsignedInt unsignedVariable)
404 {
405 if(signedVariable < 0) return false;
406 if( unsignedVariable > static_cast< size_t >(itk::NumericTraits<TSignedInt>::max()) ) return false;
407 return signedVariable == static_cast< TSignedInt >(unsignedVariable);
408 }
409 };
410
411 struct AlmostEqualsUnsignedVsSigned
412 {
413 template <typename TUnsignedInt, typename TSignedInt>
AlmostEqualsFunctionAlmostEqualsUnsignedVsSigned414 static bool AlmostEqualsFunction(TUnsignedInt unsignedVariable, TSignedInt signedVariable)
415 {
416 return AlmostEqualsSignedVsUnsigned::AlmostEqualsFunction(signedVariable, unsignedVariable);
417 }
418 };
419
420 struct AlmostEqualsPlainOldEquals
421 {
422 template <typename TIntegerType1, typename TIntegerType2>
AlmostEqualsFunctionAlmostEqualsPlainOldEquals423 static bool AlmostEqualsFunction(TIntegerType1 x1, TIntegerType2 x2)
424 {
425 return x1 == x2;
426 }
427 };
428 // end of structs that choose the specific AlmostEquals function
429
430 // Selector structs, these select the correct case based on its types
431 // input1 is int? input 1 is signed? input2 is int? input 2 is signed?
432 template<bool TInput1IsIntger, bool TInput1IsSigned, bool TInput2IsInteger, bool TInput2IsSigned>
433 struct AlmostEqualsFunctionSelector
434 { // default case
435 using SelectedVersion = AlmostEqualsPlainOldEquals;
436 };
437
438 /// \cond HIDE_SPECIALIZATION_DOCUMENTATION
439 template<>
440 struct AlmostEqualsFunctionSelector < false, true, false, true>
441 // floating type v floating type
442 {
443 using SelectedVersion = AlmostEqualsFloatVsFloat;
444 };
445
446 template<>
447 struct AlmostEqualsFunctionSelector <false, true, true, true>
448 // float vs signed int
449 {
450 using SelectedVersion = AlmostEqualsFloatVsInteger;
451 };
452
453 template<>
454 struct AlmostEqualsFunctionSelector <false, true, true,false>
455 // float vs unsigned int
456 {
457 using SelectedVersion = AlmostEqualsFloatVsInteger;
458 };
459
460 template<>
461 struct AlmostEqualsFunctionSelector <true, false, false, true>
462 // unsigned int vs float
463 {
464 using SelectedVersion = AlmostEqualsIntegerVsFloat;
465 };
466
467 template<>
468 struct AlmostEqualsFunctionSelector <true, true, false, true>
469 // signed int vs float
470 {
471 using SelectedVersion = AlmostEqualsIntegerVsFloat;
472 };
473
474 template<>
475 struct AlmostEqualsFunctionSelector<true, true, true, false>
476 // signed vs unsigned
477 {
478 using SelectedVersion = AlmostEqualsSignedVsUnsigned;
479 };
480
481 template<>
482 struct AlmostEqualsFunctionSelector<true, false, true, true>
483 // unsigned vs signed
484 {
485 using SelectedVersion = AlmostEqualsUnsignedVsSigned;
486 };
487
488 template<>
489 struct AlmostEqualsFunctionSelector<true, true, true, true>
490 // signed vs signed
491 {
492 using SelectedVersion = AlmostEqualsPlainOldEquals;
493 };
494
495 template<>
496 struct AlmostEqualsFunctionSelector<true, false, true, false>
497 // unsigned vs unsigned
498 {
499 using SelectedVersion = AlmostEqualsPlainOldEquals;
500 };
501 // end of AlmostEqualsFunctionSelector structs
502
503 // The implementor tells the selector what to do
504 template<typename TInputType1, typename TInputType2>
505 struct AlmostEqualsScalarImplementer
506 {
507 static constexpr bool TInputType1IsInteger = itk::NumericTraits<TInputType1>::IsInteger;
508 static constexpr bool TInputType1IsSigned = itk::NumericTraits<TInputType1>::IsSigned;
509 static constexpr bool TInputType2IsInteger = itk::NumericTraits<TInputType2>::IsInteger;
510 static constexpr bool TInputType2IsSigned = itk::NumericTraits<TInputType2>::IsSigned;
511
512 using SelectedVersion = typename AlmostEqualsFunctionSelector< TInputType1IsInteger, TInputType1IsSigned,
513 TInputType2IsInteger, TInputType2IsSigned >::SelectedVersion;
514 };
515
516 // The AlmostEqualsScalarComparer returns the result of an
517 // approximate comparison between two scalar values of
518 // potentially different data types.
519 template <typename TScalarType1, typename TScalarType2>
520 inline bool
521 AlmostEqualsScalarComparer( TScalarType1 x1, TScalarType2 x2 )
522 {
523 return AlmostEqualsScalarImplementer<TScalarType1, TScalarType2>::SelectedVersion:: template AlmostEqualsFunction<TScalarType1, TScalarType2>(x1, x2);
524 }
525
526 // The following structs are used to evaluate approximate comparisons between
527 // complex and scalar values of potentially different types.
528
529 // Comparisons between scalar types use the AlmostEqualsScalarComparer function.
530 struct AlmostEqualsScalarVsScalar
531 {
532 template <typename TScalarType1, typename TScalarType2>
533 static bool
534 AlmostEqualsFunction(TScalarType1 x1, TScalarType2 x2)
535 {
536 return AlmostEqualsScalarComparer(x1, x2);
537 }
538 };
539
540 // Comparisons between two complex values compare the real and imaginary components
541 // separately with the AlmostEqualsScalarComparer function.
542 struct AlmostEqualsComplexVsComplex
543 {
544 template <typename TComplexType1, typename TComplexType2>
545 static bool
546 AlmostEqualsFunction(TComplexType1 x1, TComplexType2 x2)
547 {
548 return AlmostEqualsScalarComparer(x1.real(), x2.real()) && AlmostEqualsScalarComparer( x1.imag(), x2.imag() );
549 }
550 };
551
552 // Comparisons between complex and scalar values first check to see if the imaginary component
553 // of the complex value is approximately 0. Then a ScalarComparison is done between the real
554 // part of the complex number and the scalar value.
555 struct AlmostEqualsScalarVsComplex
556 {
557 template <typename TScalarType, typename TComplexType>
558 static bool
559 AlmostEqualsFunction(TScalarType scalarVariable, TComplexType complexVariable)
560 {
561 if( !AlmostEqualsScalarComparer( complexVariable.imag(), itk::NumericTraits< typename itk::NumericTraits< TComplexType >::ValueType >::ZeroValue() ) )
562 {
563 return false;
564 }
565 return AlmostEqualsScalarComparer(scalarVariable, complexVariable.real());
566 }
567 };
568
569 struct AlmostEqualsComplexVsScalar
570 {
571 template <typename TComplexType, typename TScalarType>
572 static bool
573 AlmostEqualsFunction(TComplexType complexVariable, TScalarType scalarVariable)
574 {
575 return AlmostEqualsScalarVsComplex::AlmostEqualsFunction(scalarVariable, complexVariable);
576 }
577 };
578
579 // The AlmostEqualsComplexChooser structs choose the correct case
580 // from the input parameter types' IsComplex property
581 // The default case is scalar vs scalar
582 template < bool T1IsComplex, bool T2IsComplex > //Default is false, false
583 struct AlmostEqualsComplexChooser
584 {
585 using ChosenVersion = AlmostEqualsScalarVsScalar;
586 };
587
588 template <>
589 struct AlmostEqualsComplexChooser< true, true >
590 {
591 using ChosenVersion = AlmostEqualsComplexVsComplex;
592 };
593
594 template <>
595 struct AlmostEqualsComplexChooser< false, true >
596 {
597 using ChosenVersion = AlmostEqualsScalarVsComplex;
598 };
599
600 template <>
601 struct AlmostEqualsComplexChooser< true, false>
602 {
603 using ChosenVersion = AlmostEqualsComplexVsScalar;
604 };
605 // End of AlmostEqualsComplexChooser structs.
606
607 // The AlmostEqualsComplexImplementer determines which of the input
608 // parameters are complex and which are real, and sends that information
609 // to the AlmostEqualsComplexChooser structs to determine the proper
610 // type of evaluation.
611 template <typename T1, typename T2>
612 struct AlmostEqualsComplexImplementer
613 {
614 static constexpr bool T1IsComplex = NumericTraits< T1 >::IsComplex;
615 static constexpr bool T2IsComplex = NumericTraits< T2 >::IsComplex;
616
617 using ChosenVersion = typename AlmostEqualsComplexChooser< T1IsComplex, T2IsComplex >::ChosenVersion;
618 };
619 /// \endcond
620
621 } // end namespace Detail
622
623 /** \brief Provide consistent equality checks between values of potentially different scalar or complex types
624 *
625 * template< typename T1, typename T2 >
626 * AlmostEquals( T1 x1, T2 x2 )
627 *
628 * template< typename T1, typename T2 >
629 * NotAlmostEquals( T1 x1, T2 x2 )
630 *
631 * This function compares two scalar or complex values of potentially different types.
632 * For maximum extensibility the function is implemented through a series of templated
633 * structs which direct the AlmostEquals() call to the correct function by evaluating
634 * the parameter's types.
635 *
636 * Overall algorithm:
637 * If both values are complex...
638 * separate values into real and imaginary components and compare them separately
639 *
640 * If one of the values is complex..
641 * see if the imaginary part of the complex value is approximately 0 ...
642 * compare real part of complex value with scalar value
643 *
644 * If both values are scalars...
645 *
646 * To compare two floating point types...
647 * use FloatAlmostEqual.
648 *
649 * To compare a floating point and an integer type...
650 * Use static_cast<FloatingPointType>(integerValue) and call FloatAlmostEqual
651 *
652 * To compare signed and unsigned integers...
653 * Check for negative value or overflow, then cast and use ==
654 *
655 * To compare two signed or two unsigned integers ...
656 * Use ==
657 *
658 * To compare anything else ...
659 * Use ==
660 *
661 * \param x1 first scalar value to compare
662 * \param x2 second scalar value to compare
663 */
664
665 // The AlmostEquals function
666 template <typename T1, typename T2>
667 inline bool
668 AlmostEquals( T1 x1, T2 x2 )
669 {
670 return Detail::AlmostEqualsComplexImplementer<T1,T2>::ChosenVersion::AlmostEqualsFunction(x1, x2);
671 }
672
673 // The NotAlmostEquals function
674 template <typename T1, typename T2>
675 inline bool
676 NotAlmostEquals( T1 x1, T2 x2 )
677 {
678 return ! AlmostEquals( x1, x2 );
679 }
680
681
682 /** \brief Return the result of an exact comparison between two scalar values of potetially different types.
683 *
684 * template <typename TInput1, typename TInput2>
685 * inline bool
686 * ExactlyEquals( const TInput & x1, const TInput & x2 )
687 *
688 * template <typename TInput1, typename TInput2>
689 * inline bool
690 * NotExactlyEquals( const TInput & x1, const TInput & x2 )
691 *
692 * These functions complement the EqualsComparison functions and determine if two scalar
693 * values are exactly equal using the compilers casting rules when comparing two different types.
694 * While this is also easily accomplished by using the equality operators,
695 * use of this function demonstrates the intent of an exact equality check and thus improves
696 * readability and clarity of code. In addition, this function suppresses float-equal warnings
697 * produced when using Clang.
698 *
699 * \param x1 first floating point value to compare
700 * \param x2 second floating point value to compare
701 */
702
703 // The ExactlyEquals function
704 template <typename TInput1, typename TInput2>
705 inline bool
706 ExactlyEquals( const TInput1 & x1, const TInput2 & x2 )
707 {
708 CLANG_PRAGMA_PUSH
709 CLANG_SUPPRESS_Wfloat_equal
710 return x1 == x2;
711 CLANG_PRAGMA_POP
712 }
713
714 //The NotExactlyEquals function
715 template <typename TInput1, typename TInput2>
716 inline bool
717 NotExactlyEquals( const TInput1 & x1, const TInput2 & x2 )
718 {
719 return !ExactlyEquals(x1, x2);
720 }
721
722
723 /** Return whether the number is a prime number or not.
724 *
725 * \note Negative numbers cannot be prime.
726 */
727 ITKCommon_EXPORT bool IsPrime( unsigned short n );
728 ITKCommon_EXPORT bool IsPrime( unsigned int n );
729 ITKCommon_EXPORT bool IsPrime( unsigned long n );
730 ITKCommon_EXPORT bool IsPrime( unsigned long long n );
731
732
733 /** Return the greatest factor of the decomposition in prime numbers. */
734 ITKCommon_EXPORT unsigned short GreatestPrimeFactor( unsigned short n );
735 ITKCommon_EXPORT unsigned int GreatestPrimeFactor( unsigned int n );
736 ITKCommon_EXPORT unsigned long GreatestPrimeFactor( unsigned long n );
737 ITKCommon_EXPORT unsigned long long GreatestPrimeFactor( unsigned long long n );
738
739
740 /*==========================================
741 * Alias the vnl_math functions in the itk::Math
742 * namespace. If possible, use the std:: equivalents
743 */
744 using std::isnan;
745 using std::isinf;
746 using std::isfinite;
747 using std::isnormal;
748 using std::cbrt;
749 using std::hypot;
750 using vnl_math::angle_0_to_2pi;
751 using vnl_math::angle_minuspi_to_pi;
752 using vnl_math::rnd_halfinttoeven;
753 using vnl_math::rnd_halfintup;
754 using vnl_math::rnd;
755 using vnl_math::floor;
756 using vnl_math::ceil;
757 using vnl_math::sgn;
758 using vnl_math::sgn0;
759 using vnl_math::remainder_truncated;
760 using vnl_math::remainder_floored;
761 using vnl_math::abs;
762 using vnl_math::sqr;
763 using vnl_math::cube;
764 using vnl_math::squared_magnitude;
765
766 } // end namespace Math
767 } // end namespace itk
768
769 #endif // end of itkMath.h
770