CGAL内置的边塌陷算法代码解析

这个algorithm.run()就是实际的边塌陷算法具体实现

int edge_collapse(TM& tmesh,
                  const ShouldStop& should_stop,
                  // optional mesh information policies
                  const GT& traits,
                  const VertexIndexMap& vim, // defaults to get(vertex_index, tmesh)
                  const VertexPointMap& vpm, // defaults to get(vertex_point, tmesh)
                  const HalfedgeIndexMap& him, // defaults to get(edge_index, tmesh)
                  const EdgeIsConstrainedMap& ecm, // defaults to No_constrained_edge_map<TM>()
                  // optional strategy policies - defaults to LindstomTurk
                  const GetCost& get_cost,
                  const GetPlacement& get_placement,
                  const ShouldIgnore& should_ignore,
                  Visitor visitor)
{
  typedef EdgeCollapse<TM, GT, ShouldStop,
                       VertexIndexMap, VertexPointMap, HalfedgeIndexMap, EdgeIsConstrainedMap,
                       GetCost, GetPlacement, ShouldIgnore, Visitor> Algorithm;

  Algorithm algorithm(tmesh, traits, should_stop, vim, vpm, him, ecm, get_cost, get_placement, should_ignore, visitor);

  return algorithm.run();
}

run()实现函数，重点关注collect()和loop()

template<class TM, class GT, class SP, class VIM, class VPM,class HIM, class ECM, class CF, class PF, class SI, class V>
int
EdgeCollapse<TM,GT,SP,VIM,VPM,HIM,ECM,CF,PF,SI,V>::
run()
{
  CGAL_expensive_precondition(is_valid_polygon_mesh(m_tm) && CGAL::is_triangle_mesh(m_tm));

  m_visitor.OnStarted(m_tm);

  // this is similar to the visitor, but for the cost/stop/placement oracles
  internal::Oracles_initializer<Self>(*this)();

  // First collect all candidate edges in a PQ
  collect();

  // Then proceed to collapse each edge in turn
  loop();

  CGAL_SMS_TRACE(0, "Finished: " << (m_initial_edge_count - m_current_edge_count) << " edges removed.");

  int r = int(m_initial_edge_count - m_current_edge_count);

  m_visitor.OnFinished(m_tm);

  return r;
}

1、数据结构 

它使用C++模板编程来定义多种类型，这些类型用于处理网格数据结构中的不同元素和操作。下面是每个类型定义的解释：

1. **基本网格和几何类型**：
   – `Triangle_mesh`: 表示三角形网格的类型。
   – `Geom_traits`: 几何特征类，包含用于处理网格的几何计算的函数和类型。
   – `Should_ignore`: 决定是否应该忽略某些边的策略类型。
   – `Vertex_index_map`: 顶点索引映射类型，用于访问顶点的索引。
   – `Vertex_point_map`: 顶点点映射类型，用于访问顶点的位置。
   – `Halfedge_index_map`: 半边索引映射类型，用于访问半边的索引。
   – `Edge_is_constrained_map`: 表示边是否受约束的映射类型。
   – `Get_cost`: 计算边塌陷成本的函数类型。
   – `Get_placement`: 在边塌陷时确定新顶点位置的函数类型。
   – `Should_stop`: 决定何时停止边塌陷算法的策略类型。
   – `Visitor`: 用于在算法的不同阶段进行通知和记录的访问者类型。

2. **边塌陷算法的专用类型**：
   – `Profile`: 存储有关特定边的信息（如顶点和相邻面）的类型。
   – `Graph_traits`: 表示网格的图形特征的类型。

3. **图形元素描述符和迭代器**：
   – `vertex_descriptor`: 顶点的描述符。
   – `vertex_iterator`: 遍历顶点的迭代器。
   – `halfedge_descriptor`: 半边的描述符。
   – `halfedge_iterator`: 遍历半边的迭代器。
   – `edge_descriptor`: 边的描述符。
   – `edge_iterator`: 遍历边的迭代器。
   – `out_edge_iterator`: 遍历某一顶点发出的边的迭代器。

