include/geometry/euclidean.hpp
 author Dmitriy Morozov Mon, 05 Jul 2010 23:51:26 -0700 branch dev changeset 211 347b3461965a parent 19 efa14432761a permissions -rw-r--r--
Bug in the tutorial: simplices = [] -> simplices = Filtration()
```
/* --- Point --- */
template<class NumberType_>
Kernel<NumberType_>::
Kernel(DimensionType dimension): dimension_(dimension), origin_(dimension)
{
for (unsigned int i = 0; i < dimension; ++i)
origin_(i) = NumberType(0);
}

template<class NumberType_>
Kernel<NumberType_>::Point::
Point(const Point& p, const NumberType& pp): Parent(p.size() + 1)
{
using boost::numeric::ublas::vector_range;
using boost::numeric::ublas::range;

vector_range<VectorType> vr(*this, range(0, size() - 1));
vr = p;
(*this)(size() - 1) = pp;
}

template<class NumberType_>
typename Kernel<NumberType_>::NumberType
Kernel<NumberType_>::Point::
squared_distance(const Point& p) const
{
return boost::numeric::ublas::inner_prod(*this - p, *this - p);
}

template<class NumberType_>
std::istream&
operator>>(std::istream& in, typename Kernel<NumberType_>::Point& p)
{
for (unsigned int i = 0; i < p.size(); ++i)
in >> p[i];
return in;
}

template<class NumberType_>
std::ostream&
operator<<(std::ostream& out, typename Kernel<NumberType_>::Point& p)
{
for (unsigned int i = 0; p.size(); ++i)
out << p[i];
return out;
}

/* --- Sphere --- */
template<class NumberType_>
typename Kernel<NumberType_>::NumberType
Kernel<NumberType_>::Sphere::
side_of(const Point& p) const
{ return p.squared_distance(center_) - squared_radius(); }

/* --- Kernel --- */
template<class NumberType_>
typename Kernel<NumberType_>::Sphere
Kernel<NumberType_>::
circumsphere(const PointContainer& points) const
{
if (points.size() == 0) return Sphere(origin(), NumberType(0));		// FIXME: should this be an assertion instead?

using boost::numeric::ublas::inner_prod;
using boost::numeric::ublas::vector_range;
using boost::numeric::ublas::range;

std::vector<Point> basis;
basis.reserve(points.size() - 1);
for (unsigned int i = 1; i < points.size(); ++i)
{
Point pt = *points[i] - *points[0];
basis.push_back(Point(pt, inner_prod(pt, pt)));
for (unsigned int j = 0; j < i - 1; ++j)
{
basis[i-1] -= basis[j] * inner_prod(basis[i-1], basis[j]) / inner_prod(basis[j], basis[j]);
normalize(basis[i-1]);
}
}
Point clifted(origin(), NumberType(1));
for (unsigned int j = 0; j < basis.size(); ++j)
{
clifted -= basis[j] * inner_prod(clifted, basis[j]) / inner_prod(basis[j], basis[j]);
normalize(clifted);
}

Point center(vector_range<VectorType>(clifted, range(0, dimension())));
center /= NumberType_(-2) * clifted(dimension());
center += *points[0];

normalize(center);

}

template<class NumberType_>
typename Kernel<NumberType_>::NumberType
Kernel<NumberType_>::
{
}

template<class NumberType_>
typename Kernel<NumberType_>::Point
Kernel<NumberType_>::
circumcenter(const PointContainer& points) const
{
return circumsphere(points).center();
}

template<class NumberType_>
typename Kernel<NumberType_>::NumberType
Kernel<NumberType_>::
side_of_circumsphere(const PointContainer& points, const Point& p) const
{
Sphere s = circumsphere(points);
return s.side_of(p);
}

/* --- Kernel Private --- */
template<class NumberType_>
typename Kernel<NumberType_>::NumberType&
Kernel<NumberType_>::
normalize(NumberType& n) const
{
return number_traits<typename Point::NumberType>::normalize(n);
}

template<class NumberType_>
typename Kernel<NumberType_>::Point&
Kernel<NumberType_>::
normalize(Point& p) const
{
for (unsigned int i = 0; i < p.size(); ++i)
number_traits<typename Point::NumberType>::normalize(p(i));
return p;
}
```