vtkMathUtilities.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
// SPDX-License-Identifier: BSD-3-Clause

#ifndef vtkMathUtilities_h
#define vtkMathUtilities_h

#include "vtkABINamespace.h"
#include "vtkAssume.h"

#include <algorithm>
#include <cmath>
#include <limits>
#include <tuple>
#include <type_traits>

namespace vtkMathUtilities
{
VTK_ABI_NAMESPACE_BEGIN

template <class A>
bool FuzzyCompare(A a, A b, A epsilon = std::numeric_limits<A>::epsilon())
{
  return fabs(a - b) < epsilon;
}

template <class A>
A SafeDivision(A a, A b)
{
  // Avoid overflow
  if ((b < static_cast<A>(1)) && (a > b * std::numeric_limits<A>::max()))
  {
    return std::numeric_limits<A>::max();
  }

  // Avoid underflow
  if ((a == static_cast<A>(0)) ||
    ((b > static_cast<A>(1)) && (a < b * std::numeric_limits<A>::min())))
  {
    return static_cast<A>(0);
  }

  // safe to do the division
  return (a / b);
}

template <class A>
bool NearlyEqual(A a, A b, A tol = std::numeric_limits<A>::epsilon())
{
  A absdiff = fabs(a - b);
  A d1 = vtkMathUtilities::SafeDivision<A>(absdiff, fabs(a));
  A d2 = vtkMathUtilities::SafeDivision<A>(absdiff, fabs(b));

  return ((d1 <= tol) || (d2 <= tol));
}

template <class A>
inline void UpdateRange(A& min, A& max, const A& value)
{
  if constexpr (std::is_floating_point_v<A>)
  {
    if (VTK_UNLIKELY(std::isnan(value)))
    {
      return;
    }
  }
  if (value < min)
  {
    min = value;
    max = std::max(max, value);
  }
  else if (value > max)
  {
    min = std::min(min, value);
    max = value;
  }
}

template <class A>
inline void UpdateRangeFinite(A& min, A& max, const A& value)
{
  if constexpr (std::is_floating_point_v<A>)
  {
    if (VTK_UNLIKELY(std::isinf(value) || std::isnan(value)))
    {
      return;
    }
  }
  if (value < min)
  {
    min = value;
    max = std::max(max, value);
  }
  else if (value > max)
  {
    min = std::min(min, value);
    max = value;
  }
}

template <class A>
inline A SafeCastFromDouble(double value)
{
  if constexpr (std::is_floating_point_v<A>)
  {
    return static_cast<A>(value);
  }
  else
  {
    constexpr double typeMin = static_cast<double>(std::numeric_limits<A>::lowest());
    constexpr double typeMax = static_cast<double>(std::numeric_limits<A>::max());
    // The comparisons are done in a specific way to support 64-bit ints.
    // Explanation: casting the int64 max (equal to 2^63-1) to double gives
    // 2^63 because double cannot exactly represent 2^63-1. So the value must
    // be strictly less than typeMax to avoid overflow when casting to int64.
    if (VTK_LIKELY(value >= typeMin))
    {
      if (VTK_LIKELY(value < typeMax))
      {
        return static_cast<A>(value);
      }
      return std::numeric_limits<A>::max();
    }
    if (VTK_UNLIKELY(std::isnan(value)))
    {
      return 0;
    }
    return std::numeric_limits<A>::lowest();
  }
}

VTK_ABI_NAMESPACE_END
} // End vtkMathUtilities namespace.

#endif // vtkMathUtilities_h
// VTK-HeaderTest-Exclude: vtkMathUtilities.h