4. **其他基本类型**：
   – `size_type`: 表示大小或数量的类型。
   – `Point`: 表示点的类型。
   – `FT`: 几何特征中用于表示数值的类型（通常是浮点类型）。
   – `Vector`: 表示向量的类型。
   – `Equal_3`: 用于比较三维点或向量是否相等的函数类型。

5. **可选类型**：
   – `Cost_type`: 表示边塌陷成本的可选类型。
   – `Placement_type`: 表示边塌陷后新顶点位置的可选类型。

2、collect()

这个collect()方法是CGAL中EdgeCollapse类的一个成员函数，用于收集网格中所有可以进行塌陷操作的边，并将它们添加到优先队列中

• 方法首先初始化一些局部变量来跟踪边的数量和边数据。
• 它遍历网格中的所有边，并对每条边创建一个Profile对象。这个Profile对象包含关于边的信息，如它的顶点和相邻面。
• 对于每条边，如果它不是受约束的边（即可以塌陷的边），它会计算这条边的塌陷成本，并将其添加到优先队列中。
• 如果边的两个端点位置相同（即边的长度为零），则将这条边加入到一个特殊的集合中以后续处理。
• 最后，对于长度为零的边，它会执行一些特殊的处理，包括从优先队列中移除相关边，并执行塌陷操作。

