songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Revive the function In_smallest_orthogonal_sphere so that we can compare it with the new function

This commit is contained in:
Andreas Fabri 2016-11-23 13:05:14 +01:00
parent e5e0cdd5f9
commit d1184c91e2
3 changed files with 179 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

31
Skin_surface_3/include/CGAL/Compute_anchor_3.h
View File

@ -23,6 +23,7 @@
#include <CGAL/Regular_triangulation_3.h>
#include <CGAL/Triangulation_simplex_3.h>
#include <CGAL/In_smallest_orthogonal_sphere_3.h>
namespace CGAL {
@ -33,7 +34,9 @@ public:
typedef RegularTriangulation_3 Regular_triangulation;
typedef typename Regular_triangulation::Geom_traits Geom_traits;
typedef typename Geom_traits::Weighted_point_3 Weighted_point;
typedef In_smallest_orthogonal_sphere_3<Geom_traits> In_smallest_orthogonal_sphere_3;
typedef typename RegularTriangulation_3::Vertex_handle Vertex_handle;
typedef typename RegularTriangulation_3::Cell_handle Cell_handle;
typedef typename RegularTriangulation_3::Facet Facet;
@ -144,16 +147,22 @@ private:
// Test whether the anchor of edge (wp1,wp2) and wp2 are equal
Sign test_anchor(Weighted_point &wp1, Weighted_point &wp2) {
In_smallest_orthogonal_sphere_3 iso;
assert(iso(wp1,wp2) == enum_cast<Sign>(reg.geom_traits().power_side_of_bounded_power_sphere_3_object()(wp1, wp2)));
return
enum_cast<Sign>(reg.geom_traits().power_side_of_bounded_power_sphere_3_object()(wp1, wp2));
}
Sign test_anchor(Weighted_point const& wp1, Weighted_point const& wp2,
Weighted_point const& wp3) {
In_smallest_orthogonal_sphere_3 iso;
assert(iso(wp1,wp2,wp3) == enum_cast<Sign>(reg.geom_traits().power_side_of_bounded_power_sphere_3_object()(wp1, wp2,wp3)));
return
enum_cast<Sign>(reg.geom_traits().power_side_of_bounded_power_sphere_3_object()(wp1, wp2, wp3));
}
Sign test_anchor(Weighted_point const& wp1, Weighted_point const& wp2,
Weighted_point const& wp3, Weighted_point const& wp4) {
In_smallest_orthogonal_sphere_3 iso;
assert(iso(wp1,wp2,wp3,wp4) == enum_cast<Sign>(reg.geom_traits().power_side_of_bounded_power_sphere_3_object()(wp1, wp2,wp3,wp4)));
return
enum_cast<Sign>(reg.geom_traits().power_side_of_bounded_power_sphere_3_object()(wp1, wp2, wp3, wp4));
}
@ -194,7 +203,8 @@ private:
default:
CGAL_error();
}
In_smallest_orthogonal_sphere_3 iso;
assert(iso(wp1,wp2,wp3) == enum_cast<Sign>(reg.geom_traits().power_side_of_bounded_power_sphere_3_object()(wp1, wp2,wp3)));
return
enum_cast<Sign>(reg.geom_traits().power_side_of_bounded_power_sphere_3_object()(wp1, wp2, wp3));
}
@ -203,7 +213,12 @@ private:
// Test whether the anchor of ch and anchor of the facet (ch,i) are equal
Sign test_anchor(Cell_handle ch, int i) {
CGAL_assertion(!reg.is_infinite(ch));
Weighted_point wp1 = ch->vertex((i+1)&3)->point(),
wp2 = ch->vertex((i+2)&3)->point(),
wp3 = ch->vertex((i+3)&3)->point(),
wp4 = ch->vertex(i)->point();
In_smallest_orthogonal_sphere_3 iso;
assert(iso(wp1,wp2,wp3,wp4) == enum_cast<Sign>(reg.geom_traits().power_side_of_bounded_power_sphere_3_object()(wp1, wp2,wp3,wp4)));
return enum_cast<Sign>(reg.geom_traits().power_side_of_bounded_power_sphere_3_object()(
ch->vertex((i+1)&3)->point(),
ch->vertex((i+2)&3)->point(),
@ -215,6 +230,14 @@ private:
CGAL_assertion(!reg.is_infinite(ch2));
int index = ch2->index(ch);
Weighted_point wp1 = ch->vertex(0)->point(),
wp2 = ch->vertex(1)->point(),
wp3 = ch->vertex(2)->point(),
wp4 = ch->vertex(3)->point(),
wp5 = ch2->vertex(index)->point();
In_smallest_orthogonal_sphere_3 iso;
assert(iso(wp1,wp2,wp3,wp4,wp5) == enum_cast<Sign>(reg.geom_traits().power_side_of_bounded_power_sphere_3_object()(wp1, wp2,wp3,wp4,wp5)));
return enum_cast<Sign>(
reg.geom_traits().power_side_of_bounded_power_sphere_3_object()(
ch->vertex(0)->point(),
@ -227,6 +250,8 @@ private:
Weighted_point const& wp1, Weighted_point const& wp2,
Weighted_point const& wp3, Weighted_point const& wp4,
Weighted_point const& wp5) {
In_smallest_orthogonal_sphere_3 iso;
assert(iso(wp1,wp2,wp3,wp4,wp5) == enum_cast<Sign>(reg.geom_traits().power_side_of_bounded_power_sphere_3_object()(wp1, wp2,wp3,wp4,wp5)));
return enum_cast<Sign>(
reg.geom_traits().power_side_of_bounded_power_sphere_3_object()(wp1, wp2, wp3, wp4, wp5));
}

149
Skin_surface_3/include/CGAL/In_smallest_orthogonal_sphere_3.h Normal file
View File

@ -0,0 +1,149 @@
#ifndef CGAL_IN_SMALLEST_ORTHOGONAL_SPHEREC3
#define CGAL_IN_SMALLEST_ORTHOGONAL_SPHEREC3
namespace CGAL {
//return the sign of the power test of weighted point (sx,sy,sz,sw)
//with respect to the smallest sphere orthogonal to
//p,q,r
template< class FT >
CGAL_MEDIUM_INLINE
Sign
in_smallest_orthogonal_sphereC3(
const FT &px, const FT &py, const FT &pz, const FT &pw,
const FT &qx, const FT &qy, const FT &qz, const FT &qw,
const FT &rx, const FT &ry, const FT &rz, const FT &rw,
const FT &sx, const FT &sy, const FT &sz, const FT &sw)
{
// Translate p to origin and compute determinants
FT qpx = qx-px;
FT qpy = qy-py;
FT qpz = qz-pz;
FT rpx = rx-px;
FT rpy = ry-py;
FT rpz = rz-pz;
FT qq = CGAL_NTS square(qpx) + CGAL_NTS square(qpy) + CGAL_NTS square(qpz);
FT rr = CGAL_NTS square(rpx) + CGAL_NTS square(rpy) + CGAL_NTS square(rpz);
FT qr = qpx*rpx + qpy*rpy + qpz*rpz;
FT qpw = qq - qw + pw ;
FT rpw = rr - rw + pw ;
FT den = determinant(qq,qr,
qr,rr);
FT detq = determinant(qpw,qr,
rpw,rr);
FT detr = determinant(qq,qpw,
qr,rpw);
// Smallest smallest orthogonal sphere center
// c = detq/2*den q + detr/2*den r (origin at p)
// square radius c^2 - pw
FT spx = sx-px;
FT spy = sy-py;
FT spz = sz-pz;
FT ss = CGAL_NTS square(spx) + CGAL_NTS square(spy) + CGAL_NTS square(spz);
FT sq = spx*qpx + spy*qpy + spz*qpz;
FT sr = spx*rpx + spy*rpy + spz*rpz;
CGAL_triangulation_assertion( ! CGAL_NTS is_zero(den) );
// return sign of (c- s)^2 - (c^2 - pw) - sw note that den >= 0 -
return CGAL_NTS sign( den*(ss - sw + pw)- detq*sq - detr*sr);
}
// return the sign of the power test of weighted point (rx,ry,rz,rw)
// with respect to the smallest sphere orthogoanal to
// p,q
template< class FT >
Sign
in_smallest_orthogonal_sphereC3(
const FT &px, const FT &py, const FT &pz, const FT &pw,
const FT &qx, const FT &qy, const FT &qz, const FT &qw,
const FT &rx, const FT &ry, const FT &rz, const FT &rw)
{
FT FT2(2);
FT FT4(4);
FT dpx = px - qx;
FT dpy = py - qy;
FT dpz = pz - qz;
FT dpw = pw - qw;
FT dp2 = CGAL_NTS square(dpx) + CGAL_NTS square(dpy) + CGAL_NTS square(dpz);
FT drx = rx - (px + qx)/FT2;
FT dry = ry - (py + qy)/FT2;
FT drz = rz - (pz + qz)/FT2;
FT drw = rw - (pw + qw)/FT2;
FT dr2 = CGAL_NTS square(drx) + CGAL_NTS square(dry) + CGAL_NTS square(drz);
FT dpr = dpx*drx + dpy*dry +dpz*drz;
return CGAL_NTS sign (dr2 - dp2/FT4 + dpr*dpw/dp2 - drw );
}
template< typename K >
class In_smallest_orthogonal_sphere_3
{
public:
typedef typename K::Weighted_point_3 Weighted_point_3;
typedef typename K::Sign Sign;
typedef Sign result_type;
Sign operator() ( const Weighted_point_3 & p,
const Weighted_point_3 & q,
const Weighted_point_3 & r,
const Weighted_point_3 & s,
const Weighted_point_3 & t) const
{
K traits;
typename K::Power_side_of_oriented_power_sphere_3 power_test_3 = traits.power_side_of_oriented_power_sphere_3_object();
typename K::Orientation_3 orientation = traits.orientation_3_object();
typename K::Orientation o = orientation(p,q,r,s);
typename K::Oriented_side os = power_test_3(p,q,r,s,t);
CGAL_triangulation_assertion( o != COPLANAR);
// the minus sign below is due to the fact that power_test_3
// return in fact minus the 5x5 determinant of lifted (p,q,r,s.t)
return - o * os;
}
Sign operator() ( const Weighted_point_3 & p,
const Weighted_point_3 & q,
const Weighted_point_3 & r,
const Weighted_point_3 & s) const
{
return in_smallest_orthogonal_sphereC3(
p.x(), p.y(), p.z(), p.weight(),
q.x(), q.y(), q.z(), q.weight(),
r.x(), r.y(), r.z(), r.weight(),
s.x(), s.y(), s.z(), s.weight());
}
Sign operator() ( const Weighted_point_3 & p,
const Weighted_point_3 & q,
const Weighted_point_3 & r) const
{
return in_smallest_orthogonal_sphereC3(
p.x(), p.y(), p.z(), p.weight(),
q.x(), q.y(), q.z(), q.weight(),
r.x(), r.y(), r.z(), r.weight());
}
Sign operator() ( const Weighted_point_3 & p,
const Weighted_point_3 & q) const
{
return CGAL_NTS sign( CGAL_NTS square(p.x()-q.x()) +
CGAL_NTS square(p.y()-q.y()) +
CGAL_NTS square(p.z()-q.z()) +
p.weight() - q.weight());
}
};
}
#endif

3
Skin_surface_3/test/Skin_surface_3/anchor_test.cpp
View File

@ -6,7 +6,8 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <cassert>
#include <CGAL/Regular_triangulation_3.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
// typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Simple_cartesian<CGAL::Gmpq> K;
typedef CGAL::Regular_triangulation_3<K> Regular;
#include <CGAL/Compute_anchor_3.h>