Dummy infimaximal box added. Shell selection on ray shotting corrected.

This commit is contained in:
Miguel Granados 2002-05-02 19:42:23 +00:00
parent d6f15d98cf
commit 9fc9087a7c
7 changed files with 247 additions and 114 deletions

View File

@ -1,10 +1,6 @@
Nef_3 Package: TODO
---------------------------------------------------------------------
* Create testing models for volume synthesizing.
* Synthesizing volumes from SM representation (finish tests)
- Initializate all Nef_polyhedron structures with a bounding box of
fixed corners.
- Simplification of the SNC structure
- Simplify during the construction from polyhedral surfaces
(based on the simplification of SM)

View File

@ -47,7 +47,7 @@
#define _DEBUG 43
#include <CGAL/Nef_3/debug.h>
#define IMMN 999999999
#define IMMN 12345
CGAL_BEGIN_NAMESPACE
@ -199,6 +199,7 @@ public:
USING(Point_3);
USING(Vector_3);
USING(Direction_3);
USING(Segment_3);
USING(Line_3);
USING(Plane_3);
@ -322,77 +323,189 @@ public:
return c;
}
Halffacet_handle get_visible_facet( Vertex_handle v, Sphere_point p) const {
TRACEN( "Locating " << p <<" in " << point(v));
Halffacet_handle get_visible_facet( const Vertex_handle v,
const Segment_3& ray) const {
Halffacet_handle f_visible;
CGAL_assertion( ray.target() == point(v));
Sphere_point sp(ray.source() - point(v));
TRACEN( "Locating " << sp <<" in " << point(v));
SM_point_locator L(v);
SObject_handle o = L.locate(p);
SObject_handle o = L.locate(sp);
SFace_const_handle sf;
CGAL_assertion( assign( sf, o));
/* the ray must not intersect the objects incident to v */
assign( sf, o);
SFace_cycle_const_iterator fc;
CGAL_forall_sface_cycles_of( fc, sf) {
SFace_cycle_const_iterator fc = sf->sface_cycles_begin(),
fce = sf->sface_cycles_end();
if( is_empty_range( fc, fce)) { // TO TEST: condition never satisfied
TRACEN( "no adjacent facet found.");
f_visible = Halffacet_handle();
}
else {
SHalfedge_handle se;
if ( assign( se, fc) ) {
TRACEN( "adjacent facet found.");
return facet(twin(se));
SHalfloop_handle sl;
if ( assign( se, fc)) {
TRACEN( "adjacent facet found (SEdges cycle).");
f_visible = facet(twin(se));
}
else if ( assign( sl, fc)) {
TRACEN( "adjacent facet found (SHalfloop cycle).");
f_visible = facet(twin(sl));
}
else
CGAL_assertion_msg(0, "Damn, wrong handle.");
}
return f_visible;
}
typedef typename Kernel::Direction_2 Direction_2;
typedef typename Kernel::Vector_2 Vector_2;
typedef typename Kernel::Point_2 Point_2;
bool strictly_ordered_ccw(const Direction_2& d1,
const Direction_2& d2,
const Direction_2& d3) const
/*{\Mop returns |true| iff |d2| is in the interior of the
counterclockwise angular sector between |d1| and |d3|.}*/
{
if ( d1 < d2 ) return ( d2 < d3 )||( d3 <= d1 );
if ( d1 > d2 ) return ( d2 < d3 )&&( d3 <= d1 );
return false;
}
bool strictly_ordered_ccw_3( const Direction_3& a,
const Direction_3& b, const Direction_3& c, const Plane_3& h) const {
TRACEN("--> facet plane " << h);
CGAL_assertion( !h.is_degenerate());
Point_2 (*t)(Point_3);
Vector_3 hv( h.orthogonal_vector()),
hxy( 0, 0, 1), hyz( 1, 0, 0), hxz( 0, 1, 0);
if( !is_zero( hv * hxz) )
/* the plane is not perpendicular to the XZ plane */
t = &point_3_get_x_z_point_2< Point_2, Point_3>;
else if( !is_zero( hv * hyz) )
/* the plane is not perpendicular to the XZ plane */
t = &point_3_get_y_z_point_2< Point_2, Point_3>;
else {
CGAL_assertion( !is_zero( hv * hxy) );
/* the plane is not perpendicular to the XY plane */
t = &point_3_get_x_y_point_2< Point_2, Point_3>;
}
Direction_2 da( t( ORIGIN + Vector_3(a) ) - ORIGIN),
db( t( ORIGIN + Vector_3(b)) - ORIGIN),
dc( t( ORIGIN + Vector_3(c)) - ORIGIN);
TRACE( "CCW ordered " << a << ", " << b << ", " << c << ": ");
TRACE( "CCW ordered " << da << ", " << db << ", " << dc << ": ");
TRACEN( strictly_ordered_ccw( da, db, dc));
return strictly_ordered_ccw( da, db, dc);
}
Halffacet_handle get_visible_facet( const Halfedge_handle e,
const Segment_3& ray) const {
CGAL_assertion( segment(e).has_on(ray.target()));
SM_decorator D;
if( D.is_isolated(e))
return Halffacet_handle();
Direction_3 d_r = -ray.direction(), d_e = segment(e).direction();
Plane_3 h_ref( segment(e).source(), d_e);
SHalfedge_handle se = D.first_out_edge(e);
Plane_3 h_se = plane(facet(twin(se)));
Vector_3 v_se = cross_product( Vector_3(d_e),
Vector_3(h_se.orthogonal_direction()));
Direction_3 d_se = Direction_3(v_se);
SHalfedge_around_svertex_circulator sc(D.out_edges(e)), sce(sc);
CGAL_For_all(sc, sce) {
Plane_3 h_sc = plane(facet(twin(sc)));
Vector_3 v_sc = cross_product( Vector_3(d_e),
Vector_3(h_sc.orthogonal_direction()));
Direction_3 d_sc = Direction_3(v_sc);
if ( strictly_ordered_ccw_3(d_se, d_sc, d_r, h_ref)) {
se = sc;
d_se = d_sc;
}
}
// TO TEST: code never reached
TRACEN( "no adjacent facet found.");
return Halffacet_handle();
TRACEN(" --> CCW order end.");
return facet(twin(se));
}
Halffacet_handle get_visible_facet( const Halffacet_handle f,
const Segment_3& ray) const {
Halffacet_handle f_visible;
Plane_3 h = plane(f);
CGAL_assertion( h.has_on(ray.target()));
CGAL_assertion( !is_zero(Vector_3(ray.direction()) *
Vector_3(h.orthogonal_vector())));
/* is imposible to reach the interior or f using a ray coplanar with f */
if( is_negative(Vector_3(ray.direction()) *
Vector_3(h.orthogonal_vector())))
f_visible = f;
else
f_visible = twin(f);
CGAL_assertion( is_negative( Vector_3(ray.direction()) *
Vector_3(plane(f_visible).
orthogonal_vector())));
return f_visible;
}
Halffacet_handle get_facet_below( Vertex_handle vi) const {
Segment_3 s(point(vi), Point_3( 0, 0, -IMMN));
/* TODO: replace IMMN constant for a real infimaximal number */
TRACEN( "Shoting ray " << s);
Halffacet_handle f_below;
Segment_3 ray(point(vi), Point_3( 0, 0, -IMMN));
/* TODO: replace IMMN constant for a real infimaximal number */
Object_handle o = ray_shot(ray);
Vertex_handle v;
Halfedge_handle e;
Halffacet_handle f;
if( assign(v, o)) {
f_below = get_visible_facet(v, ray);
if( f_below == Halffacet_handle())
f_below = get_facet_below(v);
}
else if( assign(e, o)) {
f_below = get_visible_facet(e, ray);
if( f_below == Halffacet_handle())
f_below = get_facet_below(vertex(e));
}
else if( assign(f, o)) {
f_below = get_visible_facet(f, ray);
CGAL_assertion( f_below != Halffacet_handle());
}
return f_below;
}
Object_handle ray_shot( Segment_3& ray) const {
TRACEN( "Shoting ray " << ray);
Object_handle o;
Vertex_handle v;
CGAL_forall_vertices( v, *sncp()) {
if ( contains_internally( s, point(v))) {
if ( ray.source() != point(v) && ray.has_on(point(v)) ) {
TRACEN("ray hit vertex case");
Sphere_point sp(point(vi)-point(v));
Halffacet_handle f_visible = get_visible_facet(v, sp);
if ( f_visible != Halffacet_handle()) {
shorten( s, point(v));
f_below = f_visible;
}
shorten( ray, point(v));
o = Object_handle(v);
}
}
Halfedge_handle e;
CGAL_forall_edges( e, *sncp()) {
Point_3 q;
if ( does_intersect_internally( s, e, q) ) {
// TO TEST: code never reached
if ( does_ray_intersect_internally( ray, e, q) ) {
TRACEN("ray hit edge case");
v = vertex(e);
Sphere_point sp(point(vi)-point(v));
Halffacet_handle f_visible = get_visible_facet(v, sp);
if ( f_visible != Halffacet_handle()) {
shorten( s, q);
f_below = f_visible;
}
shorten( ray, q);
o = Object_handle(e);
}
}
Halffacet_handle f;
CGAL_forall_halffacets( f, *sncp()) {
Point_3 q;
if ( does_intersect_internally( s, f, q) ) {
if ( does_ray_intersect_internally( ray, f, q) ) {
TRACEN("ray hit facet case");
Halffacet_cycle_iterator fc(f->facet_cycles_begin());
CGAL_assertion( fc != f->facet_cycles_end());
SHalfedge_handle se;
CGAL_assertion( assign( se, fc) );
v = vertex(se);
Sphere_point sp(point(vi)-point(v));
Halffacet_handle f_visible = get_visible_facet(v, sp);
if ( f_visible != Halffacet_handle()) {
shorten( s, q);
f_below = f_visible;
}
shorten( ray, q);
o = Object_handle(f);
}
}
return f_below;
return o;
}
void shorten(Segment_3& s, const Point_3& p) const {
@ -407,19 +520,22 @@ public:
return (r1 == opposite(r2));
}
bool does_intersect_internally( const Segment_3& ray,
const Halfedge_handle e,
Point_3& p) const {
CGAL_assertion( !ray.is_degenerate());
return does_intersect_internally( ray, segment(e), p);
bool does_ray_intersect_internally( const Segment_3& ray,
const Halfedge_handle e,
Point_3& p) const {
return does_ray_intersect_internally( ray, segment(e), p);
}
#ifdef LINE3_LINE3_INTERSECTION
bool does_intersect_internally( const Segment_3& ray,
const Segment_3& s,
bool does_ray_intersect_internally( const Segment_3& ray,
const Segment_3& s,
Point_3& p) const {
CGAL_assertion( !ray.is_degenerate());
if ( s.is_degenerate())
return false;
if ( s.has_on(ray.source()) )
return false;
Object o = intersection(Line_3(ray), Line_3(s));
if ( !assign(p, o))
return false;
@ -428,13 +544,16 @@ public:
#else // LINE3_LINE3_INTERSECTION
bool does_intersect_internally( const Segment_3& ray,
const Segment_3& s,
Point_3& p) const {
bool does_ray_intersect_internally( const Segment_3& ray,
const Segment_3& s,
Point_3& p) const {
CGAL_assertion( !ray.is_degenerate());
if ( s.is_degenerate())
return false;
/* the segment is degenerate so there is not internal intersection */
if ( s.has_on(ray.source()) )
return false;
/* the segment contains the ray source */
if ( orientation( ray.source(), ray.target(), s.source(), s.target())
!= COPLANAR)
return false;
@ -475,19 +594,36 @@ public:
}
#endif // LINE3_LINE3_INTERSECTION
bool does_intersect_internally( const Segment_3& ray,
const Halffacet_handle f,
Point_3& p) const {
bool does_ray_intersect_internally( const Segment_3& ray,
const Halffacet_handle f,
Point_3& p) const {
// TRACEN("-> Intersection face - ray");
Plane_3 h( plane(f));
// TRACEN("-> facet plane " << h);
// TRACEN("-> a point on " << h.point());
// TRACEN("-> ray segment " << ray);
CGAL_assertion( !h.is_degenerate());
CGAL_assertion( !ray.is_degenerate());
Plane_3 h(plane(f));
Object o = intersection( h, Line_3(ray));
if ( !assign( p, o) )
if( h.has_on( ray.source()))
return false;
Object o = intersection( h, ray);
Segment_3 s;
if ( assign( s, o) ) {
CGAL_assertion( s == ray );
// TRACEN( "-> ray belongs to facet's plane." << p );
return false;
}
else if( !assign( p, o))
return false;
// TRACEN( "-> intersection point " << p );
Oriented_side os1 = h.oriented_side(ray.source());
Oriented_side os2 = h.oriented_side(ray.target());
if (os1 != opposite(os2))
// TRACEN( "-> endpoint plane side " << os1 << " " << os2);
CGAL_assertion( h.has_on(p));
CGAL_assertion( ray.has_on(p));
if (os1 == os2)
return false;
// TRACEN( "Point in facet result = "<<locate_point_in_halffacet(p, f));
TRACEN( "-> point in facet? "<<locate_point_in_halffacet(p, f));
return (locate_point_in_halffacet( p, f) == CGAL::ON_BOUNDED_SIDE);
}
@ -567,8 +703,9 @@ create_box_corner(int x, int y, int z,
bool boundary=true) const
{
CGAL_nef3_assertion(x*y*z != 0);
CGAL_nef3_assertion(x==y==z==1);
Vertex_handle v = sncp()->new_vertex();
int R=3; point(v) = Point_3(x*R,y*R,z*R);
int R=IMMN; point(v) = Point_3(x*R,y*R,z*R);
Sphere_point px(-x,0,0), py(0,-y,0), pz(0,0,-z);
std::list<Sphere_segment> L;
L.push_back(Sphere_segment(px,py));
@ -843,24 +980,20 @@ create_volumes() const
V.increment_shell_number();
}
Volume_handle outer_volume = sncp()->new_volume();
for (unsigned i = 0; i < MinimalVertex.size(); ++i) {
for( unsigned i = 0; i < MinimalVertex.size(); ++i) {
Vertex_handle v = MinimalVertex[i];
TRACEN("Shell #"<<i<<" minimal vertex: "<<point(v));
if ( true) { // TODO: if v is not a bounding box vertex
bool create_shell = false;
SM_point_locator D(v);
SObject_handle o = D.locate(Sphere_point(-1,0,0));
SFace_const_handle sfc;
if( !assign(sfc, o) || Shell[sfc] != i) {
// TO TEST: !assign(sfc, o) case
SFace_handle f = EntrySFace[i];
CGAL_assertion( Shell[EntrySFace[i]] == i );
if( Closed[f] ) {
Volume_handle c = sncp()->new_volume();
link_as_outer_shell(f, c );
TRACE( "Shell #" << i <<" linked as outer shell");
TRACEN( "(sface" << (assign(sfc,o)?"":" not") << " hit case)");
}
TRACEN( "Shell #" << i << " minimal vertex: " << point(v));
SM_point_locator D(v);
SObject_handle o = D.locate(Sphere_point(-1,0,0));
SFace_const_handle sfc;
if( !assign(sfc, o) || Shell[sfc] != i) { // TO TEST: !assign(sfc, o) case
SFace_handle f = EntrySFace[i];
CGAL_assertion( Shell[EntrySFace[i]] == i );
if( Closed[f] ) {
Volume_handle c = sncp()->new_volume();
link_as_outer_shell(f, c );
TRACE( "Shell #" << i <<" linked as outer shell");
TRACEN( "(sface" << (assign(sfc,o)?"":" not") << " hit case)");
}
}
}

View File

@ -86,6 +86,7 @@ public:
USING(Segment_3);
USING(Line_3);
USING(Vector_3);
USING(Direction_3);
USING(Sphere_point);
USING(Sphere_segment);
USING(Sphere_circle);
@ -102,6 +103,8 @@ public:
Vertex_handle vertex( Halfedge_handle e) const
{ return e->center_vertex_; }
Vertex_const_handle vertex( Halfedge_const_handle e) const
{ return e->center_vertex_; }
Halfedge_handle twin( Halfedge_handle e) const
{ return e->twin_; }
Vertex_handle source( Halfedge_handle e) const
@ -176,6 +179,8 @@ public:
{ return f->twin_; }
Volume_handle volume(Halffacet_handle f) const
{ return f->volume_; }
Volume_const_handle volume(Halffacet_const_handle f) const
{ return f->volume_; }
/* Halffacet queries */
SFace_handle adjacent_sface(Halffacet_handle f) const {
@ -197,6 +202,8 @@ public:
// attributes::
Point_3& point(Vertex_handle v) const
{ return v->point(); }
const Point_3& point(Vertex_const_handle v) const
{ return v->point(); }
Sphere_point tmp_point(Halfedge_handle e) const
{ return e->tmp_point(); }

View File

@ -97,6 +97,8 @@ public:
/*{\Mtypemember supporting facets.}*/
typedef typename Kernel::Vector_3 Vector_3;
/*{\Mtypemember normal vectors.}*/
typedef typename Kernel::Direction_3 Direction_3;
/*{\Mtypemember normal directions.}*/
typedef typename Kernel::Segment_3 Segment_3;
/*{\Mtypemember segments in space.}*/
typedef typename Kernel::Line_3 Line_3;

View File

@ -69,6 +69,7 @@ Bounded_side bounded_side_3(IteratorForward first,
typedef typename R::Plane_3 Plane_3;
if(plane == Plane_3()) {
// TO TEST: code never tested
IteratorForward p(first);
Point_3 p0(*(p++));
CGAL_assertion(p != last);
@ -80,15 +81,18 @@ Bounded_side bounded_side_3(IteratorForward first,
we don't need to care about the plane orientation */
}
CGAL_assertion(!plane.is_degenerate());
TRACEN(plane);
Point_2 (*t)(Point_3);
Vector_3 pv(plane.orthogonal_vector()),
pxy(0,0,1), pyz(1,0,0), pxz(0,1,0);
Vector_3 pv(plane.orthogonal_vector()), pxy(0,0,1), pyz(1,0,0), pxz(0,1,0);
TRACEN("pv*pxz: "<<pv*pxz);
TRACEN("pv*pyz: "<<pv*pyz);
TRACEN("pv*pxy: "<<pv*pxy);
if( !is_zero(pv*pxz) )
/* the plane is not perpendicular to the XZ plane */
t = &point_3_get_y_z_point_2< Point_2, Point_3>;
t = &point_3_get_x_z_point_2< Point_2, Point_3>;
else if( !is_zero(pv*pyz) )
/* the plane is not perpendicular to the XZ plane */
t = &point_3_get_x_z_point_2< Point_2, Point_3>;
t = &point_3_get_y_z_point_2< Point_2, Point_3>;
else {
CGAL_assertion( !is_zero(pv*pxy) );
/* the plane is not perpendicular to the XY plane */
@ -103,7 +107,7 @@ Bounded_side bounded_side_3(IteratorForward first,
}
Bounded_side side = bounded_side_2( points.begin(), points.end(), t(point));
points.clear();
return side;
return side;
}
CGAL_END_NAMESPACE
@ -147,9 +151,9 @@ struct Project_XZ {
}
};
template <class ForwardIterator, class R>
Bounded_side bounded_side_3(ForwardIterator first,
ForwardIterator last,
template <class IC, class R>
Bounded_side bounded_side_3(IC first,
IC last,
const Point_3<R>& point,
Plane_3<R> plane = Plane_3<R>()) {
@ -159,21 +163,19 @@ Bounded_side bounded_side_3(ForwardIterator first,
typedef typename R::Direction_3 Direction_3;
typedef typename R::Plane_3 Plane_3;
CGAL_assertion( !CGAL::is_empty_range( first, last));
if(plane == Plane_3()) {
// CGAL_assertion(last-first >= 3);
/* we need at least 3 points to discover the original plane */
ForwardIterator p(first);
Point_3 p0(*(p++)), p1(*(p++)), p2(*(p++));
plane = Plane_3(p0, p1, p2);
/* since we just need to project the points to a non-perpendicular plane
we don't need to care about the plane orientation */
Vector_3 hv;
normal_vector_newell_3( first, last, hv);
plane = Plane_3( *first, Direction_3(hv));
}
CGAL_assertion(!plane.is_degenerate());
Direction_3 pd(plane.orthogonal_vector()), pyz(1,0,0), pxz(0,1,0);
if(pd == pyz || pd == -pyz) {
/* the plane is parallel to the YZ plane */
typedef Project_YZ< Point_2, Point_3> Project_YZ;
typedef Iterator_project< ForwardIterator, Project_YZ> Iterator_YZ;
typedef Iterator_project< IC, Project_YZ> Iterator_YZ;
Project_YZ project;
Point_2 p = project(point);
Iterator_YZ pfirst(first), plast(last);
@ -182,7 +184,7 @@ Bounded_side bounded_side_3(ForwardIterator first,
else if(pd == pxz || pd ==- pxz) {
/* the plane is parallel to the XZ plane */
typedef Project_XZ< Point_2, Point_3> Project_XZ;
typedef Iterator_project< ForwardIterator, Project_XZ> Iterator_XZ;
typedef Iterator_project< IC, Project_XZ> Iterator_XZ;
Project_XZ project;
Point_2 p = project(point);
Iterator_XZ pfirst(first), plast(last);
@ -192,7 +194,7 @@ Bounded_side bounded_side_3(ForwardIterator first,
CGAL_assertion(cross_product(pd.vector(), Vector_3(0,0,1)) == NULL_VECTOR);
/* the plane is not perpendicular to the XY plane */
typedef Project_XY< Point_2, Point_3> Project_XY;
typedef Iterator_project< ForwardIterator, Project_XY> Iterator_XY;
typedef Iterator_project< IC, Project_XY> Iterator_XY;
Project_XY project;
Point_2 p = project(point);
Iterator_XY pfirst(first), plast(last);

View File

@ -267,17 +267,6 @@ void polyhedron_3_to_nef_3(Polyhedron_& P, SNC_structure& S)
//CGAL::OGL::start_viewer();
#endif
// create nef faces from spheres
TRACEN("pair_up_halfedges()...");
C.pair_up_halfedges();
TRACEN("link_shalfedges()...");
C.link_shalfedges_to_facet_cycles();
TRACEN("categorize_facet_cycles_and_create_facets()...");
C.categorize_facet_cycles_and_create_facets();
TRACEN("create_volumes()");
C.create_volumes();
TRACEN("done");
return;
}

View File

@ -219,7 +219,8 @@ protected:
CGAL_NTS abs(y) ==
CGAL_NTS abs(z) == 1);
Vertex_handle v = ews().new_vertex();
ews().point(v) = Point_3(x, y, z);
ews().point(v) = Point_3(x*IMMN, y*IMMN, z*IMMN);
// TODO: to replace the IMMN constant by a real infimaximal number
SM_decorator D(v);
Sphere_point sp[] = { Sphere_point(-x, 0, 0),
Sphere_point(0, -y, 0),
@ -311,8 +312,10 @@ public:
typedef Polyhedron_3< Kernel> Polyhedron;
Nef_polyhedron_3( Polyhedron& P) {
initialize_simple_cube_vertices(EMPTY);
polyhedron_3_to_nef_3< Polyhedron, SNC_structure, SNC_constructor>
( P, ews() );
build_external_structure();
}
void dump() { SNC_io_parser::dump( ews()); }
@ -663,6 +666,7 @@ template <typename T>
Nef_polyhedron_3<T>::
Nef_polyhedron_3(const Plane_3& h, Boundary b) : Base(Nef_rep()) {
TRACEN("construction from plane "<<h);
initialize_simple_cube_vertices(space);
//add_box_corners(h, b);
build_external_structure();
}