template<class TM, class GT, class SP, class VIM, class VPM,class HIM, class ECM, class CF, class PF, class SI, class V>
void
EdgeCollapse<TM,GT,SP,VIM,VPM,HIM,ECM,CF,PF,SI,V>::
collect()
{
  CGAL_SMS_TRACE(0, "collecting edges...");

  // loop over all the _undirected_ edges in the surface putting them in the PQ
  //获取网格m_tm中边的总数，并设置初始边数（m_initial_edge_count）和当前边数（m_current_edge_count）。
  const size_type ne = num_edges(m_tm); // if the mesh has garbage, you might have "ne > edges(tm).size()"
  m_initial_edge_count = m_current_edge_count = size_type(edges(m_tm).size());

  //分配一个新的Edge_data数组，用于存储每条边的数据（如成本）。
  //初始化优先队列mPQ，该队列将根据成本比较函数Compare_cost对边进行排序。
  m_edge_data.reset(new Edge_data[ne]);
  mPQ.reset(new PQ(ne, Compare_cost(this), edge_id(this)));

  //初始化用于调试的计数器，分别记录插入和未插入优先队列的边的数量。
  CGAL_assertion_code(size_type num_inserted = 0);
  CGAL_assertion_code(size_type num_not_inserted = 0);

  //初始化一个集合，用于存储长度为零的边
  std::set<halfedge_descriptor> zero_length_edges;

  //遍历网格中的每一条边。
  //如果一条边是受约束的（如边界边或不可塌陷的边），则跳过该边并增加未插入计数。
  for(edge_descriptor e : edges(m_tm))
  {
    //e是一个edge_descriptor，它是对应于网格m_tm中一条边的唯一标识符。
    //halfedge(e, m_tm)函数调用返回与边e关联的一个半边的描述符。
    const halfedge_descriptor h = halfedge(e, m_tm);

    if(is_constrained(h))
    {
      CGAL_assertion_code(++num_not_inserted);
      continue; // no not insert constrainted edges
    }


    //为每条可处理的边创建一个Profile对象，该对象包含关于边的详细信息。
    //计算边的塌陷成本，并将其加入优先队列中。
   //如果边的两端点相同（长度为零），则将其添加到zero_length_edges集合中。
    const Profile profile = create_profile(h);
    if(!m_traits.equal_3_object()(profile.p0(), profile.p1()))
    {
      Edge_data& data = get_data(h);

      data.cost() = cost(profile);
      insert_in_PQ(h, data);

      m_visitor.OnCollected(profile, data.cost());

      CGAL_assertion_code(++num_inserted);
    }
    else
    {
      zero_length_edges.insert(primary_edge(h));
      CGAL_assertion_code(++num_not_inserted);
    }

    //调试输出，显示当前处理的边的信息。
    CGAL_SMS_TRACE(2, edge_to_string(h));
  }

  //断言以确保所有的边都被处理（无论是插入还是未插入）。
  CGAL_assertion(num_inserted + num_not_inserted == m_initial_edge_count);

  //遍历zero_length_edges集合中的边，并对每条边进行额外处理。
  for(halfedge_descriptor hd : zero_length_edges)
  {
    //对于每条长度为零的边，创建Profile对象，并检查其塌陷是否在拓扑上有效。如果不有效，则跳过。
    const Profile profile = create_profile(hd);

    if(!is_collapse_topologically_valid(profile))
      continue;

    // edges of length 0 removed no longer need to be treated
    if(profile.left_face_exists())
    {
      halfedge_descriptor h_to_remove = is_constrained(profile.vL_v0()) ?
                                          primary_edge(profile.v1_vL()) :
                                          primary_edge(profile.vL_v0());
      zero_length_edges.erase(h_to_remove);
      Edge_data& data = get_data(h_to_remove);

      if(data.is_in_PQ())
      {
        CGAL_SMS_TRACE(2, "Removing E" << get_edge_id(h_to_remove) << " from PQ");
        remove_from_PQ(h_to_remove, data);
      }

      --m_current_edge_count;
    }

    if(profile.right_face_exists())
    {
      halfedge_descriptor h_to_remove = is_constrained(profile.vR_v1()) ?
                                          primary_edge(profile.v0_vR()) :
                                          primary_edge(profile.vR_v1());
      zero_length_edges.erase(h_to_remove);
      Edge_data& data = get_data(h_to_remove);

      if(data.is_in_PQ())
      {
        CGAL_SMS_TRACE(2, "Removing E" << get_edge_id(h_to_remove) << " from PQ");
        remove_from_PQ(h_to_remove, data);
      }

      --m_current_edge_count;
    }

    --m_current_edge_count;

    //执行实际的塌陷操作，将顶点折叠到长度为零的边的一个端点上。
    //更新顶点位置，并调用访问者和Oracle更新器以处理塌陷后的变化。
    //the placement is trivial, it's always the point itself
    Placement_type placement = profile.p0();
    vertex_descriptor v = halfedge_collapse_bk_compatibility(profile.v0_v1(), m_ecm);
    put(m_vpm, v, *placement);

    m_visitor.OnCollapsed(profile, v);
    internal::After_collapse_oracles_updater<Self>(*this)(profile, v);
  }

  CGAL_SMS_TRACE(0, "Initial edge count: " << m_initial_edge_count);
}

 3、is_collapse_topologically_valid()

这个collect中的方法，用于判断边长为0的边进行塌陷是否满足拓扑结构

  template<class TM, class GT, class SP, class VIM, class VPM, class HIM, class ECM, class CF, class PF, class SI, class V>
bool
  EdgeCollapse<TM,GT,SP,VIM,VPM,HIM,ECM,CF,PF,SI,V>::
