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

simplified and more understandable predicate also computes the right answer

This commit is contained in:
Daniel Russel 2006-09-19 13:24:12 +00:00
parent cdcc59cc35
commit 2f3dfd51a6
18 changed files with 535 additions and 490 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/data/intersect.spheres
View File

@ -1,2 +1,2 @@
0 0 0 1.2 1
.2 1 1 1 1
0 0 0 1 1
0 1 1 1 1

4
Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/interactive_ray_shooting.cpp
View File

@ -1,7 +1,7 @@
#define CGAL_CHECK_EXPENSIVE
#define CGAL_CHECK_EXACTNESS
#include <CGAL/IO/Qt_debug_viewer_2.h>
#include <CGAL/IO/Qt_multithreaded_examiner_viewer_2.h>
#include <CGAL/Arrangement_of_spheres_3/Slice_arrangement.h>
#include <CGAL/Arrangement_of_spheres_traits_3.h>
#include <CGAL/Arrangement_of_spheres_3/Slice.h>
@ -157,7 +157,7 @@ int main(int argc, char *argv[]){
Do_work dw(atof(argv[1]), argv[2]);
Qt_debug_viewer_2<Do_work> qtd(dw, argc, argv);
Qt_multithreaded_examiner_viewer_2<Do_work> qtd(dw, argc, argv);
return qtd();
}

4
Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/makefile
View File

@ -37,7 +37,9 @@ LDFLAGS = \
# target entries
#---------------------------------------------------------------------#
all: mylibs interactive_point_insertion slice_kds show_slice show_circles interactive_ray_shooting interactive_point_location
all: mylibs interactive_point_insertion slice_kds show_slice show_circles interactive_ray_shooting
#interactive_point_location
#interactive_ray_shooting interactive_point_location random_ray_shoot random_point_locate point_locate

8
Arrangement_of_spheres_3/demo/Arrangement_of_spheres_3/show_circles.cpp
View File

@ -1,7 +1,7 @@
#define CGAL_CHECK_EXPENSIVE
#define CGAL_CHECK_EXACTNESS
//#define CGAL_CHECK_EXPENSIVE
//#define CGAL_CHECK_EXACTNESS
#include <CGAL/IO/Qt_debug_viewer_2.h>
#include <CGAL/IO/Qt_examiner_viewer_2.h>
#include <CGAL/Arrangement_of_spheres_3/Slice_arrangement.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Arrangement_of_spheres_3/coordinates.h>
@ -82,7 +82,7 @@ int main(int argc, char *argv[]){
}
app.setMainWidget( qtd);
app.setMainWidget( qtd->window());
qtd->show_everything();
qtd->show();
return app.exec();

27
Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Slice.h
View File

