code polishing; added constructive version of the predicate; added versions for different algebraic structures

This commit is contained in:
Menelaos Karavelas 2007-07-18 21:46:47 +00:00
parent c6fc1d29a8
commit ae7a08424b
1 changed files with 216 additions and 14 deletions

View File

@ -57,6 +57,37 @@ public:
s3.point());
}
inline
Sign sqrt_ext_sign(const FT& A, const FT& B,
const FT& Exp, const FT& Eyp, const FT& Erp,
const FT& Exy2, const FT& Exr, const FT& Eyr,
const FT dx, const FT& dy,
const Field_with_sqrt_tag&) const
{
FT G = CGAL::square(Exp) + CGAL::square(Eyp) - CGAL::square(Erp);
return CGAL::sign(A + B * CGAL::sqrt(G));
}
inline
Sign sqrt_ext_sign(const FT& A, const FT& B,
const FT& Exp, const FT& Eyp, const FT& Erp,
const FT& Exy2, const FT& Exr, const FT& Eyr,
const FT dx, const FT& dy,
const Integral_domain_without_division_tag&) const
{
Sign sA = CGAL::sign(A);
Sign sB = CGAL::sign(B);
if ( sA == CGAL::ZERO ) { return sB; }
if ( sB == CGAL::ZERO ) { return sA; }
if ( sA == sB ) { return sA; }
Sign s = CGAL::sign(CGAL::square(Exy2 * Exr - Erp * dy)
+ CGAL::square(Exy2 * Eyr + Erp * dx)
- CGAL::square(B));
return sA * s;
}
Orientation predicate(const Site_2& s1, const Site_2& s2,
const Site_2& s3, const Site_2& p1,
const Site_2& p2) const
@ -97,17 +128,147 @@ public:
FT Exy2 = 2 * det2x2_by_formula(xl, yl, xm, ym);
#if 0
FT A = (Exp * Exr + Eyp * Eyr) * Exy2 + (Eyp * dx - Exp * dy) * Erp;
FT B = Exy * Exy2 - Exp * dx - Eyp * dy;
FT C = CGAL::square(Exp) + CGAL::square(Eyp) - CGAL::square(Erp);
return sign_a_plus_b_x_sqrt_c(A, B, C);
#else
Sign sA = CGAL::sign((Exp * Exr + Eyp * Eyr) * Exy2
+ (Eyp * dx - Exp * dy) * Erp);
return sqrt_ext_sign(A, B, Exp, Eyp, Erp,
Exy2, Exr, Eyr, dx, dy, Method_tag());
}
FT B = Exy * Exy2 - Exp * dx - Eyp * dy;
Orientation operator()(const Site_2& s1, const Site_2& s2,
const Site_2& s3, const Site_2& p1,
const Site_2& p2) const
{
Orientation o = predicate(s1, s2, s3, p1, p2);
#ifndef NDEBUG
Orientation o_old = Base::operator()(s1, s2, s3, p1, p2);
CGAL_assertion( o == o_old );
#endif
return o;
}
};
//--------------------------------------------------------------------
template<class K, class MTag>
class Constructive_orientation8_C2
{
public:
typedef K Kernel;
typedef MTag Method_tag;
typedef typename K::Site_2 Site_2;
typedef typename K::Point_2 Point_2;
typedef typename K::Orientation Orientation;
typedef typename K::FT FT;
typedef Orientation result_type;
typedef Arity_tag<3> Arity;
typedef Site_2 argument_type;
private:
FT s1x, s1y;
FT xj, xk, yj, yk, rj, rk, pj, pk, Exp, Eyp, Erp, Exy, Exr, Eyr, A1;
FT xl_, yl_, ca1, ca2, ca3, cb1, cb2, cb3;
public:
Constructive_orientation8_C2(const Site_2& s1, const Site_2& s2,
const Site_2& s3)
{
s1x = s1.x();
s1y = s1.y();
xj = s2.x() - s1.x();
xk = s3.x() - s1.x();
yj = s2.y() - s1.y();
yk = s3.y() - s1.y();
rj = s2.weight() - s1.weight();
rk = s3.weight() - s1.weight();
pj = CGAL::square(xj) + CGAL::square(yj) - CGAL::square(rj);
pk = CGAL::square(xk) + CGAL::square(yk) - CGAL::square(rk);
Exp = det2x2_by_formula(xj, pj, xk, pk);
Eyp = det2x2_by_formula(yj, pj, yk, pk);
Erp = det2x2_by_formula(rj, pj, rk, pk);
Exy = det2x2_by_formula(xj, yj, xk, yk);
Exr = det2x2_by_formula(xj, rj, xk, rk);
Eyr = det2x2_by_formula(yj, rj, yk, rk);
A1 = Exp * Exr + Eyp * Eyr;
}
Constructive_orientation8_C2(const Site_2& s1, const Site_2& s2,
const Site_2& s3, const Site_2& p1)
{
s1x = s1.x();
s1y = s1.y();
xj = s2.x() - s1.x();
xk = s3.x() - s1.x();
yj = s2.y() - s1.y();
yk = s3.y() - s1.y();
rj = s2.weight() - s1.weight();
rk = s3.weight() - s1.weight();
pj = CGAL::square(xj) + CGAL::square(yj) - CGAL::square(rj);
pk = CGAL::square(xk) + CGAL::square(yk) - CGAL::square(rk);
Exp = det2x2_by_formula(xj, pj, xk, pk);
Eyp = det2x2_by_formula(yj, pj, yk, pk);
Erp = det2x2_by_formula(rj, pj, rk, pk);
Exy = det2x2_by_formula(xj, yj, xk, yk);
Exr = det2x2_by_formula(xj, rj, xk, rk);
Eyr = det2x2_by_formula(yj, rj, yk, rk);
A1 = Exp * Exr + Eyp * Eyr;
xl_ = p1.x() - s1.x();
yl_ = p1.y() - s1.y();
ca1 = 2 * A1 * xl_ - Eyp * Erp;
ca2 = -2 * A1 * yl_ + Exp * Erp;
ca3 = (Eyp * xl_ - Exp * yl_) * Erp;
cb1 = Exp + Exy * xl_;
cb2 = Eyp - Exy * yl_;
cb3 = -Exp * xl_ - Eyp * yl_;
}
public:
inline
Orientation operator()(const Site_2& s1, const Site_2& s2,
const Site_2& s3) const
{
return Kernel().orientation_2_object()(s1.point(), s2.point(),
s3.point());
}
inline
Sign sqrt_ext_sign(const FT& A, const FT& B, const FT& Exy2,
const FT dx, const FT& dy,
const Field_with_sqrt_tag&) const
{
FT G = CGAL::square(Exp) + CGAL::square(Eyp) - CGAL::square(Erp);
return CGAL::sign(A + B * CGAL::sqrt(G));
}
inline
Sign sqrt_ext_sign(const FT& A, const FT& B, const FT& Exy2,
const FT dx, const FT& dy,
const Integral_domain_without_division_tag&) const
{
Sign sA = CGAL::sign(A);
Sign sB = CGAL::sign(B);
if ( sA == CGAL::ZERO ) { return sB; }
@ -118,17 +279,57 @@ public:
+ CGAL::square(Exy2 * Eyr + Erp * dx)
- CGAL::square(B));
return sA * s;
#endif
}
Orientation operator()(const Site_2& s1, const Site_2& s2,
const Site_2& s3, const Site_2& p1,
const Site_2& p2) const
Orientation predicate(const Site_2& p1, const Site_2& p2) const
{
Orientation o = predicate(s1, s2, s3, p1, p2);
Orientation o_old = Base::operator()(s1, s2, s3, p1, p2);
// computes the orientation of the Voronoi vertex of s1, s2, s3 and
// the points p1 and p2
FT xl = p1.x() - s1x;
FT xm = p2.x() - s1x;
CGAL_assertion( o == o_old );
FT yl = p1.y() - s1y;
FT ym = p2.y() - s1y;
FT dx = xl - xm;
FT dy = yl - ym;
FT Exy2 = 2 * det2x2_by_formula(xl, yl, xm, ym);
FT A = A1 * Exy2 + (Eyp * dx - Exp * dy) * Erp;
FT B = Exy * Exy2 - Exp * dx - Eyp * dy;
return sqrt_ext_sign(A, B, Exy2, dx, dy, Method_tag());
}
Orientation predicate(const Site_2& p2) const
{
// computes the orientation of the Voronoi vertex of s1, s2, s3 and
// the points p1 and p2
FT xm = p2.x() - s1x;
FT ym = p2.y() - s1y;
FT dx = xl_ - xm;
FT dy = yl_ - ym;
FT Exy2 = 2 * det2x2_by_formula(xl_, yl_, xm, ym);
FT A = ca1 * xm + ca2 * ym + ca3;
FT B = cb1 * xm + cb2 * ym + cb3;
return sqrt_ext_sign(A, B, Exy2, dx, dy, Method_tag());
}
Orientation operator()(const Site_2& p1, const Site_2& p2) const
{
Orientation o = predicate(p1, p2);
return o;
}
Orientation operator()(const Site_2& p2) const
{
Orientation o = predicate(p2);
return o;
}
@ -136,6 +337,7 @@ public:
//--------------------------------------------------------------------
CGAL_APOLLONIUS_GRAPH_2_END_NAMESPACE
CGAL_END_NAMESPACE