is_collapse_topologically_valid(const Profile& profile)
{
  //初始化返回值res为true，假设边的塌陷在拓扑上是有效的，直到检测到相反的情况
  bool res = true;

  //用于输出正在测试的边的信息。
  CGAL_SMS_TRACE(3,"Testing topological collapsabilty of p_q=V" << get(m_vim,profile.v0()) << "(%" << degree(profile.v0(), m_tm) << ")"
                 << "->V" << get(m_vim, profile.v1()) << "(%" << degree(profile.v1(), m_tm) << ")");

  //输出边是否为边界边的调试信息。
  CGAL_SMS_TRACE(4, "is p_q border:" << profile.is_v0_v1_a_border());
  CGAL_SMS_TRACE(4, "is q_q border:" << profile.is_v1_v0_a_border());

  //定义外围边迭代器，用于遍历与顶点相关联的边。
  out_edge_iterator eb1, ee1;
  out_edge_iterator eb2, ee2;

  //输出边的左右面（如果存在）的相关信息。
  CGAL_SMS_TRACE(4, "  t=V" << (profile.left_face_exists() ? get(m_vim,profile.vL()) : -1)
                            << "(%" << (profile.left_face_exists() ? degree(profile.vL(), m_tm) : 0) << ")");
  CGAL_SMS_TRACE(4, "  b=V" << (profile.right_face_exists() ? get(m_vim,profile.vR()) : -1)
                            << "(%" << (profile.right_face_exists() ? degree(profile.vR(), m_tm) :0) << ")");

  //检查边的两个面是否存在。如果左面或右面不存在，则边可能是边界边或者孤立边。
  // Simple tests handling the case of non-manifold situations at a vertex or edge (pinching)
  // (even if we advertise one should not use a surface mesh with such features)
  if(profile.left_face_exists())
  {
    if(CGAL::is_border(opposite(profile.v1_vL(), m_tm), m_tm) &&
       CGAL::is_border(opposite(profile.vL_v0(), m_tm), m_tm))
      return false;

    if(profile.right_face_exists() &&
       CGAL::is_border(opposite(profile.vR_v1(), m_tm), m_tm) &&
       CGAL::is_border(opposite(profile.v0_vR(), m_tm), m_tm))
      return false;
  }
  else
  {
    if(profile.right_face_exists())
    {
      if(CGAL::is_border(opposite(profile.vR_v1(), m_tm), m_tm) &&
         CGAL::is_border(opposite(profile.v0_vR(), m_tm), m_tm))
        return false;
    }
    else
      return false;
  }


  //遍历顶点v0周围的半边，检查这些半边是否满足所谓的“链接条件”。
  //链接条件是指，对于每个与v0和v1都相邻的顶点k，必须存在一个面包含v0、v1和k。
  // The following loop checks the link condition for v0_v1.
  // Specifically, that for every vertex 'k' adjacent to both 'p and 'q', 'pkq' is a face of the mesh.
  //
  for(boost::tie(eb1,ee1) = halfedges_around_source(profile.v0(), m_tm); res && eb1 != ee1; ++eb1)
  {
    halfedge_descriptor v0_k = *eb1;

    if(v0_k != profile.v0_v1())
    {
      vertex_descriptor k = target(v0_k, m_tm);

      for(boost::tie(eb2,ee2) = halfedges_around_source(k, m_tm); res && eb2 != ee2; ++eb2)
      {
        halfedge_descriptor k_v1 = *eb2;

        if(target(k_v1, m_tm) == profile.v1())
        {
          // At this point we know p-q-k are connected and we need to determine if this triangle is a face of the mesh.
          //
          // Since the mesh is known to be triangular there are at most two faces sharing the edge p-q.
          //
          // If p->q is NOT a border edge, the top face is p->q->t where t is target(next(p->q))
          // If q->p is NOT a border edge, the bottom face is q->p->b where b is target(next(q->p))
          //
          // If k is either t or b then p-q-k *might* be a face of the mesh. It won't be if k==t but p->q is border
          // or k==b but q->b is a border (because in that case even though there exists triangles p->q->t (or q->p->b)
          // they are holes, not faces)
          //
         //检查顶点k是否与v0和v1共享一个面。
          bool is_face = (profile.vL() == k && profile.left_face_exists()) ||
                         (profile.vR() == k && profile.right_face_exists());

          CGAL_assertion_code(
            if(is_face)
            {
              // Is k_v1 the halfedge bounding the face 'k-v1-v0'?
              if(!is_border(k_v1, m_tm) && target(next(k_v1, m_tm), m_tm) == profile.v0())
              {
                CGAL_assertion(target(k_v1, m_tm) == profile.v1());
                CGAL_assertion(target(next(k_v1, m_tm), m_tm) == profile.v0());
                CGAL_assertion(target(next(next(k_v1, m_tm), m_tm), m_tm) == k);
              }
              else // or is it the opposite?
              {
                halfedge_descriptor v1_k = opposite(k_v1, m_tm);
                CGAL_assertion(!is_border(v1_k, m_tm));
                CGAL_assertion(target(v1_k, m_tm) == k);
                CGAL_assertion(target(next(v1_k, m_tm), m_tm) == profile.v0());
                CGAL_assertion(target(next(next(v1_k, m_tm), m_tm), m_tm) == profile.v1());
              }
            }
          );

          if(!is_face)
          {
            CGAL_SMS_TRACE(3,"  k=V" << get(m_vim,k) << " IS NOT in a face with p-q. NON-COLLAPSABLE edge.");
            res = false;
            break;
          }
          else
          {
            CGAL_SMS_TRACE(4,"  k=V" << get(m_vim,k) << " is in a face with p-q");
          }
        }
      }
    }
  }

  if(res)
  {
    //检查边是否与受约束的边相邻。如果是，这可能会阻止边的塌陷。
    /// ensure two constrained edges cannot get merged
    if(is_edge_adjacent_to_a_constrained_edge(profile, m_ecm))
      return false;

    //检查边是否是边界边，以及是否满足其他拓扑约束条件。
    if(profile.is_v0_v1_a_border())
    {
      if(is_open_triangle(profile.v0_v1()))
      {
        res = false;
        CGAL_SMS_TRACE(3,"  p-q belongs to an open triangle. NON-COLLAPSABLE edge.");
      }
    }
    else if(profile.is_v1_v0_a_border())
    {
      if(is_open_triangle(profile.v1_v0()))
      {
        res = false;
        CGAL_SMS_TRACE(3,"  p-q belongs to an open triangle. NON-COLLAPSABLE edge.");
      }
    }
    else
    {
      if(is_border(profile.v0(), m_tm) && is_border(profile.v1(), m_tm))
      {
        res = false;
        CGAL_SMS_TRACE(3,"  both p and q are boundary vertices but p-q is not. NON-COLLAPSABLE edge.");
      }
      else
      {
        bool tetra = is_tetrahedron(profile.v0_v1());

        //CGAL_assertion(tetra == m_tm.is_tetrahedron(profile.v0_v1()));

        if(tetra)
        {
          res = false;
          CGAL_SMS_TRACE(3,"  p-q belongs to a tetrahedron. NON-COLLAPSABLE edge.");
        }

        if(next(profile.v0_v1(), m_tm) == opposite(prev(profile.v1_v0(), m_tm), m_tm) &&
           prev(profile.v0_v1(), m_tm) == opposite(next(profile.v1_v0(), m_tm), m_tm))
        {
          CGAL_SMS_TRACE(3,"  degenerate volume.");
          return false;
        }
      }
    }
  }

  return res;
}

 2、loop()

