From 8f4555ea5d23edc3eedf0d73f060d6ca58bfc462 Mon Sep 17 00:00:00 2001 From: Nico Kruithof Date: Fri, 22 Dec 2006 13:07:26 +0000 Subject: [PATCH] Added triangulated_mixed_complex back. Next, do the locate in the triangulated mixed complex --- .../NGHK_skin_surface_simple.cpp | 2 +- .../NGHK_skin_surface_subdiv.cpp | 2 +- .../examples/Skin_surface_3/makefile | 2 +- .../Skin_surface_3/skin_surface_writer.h | 6 +- Skin_surface_3/include/CGAL/Skin_surface_3.h | 187 ++++- .../Triangulated_mixed_complex_observer_3.h | 136 ++-- .../CGAL/triangulate_mixed_complex_3.h | 678 ++++++++++++++---- .../test/Skin_surface_3/surface_test.C | 88 +-- 8 files changed, 814 insertions(+), 287 deletions(-) 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 b070fbfadae..6fe578b0909 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 @@ -46,7 +46,7 @@ int main(int argc, char *argv[]) { std::cout << "Is closed: " << (p.is_closed() ? "Yes" : "No") << std::endl; std::ofstream out("mesh.off"); - // // write_polyhedron_with_normals(p, skin_surface, out); + //write_polyhedron_with_normals(p, skin_surface, out); out << p; return 0; diff --git a/Skin_surface_3/examples/Skin_surface_3/NGHK_skin_surface_subdiv.cpp b/Skin_surface_3/examples/Skin_surface_3/NGHK_skin_surface_subdiv.cpp index 69f7be8cd76..e6488cc2806 100644 --- a/Skin_surface_3/examples/Skin_surface_3/NGHK_skin_surface_subdiv.cpp +++ b/Skin_surface_3/examples/Skin_surface_3/NGHK_skin_surface_subdiv.cpp @@ -1,5 +1,5 @@ // examples/Skin_surface_3/NGHK_skin_surface_subdiv.C -//#define CGAL_PROFILE +#define CGAL_PROFILE //#define CGAL_NO_ASSERTIONS #include diff --git a/Skin_surface_3/examples/Skin_surface_3/makefile b/Skin_surface_3/examples/Skin_surface_3/makefile index 2589a43af1a..366e7352f52 100644 --- a/Skin_surface_3/examples/Skin_surface_3/makefile +++ b/Skin_surface_3/examples/Skin_surface_3/makefile @@ -13,7 +13,7 @@ include $(CGAL_MAKEFILE) # compiler flags #---------------------------------------------------------------------# -CXXFLAGS = -g \ +CXXFLAGS = \ -I.\ -Idsrpdb/include\ -I../../include \ diff --git a/Skin_surface_3/examples/Skin_surface_3/skin_surface_writer.h b/Skin_surface_3/examples/Skin_surface_3/skin_surface_writer.h index c30d475853c..b8298fd18e8 100644 --- a/Skin_surface_3/examples/Skin_surface_3/skin_surface_writer.h +++ b/Skin_surface_3/examples/Skin_surface_3/skin_surface_writer.h @@ -15,8 +15,8 @@ void write_polyhedron_with_normals(SkinSurface &skin, typedef typename Polyhedron::Facet_iterator Facet_iterator; typedef typename Polyhedron::Halfedge_around_facet_circulator HFC; typedef typename Polyhedron::Vertex_handle Vertex_handle; - typedef typename Polyhedron::Traits::Vector_3 Vector_3; - + typedef typename Polyhedron::Traits::Vector_3 Vector; + //typedef typename SkinSurface::Vector Vector; // Write header out << "NOFF " << p.size_of_vertices () << " " << p.size_of_facets() @@ -29,7 +29,7 @@ void write_polyhedron_with_normals(SkinSurface &skin, // Subdivision_policy *policy = get_subdivision_policy(p, skin); for (Vertex_iterator vit = p.vertices_begin(); vit != p.vertices_end(); vit ++) { - Vector_3 n = skin.normal(vit->point()); + Vector n = skin.normal(vit->point()); n = n/sqrt(n*n); out << vit->point() << " " << n diff --git a/Skin_surface_3/include/CGAL/Skin_surface_3.h b/Skin_surface_3/include/CGAL/Skin_surface_3.h index e138c35f02a..fb38583f799 100644 --- a/Skin_surface_3/include/CGAL/Skin_surface_3.h +++ b/Skin_surface_3/include/CGAL/Skin_surface_3.h @@ -105,6 +105,7 @@ private: typedef typename TMC::Finite_cells_iterator TMC_Cell_iterator; typedef typename TMC::Vertex_handle TMC_Vertex_handle; typedef typename TMC::Cell_handle TMC_Cell_handle; + typedef typename TMC::Point TMC_Point; public: template < class WP_iterator > Skin_surface_3(WP_iterator begin, WP_iterator end, @@ -170,8 +171,13 @@ public: template void subdivide_skin_surface_mesh_3(Polyhedron_3 &p) const; - Sign sign(const Bare_point &p, const Simplex &start = Simplex()) const { - return get_sign(locate_mixed(p,start), p); + Sign sign(const Bare_point &p, + const TMC_Cell_handle start = TMC_Cell_handle()) const { + if (start == TMC_Cell_handle()) { + return sign(locate_mixed(p,start), p); + } else { + return sign(start, p); + } } Sign sign(TMC_Vertex_handle vit) const { CGAL_assertion(!tmc.is_infinite(vit)); @@ -187,7 +193,7 @@ public: } CGAL_assertion(!tmc.is_infinite(ch)); - // don't use get_sign, since the point is constructed: + // don't use sign, since the point is constructed: try { CGAL_PROFILER(std::string("NGHK: calls to : ") + @@ -210,8 +216,13 @@ public: EK() ).sign(p_exact); } Vector - normal(const Bare_point &p, const Simplex &start = Simplex()) const { - return get_normal(locate_mixed(p,start), p); + normal(const Bare_point &p, + const TMC_Cell_handle start = TMC_Cell_handle()) const { + if (start == TMC_Cell_handle()) { + return get_normal(locate_mixed(p,start), p); + } else { + return get_normal(start, p); + } } template @@ -361,7 +372,7 @@ private: public: TMC_Cell_handle locate_mixed(const Bare_point &p, - const TMC_Cell_handle &start = TMC_Cell_handle()) const; + TMC_Cell_handle start = TMC_Cell_handle()) const; // exact computation of the sign on a vertex of the TMC // Sign sign(const CMCT_Vertex_handle vh) const { @@ -427,9 +438,37 @@ public: ).value(p1)); } FT + less(Cell_info &info1, + const Bare_point &p1, + Cell_info &info2, + const Bare_point &p2) const { + try + { + CGAL_PROFILER(std::string("NGHK: calls to : ") + + std::string(CGAL_PRETTY_FUNCTION)); + Protect_FPU_rounding P; + Sign result = CGAL_NTS sign(info2.second->value(p2) - + info1.second->value(p1)); + if (! is_indeterminate(result)) + return result==POSITIVE; + } + catch (Interval_nt_advanced::unsafe_comparison) {} + CGAL_PROFILER(std::string("NGHK: failures of : ") + + std::string(CGAL_PRETTY_FUNCTION)); + Protect_FPU_rounding P(CGAL_FE_TONEAREST); + + return + CGAL_NTS sign(construct_surface(info2.first, + Exact_predicates_exact_constructions_kernel() + ).value(p2) - + construct_surface(info1.first, + Exact_predicates_exact_constructions_kernel() + ).value(p1)); + } + FT value(const Bare_point &p) const { - Simplex sim = locate_mixed(p); - return value(sim,p); + TMC_Cell_handle ch = locate_mixed(p); + return value(Simplex(ch),p); } FT @@ -444,31 +483,36 @@ public: return value(ch->info(), p); } Vector - get_normal(const Simplex &mc, const Bare_point &p) const { - return construct_surface(mc).gradient(p); + get_normal(TMC_Cell_handle ch, const Bare_point &p) const { + CGAL_assertion(!tmc.is_infinite(ch)); + + return ch->info().second->gradient(p); } // Move the point in the direction of the gradient void to_surface(Bare_point &p, - const Simplex &start = Simplex()) const { + const TMC_Cell_handle &start = TMC_Cell_handle()) const { Bare_point p1 = p; - Simplex s1 = locate_mixed(p,start); - Sign sign1 = get_sign(s1, p1); + TMC_Cell_handle ch1 = start; + if (start != TMC_Cell_handle()) { + ch1 = locate_mixed(p,ch1); + } + Sign sign1 = sign(ch1, p1); - Vector n = get_normal(s1,p); - if (sign1 == POSITIVE) n = -value(s1,p)*n; + Vector n = get_normal(ch1,p); + if (sign1 == POSITIVE) n = -value(ch1,p)*n; int k=2; Bare_point p2 = p+k*n; - Simplex s2 = locate_mixed(p2, s1); - while (get_sign(s2,p2) == sign1) { + TMC_Cell_handle ch2 = locate_mixed(p2, ch1); + while (sign(ch2,p2) == sign1) { k++; p1 = p2; - s1 = s2; + ch1 = ch2; p2 = p+k*n; - s2 = locate_mixed(p2, s2); + ch2 = locate_mixed(p2, ch2); } - intersect(p1,p2, s1,s2, p); + intersect(p1,p2, ch1,ch2, p); } // void intersect(const CMCT_Vertex_handle vh1, // const CMCT_Vertex_handle vh2, @@ -500,21 +544,21 @@ public: FT sq_dist = squared_distance(p1,p2); // Use value to make the computation robust (endpoints near the surface) - if (value(s1, p1) > value(s2, p2)) std::swap(p1, p2); + if (less(s2->info(), p2, s1->info(), p1)) std::swap(p1, p2); TMC_Cell_handle sp = s1; while ((s1 != s2) && (sq_dist > 1e-8)) { p = midpoint(p1, p2); sp = locate_mixed(converter(p), sp); - if (get_sign(sp, p) == NEGATIVE) { p1 = p; s1 = sp; } + if (sign(sp, p) == NEGATIVE) { p1 = p; s1 = sp; } else { p2 = p; s2 = sp; } sq_dist *= .25; } while (sq_dist > 1e-8) { p = midpoint(p1, p2); - if (get_sign(s1, p) == NEGATIVE) { p1 = p; } + if (sign(s1, p) == NEGATIVE) { p1 = p; } else { p2 = p; } sq_dist *= .25; } @@ -581,7 +625,7 @@ public: if (nIn==1) { p1 = tet_pts[sortedV[0]]; obj = CGAL::intersection(Plane(tet_pts[sortedV[1]], - tet_pts[sortedV[3]], + tet_pts[sortedV[2]], tet_pts[sortedV[3]]), Line(p1, p)); if ( !assign(p2, obj) ) { @@ -614,6 +658,7 @@ public: CGAL_assertion_msg(false,"intersection: no intersection."); } } else { + std::cout << "nIn == " << nIn << std::endl; CGAL_assertion(false); } @@ -706,6 +751,9 @@ public: FT shrink_factor() const { return gt.get_shrink(); } + const TMC &triangulated_mixed_complex() const { + return tmc; + } private: void construct_bounding_box(Regular ®ular); @@ -772,12 +820,86 @@ construct_bounding_box(Regular ®ular) template typename Skin_surface_3::TMC_Cell_handle Skin_surface_3:: -locate_mixed(const Bare_point &p, - const TMC_Cell_handle &start) const { - Cartesian_converter converter; - - // NGHK: add a try ... catch? - return tmc.locate(converter(p)); +locate_mixed(const Bare_point &p0, + TMC_Cell_handle start) const { + typedef Exact_predicates_exact_constructions_kernel EK; + Cartesian_converter converter; + Skin_surface_traits_3 exact_traits(shrink_factor()); + + typename EK::Point_3 p = converter(p0); + + Protect_FPU_rounding P(CGAL_FE_TONEAREST); + + typename EK::Point_3 e_pts[4]; + const typename EK::Point_3 *pts[4]; + + // Make sure we continue from here with a finite cell. + if ( start == TMC_Cell_handle() ) + start = tmc.infinite_cell(); + + int ind_inf; + if (start->has_vertex(tmc.infinite_vertex(), ind_inf) ) + start = start->neighbor(ind_inf); + + CGAL_triangulation_precondition(start != TMC_Cell_handle()); + CGAL_triangulation_precondition(!start->has_vertex(tmc.infinite_vertex())); + + // We implement the remembering visibility/stochastic walk. + + // Remembers the previous cell to avoid useless orientation tests. + TMC_Cell_handle previous = TMC_Cell_handle(); + TMC_Cell_handle c = start; + + // Stores the results of the 4 orientation tests. It will be used + // at the end to decide if p lies on a face/edge/vertex/interior. + Orientation o[4]; + + // Now treat the cell c. + try_next_cell: + + // We know that the 4 vertices of c are positively oriented. + // So, in order to test if p is seen outside from one of c's facets, + // we just replace the corresponding point by p in the orientation + // test. We do this using the array below. + for (int j=0; j<4; j++) { + e_pts[j] = get_anchor_point(c->vertex(j)->info(), exact_traits); + pts[j] = &e_pts[j]; + } + + // For the remembering stochastic walk, + // we need to start trying with a random index : + int i = rng.template get_bits<2>(); + // For the remembering visibility walk (Delaunay only), we don't : + // int i = 0; + + for (int j=0; j != 4; ++j, i = (i+1)&3) { + TMC_Cell_handle next = c->neighbor(i); + if (previous == next) { + o[i] = POSITIVE; + continue; + } + // We temporarily put p at i's place in pts. + const typename EK::Point_3* backup = pts[i]; + pts[i] = &p; + try { + o[i] = orientation(*pts[0], *pts[1], *pts[2], *pts[3]); + } catch (Interval_nt_advanced::unsafe_comparison) { + } + if ( o[i] != NEGATIVE ) { + pts[i] = backup; + continue; + } + if ( next->has_vertex(tmc.infinite_vertex()) ) { + // We are outside the convex hull. + return next; + } + previous = c; + c = next; + goto try_next_cell; + } + + return c; +} // Simplex prev, s; // if (start == Simplex()) { @@ -957,15 +1079,12 @@ locate_mixed(const Bare_point &p, // // std::cout << "]"; // return s; -} +// } template template void Skin_surface_3::mesh_skin_surface_3(Polyhedron_3 &p) const { - std::cout << "Mesh_Skin_Surface_3" << std::endl; - std::cout << " TODO" << std::endl; - typedef Polyhedron_3 Polyhedron; typedef Marching_tetrahedra_traits_skin_surface_3< diff --git a/Skin_surface_3/include/CGAL/Triangulated_mixed_complex_observer_3.h b/Skin_surface_3/include/CGAL/Triangulated_mixed_complex_observer_3.h index 55158502050..b493bb23edb 100644 --- a/Skin_surface_3/include/CGAL/Triangulated_mixed_complex_observer_3.h +++ b/Skin_surface_3/include/CGAL/Triangulated_mixed_complex_observer_3.h @@ -95,78 +95,76 @@ public: Rt_Cell_handle ch; switch (s.dimension()) { - case 0: - { - vh = s; - create_sphere(r2s_converter(vh->point().point()), - -r2s_converter(vh->point().weight()), - shrink, - 1); - break; - } - case 1: - { - e = s; - Surface_weighted_point p0 = - r2s_converter(e.first->vertex(e.second)->point()); - Surface_weighted_point p1 = - r2s_converter(e.first->vertex(e.third)->point()); - - create_hyperboloid - (typename Surface_regular_traits:: - Construct_weighted_circumcenter_3()(p0,p1), - typename Surface_regular_traits:: - Compute_squared_radius_smallest_orthogonal_sphere_3()(p0,p1), - p0 - p1, - shrink, - 1); - break; - } - case 2: - { - f = s; - Surface_weighted_point p0 = - r2s_converter(f.first->vertex((f.second+1)&3)->point()); - Surface_weighted_point p1 = - r2s_converter(f.first->vertex((f.second+2)&3)->point()); - Surface_weighted_point p2 = - r2s_converter(f.first->vertex((f.second+3)&3)->point()); - - create_hyperboloid - (typename Surface_regular_traits:: - Construct_weighted_circumcenter_3()(p0,p1,p2), - typename Surface_regular_traits:: - Compute_squared_radius_smallest_orthogonal_sphere_3()(p0,p1,p2), - typename Surface_regular_traits:: - Construct_orthogonal_vector_3()(p0,p1,p2), - 1-shrink, - -1); - break; - } - case 3: - { - ch = s; - create_sphere - (typename Surface_regular_traits:: - Construct_weighted_circumcenter_3() - (r2s_converter(ch->vertex(0)->point()), - r2s_converter(ch->vertex(1)->point()), - r2s_converter(ch->vertex(2)->point()), - r2s_converter(ch->vertex(3)->point())), - typename Surface_regular_traits:: - Compute_squared_radius_smallest_orthogonal_sphere_3() - (r2s_converter(ch->vertex(0)->point()), - r2s_converter(ch->vertex(1)->point()), - r2s_converter(ch->vertex(2)->point()), - r2s_converter(ch->vertex(3)->point())), - 1-shrink, - -1); - } + case 0: + { + vh = s; + create_sphere(r2s_converter(vh->point().point()), + -r2s_converter(vh->point().weight()), + shrink, + 1); + break; + } + case 1: + { + e = s; + Surface_weighted_point p0 = + r2s_converter(e.first->vertex(e.second)->point()); + Surface_weighted_point p1 = + r2s_converter(e.first->vertex(e.third)->point()); + + create_hyperboloid + (typename Surface_regular_traits:: + Construct_weighted_circumcenter_3()(p0,p1), + typename Surface_regular_traits:: + Compute_squared_radius_smallest_orthogonal_sphere_3()(p0,p1), + p0 - p1, + shrink, + 1); + break; + } + case 2: + { + f = s; + Surface_weighted_point p0 = + r2s_converter(f.first->vertex((f.second+1)&3)->point()); + Surface_weighted_point p1 = + r2s_converter(f.first->vertex((f.second+2)&3)->point()); + Surface_weighted_point p2 = + r2s_converter(f.first->vertex((f.second+3)&3)->point()); + + create_hyperboloid + (typename Surface_regular_traits:: + Construct_weighted_circumcenter_3()(p0,p1,p2), + typename Surface_regular_traits:: + Compute_squared_radius_smallest_orthogonal_sphere_3()(p0,p1,p2), + typename Surface_regular_traits:: + Construct_orthogonal_vector_3()(p0,p1,p2), + 1-shrink, + -1); + break; + } + case 3: + { + ch = s; + create_sphere + (typename Surface_regular_traits:: + Construct_weighted_circumcenter_3() + (r2s_converter(ch->vertex(0)->point()), + r2s_converter(ch->vertex(1)->point()), + r2s_converter(ch->vertex(2)->point()), + r2s_converter(ch->vertex(3)->point())), + typename Surface_regular_traits:: + Compute_squared_radius_smallest_orthogonal_sphere_3() + (r2s_converter(ch->vertex(0)->point()), + r2s_converter(ch->vertex(1)->point()), + r2s_converter(ch->vertex(2)->point()), + r2s_converter(ch->vertex(3)->point())), + 1-shrink, + -1); + } } } - // NGHK: uncomment: ch->info() = typename SkinSurface_3::Cell_info(s,surf); - //ch->simp = s; } FT shrink; 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 fd5c9287ea9..3fe76460b5e 100644 --- a/Skin_surface_3/include/CGAL/triangulate_mixed_complex_3.h +++ b/Skin_surface_3/include/CGAL/triangulate_mixed_complex_3.h @@ -17,9 +17,10 @@ // // Author(s) : Nico Kruithof -#ifndef CGAL_TRIANGULATE_MIXED_COMPLEX_3_H -#define CGAL_TRIANGULATE_MIXED_COMPLEX_3_H +#ifndef CGAL_TRIANGULATE_MIXED_COMPLEX_3 +#define CGAL_TRIANGULATE_MIXED_COMPLEX_3 +// #include #include #include @@ -96,6 +97,8 @@ private: typedef std::pair Symb_anchor; // You might get type differences here: + // The map that maps a Rt_Simplex to an iterator of the map + // (used as union_find_structure) struct Anchor_map_iterator_tmp; typedef std::map Anchor_map; struct Anchor_map_iterator_tmp : Anchor_map::iterator { @@ -105,12 +108,15 @@ private: : Anchor_map::iterator(it) {} }; typedef typename Anchor_map::iterator Anchor_map_iterator; -public: +public: Mixed_complex_triangulator_3(Regular const ®ular, - Triangulated_mixed_complex const &triangulated_mixed_complex, + Rt_FT const &shrink, + Triangulated_mixed_complex + &triangulated_mixed_complex, bool verbose) : regular(regular), + shrink(shrink), _tmc(triangulated_mixed_complex), triangulation_incr_builder(triangulated_mixed_complex), compute_anchor_obj(regular), @@ -119,11 +125,14 @@ public: build(); } - Mixed_complex_triangulator_3(Regular const®ular, - Triangulated_mixed_complex &triangulated_mixed_complex, + Mixed_complex_triangulator_3(Regular ®ular, + Rt_FT const &shrink, + Triangulated_mixed_complex + &triangulated_mixed_complex, Triangulated_mixed_complex_observer &observer, bool verbose) : regular(regular), + shrink(shrink), _tmc(triangulated_mixed_complex), observer(observer), triangulation_incr_builder(triangulated_mixed_complex), @@ -141,27 +150,64 @@ private: if (verbose) std::cout << "Construct vertices" << std::endl; construct_vertices(); - if (verbose) std::cout << "Construct cells" << std::endl; - construct_cells(); // mixed cells corresponding to regular vertices + // mixed cells corresponding to regular vertices + if (verbose) std::cout << "Construct 0 cells" << std::endl; + for (Rt_Finite_vertices_iterator vit = regular.finite_vertices_begin(); + vit != regular.finite_vertices_end(); vit ++) { + construct_0_cell(vit); + } + + // mixed cells corresponding to regular edges + if (verbose) std::cout << "Construct 1 cells" << std::endl; + for (Rt_Finite_edges_iterator eit = regular.finite_edges_begin(); + eit != regular.finite_edges_end(); eit ++) { + construct_1_cell(eit); + } + + // mixed cells corresponding to regular facets + if (verbose) std::cout << "Construct 2 cells" << std::endl; + for (Rt_Finite_facets_iterator fit = regular.finite_facets_begin(); + fit != regular.finite_facets_end(); fit ++) { + construct_2_cell(fit); + } + + // mixed cells corresponding to regular cells + if (verbose) std::cout << "Construct 3 cells" << std::endl; + for (Rt_Finite_cells_iterator cit = regular.finite_cells_begin(); + cit != regular.finite_cells_end(); + cit++) { + construct_3_cell(cit); + } triangulation_incr_builder.end_triangulation(); anchors.clear(); - CGAL_assertion(_tmc.is_valid()); - //remove_small_edges(); - +// { // NGHK: debug code: +// CGAL_assertion(_tmc.is_valid()); +// std::vector ch_vertices; +// _tmc.incident_vertices(_tmc.infinite_vertex(), +// std::back_inserter(ch_vertices)); +// for (typename std::vector::iterator +// vit = ch_vertices.begin(); vit != ch_vertices.end(); vit++) { +// CGAL_assertion((*vit)->sign() == POSITIVE); +// } +// } } - Tmc_Vertex_handle add_vertex(Rt_Simplex const &anchor); + Tmc_Vertex_handle add_vertex(Symb_anchor const &anchor); Tmc_Cell_handle add_cell(Tmc_Vertex_handle vh[], int orient, Rt_Simplex s); - Tmc_Vertex_handle get_vertex(Rt_Simplex &sVor); + Tmc_Vertex_handle get_vertex(Rt_Simplex &sDel, Rt_Simplex &sVor); + void construct_anchor_del(Rt_Simplex const &sDel); void construct_anchor_vor(Rt_Simplex const &sVor); void construct_anchors(); + Rt_Simplex get_anchor_del(Rt_Simplex const &sDel) { + return find_anchor(anchor_del2, sDel)->first; + } Rt_Simplex get_anchor_vor(Rt_Simplex const &sVor) { return find_anchor(anchor_vor2, sVor)->first; } @@ -183,20 +229,17 @@ private: void construct_vertices(); Tmc_Point get_orthocenter(Rt_Simplex const &s); - Tmc_Point get_anchor(Rt_Simplex const &sVor); + Tmc_Point get_anchor(Rt_Simplex const &sDel, Rt_Simplex const &sVor); template - Point construct_anchor_point(const Point ¢er_vor) { - std::cout << "still union_of_balls" << std::endl; -// typename Other_MixedComplexTraits_3::Bare_point p_del = -// orthocenter(v.first, traits); -// typename Other_MixedComplexTraits_3::Bare_point p_vor = -// orthocenter(v.second, traits); - -// return traits.construct_anchor_point_3_object()(p_del, p_vor); - return center_vor; + Point construct_anchor_point(const Point ¢er_del, + const Point ¢er_vor) { + return center_del + shrink*(center_vor-center_del); } - void construct_cells(); + void construct_0_cell(Rt_Vertex_handle rt_vh); + void construct_1_cell(const Rt_Finite_edges_iterator &eit); + void construct_2_cell(const Rt_Finite_facets_iterator &fit); + void construct_3_cell(Rt_Cell_handle rt_ch); void remove_small_edges(); bool is_collapsible(Tmc_Vertex_handle vh, @@ -207,6 +250,7 @@ private: private: Regular const ®ular; + Rt_FT const &shrink; Triangulated_mixed_complex &_tmc; Triangulated_mixed_complex_observer &observer; @@ -240,8 +284,8 @@ private: Unique_hash_map < Rt_Cell_handle, Index_c4 > index_03; - Anchor_map anchor_vor2; - std::map anchors; + Anchor_map anchor_del2, anchor_vor2; + std::map anchors; }; template < @@ -255,6 +299,47 @@ const int Mixed_complex_triangulator_3< edge_index[4][4] = {{-1,0,1,2},{0,-1,3,4},{1,3,-1,5},{2,4,5,-1}}; +template < + class RegularTriangulation_3, + class TriangulatedMixedComplex_3, + class TriangulatedMixedComplexObserver_3> +void +Mixed_complex_triangulator_3< + RegularTriangulation_3, + TriangulatedMixedComplex_3, + TriangulatedMixedComplexObserver_3>:: +construct_anchor_del(Rt_Simplex const &sDel) { + Rt_Simplex s = compute_anchor_obj.anchor_del(sDel); + anchor_del2[sDel] = Anchor_map_iterator(); + + Anchor_map_iterator it = anchor_del2.find(sDel); + Anchor_map_iterator it2 = anchor_del2.find(s); + CGAL_assertion(it != anchor_del2.end()); + CGAL_assertion(it2 != anchor_del2.end()); + it->second = it2; + + // degenerate simplices: + if (compute_anchor_obj.is_degenerate()) { + it = find_anchor(anchor_del2, it); + typename Compute_anchor::Simplex_iterator degenerate_it; + for (degenerate_it = compute_anchor_obj.equivalent_anchors_begin(); + degenerate_it != compute_anchor_obj.equivalent_anchors_end(); + degenerate_it++) { + Anchor_map_iterator tmp; + it2 = anchor_del2.find(*degenerate_it); + CGAL_assertion(it2 != anchor_del2.end()); + // Merge sets: + while (it2 != it2->second) { + tmp = it2->second; + it2->second = it->second; + it2 = tmp; + CGAL_assertion(it2 != anchor_del2.end()); + } + it2->second = it->second; + } + } +} + template < class RegularTriangulation_3, class TriangulatedMixedComplex_3, @@ -318,9 +403,26 @@ construct_anchors() { Rt_Simplex s; // Compute anchor points: + for (vit=regular.finite_vertices_begin(); + vit!=regular.finite_vertices_end(); vit++) { + construct_anchor_del(Rt_Simplex(vit)); + } + for (eit=regular.finite_edges_begin(); + eit!=regular.finite_edges_end(); eit++) { + s = Rt_Simplex(*eit); + construct_anchor_del(s); + CGAL_assertion(s.dimension() == 1); + } + for (fit=regular.finite_facets_begin(); + fit!=regular.finite_facets_end(); fit++) { + s = Rt_Simplex(*fit); + construct_anchor_del(s); + CGAL_assertion(s.dimension() == 2); + } for (cit=regular.finite_cells_begin(); cit!=regular.finite_cells_end(); cit++) { s = Rt_Simplex(cit); + construct_anchor_del(s); construct_anchor_vor(s); CGAL_assertion(s.dimension() == 3); } @@ -366,55 +468,178 @@ construct_vertices() { Rt_Vertex_handle v1, v2, v3; Rt_Edge e; Rt_Cell_handle c1, c2; - Rt_Simplex sVor; + Rt_Simplex sDel, sVor; Tmc_Vertex_handle vh; if (verbose) std::cout << "construct_anchors" << std::endl; construct_anchors(); - if (verbose) std::cout << "4 "; + if (verbose) std::cout << "9 "; // anchor dimDel=0, dimVor=3 for (cit=regular.finite_cells_begin(); cit!=regular.finite_cells_end(); cit++) { sVor = get_anchor_vor(Rt_Simplex(cit)); - if (anchors.find(sVor) == anchors.end()) { - vh = add_vertex(sVor); - anchors[sVor] = vh; - CGAL_assertion(vh == get_vertex(sVor)); + for (int i=0; i<4; i++) { + sDel = get_anchor_del(Rt_Simplex(cit->vertex(i))); + if (anchors.find(Symb_anchor(sDel,sVor)) == anchors.end()) { + vh = add_vertex(Symb_anchor(sDel,sVor)); + anchors[Symb_anchor(sDel,sVor)] = vh; + CGAL_assertion(vh == get_vertex(sDel, sVor)); + } } } - if (verbose) std::cout << "3 "; + if (verbose) std::cout << "8 "; + // anchor dimDel=1, dimVor=3 + for (cit=regular.finite_cells_begin(); cit!=regular.finite_cells_end(); cit++) { + sVor = get_anchor_vor(Rt_Simplex(cit)); + for (int i=0; i<3; i++) { + for (int j=i+1; j<4; j++) { + sDel = get_anchor_del(Rt_Simplex(Rt_Edge(cit,i,j))); + if (anchors.find(Symb_anchor(sDel,sVor)) == anchors.end()) { + vh = add_vertex(Symb_anchor(sDel,sVor)); + anchors[Symb_anchor(sDel,sVor)] = vh; + assert(vh == get_vertex(sDel, sVor)); + } + } + } + } + + if (verbose) std::cout << "7 "; // anchor dimDel=2, dimVor=3 and dimDel=0, dimVor=2 for (fit=regular.finite_facets_begin(); fit!=regular.finite_facets_end(); fit++) { + // anchor dimDel=2, dimVor=3 + c1 = fit->first; + c2 = c1->neighbor(fit->second); + + sDel = get_anchor_del(*fit); + if (!regular.is_infinite(c1)) { + sVor = get_anchor_vor(c1); + if (anchors.find(Symb_anchor(sDel,sVor)) == anchors.end()) { + vh = add_vertex(Symb_anchor(sDel,sVor)); + anchors[Symb_anchor(sDel,sVor)] = vh; + assert(vh == get_vertex(sDel, sVor)); + } + } + if (!regular.is_infinite(c2)) { + sVor = get_anchor_vor(c2); + if (anchors.find(Symb_anchor(sDel,sVor)) == anchors.end()) { + vh = add_vertex(Symb_anchor(sDel,sVor)); + anchors[Symb_anchor(sDel,sVor)] = vh; + assert(vh == get_vertex(sDel, sVor)); + } + } // anchor dimDel=0, dimVor=2 sVor = get_anchor_vor(*fit); - if (anchors.find(sVor) == anchors.end()) { - vh = add_vertex(sVor); - anchors[sVor] = vh; - assert(vh == get_vertex(sVor)); + for (int i=1; i<4; i++) { + sDel = get_anchor_del(Rt_Simplex(c1->vertex((fit->second+i)&3))); + if (anchors.find(Symb_anchor(sDel,sVor)) == anchors.end()) { + vh = add_vertex(Symb_anchor(sDel,sVor)); + anchors[Symb_anchor(sDel,sVor)] = vh; + assert(vh == get_vertex(sDel, sVor)); + } else { + vh = get_vertex(sDel, sVor); + } + } + } + + if (verbose) std::cout << "6 "; + // anchor dimDel=0, dimVor=1 + for (eit=regular.finite_edges_begin(); eit!=regular.finite_edges_end(); eit++) { + sVor = get_anchor_vor(*eit); + + v1 = eit->first->vertex(eit->second); + v2 = eit->first->vertex(eit->third); + sDel = get_anchor_del(v1); + if (anchors.find(Symb_anchor(sDel,sVor)) == anchors.end()) { + vh = add_vertex(Symb_anchor(sDel,sVor)); + anchors[Symb_anchor(sDel,sVor)] = vh; + assert(vh == get_vertex(sDel, sVor)); + } + + sDel = get_anchor_del(v2); + if (anchors.find(Symb_anchor(sDel,sVor)) == anchors.end()) { + vh = add_vertex(Symb_anchor(sDel,sVor)); + anchors[Symb_anchor(sDel,sVor)] = vh; + assert(vh == get_vertex(sDel, sVor)); + } + } + + if (verbose) std::cout << "5 "; + // anchor dimDel=3, dimVor=3 + for (cit=regular.finite_cells_begin(); cit!=regular.finite_cells_end(); cit++) { + sDel = get_anchor_del(Rt_Simplex(cit)); + sVor = get_anchor_vor(Rt_Simplex(cit)); + if (anchors.find(Symb_anchor(sDel,sVor)) == anchors.end()) { + vh = add_vertex(Symb_anchor(sDel,sVor)); + anchors[Symb_anchor(sDel,sVor)] = vh; + assert(vh == get_vertex(sDel, sVor)); + } + } + + + if (verbose) std::cout << "4 "; + // anchor dimDel=0, dimVor=0 + for (vit=regular.finite_vertices_begin(); vit!=regular.finite_vertices_end(); vit++) { + sDel = get_anchor_del(Rt_Simplex(vit)); + sVor = get_anchor_vor(Rt_Simplex(vit)); + if (anchors.find(Symb_anchor(sDel,sVor)) == anchors.end()) { + vh = add_vertex(Symb_anchor(sDel,sVor)); + anchors[Symb_anchor(sDel,sVor)] = vh; + assert(vh == get_vertex(sDel, sVor)); + } + } + + if (verbose) std::cout << "3 "; + // anchor dimDel=1, dimVor=2 + for (fit=regular.finite_facets_begin(); fit!=regular.finite_facets_end(); fit++) { + c1 = fit->first; + c2 = c1->neighbor(fit->second); + + sVor = get_anchor_vor(Rt_Simplex(*fit)); + for (int i=1; i<3; i++) { + for (int j=i+1; j<4; j++) { + e.first = c1; + e.second = (fit->second+i)&3; + e.third = (fit->second+j)&3; + sDel = get_anchor_del(Rt_Simplex(e)); + if (anchors.find(Symb_anchor(sDel,sVor)) == anchors.end()) { + vh = add_vertex(Symb_anchor(sDel,sVor)); + anchors[Symb_anchor(sDel,sVor)] = vh; + assert(vh == get_vertex(sDel, sVor)); + } + } } } if (verbose) std::cout << "2 "; - // anchor dimDel=0, dimVor=1 - for (eit=regular.finite_edges_begin(); eit!=regular.finite_edges_end(); eit++) { - sVor = get_anchor_vor(*eit); - if (anchors.find(sVor) == anchors.end()) { - vh = add_vertex(sVor); - anchors[sVor] = vh; - assert(vh == get_vertex(sVor)); + // anchor dimDel=2, dimVor=2 + for (fit=regular.finite_facets_begin(); fit!=regular.finite_facets_end(); fit++) { + c1 = fit->first; + c2 = c1->neighbor(fit->second); + + sVor = get_anchor_vor(Rt_Simplex(*fit)); + sDel = get_anchor_del(Rt_Simplex(*fit)); + if (anchors.find(Symb_anchor(sDel,sVor)) == anchors.end()) { + vh = add_vertex(Symb_anchor(sDel,sVor)); + anchors[Symb_anchor(sDel,sVor)] = vh; + assert(vh == get_vertex(sDel, sVor)); } } - if (verbose) std::cout << "1 "; - // anchor dimDel=0, dimVor=0 - for (vit=regular.finite_vertices_begin(); vit!=regular.finite_vertices_end(); vit++) { - sVor = get_anchor_vor(Rt_Simplex(vit)); - if (anchors.find(sVor) == anchors.end()) { - vh = add_vertex(sVor); - anchors[sVor] = vh; - assert(vh == get_vertex(sVor)); + if (verbose) std::cout << "1" << std::endl; + // anchor dimDel=1, dimVor=1 + for (eit=regular.finite_edges_begin(); eit!=regular.finite_edges_end(); eit++) { + v1 = eit->first->vertex(eit->second); + v2 = eit->first->vertex(eit->third); + + sVor = get_anchor_vor(Rt_Simplex(*eit)); + sDel = get_anchor_del(Rt_Simplex(*eit)); + + if (anchors.find(Symb_anchor(sDel,sVor)) == anchors.end()) { + vh = add_vertex(Symb_anchor(sDel,sVor)); + anchors[Symb_anchor(sDel,sVor)] = vh; + assert(vh == get_vertex(sDel, sVor)); } } } @@ -430,49 +655,253 @@ Mixed_complex_triangulator_3< RegularTriangulation_3, TriangulatedMixedComplex_3, TriangulatedMixedComplexObserver_3>:: -construct_cells() { - Rt_Simplex sVor_v, sVor_e, sVor_f, sVor_c; +construct_0_cell(Rt_Vertex_handle rt_vh) { + Rt_Simplex sDel_v, sVor_v, sVor_e, sVor_f, sVor_c; Tmc_Vertex_handle vh[4]; - for (Rt_Finite_vertices_iterator vit=regular.finite_vertices_begin(); - vit!=regular.finite_vertices_end(); vit++) { + Rt_Simplex simplex(rt_vh); + sDel_v = get_anchor_del(Rt_Simplex(rt_vh)); + sVor_v = get_anchor_vor(Rt_Simplex(rt_vh)); + vh[0] = get_vertex(sDel_v,sVor_v); - Rt_Simplex simplex(vit); - sVor_v = get_anchor_vor(Rt_Simplex(vit)); - vh[0] = get_vertex(sVor_v); + std::list adj_cells; + typename std::list::iterator adj_cell; + regular.incident_cells(rt_vh, std::back_inserter(adj_cells)); - std::list adj_cells; - typename std::list::iterator adj_cell; - regular.incident_cells(vit, std::back_inserter(adj_cells)); - - // Construct cells: - for (adj_cell = adj_cells.begin(); - adj_cell != adj_cells.end(); - adj_cell ++) { - if (!regular.is_infinite(*adj_cell)) { - sVor_c = get_anchor_vor(Rt_Simplex(*adj_cell)); - vh[3] = get_vertex(sVor_c); - int index = (*adj_cell)->index(vit); - for (int i=1; i<4; i++) { - sVor_f = get_anchor_vor( - Rt_Simplex(Rt_Facet(*adj_cell,(index+i)&3))); - vh[2] = get_vertex(sVor_f); + // Construct cells: + for (adj_cell = adj_cells.begin(); + adj_cell != adj_cells.end(); + adj_cell ++) { + if (!regular.is_infinite(*adj_cell)) { + sVor_c = get_anchor_vor(Rt_Simplex(*adj_cell)); + vh[3] = get_vertex(sDel_v,sVor_c); + int index = (*adj_cell)->index(rt_vh); + for (int i=1; i<4; i++) { + sVor_f = get_anchor_vor(Rt_Simplex(Rt_Facet(*adj_cell,(index+i)&3))); + vh[2] = get_vertex(sDel_v,sVor_f); - for (int j=1; j<4; j++) { - if (j!=i) { - sVor_e = get_anchor_vor( - Rt_Simplex(Rt_Edge(*adj_cell,index,(index+j)&3))); - vh[1] = get_vertex(sVor_e); - if ((vh[0] != vh[1]) && (vh[1] != vh[2]) && (vh[2] != vh[3])) { - CGAL_assertion(sVor_v != sVor_e); - CGAL_assertion(sVor_e != sVor_f); - CGAL_assertion(sVor_f != sVor_c); - Tmc_Cell_handle ch = - add_cell(vh,(index + (j==(i%3+1)? 1:0))&1,simplex); + for (int j=1; j<4; j++) { + if (j!=i) { + sVor_e = get_anchor_vor( + Rt_Simplex(Rt_Edge(*adj_cell,index,(index+j)&3))); + vh[1] = get_vertex(sDel_v,sVor_e); + if ((vh[0] != vh[1]) && (vh[1] != vh[2]) && (vh[2] != vh[3])) { + CGAL_assertion(sVor_v != sVor_e); + CGAL_assertion(sVor_e != sVor_f); + CGAL_assertion(sVor_f != sVor_c); + Tmc_Cell_handle ch = + add_cell(vh,(index + (j==(i%3+1)? 1:0))&1,simplex); + } + } + } + } + } + } +} + +// Constructs 1-cells of the mixed complex corresponding to edges +// of the regular triangulation +template < + class RegularTriangulation_3, + class TriangulatedMixedComplex_3, + class TriangulatedMixedComplexObserver_3> +void +Mixed_complex_triangulator_3< + RegularTriangulation_3, + TriangulatedMixedComplex_3, + TriangulatedMixedComplexObserver_3>:: +construct_1_cell(const Rt_Finite_edges_iterator &e) { + Rt_Simplex sDel_v, sDel_e, sVor_e, sVor_f, sVor_c; + Tmc_Vertex_handle vh[4]; + Rt_Vertex_handle v[2]; + Tmc_Cell_handle ch; + + Rt_Simplex mixed_cell_simplex(*e); + sDel_e = get_anchor_del(Rt_Simplex(*e)); + sVor_e = get_anchor_vor(Rt_Simplex(*e)); + + v[0] = e->first->vertex(e->second); + v[1] = e->first->vertex(e->third); + + // Construct cells on the side of v[vi]: + for (int vi=0; vi<2; vi++) { + sDel_v = get_anchor_del(Rt_Simplex(v[vi])); + if (!(sDel_v == sDel_e)) { + Rt_Cell_circulator ccir, cstart; + ccir = cstart = regular.incident_cells(*e); + do { + if (!regular.is_infinite(ccir)) { + int index0 = ccir->index(v[vi]); + int index1 = ccir->index(v[1-vi]); + + sVor_c = get_anchor_vor(Rt_Simplex(ccir)); + + for (int fi=1; fi<4; fi++) { + if (((index0+fi)&3) != index1) { + sVor_f = + get_anchor_vor(Rt_Simplex(Rt_Facet(ccir,(index0+fi)&3))); + if ((sVor_c != sVor_f) && (sVor_f != sVor_e)) { + vh[0] = get_vertex(sDel_v, sVor_e); + vh[1] = get_vertex(sDel_e, sVor_e); + vh[2] = get_vertex(sDel_e, sVor_f); + vh[3] = get_vertex(sDel_e, sVor_c); + int orient; + if (((4+index1-index0)&3) == 1) { + orient = (index1 + (fi==2))&1; + } else { + orient = (index1 + (fi==1))&1; + } + // vh: dimension are (01,11,12,13) + ch = add_cell(vh,orient,mixed_cell_simplex); + + vh[1] = get_vertex(sDel_v, sVor_f); + // vh: dimension are (01,02,12,13) + ch = add_cell(vh,1-orient,mixed_cell_simplex); + + vh[2] = get_vertex(sDel_v, sVor_c); + // vh: dimension are (01,02,03,13) + ch = add_cell(vh,orient,mixed_cell_simplex); } } } } + ccir ++; + } while (ccir != cstart); + } + } +} + + +// Constructs 2-cells of the mixed complex corresponding to facets +// of the regular triangulation +template < + class RegularTriangulation_3, + class TriangulatedMixedComplex_3, + class TriangulatedMixedComplexObserver_3> +void +Mixed_complex_triangulator_3< + RegularTriangulation_3, + TriangulatedMixedComplex_3, + TriangulatedMixedComplexObserver_3>:: +construct_2_cell(const Rt_Finite_facets_iterator &fit) { + Rt_Simplex sDel_v, sDel_e, sDel_f, sVor_f, sVor_c; + Tmc_Vertex_handle vh[4]; // Implicit function over vLabels is increasing ... + Rt_Cell_handle rt_ch; + int index; + + rt_ch = fit->first; + index = fit->second; + Rt_Simplex simplex(*fit); + sDel_f = get_anchor_del(Rt_Simplex(*fit)); + sVor_f = get_anchor_vor(Rt_Simplex(*fit)); + + for (int i=0; i<2; i++) { // Do this twice + if (!regular.is_infinite(rt_ch)) { + sVor_c = get_anchor_vor(Rt_Simplex(rt_ch)); + + vh[3] = get_vertex(sDel_f, sVor_c); + Tmc_Vertex_handle vh2 = get_vertex(sDel_f, sVor_f); + if (vh2 != vh[3]) { + // Facet and cell do not coincide .. + for (int vi=1; vi<4; vi++) { + sDel_v = get_anchor_del(Rt_Simplex(rt_ch->vertex((index+vi)&3))); + //index_02[rt_ch].V[index][(index+vi)&3]; + vh[0] = get_vertex(sDel_v, sVor_f); + for (int ei=1; ei<4; ei++) { + if (vi != ei) { + vh[2] = vh2; + int index0 = (index+vi)&3; + int index1 = (index+ei)&3; + int fi = (6+index-vi-ei)&3;//6-index-index0-index1; + sDel_e = + get_anchor_del(Rt_Simplex(Rt_Edge(rt_ch, index0, index1))); + vh[1] = get_vertex(sDel_e, sVor_f); + //index_12[rt_ch].V[index][(6+index-vi-ei)&3]; + if ((vh[0] != vh[1]) && (vh[1] != vh[2])) { + // index0: v0 + // index1: v1 + // index0+fi&3 == facet + int orient; + + if (((4+index1-index0)&3) == 3) { + orient = (index1 + (((4+index0-fi)&3)==2))&1; + } else { + orient = (index1 + (((4+index0-fi)&3)==1))&1; + } + + add_cell(vh,orient,simplex); + + vh[2] = get_vertex(sDel_e, sVor_c); + add_cell(vh,1-orient,simplex); + + vh[1] = get_vertex(sDel_v, sVor_c); + add_cell(vh,orient,simplex); + } + } + } + } + } + } + // swap to the other cell + Rt_Cell_handle ch_old = rt_ch; + rt_ch = rt_ch->neighbor(index); + index = rt_ch->index(ch_old); + } + + CGAL_assertion(rt_ch == fit->first); + CGAL_assertion(index == fit->second); +} + + +// Constructs 3-cells of the mixed complex corresponding to cells +// of the regular triangulation +template < + class RegularTriangulation_3, + class TriangulatedMixedComplex_3, + class TriangulatedMixedComplexObserver_3> +void +Mixed_complex_triangulator_3< + RegularTriangulation_3, + TriangulatedMixedComplex_3, + TriangulatedMixedComplexObserver_3>:: +construct_3_cell(Rt_Cell_handle rt_ch) { + Rt_Simplex sDel_v, sDel_e, sDel_f, sDel_c, sVor_c; + Tmc_Vertex_handle vh[4]; + Tmc_Cell_handle ch; + + // construct the tetrahedron: + // C[ch], C[Facet(ch,fi)], C[Edge(ch,ei,vi)], C[ch->vertex(vi)] + sDel_c = get_anchor_del(Rt_Simplex(rt_ch)); + sVor_c = get_anchor_vor(Rt_Simplex(rt_ch)); + Rt_Simplex simplex = Rt_Simplex(rt_ch); + vh[0] = get_vertex(sDel_c, sVor_c); + for (int fi=0; fi<4; fi++) { + sDel_f = get_anchor_del(Rt_Simplex(Rt_Facet(rt_ch, fi))); + vh[1] = get_vertex(sDel_f, sVor_c); + if (vh[0] != vh[1]) { + for (int vi=1; vi<4; vi++) { + int index0 = (fi+vi)&3; + sDel_v = get_anchor_del(Rt_Simplex(rt_ch->vertex(index0))); + for (int ei=1; ei<4; ei++) { + int index1 = (fi+ei)&3; + if (vi != ei) { + sDel_e = get_anchor_del(Rt_Simplex(Rt_Edge(rt_ch, index0, index1))); + vh[2] = get_vertex(sDel_e, sVor_c); + // index_13[rt_ch].V[edge_index[index0][index1]]; + vh[3] = get_vertex(sDel_v, sVor_c); + // index_03[rt_cit].V[index0]; + if ((vh[1] != vh[2]) && (vh[2] != vh[3])) { + int orient; + + if (((4+index1-index0)&3) == 1) { + orient = (index1 + (vi==2))&1; + } else { + orient = (index1 + (vi==3))&1; + } + ch = add_cell(vh, orient, simplex); + } + } + } } } } @@ -491,12 +920,12 @@ Mixed_complex_triangulator_3< RegularTriangulation_3, TriangulatedMixedComplex_3, TriangulatedMixedComplexObserver_3>:: -add_vertex (Rt_Simplex const &anchor) +add_vertex (Symb_anchor const &anchor) { Tmc_Vertex_handle vh; vh = triangulation_incr_builder.add_vertex(); - vh->point() = get_anchor(anchor); - observer.after_vertex_insertion(anchor, anchor, vh); + vh->point() = get_anchor(anchor.first, anchor.second); + observer.after_vertex_insertion(anchor.first, anchor.second, vh); return vh; } @@ -513,11 +942,14 @@ typename Mixed_complex_triangulator_3< Mixed_complex_triangulator_3< RegularTriangulation_3, TriangulatedMixedComplex_3, - TriangulatedMixedComplexObserver_3>::get_vertex (Rt_Simplex &sVor) + TriangulatedMixedComplexObserver_3>::get_vertex ( + Rt_Simplex &sDel, Rt_Simplex &sVor) { + Rt_Simplex sDel2 = get_anchor_del(sDel); Rt_Simplex sVor2 = get_anchor_vor(sVor); + CGAL_assertion(sDel == sDel2); CGAL_assertion(sVor == sVor2); - Tmc_Vertex_handle vh = anchors[sVor2]; + Tmc_Vertex_handle vh = anchors[Symb_anchor(sDel2,sVor2)]; CGAL_assertion(vh != Tmc_Vertex_handle()); return vh; } @@ -622,9 +1054,12 @@ Mixed_complex_triangulator_3< RegularTriangulation_3, TriangulatedMixedComplex_3, TriangulatedMixedComplexObserver_3>:: -get_anchor(Rt_Simplex const &sVor) +get_anchor(Rt_Simplex const &sDel, Rt_Simplex const &sVor) { - return get_orthocenter(sVor); + Tmc_Point dfoc = get_orthocenter(sDel); + Tmc_Point vfoc = get_orthocenter(sVor); + + return construct_anchor_point(dfoc, vfoc); } template < @@ -651,17 +1086,14 @@ remove_small_edges() // NGHK: This may intrudoce rounding errors, since the quadratic surface // may change: Tmc_Vertex_handle vh, vh_collapse_to; - Tmc_Finite_vertices_iterator vit = _tmc.finite_vertices_begin(); - int nCollapsed=0; - while (vit != _tmc.finite_vertices_end()) { + for (Tmc_Finite_vertices_iterator vit = _tmc.finite_vertices_begin(); + vit != _tmc.finite_vertices_end(); ) { vh = vit; vit++; if (is_collapsible(vh, vh_collapse_to,sq_length)) { - nCollapsed ++; do_collapse(vh,vh_collapse_to); } } - std::cout << "Collapsed: " << nCollapsed << std::endl; } template < @@ -696,8 +1128,8 @@ is_collapsible(Tmc_Vertex_handle vh, it = incident_vertices.begin(); it != incident_vertices.end(); it++) { if ((_tmc.geom_traits().compute_squared_distance_3_object()(vh->point(), - (*it)->point()) - < sq_length) && + (*it)->point()) + < sq_length) && (vh->cell()->surf == (*it)->cell()->surf) && (vh->sign() == (*it)->sign())) { bool ok = true; @@ -774,37 +1206,37 @@ template < class TriangulatedMixedComplex_3, class TriangulatedMixedComplexObserver_3> void -triangulate_mixed_complex_3 -(RegularTriangulation_3 const &rt, - typename RegularTriangulation_3::Geom_traits::FT shrink, - TriangulatedMixedComplex_3 &tmc, - TriangulatedMixedComplexObserver_3 &observer, - bool verbose) +triangulate_mixed_complex_3(RegularTriangulation_3 &rt, + typename RegularTriangulation_3::Geom_traits::FT + const & shrink_factor, + TriangulatedMixedComplex_3 &tmc, + TriangulatedMixedComplexObserver_3 &observer, + bool verbose) { typedef Mixed_complex_triangulator_3< RegularTriangulation_3, TriangulatedMixedComplex_3, TriangulatedMixedComplexObserver_3> Mixed_complex_triangulator; - Mixed_complex_triangulator(rt, tmc, observer, verbose); + Mixed_complex_triangulator(rt, shrink_factor, tmc, observer, verbose); } -// template < -// class RegularTriangulation_3, -// class TriangulatedMixedComplex_3> -// void -// triangulate_mixed_complex_3 -// (RegularTriangulation_3 const ®ular, -// typename RegularTriangulation_3::Geom_Traits::FT shrink, -// TriangulatedMixedComplex_3 &tmc, -// bool verbose) -// { -// Triangulated_mixed_complex_observer_3< -// TriangulatedMixedComplex_3, const RegularTriangulation_3> -// observer(1); -// triangulate_mixed_complex_3(regular, shrink, tmc, observer, verbose); -// } +template < + class RegularTriangulation_3, + class TriangulatedMixedComplex_3> +void +triangulate_mixed_complex_3(RegularTriangulation_3 const ®ular, + typename RegularTriangulation_3::Geom_traits::FT + const &shrink_factor, + TriangulatedMixedComplex_3 &tmc, + bool verbose) +{ + Triangulated_mixed_complex_observer_3< + TriangulatedMixedComplex_3, const RegularTriangulation_3> + observer(shrink_factor); + triangulate_mixed_complex_3(regular, shrink_factor, tmc, observer, verbose); +} CGAL_END_NAMESPACE -#endif // CGAL_TRIANGULATE_MIXED_COMPLEX_3_H +#endif // CGAL_TRIANGULATE_MIXED_COMPLEX_H diff --git a/Skin_surface_3/test/Skin_surface_3/surface_test.C b/Skin_surface_3/test/Skin_surface_3/surface_test.C index 693ac26f7b8..777ebf2faba 100644 --- a/Skin_surface_3/test/Skin_surface_3/surface_test.C +++ b/Skin_surface_3/test/Skin_surface_3/surface_test.C @@ -10,15 +10,15 @@ #include typedef CGAL::Exact_predicates_inexact_constructions_kernel Inexact_K; -typedef CGAL::Skin_surface_traits_3 Traits; +typedef CGAL::Exact_predicates_exact_constructions_kernel Exact_K; +typedef CGAL::Skin_surface_traits_3 Traits; +//typedef CGAL::Skin_surface_traits_3 Traits; typedef CGAL::Skin_surface_3 Skin_surface; -typedef Skin_surface::Simplex Simplex; -typedef Skin_surface::FT FT; -typedef Skin_surface::Weighted_point Weighted_point; +typedef Skin_surface::Simplex Simplex; +typedef Skin_surface::FT FT; +typedef Skin_surface::Weighted_point Weighted_point; typedef Weighted_point::Point Bare_point; typedef CGAL::Polyhedron_3 Polyhedron; -typedef CGAL::Exact_predicates_exact_constructions_kernel Exact_K; -typedef CGAL::Skin_surface_quadratic_surface_3 Quadratic_surface; CGAL::Cartesian_converter e2i_converter; CGAL::Cartesian_converter i2e_converter; @@ -51,58 +51,36 @@ public: Skin_surface skin_surface(l.begin(), l.end(), s); - TMC tmc; - CGAL::Triangulated_mixed_complex_observer_3 - observer(skin_surface.shrink_factor()); - triangulate_mixed_complex_3(skin_surface.get_regular_triangulation(), - skin_surface.shrink_factor(), - tmc, - observer, - false); -// CGAL::triangulate_mixed_complex_3(skin_surface.get_regular_triangulation(), -// skin_surface.shrink_factor(), -// tmc, -// false // verbose -// ); + const TMC &tmc = skin_surface.triangulated_mixed_complex(); +// CGAL::Triangulated_mixed_complex_observer_3 +// observer(skin_surface.shrink_factor()); +// triangulate_mixed_complex_3(skin_surface.get_regular_triangulation(), +// skin_surface.shrink_factor(), +// tmc, +// observer, +// false); for (TMC_Finite_vertices_iterator vit = tmc.finite_vertices_begin(); vit != tmc.finite_vertices_end(); vit++) { + if (tmc.is_infinite(vit->cell())) { std::cerr << "ERROR: is_infinite (main)" << std::endl; } - Quadratic_surface::FT - val = vit->cell()->info().second->value(vit->point()); + Exact_K::FT val = vit->cell()->info().second->value(vit->point()); std::list cells; tmc.incident_cells(vit, std::back_inserter(cells)); for (std::list::iterator cell = cells.begin(); cell != cells.end(); cell++) { if (!tmc.is_infinite(*cell)) { - Quadratic_surface::FT val2 = (*cell)->info().second->value(vit->point()); - CGAL_assertion(val == val2); + Exact_K::FT val2 = (*cell)->info().second->value(vit->point()); + // NGHK: Make exact: + //CGAL_assertion(val == val2); + CGAL_assertion(std::abs(CGAL::to_double(val-val2)) < 1e-8); } } } -// for (TMC_Finite_cells_iterator cit = tmc.finite_cells_begin(); -// cit != tmc.finite_cells_end(); cit++) { -// Bare_point baryc = -// e2i_converter(cit->vertex(0)->point() + -// (cit->vertex(1)->point()-cit->vertex(0)->point())/4 + -// (cit->vertex(2)->point()-cit->vertex(0)->point())/4 + -// (cit->vertex(3)->point()-cit->vertex(0)->point())/4); -// if (tmc.tetrahedron(cit).has_on_bounded_side(i2e_converter(baryc))) { -// Quadratic_surface::FT val1 = cit->info()->value(i2e_converter(baryc)); -// Simplex s = skin_surface.locate_mixed(baryc); -// Quadratic_surface::FT val2 = -// skin_surface.construct_surface(s, Exact_K()).value(i2e_converter(baryc)); -// // std::cout << val1 << " " << val2 << " " << val2-val1 << std::endl; -// CGAL_assertion(val1==val2); -// } else { -// std::cout << "Barycenter on unbounded side, due to rounding errors\n"; -// } -// } - } private: double s; @@ -110,21 +88,21 @@ private: int main(int argc, char *argv[]) { std::vector filenames; - filenames.push_back("data/degenerate.cin"); - filenames.push_back("data/test1.cin"); - filenames.push_back("data/test2.cin"); - filenames.push_back("data/test3.cin"); - filenames.push_back("data/test4.cin"); - filenames.push_back("data/test5.cin"); - filenames.push_back("data/test6.cin"); - filenames.push_back("data/test7.cin"); - filenames.push_back("data/test8.cin"); - filenames.push_back("data/test9.cin"); - filenames.push_back("data/test10.cin"); - filenames.push_back("data/test11.cin"); + filenames.push_back("data/test1.cin"); + filenames.push_back("data/test2.cin"); + filenames.push_back("data/test3.cin"); + filenames.push_back("data/test4.cin"); + filenames.push_back("data/test5.cin"); + filenames.push_back("data/test6.cin"); + filenames.push_back("data/test7.cin"); + filenames.push_back("data/test8.cin"); + filenames.push_back("data/test9.cin"); + filenames.push_back("data/test10.cin"); + filenames.push_back("data/test11.cin"); + filenames.push_back("data/degenerate.cin"); - std::for_each(filenames.begin(), filenames.end(), Test_file(.85)); std::for_each(filenames.begin(), filenames.end(), Test_file(.5)); + std::for_each(filenames.begin(), filenames.end(), Test_file(.85)); return 0; }