@ -206,6 +206,10 @@ struct Slice {
void audit_edge_collapse_event(Halfedge_handle h, Event_key ek);
void process_equal_centers_two_sphere_event(T::Key k, T::Key l, bool f);
void process_aligned_centers_two_sphere_event(T::Key k, T::Key l, bool f,
CGAL::Comparison_result c[2]);
// Functions to update certificates ----------------------------------
void clean_edge(Halfedge_handle h);
@ -342,10 +346,10 @@ struct Slice {
*/
Face_handle locate_point(const T::Sphere_point_3 &ep, T::Key ind);
Face_handle locate_point(const T::Sphere_point_3 &ep);
template <class It>
Face_handle locate_point(It b, It e, const T::Sphere_point_3 &ep, T::Key ind);
Face_handle locate_point(It b, It e, const T::Sphere_point_3 &ep);
Halfedge_handle shoot_rule(const T::Sphere_point_3& source,
Face_handle f,
@ -364,29 +368,24 @@ struct Slice {
//enum Location {L_BIT=1, R_BIT=2, T_BIT=4, B_BIT=8, IN_BIT=16, };
int sphere_location(const T::Sphere_point_3 &sp,
T::Key locate_point,
T::Key sphere) const ;
bool behind_arc(const T::Sphere_point_3 &ep, T::Key ind,
bool behind_arc(const T::Sphere_point_3 &ep,
Sds::Curve arc, int location) const;
void point_sphere_orientation(const T::Sphere_point_3 &time,
T::Key point,
T::Key sphere,
std::vector<int> &locations) const;
bool locate_point_check_face(const T::Sphere_point_3 &z,
Face_const_handle it,
T::Key ind,
Face_const_handle it,
std::vector<int> &locations) const ;
bool locate_point_check_face_arcs(const T::Sphere_point_3 &z,
T::Key ind,
Face_const_handle f,
std::vector<int> &locations) const ;
bool locate_point_check_face_vertices(const T::Sphere_point_3 &ep,
T::Key index,
Face_const_handle it) const;
CGAL::Comparison_result debug_rule_shoot_answer(const T::Sphere_point_3 &z,
@ -413,9 +412,7 @@ struct Slice {
CGAL::Comparison_result rule_shoot_compare_SS(const T::Sphere_point_3 &z,
Sds::Curve srule,
Sds::Curve arc0,
Sds::Point pt,
Sds::Curve arc1) const ;
Sds::Point pt) const ;
// return comparison of separator to intersection point on the C coordinate
@ -430,10 +427,10 @@ struct Slice {
Sds::Curve rule,
Sds::Point pt,
CGAL::Comparison_result &ret) const ;
bool rule_shoot_compare_if_rational_arc(const T::Sphere_point_3 &z,
/*bool rule_shoot_compare_if_rational_arc(const T::Sphere_point_3 &z,
Sds::Curve rule,
Sds::Curve a,
CGAL::Comparison_result &ret) const;
CGAL::Comparison_result &ret) const;*/
@ -441,7 +438,7 @@ struct Slice {
// Debug Functions ---------------------------------------------------
Face_handle locate_point(const T::Sphere_point_3 &ep);
//Face_handle locate_point(const T::Sphere_point_3 &ep);
Halfedge_handle shoot_rule(const T::Sphere_point_3 &source,
Face_handle f,

9
Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Slice_data_structure.h
View File

@ -196,7 +196,9 @@ public:
void exchange_vertices(Halfedge_handle h, Halfedge_handle p);
Face_handle intersect(Halfedge_handle ha, Halfedge_handle hb);
std::pair<Halfedge_handle,Halfedge_handle>
intersect(Halfedge_handle ha, Halfedge_handle hb,
Halfedge_handle hc, Halfedge_handle hd);
std::pair<Halfedge_handle,Halfedge_handle> unintersect(Face_handle f);
@ -252,6 +254,11 @@ public:
// a halfedge on the rule pointing to the extremal vertex
Halfedge_handle rule_halfedge(Curve::Key k, int i) const;
// a halfedge on the circle pointing to the extremal vertex
Halfedge_handle extremum_halfedge(Curve::Key k, int i) const;
void set_extremum_halfedge(Curve::Key k, int i, Halfedge_handle h);
template <class Out>
void find_halfedges(Curve c, Out o) const {

2
Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_3/Sphere_line_intersection.h
View File

@ -177,7 +177,7 @@ public:
}
CGAL::Comparison_result compare_on_line(const Point_3 &pt) const {
CGAL_precondition(line().has_on(pt));
CGAL_exactness_precondition(line().has_on(pt));
for (unsigned int i=0; i< 3; ++i){
CGAL::Comparison_result c= compare(pt,Coordinate_index(i));
if (c != CGAL::EQUAL) {

40
Arrangement_of_spheres_3/include/CGAL/Arrangement_of_spheres_traits_3.h
View File

@ -139,11 +139,15 @@ struct Arrangement_of_spheres_traits_3 {
bool sphere_intersects_rule(Key sphere, Key rule_sphere,
Coordinate_index C) const;
bool sphere_intersects_sweep(Key sphere, Sphere_point_3 ep) const;
bool sphere_intersects_rule(Key sphere, const Sphere_point_3& sp,
Coordinate_index C) const;
bool sphere_intersects_sweep(Key sphere, const Sphere_point_3& ep) const;
CGAL::Comparison_result compare_equipower_point_to_rule(Key sphere0,
Key sphere1,
Key rule_sphere,
const Sphere_point_3& sp,
Coordinate_index C) const;
CGAL::Sign sign_of_separating_plane_normal_c(Key sphere0, Key sphere1,
@ -157,8 +161,8 @@ struct Arrangement_of_spheres_traits_3 {
CGAL::Oriented_side oriented_side_of_equipower_plane(Key sphere_0, Key sphere_1,
const Sphere_point_3 &s) const;
CGAL::Oriented_side oriented_side_of_center_plane(Key sphere_0, Key sphere_1,
Key sphere_center) const ;
CGAL::Oriented_side oriented_side_of_separating_plane(Key sphere_0, Key sphere_1,
const Sphere_point_3& sp) const ;
CGAL::Comparison_result compare_sphere_centers_c(Key a, Key b, Coordinate_index C) const;
@ -170,20 +174,32 @@ struct Arrangement_of_spheres_traits_3 {
const Sphere_point_3& d,
Coordinate_index C) const;
CGAL::Bounded_side bounded_side_of_sphere_on_equipower_plane_rule_line(Key sphere0,
Key Sphere1,
Key rule,
Coordinate_index C,
const Sphere_point_3 &sp) const;
// return true if the interval defined by the sphere in C contains d.C
bool is_over_circle_c(Key c, const Sphere_point_3& d,
Coordinate_index C) const;
CGAL::Comparison_result compare_sphere_sphere_at_sweep(Key sphere0,
Key Sphere1,
const Sphere_point_3 &sw,
const Sphere_point_3 &ep,
Coordinate_index C) const;
CGAL::Bounded_side bounded_side_of_sphere_projected(Key sphere,
const Sphere_point_3 &ep,
Key plane,
Coordinate_index C) const;
CGAL::Bounded_side bounded_side_of_sphere(Key sphere,
Key sphere_x,
Key sphere_y,
const Sphere_point_3 &z) const;
CGAL::Bounded_side bounded_side_of_sphere(Key sphere,
Key x,
Key y,
const Sphere_point_3 &z) const;
CGAL::Comparison_result compare_depths(const Sphere_point_3 &a,
const Sphere_point_3 &b) const;
CGAL::Oriented_side oriented_side(const Plane_3 &p,
const Sphere_point_3 &pt) const;
private:
// HDS hds_;

298
Arrangement_of_spheres_3/src/CGAL/Arrangement_of_spheres_traits_3.C
View File

@ -23,7 +23,7 @@ bool Arrangement_of_spheres_traits_3::intersects(Key a,
Point_3 eqpt= table_->equipower_point(a,b);
std::cout << eqpt << std::endl;
CGAL::Bounded_side bs=table_->sphere(a).bounded_side(eqpt);
CGAL::Bounded_side obs=table_->sphere(b).bounded_side(eqpt);
//CGAL::Bounded_side obs=table_->sphere(b).bounded_side(eqpt);
return bs != CGAL::ON_UNBOUNDED_SIDE;
} catch (Sphere_3_table::Equal_centers_exception) {
std::cout << "Equal centers for " << a << " and " << b << std::endl;
@ -56,7 +56,7 @@ bool Arrangement_of_spheres_traits_3::intersects(Key a,
bool Arrangement_of_spheres_traits_3::sphere_intersects_sweep(Key sphere,
Sphere_point_3 ep) const {
const Sphere_point_3& ep) const {
Event_pair be= sphere_events(sphere);
if (be.first.compare(ep, sweep_coordinate()) != CGAL::LARGER
&& be.second.compare(ep, sweep_coordinate()) != CGAL::SMALLER) return true;
@ -77,6 +77,24 @@ bool Arrangement_of_spheres_traits_3::sphere_intersects_rule(Key s,
}
bool Arrangement_of_spheres_traits_3::sphere_intersects_rule(Key s,
const Sphere_point_3 &sp,
Coordinate_index C) const {
FT v[3]={0,0,0};
CGAL::Comparison_result cr=sp.compare(table_->center(s), C);
if (cr == CGAL::SMALLER) {
v[C.index()]=-1;
} else if (cr == CGAL::LARGER) {
v[C.index()]=1;
} else {
return true;
}
Sphere_point_3 ep(table_->sphere(s), Line_3(table_->center(s),
Vector_3(v[0], v[1], v[2])));
return ep.compare(sp, C) != cr;
}
CGAL::Comparison_result
Arrangement_of_spheres_traits_3::compare_sphere_center_c(Key a,
const Sphere_point_3& d,
@ -95,12 +113,11 @@ Arrangement_of_spheres_traits_3::compare_sphere_center_c(Key a,
CGAL::Comparison_result
Arrangement_of_spheres_traits_3::compare_equipower_point_to_rule(Key a,
Key b,
Key c,
const Sphere_point_3 &sp,
Coordinate_index C) const{
// NOTE redo this with computing it directly
Point_3 pt= table_->equipower_point(a,b);
return CGAL::compare(pt[C.index()], table_->center(c)[C.index()]);
return CGAL::Comparison_result(-sp.compare(pt, C));
}
@ -129,19 +146,16 @@ CGAL::Sign Arrangement_of_spheres_traits_3::sign_of_equipower_plane_normal_c(Key
return sn;
}
CGAL::Oriented_side Arrangement_of_spheres_traits_3::oriented_side_of_equipower_plane(Key a, Key b,
const Sphere_point_3 &s) const {
CGAL_assertion(s.is_valid());
Plane_3 p= table_->equipower_plane(a,b);
CGAL::Oriented_side
Arrangement_of_spheres_traits_3::oriented_side(const Plane_3 &p,
const Sphere_point_3 &s) const{
CGAL::Object o= di_(p, s.line());
Point_3 pt;
Line_3 l;
if (CGAL::assign(pt, o)){
CGAL::Comparison_result c= s.compare_on_line(pt);
FT d= p.orthogonal_vector() * s.line().to_vector();
if (d > 0) return CGAL::enum_cast<CGAL::Oriented_side>(c);
CGAL::Sign sn= CGAL::sign(p.orthogonal_vector() * s.line().to_vector());
if (sn ==CGAL::POSITIVE) return CGAL::enum_cast<CGAL::Oriented_side>(c);
else return CGAL::enum_cast<CGAL::Oriented_side>(-c);
} else if (CGAL::assign( l, o)){
return CGAL::ON_ORIENTED_BOUNDARY;
@ -150,15 +164,19 @@ CGAL::Oriented_side Arrangement_of_spheres_traits_3::oriented_side_of_equipower_
}
}
CGAL::Oriented_side Arrangement_of_spheres_traits_3::oriented_side_of_equipower_plane(Key a, Key b,
const Sphere_point_3 &s) const {
CGAL_assertion(s.is_valid());
Plane_3 p= table_->equipower_plane(a,b);
return oriented_side(p, s);
}
CGAL::Oriented_side Arrangement_of_spheres_traits_3::oriented_side_of_center_plane(Key a, Key b,
Key sphere_center) const {
Vector_3 d(table_->center(b)-table_->center(a));
Line_2 l(Point_2(table_->center(a)[plane_coordinate(0).index()],
table_->center(a)[plane_coordinate(1).index()]),
Vector_2(d[plane_coordinate(0).index()], d[plane_coordinate(1).index()]));
return l.oriented_side(Point_2(table_->center(sphere_center)[plane_coordinate(0).index()],
table_->center(sphere_center)[plane_coordinate(1).index()]));
CGAL::Oriented_side
Arrangement_of_spheres_traits_3::oriented_side_of_separating_plane(Key a, Key b,
const Sphere_point_3 &ep) const {
Plane_3 sp= table_->separating_plane(a,b);
return oriented_side(sp, ep);
}
CGAL::Comparison_result Arrangement_of_spheres_traits_3::compare_sphere_centers_c(Key a, Key b, Coordinate_index C) const {
@ -167,56 +185,192 @@ CGAL::Comparison_result Arrangement_of_spheres_traits_3::compare_sphere_centers_
CGAL::Comparison_result
Arrangement_of_spheres_traits_3::compare_sphere_sphere_at_sweep(Key sphere0,
Key sphere1,
const Sphere_point_3 &sw,
const Sphere_point_3 &ep,
Coordinate_index C) const {
{
std::cout << "Compare SS at sweep " << sphere0 << " " << sphere1
<< " sweep " << sw << " line " << ep.line() << " coord " << C
<< std::endl;
bool i0= sphere_intersects_rule(sphere0, ep, C);
bool i1= sphere_intersects_rule(sphere1, ep, C);
//std::cout << "Comparisons are " << hi << " " << hp << " and " << ohi << " " << ohp << std::endl;
if (!i0) {
CGAL::Comparison_result sc0= compare_sphere_center_c(sphere0, ep, C);
std::cout << "Returning " << sc0 << " since sphere 0 misses" << std::endl;
return sc0;
} else if (!i1) {
CGAL::Comparison_result sc1= compare_sphere_center_c(sphere1, ep, C);
std::cout << "Returning " << sc1 << " since sphere 0 misses" << std::endl;
return sc1;
}
}
Plane_3 eqp= table_->equipower_plane(sphere0, sphere1);
std::cout << "EQP is " << eqp << std::endl;
FT v[3]={0,0,0};
v[other_plane_coordinate(C).index()]=1;
Plane_3 rule_plane(ep.line().point(),
ep.line().point()+ Vector_3(v[0], v[1], v[2]),
ep.line().point()+ ep.line().to_vector());
if (rule_plane.orthogonal_vector()[C.index()] < 0) {
rule_plane=rule_plane.opposite();
}
std::cout << "rule_plane is " << rule_plane << std::endl;
CGAL_assertion(rule_plane != Plane_3());
Line_3 l;
CGAL::Object o= di_(eqp, rule_plane);
if (!CGAL::assign(l,o)) {
CGAL::Comparison_result cr= compare_equipower_point_to_rule(sphere0, sphere1,
ep, C);
std::cout << "Returning " << cr << " since sphere the planes miss" << std::endl;
return cr;
}
std::cout << "The intersection line is " << l << std::endl;
Sphere_point_3 switch_point_0(table_->sphere(sphere0), l);
if (!switch_point_0.is_valid()){
CGAL::Comparison_result cr= compare_equipower_point_to_rule(sphere0,sphere1,
ep, C);
std::cout << "Returning " << cr << " line misses" << std::endl;
return cr;
}
//if (switch_point_0 > switch_point_1) {
// std::swap(switch_point_0, switch_point_1);
//}
Plane_3 sepp= table_->separating_plane(sphere0, sphere1);
std::cout << "separating_plane is " << sepp << std::endl;
Sphere_point_3 frontp= intersection_2_events(sphere0, sphere1).first;
std::cout << "The front point is " << frontp << std::endl;
CGAL::Oriented_side fpos= oriented_side(rule_plane, frontp);
bool flip=false;
if (fpos == CGAL::ON_ORIENTED_BOUNDARY) {
std::cout << "Front point is on boundary using normals ";
fpos = CGAL::Oriented_side(CGAL::sign(sepp.orthogonal_vector()
* rule_plane.orthogonal_vector()));
std::cout << "got " << fpos << std::endl;
}
std::cout << "Oriented_side is " << fpos << std::endl;
CGAL::Comparison_result spcr0= switch_point_0.compare(sw, sweep_coordinate());
Sphere_point_3 switch_point_1(table_->sphere(sphere0), l.opposite());
std::cout << "switch points are " << switch_point_0
<< " and " << switch_point_1 << std::endl;
CGAL::Comparison_result spcr1= switch_point_1.compare(sw, sweep_coordinate());
std::cout << "spcr is " << spcr0 << spcr1 << std::endl;
CGAL::Oriented_side os_0= oriented_side(sepp, switch_point_0);
CGAL::Oriented_side os_1= oriented_side(sepp, switch_point_1);
std::cout << "os are " << os_0 << os_1 << std::endl;
if (os_0 == CGAL::ON_POSITIVE_SIDE && os_1 == CGAL::ON_POSITIVE_SIDE) {
flip = (spcr0 != spcr1);
} else if (os_0 == CGAL::ON_POSITIVE_SIDE) {
flip = (spcr0 == CGAL::SMALLER);
} else if (os_1 == CGAL::ON_POSITIVE_SIDE) {
flip = (spcr1 == CGAL::LARGER);
}
if (!flip) {
return CGAL::Comparison_result(fpos);
} else {
return CGAL::Comparison_result(-fpos);
}
}
CGAL::Bounded_side
Arrangement_of_spheres_traits_3::bounded_side_of_sphere_on_equipower_plane_rule_line(Key sphere0,
Key sphere1,
Key rule,
Coordinate_index C,
const Sphere_point_3& sp) const {
Plane_3 eqp= table_->equipower_plane(sphere0, sphere1);
FT v[3], p[3];
v[C.index()]=1;
v[other_plane_coordinate(C).index()]=0;
v[sweep_coordinate().index()] = 0;
p[plane_coordinate(0).index()]=table_->center(rule)[plane_coordinate(0).index()];
p[plane_coordinate(1).index()]=table_->center(rule)[plane_coordinate(1).index()];
p[sweep_coordinate().index()]= 0;
Plane_3 pp(Point_3(p[0],p[1],p[2]),
Vector_3(v[0], v[1], v[2]));
Line_3 l;
CGAL::Object o= di_(eqp, pp);
CGAL_assertion(CGAL::assign(l,o));
CGAL::assign(l,o);
std::cout << "The intersection line is " << l << std::endl;
//CGAL::Comparison_result dir= CGAL::compare(ep.simple_coordinate(1-C), spheres_[h->curve().key()].table_->center(1-C));*/
Sphere_point_3 p0(table_->sphere(sphere0), l);
Arrangement_of_spheres_traits_3
::bounded_side_of_sphere_projected(Key sphere,
const Sphere_point_3 &ep,
Key plane,
Coordinate_index C) const{
FT p[3],v[3];
for (unsigned int i=0; i< 3; ++i) {
p[i]= ep.line().point()[i];
v[i]= ep.line().to_vector()[i];
}
v[other_plane_coordinate(C).index()] = 0;
p[other_plane_coordinate(C).index()] = table_->center(plane)[other_plane_coordinate(C).index()];
if (!p0.is_valid()){
Line_3 l(Point_3(p[0],p[1], p[2]),Vector_3(v[0], v[1], v[2]));
std::cout << "Line is " << l << std::endl;
Sphere_point_3 p0(table_->sphere(sphere), l);
if (!p0.is_valid()) {
std::cout << "Line misses" << std::endl;
return CGAL::ON_UNBOUNDED_SIDE;
} else {
Sphere_point_3 p1(table_->sphere(sphere0), l.opposite());
if (compare_depths(p0, p1)== CGAL::LARGER) {
std::swap(p0, p1);
}
//CGAL_assertion(p0 < p1);
if (compare_depths(sp,p0)==CGAL::EQUAL
|| compare_depths(sp, p1) == CGAL::EQUAL) {
Sphere_point_3 p1(table_->sphere(sphere), l.opposite());
std::cout << "Points are " << p0 << " " << p1 << std::endl;
CGAL::Comparison_result c0= compare_depths(ep, p0), c1= compare_depths(ep, p1);
std::cout << "Depths are " << c0 << " " << c1 << std::endl;
if (c0==CGAL::EQUAL || c1 == CGAL::EQUAL) {
return CGAL::ON_BOUNDARY;
} else if (c0 != c1) {
return CGAL::ON_BOUNDED_SIDE;
} else {
return CGAL::ON_UNBOUNDED_SIDE;
}
if (compare_depths(sp, p0) == CGAL::LARGER
&& compare_depths(sp, p1) == CGAL::SMALLER) return CGAL::ON_BOUNDED_SIDE;
else return CGAL::ON_UNBOUNDED_SIDE;
}
}
// compare along C where the intersection point of the two rules
// (defined by Coord C of ep) is relative to the sphere of the arc
CGAL::Bounded_side
Arrangement_of_spheres_traits_3::bounded_side_of_sphere(Key s,
const Sphere_point_3 &z) const{
/*FT p[3];
p[plane_coordinate(0).index()]= table_->center(x)[plane_coordinate(0).index()];
p[plane_coordinate(1).index()]= table_->center(y)[plane_coordinate(1).index()];
p[sweep_coordinate().index()]=0;
FT v[3];
v[plane_coordinate(0).index()]=0;
v[plane_coordinate(1).index()]=0;
v[sweep_coordinate().index()]=1;*/
Line_3 l= z.line();//(Point_3(p[0], p[1], p[2]), Vector_3(v[0],v[1],v[2]));
Sphere_point_3 p0(table_->sphere(s), l);
if (!p0.is_valid()) {
return CGAL::ON_UNBOUNDED_SIDE;
} else {
Sphere_point_3 p1(table_->sphere(s),
l.opposite());
CGAL::Comparison_result c0= compare_depths(p0, z);
CGAL::Comparison_result c1= compare_depths(p1, z);
if (c0 == CGAL::EQUAL || c1== CGAL::EQUAL) {
//std::cout << "Equal" << std::endl;
return CGAL::ON_BOUNDARY;
}
//CGAL_assertion(p0 <= p1);
if (c0== CGAL::SMALLER && c1 == CGAL::LARGER){
return CGAL::ON_BOUNDED_SIDE;
} else {
return CGAL::ON_UNBOUNDED_SIDE;
}
}
}
CGAL::Bounded_side
Arrangement_of_spheres_traits_3::bounded_side_of_sphere(Key s,
Key x,
@ -230,8 +384,6 @@ Arrangement_of_spheres_traits_3::bounded_side_of_sphere(Key s,
v[plane_coordinate(0).index()]=0;
v[plane_coordinate(1).index()]=0;
v[sweep_coordinate().index()]=1;
/*std::cout << "Line through " << CGAL::to_double(p[0])
<< " " << CGAL::to_double(p[1]) << std::endl;*/
Line_3 l(Point_3(p[0], p[1], p[2]), Vector_3(v[0],v[1],v[2]));
Sphere_point_3 p0(table_->sphere(s), l);
@ -258,6 +410,32 @@ Arrangement_of_spheres_traits_3::bounded_side_of_sphere(Key s,
bool
Arrangement_of_spheres_traits_3::is_over_circle_c(Key s,
const Sphere_point_3 &z,
Coordinate_index C) const {
Line_3 l= z.line();
FT p[3], v[3];
for (unsigned int i=0; i< 3; ++i) {
p[i]= l.point()[i];
v[i]= l.to_vector()[i];
}
Coordinate_index oc= other_plane_coordinate(C);
p[oc.index()]= table_->center(s)[oc.index()];
v[oc.index()]= 0;
Line_3 nl=Line_3(Point_3(p[0], p[1], p[2]),Vector_3(v[0], v[1], v[2]));
Sphere_point_3 sp(table_->sphere(s), nl);
if (!sp.is_valid()) return false;
else {
Sphere_point_3 spo(table_->sphere(s), nl.opposite());
CGAL::Comparison_result c= sp.compare(z, sweep_coordinate());
CGAL::Comparison_result co= spo.compare(z, sweep_coordinate());
if (co == c) return true;
else return false;
}
}
/*

1
Arrangement_of_spheres_3/src/CGAL/Combinatorial_stuff.C
View File

@ -344,6 +344,7 @@ Combinatorial_curve::constant_coordinate() const {
r= plane_coordinate(1);
else r= plane_coordinate(0);
//if (!is_finite()) r=1-r;
//std::cout << "Constant coordinate of " << *this << " is " << r << std::endl;
return r;
}

10
Arrangement_of_spheres_3/src/CGAL/Slice_arrangement.C
View File

@ -544,11 +544,13 @@ Slice_arrangement::Point Slice_arrangement::point(CArr::Vertex_const_handle h) c
}
Slice_arrangement::Curve Slice_arrangement::curve(CArr::Halfedge_const_handle h) const {
CGAL_assertion(std::distance(arr_.originating_curves_begin(h),
arr_.originating_curves_end(h))==1);
CGAL_assertion(map_.find(arr_.originating_curves_begin(h)) != map_.end());
return map_.find(arr_.originating_curves_begin(h))->second;
Slice_arrangement::Curve c= map_.find(arr_.originating_curves_begin(h))->second;
if (std::distance(arr_.originating_curves_begin(h),
arr_.originating_curves_end(h))!=1) {
CGAL_assertion(c.is_rule());
}
return c;
//return boost::apply_visitor(lookup_curve_data_, h->curve());
/*if (Circular_k::Circular_arc_2& a=boost::get<Circular_k::Circular_arc_2>(h->curve())) {

30
Arrangement_of_spheres_3/src/CGAL/Slice_data_structure.C
View File

@ -36,6 +36,24 @@ Slice_data_structure::rule_halfedge(Curve::Key k, int i) const {
return h;
}
Slice_data_structure::Halfedge_handle
Slice_data_structure::extremum_halfedge(Curve::Key k, int i) const {
Halfedge_handle h= halfedges_[k.input_index()][i];
CGAL_assertion(h->vertex()->point().is_sphere_extremum());
CGAL_assertion(h->curve().key() == k);
CGAL_assertion(!h->curve().is_rule());
return h;
}
void
Slice_data_structure::set_extremum_halfedge(Curve::Key k, int i,
Halfedge_handle h){
CGAL_precondition(h->curve().key() ==k);
CGAL_precondition(h->curve().is_arc());
CGAL_precondition(h->vertex()->point().is_sphere_extremum());
halfedges_[k.input_index()][i];
}
Slice_data_structure::Halfedge_handle
Slice_data_structure::next_edge_on_curve(Halfedge_handle h) const {
if (h->curve().is_rule()) {
@ -271,7 +289,7 @@ Slice_data_structure::exchange_spheres(Curve::Key k, Curve::Key l) {
if (kh->vertex()->point().is_sphere_extremum()) {
// relabel rule
Halfedge_handle r= kh->opposite()->prev()->opposite();
kes[r->curve().rule_index()]= r;
//kes[r->curve().rule_index()]= r;
CGAL_assertion(r->curve().is_rule());
if (r->curve().key() == k) {
Curve nr= r->curve();
@ -646,6 +664,7 @@ bool Slice_data_structure::has_vertex(Face_const_handle fh,
Halfedge_const_handle h= vh->halfedge();
do {
if (h->face() == fh) return true;
//if (h->opposite()->face() == fh) return true;
h= h->opposite()->prev();
CGAL_assertion(h->vertex() == vh);
} while (h != vh->halfedge());
@ -1203,8 +1222,10 @@ void Slice_data_structure::exchange_vertices(Halfedge_handle h,
Slice_data_structure::Face_handle
Slice_data_structure::intersect(Halfedge_handle ha, Halfedge_handle hb) {
std::pair<Slice_data_structure::Halfedge_handle,
Slice_data_structure::Halfedge_handle>
Slice_data_structure::intersect(Halfedge_handle ha, Halfedge_handle hb,
Halfedge_handle hc, Halfedge_handle hd) {
CGAL_precondition(ha->face() == hb->face());
Halfedge_handle ha_prev= ha->prev();
Halfedge_handle ha_op_next= ha->opposite()->next();
@ -1302,7 +1323,8 @@ Slice_data_structure::intersect(Halfedge_handle ha, Halfedge_handle hb) {
write_face(hb->face()->halfedge(), std::cout);
std::cout << std::endl;
return fn;
CGAL_assertion(0);
return std::pair<Halfedge_handle, Halfedge_handle>();;
}

132
Arrangement_of_spheres_3/src/CGAL/Slice_kds.C
View File

@ -212,43 +212,55 @@ void Slice::process_one_sphere_event(T::Key k) {
}
}
void Slice::process_two_sphere_event(T::Key k, T::Key l, bool f) {
std::cout << "Two sphere kds " << k << " " << l << " " << f << std::endl;
void Slice::process_equal_centers_two_sphere_event(T::Key k, T::Key l, bool f) {
Intersection_2 i2(k,l);
if (!f) {
intersections_2_.erase(i2);
} else {
intersections_2_[i2].first=Event_key();
if (t_.compare_sphere_centers_c(k,l, sweep_coordinate()) == CGAL::LARGER) {
std::swap(k,l);
}
// check if centers align
CGAL::Comparison_result c[2];
c[0]= t_.compare_sphere_centers_c(k,l, plane_coordinate(0));
c[1] t_.compare_sphere_centers_c(k,l, plane_coordinate(1));
if (c[0]== CGAL::EQUAL && c[1] == CGAL::EQUAL) {
if (t_.compare_sphere_centers_c(k,l, sweep_coordinate()) == CGAL::LARGER) {
std::swap(k,l);
}
sds_.exchange_spheres(k,l);
sim_->delete_event(intersections_2_[i2].second);
intersections_2_.erase(i2);
for (unsigned int i=0; i< 4; ++i) {
Halfedge_handle rh= sds_.rule_halfedge(k, i);
// now outside
clean_edge(rh);
check_edge_collapse(rh);
}
for (unsigned int i=0; i< 4; ++i) {
Halfedge_handle rh= sds_.rule_halfedge(l, i);
// now outside
clean_edge(rh);
if (!(rh->vertex()->point().is_sphere_rule()
&& rh->vertex()->point().sphere_key() == k)) {
sds_.exchange_spheres(k,l);
sim_->delete_event(intersections_2_[i2].second);
intersections_2_.erase(i2);
{
Halfedge_handle kc= sds_.a_halfedge(k), ke=kc;
do {
kc->curve().set_key(l);
kc= sds_.next_edge_on_curve(kc);
} while (kc != ke);
}
{
Halfedge_handle lc= sds_.a_halfedge(l), le=lc;
do {
lc->curve().set_key(k);
lc= sds_.next_edge_on_curve(lc);
} while (lc != le);
}
Halfedge_handle keh[4];
for (unsigned int i=0; i< 4; ++i) {
Halfedge_handle rh= sds_.rule_halfedge(k, i);
keh[i]= rh;
sds_.set_extremum_halfedge(k,i,sds_.extremum_halfedge(l,i));
// now outside
clean_edge(rh);
check_edge_collapse(rh);
}
for (unsigned int i=0; i< 4; ++i) {
Halfedge_handle rh= sds_.rule_halfedge(l, i);
// now outside
sds_.set_extremum_halfedge(l,i,keh[i]);
clean_edge(rh);
if (!(rh->vertex()->point().is_sphere_rule()
&& rh->vertex()->point().sphere_key() == k)) {
check_edge_collapse(rh);
}
}
} else if (c[0]== CGAL::EQUAL || c[1]== CGAL::EQUAL) {
// figure out which extremums
}
}
void Slice::process_aligned_centers_two_sphere_event(T::Key k, T::Key l, bool f,
CGAL::Comparison_result c[2]) {
// figure out which extremums
Intersection_2 i2(k,l);
const int equal_dir= (c[0]== CGAL::EQUAL)? 0:1;
const int other_dir= 1-equal_dir;
// make sure that k is smaller
@ -272,25 +284,38 @@ void Slice::process_two_sphere_event(T::Key k, T::Key l, bool f) {
// name them for who connects later
Halfedge_handle rule_vertex_l;
if (extr_k->next()->key() != k) {
if (extr_k->next()->curve().key() != k) {
clean_edge(extr_k->next());
std::pair<Halfedge_handle, Halfedge_handle> hr
= sds_.remove_rule(extra_k->next());
= sds_.remove_rule(extr_k->next());
rule_vertex_l= hr.first;
CGAL_assertion(hr.second == extr_k);
}
Halfedge_handle rule_vertex_k;
if (extr_l->next()->key() != l) {
if (extr_l->next()->curve().key() != l) {
clean_edge(extr_l->next());
std::pair<Halfedge_handle, Halfedge_handle> hr
= sds_.remove_rule(extra_l->next());
= sds_.remove_rule(extr_l->next());
rule_vertex_k= hr.first;
CGAL_assertion(hr.second == extr_l);
}
// maybe replace by two pinch operations
std::pair<Halfedge_handle, Halfedge_handle> hp=
intersect_2(extr_l->next(), extr_l->opposite(),
extr_k->next(), extr_k->opposite());
clean_edge(extr_l->opposite());
clean_edge(extr_k->opposite());
CGAL_assertion(0);
std::pair<Halfedge_handle, Halfedge_handle> hp;
/*=
sds_.intersect_2(extr_l->next(), extr_l->opposite(),
extr_k->next(), extr_k->opposite());*/
// clean edges and create certificates
CGAL_postcondition(hp.first->face() == hp.second->face());
CGAL_postcondition(hp.first->prev()->prev()== hp.second);
CGAL_postcondition(hp.second->prev()->prev()== hp.first);
check_edge_collapse(hp.first->opposite()->prev());
check_edge_collapse(hp.first->opposite()->prev()->opposite()->prev());
check_edge_collapse(hp.second->opposite()->prev());
check_edge_collapse(hp.second->opposite()->prev()->opposite()->prev());
CGAL_assertion(hp.first->curve().key() == k);
CGAL_assertion(hp.second->curve().key() == l);
@ -310,11 +335,11 @@ void Slice::process_two_sphere_event(T::Key k, T::Key l, bool f) {
Sds::Curve::make_rule(l, rl));
// set halfedges for spheres
sds_.set_extremum_edge(k, hp.first->prev());
sds_.set_extremum_edge(l, hp.second->prev());
sds_.set_extremum_halfedge(k, rk, hp.first->prev());
sds_.set_extremum_halfedge(l, rl, hp.second->prev());
} else {
// first figure out which case we are in
const T::Event_point_3 ep& = sim_->current_time();
const T::Event_point_3 &ep = sim_->current_time();
CGAL::Comparison_result c[2];
c[0]= t_.compare_sphere_center_c(k, ep, T::Coordinate_index(other_dir));
c[1]= t_.compare_sphere_center_c(l, ep, T::Coordinate_index(other_dir));
@ -327,6 +352,26 @@ void Slice::process_two_sphere_event(T::Key k, T::Key l, bool f) {
// then handle hard case
// handle easy case by accident
}
}
void Slice::process_two_sphere_event(T::Key k, T::Key l, bool f) {
#if 0
std::cout << "Two sphere kds " << k << " " << l << " " << f << std::endl;
Intersection_2 i2(k,l);
if (!f) {
intersections_2_.erase(i2);
} else {
intersections_2_[i2].first=Event_key();
}
// check if centers align
CGAL::Comparison_result c[2];
c[0]= t_.compare_sphere_centers_c(k,l, plane_coordinate(0));
c[1] t_.compare_sphere_centers_c(k,l, plane_coordinate(1));
if (c[0]== CGAL::EQUAL && c[1] == CGAL::EQUAL) {
process_equal_centers_two_sphere_event(k,l,f);
} else if (c[0]== CGAL::EQUAL || c[1]== CGAL::EQUAL) {
process_aligned_centers_two_sphere_event(k,l,f, c);
} else if (f) {
typedef std::pair<Halfedge_handle, Halfedge_handle> EP;
@ -380,6 +425,7 @@ void Slice::process_two_sphere_event(T::Key k, T::Key l, bool f) {
//CGAL_assertion(!shared_faces.empty());
unintersect_spheres(sim_->current_time(), f);
}
#endif
}
void Slice::process_three_sphere_event(T::Key k, T::Key l, T::Key m, bool f) {

2
Arrangement_of_spheres_3/src/CGAL/Slice_modifiers.C
View File

@ -389,7 +389,7 @@ Slice::Face_handle Slice::insert_sphere(const T::Sphere_point_3 &ep,
// the halfedge points to the vertex and is on this face
Vertex_handle vhs[4];
try {
f= locate_point(ep, k);
f= locate_point(ep);
} catch (On_edge_exception e) {
std::cout << "Point hit edge ";
sds_.write( e.halfedge_handle(), std::cout) << std::endl;

29
Arrangement_of_spheres_3/src/CGAL/Slice_point_location.C
View File

@ -15,7 +15,6 @@
bool Slice::locate_point_check_face(const T::Sphere_point_3 &z,
Face_const_handle it,
T::Key index,
std::vector<int> &locations) const {
Halfedge_const_handle h= it->halfedge();
//bool finite=false;
@ -26,12 +25,12 @@ bool Slice::locate_point_check_face(const T::Sphere_point_3 &z,
if (h->curve().is_right() && h->curve().is_inside()) return false;
} else {
T::Key sphere= h->curve().key();
point_sphere_orientation(z, index, sphere, locations);
point_sphere_orientation(z, sphere, locations);
//DPRINT(std::cout << "Testing " << h->curve() <<std::endl);
if (!h->curve().is_compatible_location(locations[sphere.input_index()])) {
//DPRINT(std::cout << "Nixed by edge " << h->curve() << std::endl);
return false;
}
}
}
//finite=true;
//} else {
@ -57,7 +56,6 @@ bool Slice::locate_point_check_face(const T::Sphere_point_3 &z,
// cache checked arcs
bool Slice::locate_point_check_face_arcs(const T::Sphere_point_3 &ep,
T::Key ind,
Face_const_handle f,
std::vector<int> &locations) const {
std::set<T::Key> check_arcs;
@ -67,7 +65,7 @@ bool Slice::locate_point_check_face_arcs(const T::Sphere_point_3 &ep,
&& check_arcs.find(h->curve().key()) == check_arcs.end()){
std::cout << "Arc test " << ep << " on " << h->curve() << std::endl;
check_arcs.insert(h->curve().key());
bool ba=behind_arc(ep, ind, h->curve(),
bool ba=behind_arc(ep, h->curve(),
locations[h->curve().key().input_index()]);
if (ba) {
std::cout << "Point is behind arc " << std::endl;
@ -84,8 +82,8 @@ bool Slice::locate_point_check_face_arcs(const T::Sphere_point_3 &ep,
bool Slice::locate_point_check_face_vertices(const T::Sphere_point_3 &ep,
T::Key index,
Face_const_handle it) const {
Halfedge_const_handle h= it->halfedge();
do {
if (h->vertex()->point().is_sphere_sphere()
@ -94,17 +92,17 @@ bool Slice::locate_point_check_face_vertices(const T::Sphere_point_3 &ep,
// NOTE what about degeneracies? not sure if I need to handle
// them here
Sds::Point npt= h->vertex()->point();
if (t_.oriented_side_of_center_plane(npt.sphere_key(0),
npt.sphere_key(1),
index) == CGAL::ON_NEGATIVE_SIDE) {
CGAL_assertion(0);
if (t_.oriented_side_of_separating_plane(npt.sphere_key(0),
npt.sphere_key(1),
ep) == CGAL::ON_NEGATIVE_SIDE) {
//CGAL_assertion(0);
std::cout << "Face nixed by vertex " << npt << std::endl;
return false;
}
}
h= h->next();
} while (h != it->halfedge());
return true;
}
@ -126,18 +124,17 @@ bool Slice::locate_point_check_face_vertices(const T::Sphere_point_3 &ep,
Slice::Face_handle Slice::locate_point(const T::Sphere_point_3 & ep) {
/*Slice::Face_handle Slice::locate_point(const T::Sphere_point_3 & ep) {
t_.set_temp_sphere(ep.sphere());
return locate_point(ep, T::Key::temp_key());
}
}*/
Slice::Face_handle
Slice::locate_point(const T::Sphere_point_3 & ep,
T::Key index) {
return locate_point(sds_.faces_begin(), sds_.faces_end(), ep, index);
Slice::locate_point(const T::Sphere_point_3 & ep) {
return locate_point(sds_.faces_begin(), sds_.faces_end(), ep);
}

64
Arrangement_of_spheres_3/src/CGAL/Slice_predicates.C
View File

@ -13,12 +13,9 @@
int Slice::sphere_location(const T::Sphere_point_3& sp,
T::Key locate_point,
T::Key s) const {
int r=0;
CGAL::Bounded_side bs=t_.bounded_side_of_sphere(s,
locate_point,
locate_point,
sp);
/*T::Event_point_3 ep=sphere_start(locate_point);
@ -37,17 +34,15 @@ int Slice::sphere_location(const T::Sphere_point_3& sp,
} else {
r |= Sds::Curve::lOUT_BIT;
}
/*T::Plane_3 lrp(center(s), T::Vector_3(1, 0, 0));
T::Plane_3 tbp(center(s), T::Vector_3(0, 1, 0));
CGAL::Comparison_result xo= oriented_side(lrp, ep);
CGAL::Oriented_side yo= oriented_side(tbp, ep);*/
//CGAL::Comparison_result xo= compare_event_point_to_rule(ep, s, 0);
//CGAL::Comparison_result yo= compare_event_point_to_rule(ep, s, 1);
CGAL::Comparison_result xo= t_.compare_sphere_centers_c(s, locate_point,
/*
NOTE can make it more specific for intersection point location
*/
CGAL::Comparison_result xo= t_.compare_sphere_center_c(s, sp,
plane_coordinate(0));
CGAL::Comparison_result yo= t_.compare_sphere_centers_c(s, locate_point,
plane_coordinate(1));
CGAL::Comparison_result yo= t_.compare_sphere_center_c(s, sp,
plane_coordinate(1));
if (xo != CGAL::LARGER) {
r |= Sds::Curve::lR_BIT;
}
@ -72,39 +67,19 @@ int Slice::sphere_location(const T::Sphere_point_3& sp,
bool Slice::behind_arc(const T::Sphere_point_3 &ep,
T::Key ind, Sds::Curve arc,
Sds::Curve arc,
int location) const{
CGAL_assertion(arc.is_compatible_location(location));
CGAL_precondition(!arc.is_inside());
CGAL_precondition(arc.is_compatible_location(location));
T::Coordinate_index C=arc.is_weakly_incompatible(location);
if (C.is_valid()) {
CGAL::Bounded_side bs;
if (C== plane_coordinate(0)) {
bs= t_.bounded_side_of_sphere(arc.key(),
arc.key(),
ind,
ep);
} else {
bs= t_.bounded_side_of_sphere(arc.key(),
ind,
arc.key(),
ep);
}
/*NT v[2];
v[C]= t_.center_c(arc.key(), C);
v[1-C]= t_.center_c(ind, 1-C);
T::Line_3 l(T::Point_3(v[0], v[1], 0),
T::Vector_3(0,0,1));
T::Event_point_3 fp(sphere(arc.key()),
l);
if (!fp.is_valid()) return false;
T::Event_point_3 bp(sphere(arc.key()),
l.opposite());
if (fp <= ep && bp >=ep) {
return true;
} else {return false;}
} else return false;*/
return bs== CGAL::ON_BOUNDED_SIDE;
if (C.is_valid()) { // we know it is opposite the center
/*CGAL::Comparison_result cr = t_.compare_sphere_center_c(arc.key(),
ep, C);
if (cr == CGAL::SMALLER && arc.is_negative()
|| cr == CGAL::LARGER && !arc.is_negative()) return false;*/
bool b= t_.is_over_circle_c(arc.key(), ep, C);
return b;
} else {
return false;
}
@ -112,13 +87,12 @@ bool Slice::behind_arc(const T::Sphere_point_3 &ep,
void Slice::point_sphere_orientation(const T::Sphere_point_3 &time,
T::Key front_point,
T::Key sphere,
std::vector<int> &locations
/*,
std::vector<Sds::Curve> &edges*/) const {
if (locations[sphere.input_index()]==0){
locations[sphere.input_index()]=sphere_location(time, front_point, sphere);
locations[sphere.input_index()]=sphere_location(time, sphere);
//DPRINT(std::cout << "For sphere " << sphere << " got " << SL::decode(locations[sphere]) << std::endl);
}
}

358
Arrangement_of_spheres_3/src/CGAL/Slice_shoot_rule.C
View File

@ -11,93 +11,55 @@ Slice::Halfedge_handle Slice::shoot_rule(const T::Sphere_point_3& source,
t_.set_temp_sphere(source.sphere());
return shoot_rule(source, f, Sds::Curve(T::Key::target_key(), Sds::Curve::Part(type)));
}
}
Slice::Halfedge_handle Slice::shoot_rule(const T::Sphere_point_3 & source,
Face_handle f,
Sds::Curve rule) {
Face_handle f,
Sds::Curve rule) {
bool backwards= (rule.is_top() || rule.is_left());
DPRINT(std::cout << "Rule is " << rule << std::endl);
Halfedge_handle h= f->halfedge();
std::vector<Halfedge_handle> edges;
bool has_last=false;
CGAL::Comparison_result last;
if (rule.can_intersect(h->curve())
&& rule.can_intersect(h->prev()->curve())) {
has_last=true;
last= rule_shoot_edge_vertex(source, rule, h->prev()->curve(),
h->prev()->vertex()->point(),
h->curve());
} else {
DPRINT(std::cout << "Eliminating " << h->prev()->vertex()->point() << std::endl);
Halfedge_handle end= f->halfedge();
while (!rule.can_intersect(end->curve())) {
end=end->next();
}
while (rule.can_intersect(end->curve())) {
end=end->next();
}
Halfedge_handle h=end;
while (!rule.can_intersect(h->curve())) {
h=h->next();
}
do {
// later skip second one if first is wrong
CGAL::Comparison_result pbs, nbs;
DPRINT(std::cout << "\nTesting edge: " << h->prev()->vertex()->point() << "--"
<< h->curve() << "--"
<< h->vertex()->point() << std::endl;)
if (rule.can_intersect(h->curve())) {
if (has_last) {
DPRINT(std::cout << "Using last of " << last << std::endl);
pbs=last;
} else if (backwards) {
pbs= CGAL::LARGER;
} else {
pbs=CGAL::SMALLER;
}
/* pbs= shoot_rule_edge_vertex(ep, rule, h->prev()->curve(),
h->prev()->vertex()->point(),
h->curve());
}*/
if (rule.can_intersect(h->next()->curve())) {
nbs= rule_shoot_edge_vertex(source, rule, h->curve(),
h->vertex()->point(),
h->next()->curve());
DPRINT(std::cout << "Results are " << nbs << " " << pbs << std::endl);
has_last=true;
last = nbs;
} else {
has_last=false;
if (backwards) nbs = CGAL::SMALLER;
else nbs = CGAL::LARGER;
}
if (!backwards && pbs != CGAL::LARGER && nbs != CGAL::SMALLER
|| backwards && pbs != CGAL::SMALLER && nbs != CGAL::LARGER) {
//std::cout << "it passes" << std::endl;
edges.push_back(h);
}
DPRINT(sds_.write(h, std::cout) << " is being tested" << std::endl);
if (h->next() == end) {
DPRINT(sds_.write(h, std::cout) << " is returned by default" << std::endl);
return h;
} else {
//std::cout << "Filtered" << std::endl;
has_last=false;
CGAL::Comparison_result cr= rule_shoot_edge_vertex(source, rule, h->curve(),
h->vertex()->point(),
h->next()->curve());
DPRINT(std::cout << "Result is " << cr << std::endl);
if (cr == CGAL::EQUAL) {
DPRINT(std::cout << "On vertex " << h->vertex()->point() << std::endl);
throw On_vertex_exception(h->vertex());
}
if (backwards && cr == CGAL::SMALLER
|| !backwards && cr == CGAL::LARGER) {
DPRINT(sds_.write(h, std::cout) << " is returned " << std::endl);
return h;
}
}
h= h->next();
CGAL_assertion(h->face()== f);
} while (h != f->halfedge());
} while (h != end);
if (edges.size() == 0) {
std::cerr << "No edge found " << std::endl;
return Halfedge_handle();
} else if (edges.size() ==1) {
return edges[0];
} else {
CGAL_assertion(edges.size()== 2);
std::cout << "Both edges " << edges[0]->curve() << " and "
<< edges[1]->curve()
<< " are OK" << std::endl;
if (edges[0]->vertex()== edges[1]->opposite()->vertex()) {
throw On_vertex_exception(edges[0]->vertex());
} else {
CGAL_assertion(edges[1]->vertex()== edges[0]->opposite()->vertex());
throw On_vertex_exception(edges[1]->vertex());
}
}
CGAL_assertion(0);
return Halfedge_handle();
}
@ -150,34 +112,23 @@ CGAL::Comparison_result Slice::rule_shoot_edge_vertex(const T::Sphere_point_3 &
return rule_shoot_compare_SR(ep, rule, hp, p.rule_key(),p, arc_top);
} else {
// rule from same sphere
CGAL_assertion(hp.is_inside());
bool exact;
CGAL::Comparison_result debug_answer= debug_rule_shoot_answer(ep, rule,
Sds::Point(hp, hn),
exact);
if (hp.part() & hn.part() & rule.part()) {
// check rational should get this
CGAL_assertion(0);
/*CGAL::Comparison_result ret
= CGAL::compare(t_.center_c(hn.key(), rule.constant_coordinate()),
ep.simple_coordinate(rule.constant_coordinate()));
debug_rule_shoot_check(debug_answer, ret, exact);*/
return CGAL::SMALLER;
} else if (hp.is_top()&& hn.is_top()
|| hp.is_right() && hn.is_right()
/*hp.part() & hn.part() & Sds::Curve::T_BIT
|| hp.part() & hn.part() & Sds::Curve::R_BIT*/) {
debug_rule_shoot_check(debug_answer, CGAL::LARGER, exact);
return CGAL::LARGER;
} else {
debug_rule_shoot_check(debug_answer, CGAL::SMALLER, exact);
return CGAL::SMALLER;
}
// compare rational should have picked this up
CGAL_assertion(0);
}
}
} else {
return rule_shoot_compare_SS(ep, rule, hp, p, hn);
bool tangent= (hp.key()== hn.key());
if (tangent) {
std::pair<T::Event_point_3, T::Event_point_3> ep3
= t_.intersection_2_events(p.sphere_key(0),
p.sphere_key(1));
std::cout << "Tangent." << std::endl;
CGAL::Comparison_result cr= ep3.first.compare(ep, sweep_coordinate());
//debug_rule_shoot_check(debug_answer, cr, exact);
return cr;
} else {
return rule_shoot_compare_SS(ep, rule,p);
}
}
CGAL_assertion(0);
@ -187,9 +138,7 @@ CGAL::Comparison_result Slice::rule_shoot_edge_vertex(const T::Sphere_point_3 &
CGAL::Comparison_result Slice::rule_shoot_compare_SS(const T::Sphere_point_3 & ep,
Sds::Curve srule,
Sds::Curve arc0,
Sds::Point pt,
Sds::Curve arc1) const {
Sds::Point pt) const {
bool exact;
CGAL::Comparison_result debug_answer= debug_rule_shoot_answer(ep, srule,
pt,
@ -197,84 +146,12 @@ CGAL::Comparison_result Slice::rule_shoot_compare_SS(const T::Sphere_point_3 & e
//std::cout << "SS " << arc0 << "--" << pt << "--" << arc1 << std::endl;
T::Coordinate_index C= srule.constant_coordinate();
bool tangent= (arc0.key()==arc1.key());
if (tangent) {
CGAL_assertion(0);
std::cout << "Tangent." << std::endl;
/*if (arc0.key() == pt.sphere_key(0)) arc1= pt.sphere(1);
else arc1= pt.sphere(0);
std::cout << "Using " << arc0 << "--" << pt << "--" << arc1 << std::endl;*/
}
// NOTE maybe make these into index based predicates
bool hi= t_.sphere_intersects_rule(arc0.key(), srule.key(), C);
bool ohi= t_.sphere_intersects_rule(arc0.key(), srule.key(), C);
CGAL::Comparison_result hp= t_.compare_sphere_centers_c(arc0.key(),
srule.key(), C);
CGAL::Comparison_result ohp= t_.compare_sphere_centers_c(arc1.key(),
srule.key(), C);
//std::cout << "Comparisons are " << hi << " " << hp << " and " << ohi << " " << ohp << std::endl;
if (hi && !ohi || !ohi && !hi) {
debug_rule_shoot_check(debug_answer,ohp, exact);
return ohp;
} else if (ohi && !hi) {
debug_rule_shoot_check(debug_answer,hp, exact);
return hp;
} /*else if (!ohp && !hp) {
debug_rule_shoot_check(debug_answer,hp, exact);
CGAL_assertion(hp == ohp);
return ohp;
}*/
// now they both intersect the plane
// NOTE hmmm, need to clean
CGAL::Bounded_side bs
= t_.bounded_side_of_sphere_on_equipower_plane_rule_line(arc0.key(),
arc1.key(),
srule.key(),
T::Coordinate_index(C),
ep);
if (bs == CGAL::ON_UNBOUNDED_SIDE){
//std::cout << "Missed." << std::endl;
CGAL::Comparison_result ret= t_.compare_equipower_point_to_rule(arc0.key(),
arc1.key(),
srule.key(), C);
debug_rule_shoot_check(debug_answer,ret, exact);
return ret;
} else if (bs== CGAL::ON_BOUNDARY) {
//std::cout << "Equal one of them " << std::endl;
debug_rule_shoot_check(debug_answer,CGAL::EQUAL, exact);
return CGAL::EQUAL;
} else {
if (arc0.is_inside() || arc1.is_inside()) {
// arcs are at extremum of face if I am checking
CGAL::Sign sn= t_.sign_of_separating_plane_normal_c(pt.sphere_key(0),
pt.sphere_key(1), C);
CGAL_assertion(sn != CGAL::ZERO);
CGAL::Comparison_result cret= CGAL::enum_cast<CGAL::Comparison_result>(sn);
debug_rule_shoot_check(debug_answer,
cret,
exact);
return cret;
} else {
//bool between =true;
CGAL::Sign dir= t_.sign_of_equipower_plane_normal_c(arc0.key(), arc1.key(),
C);
//CGAL::Comparison_result dir=CGAL::compare(eqp.orthogonal_vector()[C], 0);
int cum=0;
CGAL::Oriented_side os= t_.oriented_side_of_equipower_plane(arc0.key(), arc1.key(), ep);
if (os == CGAL::ON_POSITIVE_SIDE) ++cum;
if (dir == CGAL::POSITIVE) ++cum;
if (cum%2==0) {
debug_rule_shoot_check(debug_answer,CGAL::SMALLER, exact);
return CGAL::SMALLER;
} else {
debug_rule_shoot_check(debug_answer,CGAL::LARGER, exact);
return CGAL::LARGER;
}
}
}
CGAL::Comparison_result cr= t_.compare_sphere_sphere_at_sweep(pt.sphere_key(0),
pt.sphere_key(1),
ep, ep, C);
debug_rule_shoot_check(debug_answer, cr, exact);
return cr;
}
@ -296,45 +173,32 @@ CGAL::Comparison_result Slice::rule_shoot_compare_SR(const T::Sphere_point_3 & e
//CGAL_precondition(orule.is_rule());
//CGAL_precondition(srule.constant_coordinate() != orule.constant_coordinate());
CGAL::Comparison_result c= t_.compare_sphere_center_c(arc.key(), ep,
srule.constant_coordinate());
if (srule.is_vertical()) {
CGAL::Comparison_result c= t_.compare_sphere_centers_c(srule.key(), arc.key(),
plane_coordinate(0));
if (c== CGAL::SMALLER && arc.is_right()) {
if (c== CGAL::LARGER && arc.is_right()) {
debug_rule_shoot_check(debug_answer, CGAL::LARGER, exact);
return CGAL::LARGER;
} else if (c== CGAL::LARGER && arc.is_left()) {
} else if (c== CGAL::SMALLER && arc.is_left()) {
debug_rule_shoot_check(debug_answer, CGAL::SMALLER, exact);
return CGAL::SMALLER;
}
} else {
CGAL::Comparison_result c= t_.compare_sphere_centers_c(srule.key(), arc.key(),
plane_coordinate(1));
if (c== CGAL::SMALLER && arc.is_top()) {
debug_rule_shoot_check(debug_answer, CGAL::LARGER, exact);
return CGAL::LARGER;
} else if (c== CGAL::LARGER && !arc.is_top()) {
if (c== CGAL::SMALLER && arc.is_bottom()) {
debug_rule_shoot_check(debug_answer, CGAL::SMALLER, exact);
return CGAL::SMALLER;
} else if (c== CGAL::LARGER && arc.is_top()) {
debug_rule_shoot_check(debug_answer, CGAL::LARGER, exact);
return CGAL::LARGER;
}
}
CGAL::Bounded_side intersection_side;
if (srule.constant_coordinate() == plane_coordinate(0)) {
intersection_side= t_.bounded_side_of_sphere(arc.key(),
srule.key(),
orule,
ep);
} else {
intersection_side= t_.bounded_side_of_sphere(arc.key(),
orule,
srule.key(),
ep);
}
//ep, orule,
// arc.key(),
// srule.constant_coordinate());
CGAL::Bounded_side intersection_side
=t_.bounded_side_of_sphere_projected(arc.key(),
ep,
orule,
srule.constant_coordinate());
if (intersection_side== CGAL::ON_BOUNDARY){
@ -418,7 +282,7 @@ void Slice::debug_rule_shoot_check(CGAL::Comparison_result check,
}
bool Slice::rule_shoot_compare_if_rational_arc(const T::Sphere_point_3 & ep,
/*bool Slice::rule_shoot_compare_if_rational_arc(const T::Sphere_point_3 & ep,
Sds::Curve rule,
Sds::Curve a,
CGAL::Comparison_result &ret) const {
@ -443,7 +307,7 @@ bool Slice::rule_shoot_compare_if_rational_arc(const T::Sphere_point_3 & ep,
}
}
return false;
}
}*/
bool Slice::rule_shoot_compare_if_rational(const T::Sphere_point_3 & ep,
@ -452,82 +316,20 @@ bool Slice::rule_shoot_compare_if_rational(const T::Sphere_point_3 & ep,
CGAL::Comparison_result &ret) const {
//NT coord;
if (pt.is_rule_rule()) {
ret= t_.compare_sphere_centers_c(pt.rule_key(rule.constant_coordinate()),
rule.key(),
rule.constant_coordinate());
ret= t_.compare_sphere_center_c(pt.rule_key(rule.constant_coordinate()),
ep,
rule.constant_coordinate());
return true;
} else if (pt.is_sphere_rule()
&& pt.rule_coordinate() == rule.constant_coordinate()) {
ret= t_.compare_sphere_centers_c(pt.rule_key(),
rule.key(),
rule.constant_coordinate());
ret= t_.compare_sphere_center_c(pt.rule_key(),
ep,
rule.constant_coordinate());
return true;
} else {
// need the vertex_handle
// actually, I don't think it was ever useful
return false;
}
/*if (a.is_arc() && rule_shoot_compare_if_rational_arc(ep, rule, a, ret)){
return true;
}
if (b.is_arc() && rule_shoot_compare_if_rational_arc(ep, rule, b, ret)){
return true;
}
return false;
}*/
/*std::cout << "It is rational with coordinate "
<< CGAL::to_double(coord) << std::endl;*/
/*ret= CGAL::compare(coord,
t_.center(rule.key())[rule.constant_coordinate()]);
return true;*/
}
// compare along C where the intersection point of the two rules
// (defined by Coord C of ep) is relative to the sphere of the arc
#if 0
CGAL::Bounded_side Slice::rule_shoot_source_side(const T::Sphere_point_3 & ep,
T::Key rule,
Sds::Curve orule,
T::Key arc_index,
int C) const {
return sc_.bounded_side_of_sphere_rule_rule_line(epi,
orule.index(),
arc_index,
C,
ep);
/*NT p[2];
p[C]= ep.simple_coordinate(C);
p[1-C]= center_c(orule.key(), orule.constant_coordinate());
std::cout << "Line through " << CGAL::to_double(p[0])
<< " " << CGAL::to_double(p[1]) << std::endl;
T::Line_3 l(T::Point_3(p[0], p[1], 0), T::Vector_3(0,0,1));
const T::Sphere_point_3 & p0=const T::Sphere_point_3 &(sphere(arc_index), l);
CGAL::Bounded_side intersection_side;
if (!p0.is_valid()) {
intersection_side=CGAL::ON_UNBOUNDED_SIDE;
} else if (p0 == ep) {
std::cout << "Equal" << std::endl;
intersection_side= CGAL::ON_BOUNDARY;
} else {
const T::Sphere_point_3 & p1=const T::Sphere_point_3 &(sphere(arc_index),
l.opposite());
if (p1 == ep) {
std::cout << "Equal 1" << std::endl;
intersection_side= CGAL::ON_BOUNDARY;
}
CGAL_assertion(p0 <= p1);
if (p0 < ep && p1 > ep){
intersection_side= CGAL::ON_BOUNDED_SIDE;
} else {
intersection_side= CGAL::ON_UNBOUNDED_SIDE;
}
}
return intersection_side;*/
}
#endif

3
Arrangement_of_spheres_3/src/CGAL/makefile.internal
View File

@ -15,7 +15,8 @@ include $(CGAL_MAKEFILE)
# compiler flags
#---------------------------------------------------------------------#
# -fmudflapth
CXXFLAGS = -I../../include -Wall -DCGAL_USE_GMP -DCGAL_USE_GMPXX -DCGAL_USE_CORE -DCGAL_USE_QT -DQT_THREAD_SUPPORT '-I$(CGAL_INCL_DIR)' '-I$(CGAL_INCL_DIR)/CGAL/config/$(CGAL_OS_COMPILER)' '-I$(QTDIR)/include' -g -DCGAL_CHECK_EXPENSIVE -DCGAL_CHECK_EXACTNESS
CXXFLAGS = -I../../include -Wall -DCGAL_USE_GMP -DCGAL_USE_GMPXX -DCGAL_USE_CORE -DCGAL_USE_QT -DQT_THREAD_SUPPORT '-I$(CGAL_INCL_DIR)' '-I$(CGAL_INCL_DIR)/CGAL/config/$(CGAL_OS_COMPILER)' '-I$(QTDIR)/include' -g -DCGAL_CHECK_EXPENSIVE
# -DCGAL_CHECK_EXACTNESS
# $(CGAL_LIB_CXXFLAGS) -I../include
#