这部分代码比较容易看懂，主要函数就是m_should_stop(*cost, profile, m_initial_edge_count, m_current_edge_count)

这个函数主要用于从一个优先队列中不断取出边，并根据一定的条件和成本判断是否对这些边进行折叠操作，如何可以折叠那么是否在拓扑和几何上有效。从而对模型进行简化

// 在EdgeCollapse类中定义一个模板函数。
// 类的模板参数包括TM, GT, SP, VIM, VPM, HIM, ECM, CF, PF, SI, V。
template<class TM, class GT, class SP, class VIM, class VPM, class HIM, class ECM, class CF, class PF, class SI, class V>
void
EdgeCollapse<TM,GT,SP,VIM,VPM,HIM,ECM,CF,PF,SI,V>::
loop()
{
  // 输出日志信息，表示开始折叠边
  CGAL_SMS_TRACE(0, "Collapsing edges...");

  // 断言代码，用于调试时验证条件。此处检查不可折叠计数。
  CGAL_assertion_code(size_type non_collapsable_count = 0);

  // 从优先队列中弹出并处理每条边
  boost::optional<halfedge_descriptor> opt_h;

  // 如果定义了CGAL_SURF_SIMPL_INTERMEDIATE_STEPS_PRINTING，则启用中间步骤打印
#ifdef CGAL_SURF_SIMPL_INTERMEDIATE_STEPS_PRINTING
  int i_rm = 0;
#endif

  // 当优先队列非空时循环
  while((opt_h = pop_from_PQ()))
  {
    // 输出日志信息，显示弹出的边
    CGAL_SMS_TRACE(1, "Popped " << edge_to_string(*opt_h));
    // 断言检查边是否不受约束
    CGAL_assertion(!is_constrained(*opt_h));

    // 创建边的Profile
    const Profile profile = create_profile(*opt_h);
    // 获取边的成本
    Cost_type cost = get_data(*opt_h).cost();

    // 调用访问者模式的函数
    m_visitor.OnSelected(profile, cost, m_initial_edge_count, m_current_edge_count);

    // 如果成本存在
    if(cost)
    {
      // 如果满足停止条件
      if(m_should_stop(*cost, profile, m_initial_edge_count, m_current_edge_count))
      {
        // 调用访问者模式的函数，表示达到停止条件
        m_visitor.OnStopConditionReached(profile);

        // 输出日志信息，显示达到停止条件的信息
        CGAL_SMS_TRACE(0, "Stop condition reached with initial edge count=" << m_initial_edge_count
                            << " current edge count=" << m_current_edge_count
                            << " current edge: " << edge_to_string(*opt_h));
        break;
      }

      // 如果折叠在拓扑上有效
      if(is_collapse_topologically_valid(profile))
      {
        // 获取放置位置
        Placement_type placement = get_placement(profile);

        // 如果折叠在几何上有效
        if(is_collapse_geometrically_valid(profile, placement))
        {
          // 如果定义了CGAL_SURF_SIMPL_INTERMEDIATE_STEPS_PRINTING，则打印步骤信息
#ifdef CGAL_SURF_SIMPL_INTERMEDIATE_STEPS_PRINTING
          std::cout << "step " << i_rm << " " << get(m_vpm, source(*h, m_tm))
                                       << " " << get(m_vpm, target(*h, m_tm)) << "n";
#endif
          // 如果不应该忽略此折叠
          if(m_should_ignore(profile, placement)!= boost::none){
            // 执行折叠操作
            collapse(profile, placement);
          }
          else
          {
            // 如果不能折叠，则增加不可折叠计数
            CGAL_assertion_code(++non_collapsable_count);

            // 调用访问者模式的函数，表示边不可折叠
            m_visitor.OnNonCollapsable(profile);

            // 输出日志信息，显示边不可折叠
            CGAL_SMS_TRACE(1, edge_to_string(*opt_h) << " NOT Collapsable" );
          }

          // 如果定义了CGAL_SURF_SIMPL_INTERMEDIATE_STEPS_PRINTING，则保存中间步骤的文件
#ifdef CGAL_SURF_SIMPL_INTERMEDIATE_STEPS_PRINTING
          std::stringstream sstr;
          sstr << "debug/P-";
          if(i_rm<10) sstr << "0";
          if(i_rm<100) sstr << "0";
          sstr << i_rm << ".off";
          std::ofstream out(sstr.str().c_str());
          out << m_tm;
          ++i_rm;
#endif
        }
      }
      else
      {
        // 如果不能折叠，则增加不可折叠计数
        CGAL_assertion_code(++non_collapsable_count);

        // 调用访问者模式的函数，表示边不可折叠
        m_visitor.OnNonCollapsable(profile);

        // 输出日志信息，显示边不可折叠
        CGAL_SMS_TRACE(1, edge_to_string(*opt_h) << " NOT Collapsable" );
      }
    }
    else
    {
      // 输出日志信息，显示成本不可计算
      CGAL_SMS_TRACE(1, edge_to_string(*opt_h) << " uncomputable cost." );
    }
  }
}

代码007普通

打赏 收藏 海报 链接