modulo.h

Go to the documentation of this file.

//
//
#ifndef __stir_modulo_H__
#define __stir_modulo_H__

/*
    Copyright (C) 2004- 2009, Hammersmith Imanet Ltd
    This file is part of STIR.
 
    SPDX-License-Identifier: Apache-2.0
 
    See STIR/LICENSE.txt for details
*/
 
#include "stir/common.h"
 
START_NAMESPACE_STIR



inline double
modulo(const double a, const double b)
{
  const double res = fmod(a, b);
  return res < 0 ? res + fabs(b) : res;
}


inline float
modulo(const float a, const float b)
{
  float res = static_cast<float>(modulo(static_cast<double>(a), static_cast<double>(b)));
  assert(res >= 0);
  const float abs_b = b >= 0 ? b : -b;
  if (res >= abs_b)
    res -= abs_b;
  assert(res >= 0);
  return res;
}


inline int
modulo(const int a, const int b)
{
  const int res = a % b;
  const int res2 = res < 0 ? res + (b >= 0 ? b : -b) : res;
  assert(res2 >= 0);
  assert(res2 < (b >= 0 ? b : -b));
  return res2;
}


template <int num_dimensions, typename T>
inline BasicCoordinate<num_dimensions, T>
modulo(const BasicCoordinate<num_dimensions, T>& a, const BasicCoordinate<num_dimensions, T>& b)
{
  BasicCoordinate<num_dimensions, T> result;
  for (int d = 1; d <= num_dimensions; ++d)
    result[d] = modulo(a[d], b[d]);
  return result;
}


template <typename FloatOrDouble>
inline FloatOrDouble
from_min_pi_plus_pi_to_0_2pi(const FloatOrDouble phi)
{
  static const FloatOrDouble two_pi = static_cast<FloatOrDouble>(2 * _PI);
  assert(phi >= -two_pi);
  assert(phi < two_pi);
  if (phi >= 0)
    return phi;
  FloatOrDouble res = phi + two_pi;
  // due to floating point finite precision, the following check is needed...
  if (res >= two_pi)
    res -= two_pi;
  assert(res >= 0);
  assert(res < two_pi);
  return res;
}


template <typename FloatOrDouble>
inline FloatOrDouble
to_0_2pi(const FloatOrDouble phi)
{
  return modulo(phi, static_cast<FloatOrDouble>(2 * _PI));
}

 
END_NAMESPACE_STIR
 
#endif