Merge pull request #2404 from sloriot/SM_3-fix_border_halfedge_for_border_node

Fix issues on the border of subdivided mesh

Subdivision_method_3/doc/Subdivision_method_3/Concepts/PQQMask_3.h

@@ -66,7 +66,9 @@ void vertex_node(vertex_descriptor vd, Point_3& pt);
 /*!
 
 computes the edge-point `ept` and the vertex-point `vpt`
-based on the neighborhood of the border edge `hd`.
+based on the neighborhood of the border edge of `hd`. `hd` is not
+a border halfedge (its opposite is) and `vpt` corresponds to the
+target vertex of `hd`.
 */
 void border_node(halfedge_descriptor hd, Point_3& ept, Point_3& vpt);
 

Subdivision_method_3/doc/Subdivision_method_3/Concepts/PTQMask_3.h

@@ -55,9 +55,10 @@ of the vertex `vd`.
 void vertex_node(vertex_descriptor vd, Point_3& pt);
 
 /*!
-
 computes the edge-point `ept` and the vertex-point `vpt`
-based on the neighborhood of the border edge `hd`.
+based on the neighborhood of the border edge of `hd`. `hd` is not
+a border halfedge (its opposite is) and `vpt` corresponds to the
+target vertex of `hd`.
 */
 void border_node(halfedge_descriptor hd, Point_3& ept, Point_3& vpt);
 

Subdivision_method_3/doc/Subdivision_method_3/Concepts/Sqrt3Mask_3.h

@@ -55,10 +55,12 @@ of the vertex `vd`.
 void vertex_node(vertex_descriptor vd, Point& pt);
 
 /*!
-computes the subdivided points `ept` and `vpt` based on the neighborhood
+computes the subdivided points `ept1` and `ept2` based on the neighborhood
-of the border halfedge `hd`.
+of the halfedge `hd` (whose opposite is on the border).
+Along `hd`, `ept1` comes before `ept2`.
+`vpt` is the updated point for the target vertex of `hd`.
 */
-void border_node(halfedge_descriptor hd, Point& ept, Point& vpt);
+void border_node(halfedge_descriptor hd, Point& ept1, Point& ept2, Point& vpt);
 
 /// @}
 

Subdivision_method_3/include/CGAL/Subdivision_method_3/internal/subdivision_hosts_impl_3.h

