SphereD.h

Go to the documentation of this file.

/**
 *  \file IMP/algebra/SphereD.h   \brief Simple 3D sphere class.
 *
 *  Copyright 2007-2022 IMP Inventors. All rights reserved.
 *
 */

#ifndef IMPALGEBRA_SPHERE_D_H
#define IMPALGEBRA_SPHERE_D_H

#include "algebra_macros.h"
#include "constants.h"
#include "BoundingBoxD.h"
#include "VectorD.h"
#include "utility.h"
#include "GeometricPrimitiveD.h"
#include <cereal/access.hpp>
#include <cmath>

IMPALGEBRA_BEGIN_NAMESPACE

//! Represent a sphere in D-dimensions.
/** \geometry
  */
template <int D>
class SphereD : public GeometricPrimitiveD<D> {
  friend class cereal::access;

  template<class Archive>
  void serialize(Archive &ar) {
    ar(center_, radius_);
  }

 public:
  SphereD() {
#if IMP_HAS_CHECKS >= IMP_USAGE
    radius_ = std::numeric_limits<double>::quiet_NaN();
#endif
  }
  SphereD(const VectorD<D> &center, double radius)
      : center_(center), radius_(radius) {
    IMP_USAGE_CHECK(radius >= 0, "Radius can't be negative");
  }
  double get_radius() const {
    IMP_INTERNAL_CHECK(!IMP::isnan(radius_),
                       "Attempt to use uninitialized sphere.");
    return radius_;
  }
  const VectorD<D> &get_center() const { return center_; }
  //! Return true if this sphere contains the other one
  bool get_contains(const SphereD<D> &o) const {
    double d = (get_center() - o.get_center()).get_magnitude();
    return (d + o.get_radius() < get_radius());
  }

  //! Return true if the point is in or on the surface of the sphere
  bool get_contains(const VectorD<D> &p) const {
    return ((p - center_).get_squared_magnitude() <= get_squared(radius_));
  }
  IMP_SHOWABLE_INLINE(SphereD,
  { out << "(" << spaces_io(center_) << ": " << get_radius() << ")"; });
#ifdef IMP_SWIG_WRAPPER
  static void _get_struct_size(size_t &sz, size_t &center_offset,
                               size_t &radius_offset) {
    sz = sizeof(SphereD<D>);
    center_offset = offsetof(SphereD<D>, center_);
    radius_offset = offsetof(SphereD<D>, radius_);
  }
#endif
#ifndef IMP_DOXYGEN
#ifndef SWIG
  VectorD<D> &_access_center() { return center_; }
  void _set_radius(double d) { radius_ = d; }
  void _set_center(const VectorD<D> &center){ center_ = center; }
  double &operator[](unsigned int i) {
    IMP_USAGE_CHECK(i < D + 1, "Out of range");
    if (i < D) {
      return center_[i];
    } else {
      return radius_;
    }
  }
  double operator[](unsigned int i) const {
    IMP_USAGE_CHECK(i < D + 1, "Out of range");
    if (i < D) {
      return center_[i];
    } else {
      return radius_;
    }
  }
#endif
#endif
  unsigned int get_dimension() const { return center_.get_dimension(); }

 private:
  VectorD<D> center_;
  double radius_;
};

IMP_VOLUME_GEOMETRY_METHODS_D(Sphere, sphere,
{ return PI * 4.0 * get_squared(g.get_radius()); },
{ return PI * (4.0 / 3.0) * std::pow(g.get_radius(), 3.0); },
                              return BoundingBoxD<D>(g.get_center()) +
                                     g.get_radius(););

template <unsigned int D>
inline SphereD<D> get_unit_sphere_d() {
  return SphereD<D>(get_zero_vector_d<D>(), 1.0);
}

inline SphereD<-1> get_unit_sphere_kd(unsigned int d) {
  return SphereD<-1>(get_zero_vector_kd(d), 1.0);
}

//! Return the distance between the two spheres if they are disjoint
/** If they intersect, the distances are not meaningful.
    \see SphereD
*/
template <int D>
inline double get_distance(const SphereD<D> &a, const SphereD<D> &b) {
  double d = (a.get_center() - b.get_center()).get_magnitude();
  return d - a.get_radius() - b.get_radius();
}

//! Return the power distance between the two spheres
/** The power distance is the square of the distance between the centers
    minus the sum of the square of the radii.
    \see SphereD
*/
template <int D>
inline double get_power_distance(const SphereD<D> &a, const SphereD<D> &b) {
  double d = (a.get_center() - b.get_center()).get_squared_magnitude();
  return d - square(a.get_radius()) - square(b.get_radius());
}

//! Return true if the two balls bounded by the two spheres intersect
/** \see SphereD
 */
template <int D>
inline bool get_interiors_intersect(const SphereD<D> &a, const SphereD<D> &b) {
  double sr = a.get_radius() + b.get_radius();
  for (unsigned int i = 0; i < a.get_dimension(); ++i) {
    double delta = std::abs(a.get_center()[i] - b.get_center()[i]);
    if (delta >= sr) return false;
  }
  return get_squared_distance(a.get_center(), b.get_center()) < get_squared(sr);
}

#if !defined(SWIG) && !defined(IMP_DOXYGEN)

namespace internal {
template <int D>
struct SphereSpacesIO {
  const SphereD<D> &v_;
  SphereSpacesIO(const SphereD<D> &v) : v_(v) {}
};
template <int D>
inline std::ostream &operator<<(std::ostream &out, const SphereSpacesIO<D> &s) {
  for (unsigned int i = 0; i < s.v_.get_center().get_dimension(); ++i) {
    out << s.v_.get_center()[i] << " ";
  }
  out << s.v_.get_radius();
  return out;
}
}

//! Use this before outputting to a stream with spaces delimiting
/** std::cout << spaces_io(s);
    produces "1.0 2.0 3.0 4.0" where the radius is 4.0
    \see SphereD
 */
template <int D>
inline internal::SphereSpacesIO<D> spaces_io(const SphereD<D> &v) {
  return internal::SphereSpacesIO<D>(v);
}
#endif

#ifdef IMP_DOXYGEN
//! Compute the bounding box of any geometric object
template <class Geometry>
BoundingBoxD<3> get_bounding_box(const Geometry &);
//! Compute the surface area of any volumetric object
template <class Geometry>
double get_surface_area(const Geometry &);
//! Compute the volume of any volumetric object
template <class Geometry>
double get_volume(const Geometry &);
//! Compute the area of any surface object
template <class Geometry>
double get_area(const Geometry &);

#endif
template <int D>
VectorD<D> get_vector_geometry(const SphereD<D> &s) {
  return s.get_center();
}

IMPALGEBRA_END_NAMESPACE

#endif /* IMPALGEBRA_SPHERE_D_H */