diff --git a/Skin_surface_3/examples/Skin_surface_3/NGHK_skin_surface_simple.cpp b/Skin_surface_3/examples/Skin_surface_3/NGHK_skin_surface_simple.cpp index 5c64d3f1682..ca95e1af7da 100644 --- a/Skin_surface_3/examples/Skin_surface_3/NGHK_skin_surface_simple.cpp +++ b/Skin_surface_3/examples/Skin_surface_3/NGHK_skin_surface_simple.cpp @@ -1,21 +1,25 @@ -// examples/Skin_surface_3/NGHK_skin_surface_simple.C +//#define CGAL_PROFILE + +// examples/Skin_surface_3/skin_surface_simple.C #include #include #include #include -#include +#include #include +#include #include "skin_surface_writer.h" -typedef CGAL::Exact_predicates_inexact_constructions_kernel K; -typedef CGAL::Regular_triangulation_euclidean_traits_3 Traits; -typedef CGAL::Skin_surface_3 Skin_surface_3; -typedef Skin_surface_3::RT RT; -typedef Skin_surface_3::Weighted_point Weighted_point; -typedef Weighted_point::Point Bare_point; -typedef CGAL::Skin_surface_polyhedral_items_3 Poly_items; -typedef CGAL::Polyhedron_3 Polyhedron; +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Mixed_complex_traits_3 Traits; +typedef CGAL::Skin_surface_3 Skin_surface_3; +typedef Skin_surface_3::RT RT; +typedef Skin_surface_3::Weighted_point Weighted_point; +typedef Weighted_point::Point Bare_point; +typedef CGAL::Polyhedron_3< + CGAL::Simple_cartesian, + CGAL::Skin_surface_polyhedral_items_3 > Polyhedron; int main(int argc, char *argv[]) { if (argc < 2) { @@ -28,14 +32,14 @@ int main(int argc, char *argv[]) { Weighted_point wp; std::ifstream in(argv[1]); - while (in >> wp) l.push_front(wp); + while (in >> wp) l.push_back(wp); - Skin_surface_3 skin_surface(l.begin(), l.end(), shrinkfactor, false); + Skin_surface_3 skin_surface(l.begin(), l.end(), shrinkfactor); Polyhedron p; CGAL::mesh_skin_surface_3(skin_surface, p); - CGAL::subdivide_skin_surface_mesh_3(p, skin_surface,1); +// CGAL::subdivide_skin_surface_mesh_3(p, skin_surface,1); std::cout << "Is closed: " << (p.is_closed() ? "Yes" : "No") << std::endl; diff --git a/Skin_surface_3/examples/Skin_surface_3/makefile b/Skin_surface_3/examples/Skin_surface_3/makefile index 6d10854c7da..3dc7593c80f 100644 --- a/Skin_surface_3/examples/Skin_surface_3/makefile +++ b/Skin_surface_3/examples/Skin_surface_3/makefile @@ -86,6 +86,7 @@ clean: \ skin_surface_with_surface_mesher.clean \ NGHK_skin_surface_simple.clean \ NGHK_skin_surface_subdiv.clean + rm -f $(OBJ) #---------------------------------------------------------------------# # suffix rules diff --git a/Skin_surface_3/examples/Skin_surface_3/skin_surface_simple.cpp b/Skin_surface_3/examples/Skin_surface_3/skin_surface_simple.cpp index 216210014a1..17e2792942e 100644 --- a/Skin_surface_3/examples/Skin_surface_3/skin_surface_simple.cpp +++ b/Skin_surface_3/examples/Skin_surface_3/skin_surface_simple.cpp @@ -1,5 +1,3 @@ -#define CGAL_PROFILE - // examples/Skin_surface_3/skin_surface_simple.C #include #include diff --git a/Skin_surface_3/include/CGAL/Marching_tetrahedra_traits_skin_surface_3.h b/Skin_surface_3/include/CGAL/Marching_tetrahedra_traits_skin_surface_3.h index dc7cb5a8f10..b376f503388 100644 --- a/Skin_surface_3/include/CGAL/Marching_tetrahedra_traits_skin_surface_3.h +++ b/Skin_surface_3/include/CGAL/Marching_tetrahedra_traits_skin_surface_3.h @@ -49,7 +49,7 @@ public: HDS_point intersection(Cell_iterator const ch, int i, int j) const { // Precondition: ch is not an infinite cell: their surface is not set Skin_point p; - ss_3.intersect(ch->vertex(i), ch->vertex(j), p); + ss_3.intersect(ch->vertex(i), ch->vertex(j), ch->mixed_cell(), p); return Cartesian_converter()(p); diff --git a/Skin_surface_3/include/CGAL/Skin_surface_3.h b/Skin_surface_3/include/CGAL/Skin_surface_3.h index 628fb6edf8e..004103d6dec 100644 --- a/Skin_surface_3/include/CGAL/Skin_surface_3.h +++ b/Skin_surface_3/include/CGAL/Skin_surface_3.h @@ -353,19 +353,29 @@ public: return construct_surface(vh->first, K()).sign(p); } + // Trivial caching: check wether the surface is the same as the previous: + mutable Skin_surface_quadratic_surface_3< + Simple_cartesian > previous_sign_surface; + mutable Simplex previous_sign_simplex; Sign sign(const Simplex &sim, const Bare_point &p) const { + if (previous_sign_simplex != sim) { + previous_sign_simplex = sim; + previous_sign_surface = + construct_surface(sim, + Simple_cartesian()); + } try { - CGAL_PROFILER(std::string("NGHK: calls to : ") + std::string(CGAL_PRETTY_FUNCTION)); + CGAL_PROFILER(std::string("NGHK: calls to : ") + + std::string(CGAL_PRETTY_FUNCTION)); Protect_FPU_rounding P; - Sign result = construct_surface - (sim, - Exact_predicates_inexact_constructions_kernel()).sign(p); + Sign result = previous_sign_surface.sign(p); if (! is_indeterminate(result)) return result; } catch (Interval_nt_advanced::unsafe_comparison) {} - CGAL_PROFILER(std::string("NGHK: failures of : ") + std::string(CGAL_PRETTY_FUNCTION)); + CGAL_PROFILER(std::string("NGHK: failures of : ") + + std::string(CGAL_PRETTY_FUNCTION)); Protect_FPU_rounding P(CGAL_FE_TONEAREST); return construct_surface (sim, @@ -400,6 +410,15 @@ public: Simplex s2 = vh2->first; intersect(p1,p2, s1,s2, p); } + void intersect(const CMCT_Vertex_handle vh1, + const CMCT_Vertex_handle vh2, + const Simplex &s, + Bare_point &p) const { + Bare_point p1 = mc_triangulator->location(vh1, gt); + Bare_point p2 = mc_triangulator->location(vh2, gt); + Simplex sp = s; + intersect(p1,p2, sp,sp, p); + } void intersect(Bare_point &p1, Bare_point &p2, Simplex &s1, Simplex &s2, @@ -432,20 +451,15 @@ public: p = midpoint(p1, p2); } - template - void intersect_with_transversal_segment(Point &p) const { - typedef typename Point::R Traits; - typedef typename Traits::Plane_3 Plane; - typedef typename Traits::Line_3 Line; - Mixed_complex_traits_3 point_traits(gt.get_shrink()); - Cartesian_converter converter; + void intersect_with_transversal_segment(Bare_point &p) const { + typedef typename Geometric_traits::Kernel::Plane_3 Plane; + typedef typename Geometric_traits::Kernel::Line_3 Line; - Simplex sim = locate_mixed(converter(p)); + Simplex sim = locate_mixed(p); CMCT_Cell tet = locate_tet(p, sim); // get transversal segment: - Point p1, p2; + Bare_point p1, p2; // Compute signs on vertices and sort them: int nIn = 0; @@ -460,12 +474,12 @@ public: Object obj; if (nIn==1) { - p1 = mc_triangulator->location(tet.vertex(sortedV[0]), point_traits); + p1 = mc_triangulator->location(tet.vertex(sortedV[0]), gt); obj = CGAL::intersection( Plane - (mc_triangulator->location(tet.vertex(sortedV[1]), point_traits), - mc_triangulator->location(tet.vertex(sortedV[2]), point_traits), - mc_triangulator->location(tet.vertex(sortedV[3]), point_traits)), + (mc_triangulator->location(tet.vertex(sortedV[1]), gt), + mc_triangulator->location(tet.vertex(sortedV[2]), gt), + mc_triangulator->location(tet.vertex(sortedV[3]), gt)), Line(p1, p)); if ( !assign(p2, obj) ) { CGAL_assertion_msg(false,"intersection: no intersection."); @@ -473,33 +487,33 @@ public: } else if (nIn==2) { obj = CGAL::intersection( Plane - (mc_triangulator->location(tet.vertex(sortedV[2]), point_traits), - mc_triangulator->location(tet.vertex(sortedV[3]), point_traits), + (mc_triangulator->location(tet.vertex(sortedV[2]), gt), + mc_triangulator->location(tet.vertex(sortedV[3]), gt), p), Line( - mc_triangulator->location(tet.vertex(sortedV[0]), point_traits), - mc_triangulator->location(tet.vertex(sortedV[1]), point_traits))); + mc_triangulator->location(tet.vertex(sortedV[0]), gt), + mc_triangulator->location(tet.vertex(sortedV[1]), gt))); if ( !assign(p1, obj) ) { CGAL_assertion_msg(false,"intersection: no intersection."); } obj = CGAL::intersection( Plane - (mc_triangulator->location(tet.vertex(sortedV[0]), point_traits), - mc_triangulator->location(tet.vertex(sortedV[1]), point_traits), + (mc_triangulator->location(tet.vertex(sortedV[0]), gt), + mc_triangulator->location(tet.vertex(sortedV[1]), gt), p), Line( - mc_triangulator->location(tet.vertex(sortedV[2]), point_traits), - mc_triangulator->location(tet.vertex(sortedV[3]), point_traits))); + mc_triangulator->location(tet.vertex(sortedV[2]), gt), + mc_triangulator->location(tet.vertex(sortedV[3]), gt))); if ( !assign(p2, obj) ) { CGAL_assertion_msg(false,"intersection: no intersection."); } } else if (nIn==3) { - p2 = mc_triangulator->location(tet.vertex(sortedV[3]), point_traits); + p2 = mc_triangulator->location(tet.vertex(sortedV[3]), gt); obj = CGAL::intersection( Plane( - mc_triangulator->location(tet.vertex(sortedV[0]), point_traits), - mc_triangulator->location(tet.vertex(sortedV[1]), point_traits), - mc_triangulator->location(tet.vertex(sortedV[2]), point_traits)), + mc_triangulator->location(tet.vertex(sortedV[0]), gt), + mc_triangulator->location(tet.vertex(sortedV[1]), gt), + mc_triangulator->location(tet.vertex(sortedV[2]), gt)), Line(p2, p)); if ( !assign(p1, obj) ) { CGAL_assertion_msg(false,"intersection: no intersection."); diff --git a/Skin_surface_3/include/CGAL/Skin_surface_polyhedral_items_3.h b/Skin_surface_3/include/CGAL/Skin_surface_polyhedral_items_3.h index f78e7a9042d..d95a8f5517d 100644 --- a/Skin_surface_3/include/CGAL/Skin_surface_polyhedral_items_3.h +++ b/Skin_surface_3/include/CGAL/Skin_surface_polyhedral_items_3.h @@ -28,10 +28,9 @@ CGAL_BEGIN_NAMESPACE template struct Skin_Surface_polyhedral_face : public CGAL::HalfedgeDS_face_base { typedef SkinSurface3 Skin_surface; -// typedef typename Skin_surface::Triangulated_mixed_complex TMC; -// typedef typename TMC::Cell_handle Triang_Cell_handle; + typedef typename SkinSurface3::Simplex Simplex; -// Triang_Cell_handle triang_ch; + Simplex sim; }; template < class SkinSurface3 > diff --git a/Skin_surface_3/include/CGAL/Skin_surface_refinement_traits_3.h b/Skin_surface_3/include/CGAL/Skin_surface_refinement_traits_3.h index e5ce4fd9471..3a34c57ac25 100644 --- a/Skin_surface_3/include/CGAL/Skin_surface_refinement_traits_3.h +++ b/Skin_surface_3/include/CGAL/Skin_surface_refinement_traits_3.h @@ -48,13 +48,24 @@ public: } P_point to_surface(P_vertex_handle vh) { - P_point result = vh->point(); + typename Skin_surface::Bare_point result = + Cartesian_converter()(vh->point()); ss_3.intersect_with_transversal_segment(result); - return result; + return + Cartesian_converter + ()( result ); } P_vector normal(P_vertex_handle vh) { - return ss_3.normal(vh->point()); + // Convert to and from the skin surface kernel + return + Cartesian_converter + ()( ss_3.normal + (Cartesian_converter() + (vh->point()))); } protected: diff --git a/Skin_surface_3/include/CGAL/triangulate_mixed_complex_3.h b/Skin_surface_3/include/CGAL/triangulate_mixed_complex_3.h index c50ba657e9e..99d923e7ce4 100644 --- a/Skin_surface_3/include/CGAL/triangulate_mixed_complex_3.h +++ b/Skin_surface_3/include/CGAL/triangulate_mixed_complex_3.h @@ -36,17 +36,19 @@ template < class CMCT > class CMCT_Cell { public: typedef CMCT_Cell Self; - typedef typename CMCT::Vertex_handle Vertex_handle; + typedef typename CMCT::Vertex_handle Vertex_handle; + typedef typename CMCT::Rt_Simplex Rt_Simplex; CMCT_Cell() { } - CMCT_Cell(Vertex_handle vs[]) { + CMCT_Cell(Vertex_handle vs[], Rt_Simplex &sim) : sim_(sim) { for (int i=0; i<4; i++) _vs[i] = vs[i]; } CMCT_Cell(const Vertex_handle &vh0, - const Vertex_handle &vh1, - const Vertex_handle &vh2, - const Vertex_handle &vh3) + const Vertex_handle &vh1, + const Vertex_handle &vh2, + const Vertex_handle &vh3, + Rt_Simplex &sim) : sim_(sim) { _vs[0] = vh0; _vs[1] = vh1; @@ -61,6 +63,9 @@ public: CGAL_assertion((0 <= i) && (i<4)); return _vs[i]; } + const Rt_Simplex &mixed_cell() const { + return sim_; + } bool operator==(const Self &other) const { return ((_vs[0] == other._vs[0]) && (_vs[1] == other._vs[1]) && @@ -70,6 +75,7 @@ public: bool operator!=(const Self &other) const { return !operator==(other); } private: Vertex_handle _vs[4]; + Rt_Simplex sim_; }; template < class CMCT > @@ -501,6 +507,7 @@ private: template void add_cell(Vertex_handle vh[], int orient, + Rt_Simplex &simplex, OutputIteratorCells cells) const; Vertex_handle get_vertex(Rt_Simplex &sDel, Rt_Simplex &sVor) const; @@ -941,7 +948,8 @@ construct_0_cell(Rt_Vertex_handle rt_vh, OutputIteratorCells cells) const CGAL_assertion(sVor_e != sVor_f); CGAL_assertion(sVor_f != sVor_c); - add_cell(vh,(index + (j==(i%3+1)? 1:0))&1, cells); + add_cell(vh,(index + (j==(i%3+1)? 1:0))&1, + simplex, cells); } } } @@ -998,15 +1006,15 @@ construct_1_cell(const Rt_Edge &e, orient = (index1 + (fi==1))&1; } // vh: dimension are (01,11,12,13) - add_cell(vh,orient,cells); + add_cell(vh,orient,mixed_cell_simplex,cells); vh[1] = get_vertex(sDel_v, sVor_f); // vh: dimension are (01,02,12,13) - add_cell(vh,1-orient,cells); + add_cell(vh,1-orient,mixed_cell_simplex,cells); vh[2] = get_vertex(sDel_v, sVor_c); // vh: dimension are (01,02,03,13) - add_cell(vh,orient,cells); + add_cell(vh,orient,mixed_cell_simplex,cells); } } } @@ -1071,13 +1079,13 @@ construct_2_cell(const Rt_Facet &f, orient = (index1 + (((4+index0-fi)&3)==1))&1; } - add_cell(vh,orient,cells); + add_cell(vh,orient,simplex,cells); vh[2] = get_vertex(sDel_e, sVor_c); - add_cell(vh,1-orient,cells); + add_cell(vh,1-orient,simplex,cells); vh[1] = get_vertex(sDel_v, sVor_c); - add_cell(vh,orient,cells); + add_cell(vh,orient,simplex,cells); } } } @@ -1135,7 +1143,7 @@ construct_3_cell(Rt_Cell_handle rt_ch, } else { orient = (index1 + (vi==3))&1; } - add_cell(vh, orient, cells); + add_cell(vh, orient, simplex, cells); } } } @@ -1181,7 +1189,9 @@ template template void Combinatorial_mixed_complex_triangulator_3:: -add_cell(Vertex_handle vh[], int orient, OutputIteratorCells cells) const { +add_cell(Vertex_handle vh[], int orient, + Rt_Simplex &simplex, + OutputIteratorCells cells) const { assert((orient==0) || (orient==1)); assert(*vh[0] != Vertex()); assert(*vh[1] != Vertex()); assert(*vh[2] != Vertex()); assert(*vh[3] != Vertex()); @@ -1189,9 +1199,9 @@ add_cell(Vertex_handle vh[], int orient, OutputIteratorCells cells) const { assert(*vh[1] != *vh[2]); assert(*vh[1] != *vh[3]); assert(*vh[2] != *vh[3]); if (orient) { - *cells++ = Cell(vh[0], vh[1], vh[2], vh[3]); + *cells++ = Cell(vh[0], vh[1], vh[2], vh[3], simplex); } else { - *cells++ = Cell(vh[0], vh[1], vh[3], vh[2]); + *cells++ = Cell(vh[0], vh[1], vh[3], vh[2], simplex); } }