From 130267384f7cda1d153e18c2b846acbc55209a5c Mon Sep 17 00:00:00 2001 From: Christina Vaz Date: Wed, 27 Jun 2018 21:34:48 -0400 Subject: [PATCH] fix for IDT as well as structures needed for hm integration --- .../CGAL/Heat_method_3/Heat_method_3.h | 2 + .../Intrinsic_Delaunay_Triangulation_3.h | 44 +++++++++++++++---- .../intrinsic_delaunay_triangulation_test.cpp | 2 +- 3 files changed, 39 insertions(+), 9 deletions(-) diff --git a/Heat_method_3/include/CGAL/Heat_method_3/Heat_method_3.h b/Heat_method_3/include/CGAL/Heat_method_3/Heat_method_3.h index c455d5faceb..be4d5e48d67 100644 --- a/Heat_method_3/include/CGAL/Heat_method_3/Heat_method_3.h +++ b/Heat_method_3/include/CGAL/Heat_method_3/Heat_method_3.h @@ -39,6 +39,8 @@ #include #include #include +#include + diff --git a/Heat_method_3/include/CGAL/Heat_method_3/Intrinsic_Delaunay_Triangulation_3.h b/Heat_method_3/include/CGAL/Heat_method_3/Intrinsic_Delaunay_Triangulation_3.h index cd527c6f4a6..93776cbee3f 100644 --- a/Heat_method_3/include/CGAL/Heat_method_3/Intrinsic_Delaunay_Triangulation_3.h +++ b/Heat_method_3/include/CGAL/Heat_method_3/Intrinsic_Delaunay_Triangulation_3.h @@ -44,6 +44,7 @@ #include #include #include +#include @@ -109,6 +110,15 @@ namespace Intrinsic_Delaunay_Triangulation_3 { typedef typename std::stack > edge_stack; + typedef typename std::array vertex_for_face; + typedef typename std::array Face_Vertex_Array; + + //all z coords will be 0 though + typedef typename std::array coords_for_face; + typedef typename std::array Face_Coords_Array; + + + public: Intrinsic_Delaunay_Triangulation_3(TriangleMesh tm, VertexDistanceMap vdm) @@ -178,21 +188,20 @@ namespace Intrinsic_Delaunay_Triangulation_3 { Index b_i = get(edge_id_map, edge(hd2,tm)); Index c_i = get(edge_id_map, edge(hd3,tm)); double a = edge_lengths(i,0); - double b = edge_lengths(b_i,0); + double b1 = edge_lengths(b_i,0); double c1 = edge_lengths(c_i,0); - double tan2a = CGAL::sqrt(std::abs(((-a+b+c1)*(a+b-c1))/((a+b+c1)*(-b+a+c1)))); + double tan2a = CGAL::sqrt(CGAL::abs(((c1-a+b1)*(-b1+a+c1))/((a+b1+c1)*(b1+a-c1)))); hd = opposite(hd,tm); hd2 =next(hd,tm); hd3 = next(hd2,tm); b_i = get(edge_id_map, edge(hd2,tm)); c_i = get(edge_id_map, edge(hd3,tm)); - b = edge_lengths(b_i,0); + double b2 = edge_lengths(b_i,0); double c2 = edge_lengths(c_i,0); - - double tan2d = CGAL::sqrt(std::abs(((-a+b+c2)*(a+b-c2))/((a+b+c2)*(a-b+c2)))); + double tan2d = CGAL::sqrt(CGAL::abs(((-a+b2+c2)*(a+b2-c2))/((a+b2+c2)*(a-b2+c2)))); double tan2ad = (tan2a + tan2d)/(1-tan2a*tan2d); double cosad = (1-tan2ad*tan2ad)/(1+tan2ad*tan2ad); - double new_length = CGAL::sqrt( CGAL::abs(c1*c1 + c2*c2 - 2*cosad)); + double new_length = CGAL::sqrt( CGAL::abs(b1*b1 + c2*c2 - 2*cosad)); edge_lengths(i,0) = new_length; } @@ -213,9 +222,10 @@ namespace Intrinsic_Delaunay_Triangulation_3 { stack.pop(); Index edge_i = get(edge_id_map,ed); + marked_edges(edge_i,0)=0; //if the edge itself is not locally delaunay, go back - if(!(is_edge_locally_delaunay(ed))) // && !CGAL::is_boundary(ed,tm)) need to fix this + if(!(is_edge_locally_delaunay(ed))) { if(!(is_border(ed,tm))) { @@ -223,9 +233,9 @@ namespace Intrinsic_Delaunay_Triangulation_3 { change_edge_length(edge_i,ed); halfedge_descriptor hd = (halfedge(ed, tm)); CGAL::Euler::flip_edge(hd, tm); - edge_descriptor next_edge= edge(next(hd,tm),tm); Index next_edge_i = get(edge_id_map, next_edge); + //if edge was already checked, go back and check again //for the 4 surrounding edges, since local 'geometry' changed, if(!(marked_edges(next_edge_i,0))) @@ -289,7 +299,21 @@ namespace Intrinsic_Delaunay_Triangulation_3 { } loop_over_edges(stack, mark_edges); + //now that edges are calculated, go through and for each face, calculate the vertex positions around it + BOOST_FOREACH(face_descriptor f, faces(tm)) + { + CGAL::Vertex_around_face_iterator vbegin, vend, vmiddle; + Index face_i = get(face_id_map, f); + boost::tie(vbegin, vend) = vertices_around_face(halfedge(f,tm),tm); + vertex_descriptor current = *(vbegin); + vertex_descriptor n_two = *(++vbegin); + vertex_descriptor n_three = *(++vbegin); + //fill in Face_Vertex_Array with these descriptors + //fill in face_coord_array with (0,0),(length,0)) and then (lki(costheta, sintheta)) + //use those to compute in heat method + + } } //todo:: determine which can be const @@ -299,8 +323,12 @@ namespace Intrinsic_Delaunay_Triangulation_3 { EdgeIndexMap epm; VertexDistanceMap vdm; int number_of_edges = num_edges(tm); + int number_of_faces = num_faces(tm); Eigen::VectorXd edge_lengths; Eigen::VectorXd mark_edges; + Face_Coords_Array *face_coords_array = new Face_Coords_Array[number_of_faces]; + Face_Vertex_Array *face_vertex_array = new Face_Vertex_Array[number_of_faces]; + }; } // namespace Intrinsic_Delaunay_Triangulation_3 diff --git a/Heat_method_3/test/Heat_method_3/intrinsic_delaunay_triangulation_test.cpp b/Heat_method_3/test/Heat_method_3/intrinsic_delaunay_triangulation_test.cpp index e4a01505302..1ec9bdc605a 100644 --- a/Heat_method_3/test/Heat_method_3/intrinsic_delaunay_triangulation_test.cpp +++ b/Heat_method_3/test/Heat_method_3/intrinsic_delaunay_triangulation_test.cpp @@ -27,7 +27,7 @@ int main() Mesh sm; Vertex_distance_map vertex_distance_map = get(Vertex_distance_tag(),sm); - std::ifstream in("data/camel_not_delaunay.off"); + std::ifstream in("data/brain100k.off"); in >> sm; if(!in || num_vertices(sm) == 0) { std::cerr << "Problem loading the input data" << std::endl;