@@ -97,9 +97,11 @@ void PQQ_1step(Poly& p, VertexPointMap vpm, Mask mask) {
     std::size_t i = 0;
     BOOST_FOREACH(edge_descriptor ed, edges(p)){
       if(is_border(ed,p)){
-        int v = v_index[target(ed,p)];
+        halfedge_descriptor h=halfedge(ed,p);
+        if (is_border(h,p)) h=opposite(h,p);
+        int v = v_index[target(h,p)];
         v_onborder[v] = true;
-        mask.border_node(halfedge(ed,p), edge_point_buffer[i], vertex_point_buffer[v]);
+        mask.border_node(h, edge_point_buffer[i], vertex_point_buffer[v]);
 
       }else{
         mask.edge_node(halfedge(ed,p), edge_point_buffer[i]);
@@ -211,9 +213,11 @@ void PTQ_1step(Poly& p, VertexPointMap vpm, Mask mask) {
       if(! is_border(ed,p)){
         mask.edge_node(halfedge(ed,p), edge_point_buffer[i]);
       } else{
-        int v = v_index[target(ed,p)];
+        halfedge_descriptor h = halfedge(ed,p);
+        if (is_border(h, p)) h = opposite(h,p);
+        int v = v_index[target(h,p)];
         v_onborder[v] = true;
-        mask.border_node(halfedge(ed,p), edge_point_buffer[i], vertex_point_buffer[v]);
+        mask.border_node(h, edge_point_buffer[i], vertex_point_buffer[v]);
       }
       ++i;
     }
@@ -577,6 +581,7 @@ void Sqrt3_1step(Poly& p, VertexPointMap vpm, Mask mask,
   // halfedge corresponds to THE (there can only be one) border edge.
   std::vector<halfedge_descriptor> face_halfedge_border(num_f,
                                                         boost::graph_traits<Poly>::null_halfedge());
+  std::list<std::pair<vertex_descriptor, Point> > new_positions;
 
   // compute the positions of new points
   std::size_t i = 0;
@@ -587,8 +592,9 @@ void Sqrt3_1step(Poly& p, VertexPointMap vpm, Mask mask,
       BOOST_FOREACH(halfedge_descriptor hd, halfedges_around_face(halfedge(fd, p), p)) {
         if(is_border(opposite(hd, p),p)) {
           face_halfedge_border[face_id] = hd;
-          halfedge_descriptor bhd = opposite(hd, p);
+          Point vpt;
-          mask.border_node(bhd, cpt[i], cpt[i+1]);
+          mask.border_node(hd, cpt[i+1], cpt[i], vpt);
+          new_positions.push_back(std::make_pair(target(hd,p), vpt));
           i += 2;
 
           // the border subdivision is only performed every second subdivision
@@ -607,7 +613,6 @@ void Sqrt3_1step(Poly& p, VertexPointMap vpm, Mask mask,
   }
 
   // smooth the position of existing vertices
-  std::list<std::pair<vertex_descriptor, Point> > new_positions;
   BOOST_FOREACH(vertex_descriptor vd, vertices(p)){
     Point pt;
     if(!is_border(vd, p)) {
@@ -636,21 +641,6 @@ void Sqrt3_1step(Poly& p, VertexPointMap vpm, Mask mask,
     }
   }
 
-  if(refine_border) {
-    // collect the new positions for border vertices
-    BOOST_FOREACH(halfedge_descriptor hd, halfedges(p)) {
-      // switch to the border halfedge
-      hd = opposite(hd, p);
-      if(!is_border(hd, p))
-        continue;
-
-      Point pt;
-      mask.border_node(hd, pt);
-      vertex_descriptor vd = target(hd, p);
-      new_positions.push_back(std::make_pair(vd, pt));
-    }
-  }
-
   // actually inserts the new positions in the vertex point property map
   typename std::list<std::pair<vertex_descriptor, Point> >::iterator it = new_positions.begin(),
       end = new_positions.end();

Subdivision_method_3/include/CGAL/Subdivision_method_3/subdivision_masks_3.h

@@ -585,11 +585,13 @@ public:
     pt = CGAL::ORIGIN + cv;
   }
 
-  /// computes the \f$ \sqrt{3}\f$ edge-points `ept` and `vpt` of the halfedge `hd`.
+  /// computes the \f$ \sqrt{3}\f$ edge-points `ept1` and `ept2` of the halfedge `hd`.
+  /// The updated point coordinates for the target vertex of `hd` is also computed and put in `vpt`.
   /// \attention Border subdivision only happens every second step of a <em>single</em>
   ///            successive \f$ \sqrt{3}\f$ subdivision (thus requiring a depth larger than 1).
-  void border_node(halfedge_descriptor bhd, Point& ept, Point& vpt) {
+  void border_node(halfedge_descriptor hd, Point& ept1, Point& ept2, Point& vpt) {
-    // this function takes a BORDER halfedge
+    // this function takes the opposite of a BORDER halfedge
+    halfedge_descriptor bhd = opposite(hd, *(this->pmesh));
     CGAL_precondition(is_border(bhd, *(this->pmesh)));
     vertex_descriptor prev_s = source(prev(bhd, *(this->pmesh)), *(this->pmesh));
     Vector prev_sv = get(this->vpmap, prev_s) - CGAL::ORIGIN;
@@ -604,22 +606,9 @@ public:
     Vector next_tv = get(this->vpmap, next_t) - CGAL::ORIGIN;
 
     const FT denom = 1./27.;
-    ept = CGAL::ORIGIN + denom * ( 10.*sv + 16.*tv + next_tv );
+    ept1 = CGAL::ORIGIN + denom * ( prev_sv + 16.*sv + 10.*tv );
-    vpt = CGAL::ORIGIN + denom * ( prev_sv + 16.*sv + 10.*tv );
+    ept2 = CGAL::ORIGIN + denom * ( 10.*sv + 16.*tv + next_tv );
-  }
+    vpt = CGAL::ORIGIN + 1./27. * ( 4*prev_sv + 19*sv + 4*tv );
-
-  void border_node(halfedge_descriptor bhd, Point& pt) {
-    // this function takes a BORDER halfedge
-    CGAL_precondition(is_border(bhd, *(this->pmesh)));
-
-    vertex_descriptor s = source(bhd, *(this->pmesh));
-    Vector sv = get(this->vpmap, s) - CGAL::ORIGIN;
-    vertex_descriptor c = target(bhd, *(this->pmesh));
-    Vector cv = get(this->vpmap, c) - CGAL::ORIGIN;
-    vertex_descriptor n = target(next(bhd, *(this->pmesh)), *(this->pmesh));
-    Vector nv = get(this->vpmap, n) - CGAL::ORIGIN;
-
-    pt = CGAL::ORIGIN + 1./27. * ( 4*sv + 19*cv + 4*nv );
   }
 };