diff --git a/Packages/Nef_3/TODO b/Packages/Nef_3/TODO index 9706a660445..316f87b6d4d 100644 --- a/Packages/Nef_3/TODO +++ b/Packages/Nef_3/TODO @@ -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) diff --git a/Packages/Nef_3/include/CGAL/Nef_3/SNC_constructor.h b/Packages/Nef_3/include/CGAL/Nef_3/SNC_constructor.h index 28d21f3e3f1..4ecb18155a9 100644 --- a/Packages/Nef_3/include/CGAL/Nef_3/SNC_constructor.h +++ b/Packages/Nef_3/include/CGAL/Nef_3/SNC_constructor.h @@ -47,7 +47,7 @@ #define _DEBUG 43 #include -#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 = "< point in facet? "<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 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 #"<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)"); } } } diff --git a/Packages/Nef_3/include/CGAL/Nef_3/SNC_decorator.h b/Packages/Nef_3/include/CGAL/Nef_3/SNC_decorator.h index e03abd305b7..a0a6f9ccce7 100644 --- a/Packages/Nef_3/include/CGAL/Nef_3/SNC_decorator.h +++ b/Packages/Nef_3/include/CGAL/Nef_3/SNC_decorator.h @@ -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(); } diff --git a/Packages/Nef_3/include/CGAL/Nef_3/SNC_structure.h b/Packages/Nef_3/include/CGAL/Nef_3/SNC_structure.h index 374a85c9f0b..db80285bcf8 100644 --- a/Packages/Nef_3/include/CGAL/Nef_3/SNC_structure.h +++ b/Packages/Nef_3/include/CGAL/Nef_3/SNC_structure.h @@ -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; diff --git a/Packages/Nef_3/include/CGAL/Nef_3/bounded_side_3.h b/Packages/Nef_3/include/CGAL/Nef_3/bounded_side_3.h index d195faf75b5..38bceb33ed7 100644 --- a/Packages/Nef_3/include/CGAL/Nef_3/bounded_side_3.h +++ b/Packages/Nef_3/include/CGAL/Nef_3/bounded_side_3.h @@ -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: "<; + 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 -Bounded_side bounded_side_3(ForwardIterator first, - ForwardIterator last, +template +Bounded_side bounded_side_3(IC first, + IC last, const Point_3& point, Plane_3 plane = Plane_3()) { @@ -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); diff --git a/Packages/Nef_3/include/CGAL/Nef_3/polyhedron_3_to_nef_3.h b/Packages/Nef_3/include/CGAL/Nef_3/polyhedron_3_to_nef_3.h index f032b096d33..5913f6d3b78 100644 --- a/Packages/Nef_3/include/CGAL/Nef_3/polyhedron_3_to_nef_3.h +++ b/Packages/Nef_3/include/CGAL/Nef_3/polyhedron_3_to_nef_3.h @@ -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; } diff --git a/Packages/Nef_3/include/CGAL/Nef_polyhedron_3.h b/Packages/Nef_3/include/CGAL/Nef_polyhedron_3.h index 2a504f23b81..150d7ab8d26 100644 --- a/Packages/Nef_3/include/CGAL/Nef_polyhedron_3.h +++ b/Packages/Nef_3/include/CGAL/Nef_polyhedron_3.h @@ -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 Nef_polyhedron_3:: Nef_polyhedron_3(const Plane_3& h, Boundary b) : Base(Nef_rep()) { TRACEN("construction from plane "<