clusters.h

很多二维 三维几何计算算法 C++ 类库

update_cluster(Cluster& c, iterator it, Vertex_handle va,
               Vertex_handle vb, Vertex_handle vm, bool reduction)
{
  typename Geom_traits::Compute_squared_distance_2 squared_distance =
    tr.geom_traits().compute_squared_distance_2_object();

  cluster_map.erase(it);

  c.vertices.erase(vb);
  c.vertices[vm] = reduction;

  if(vb==c.smallest_angle.first)
    c.smallest_angle.first = vm;
  if(vb==c.smallest_angle.second)
    c.smallest_angle.second = vm;

  FT l = squared_distance(va->point(),vm->point());
  if(l<c.minimum_squared_length)
    c.minimum_squared_length = l;

  if(!c.is_reduced())
    {
      typename Cluster::Vertices_map::iterator it = c.vertices.begin();
      while(it!=c.vertices.end() && c.is_reduced(it->first))
        ++it; // @todo: use std::find and an object class
      if(it==c.vertices.end())
        c.reduced = true;
    }

  if(c.is_reduced())
    c.rmin = squared_distance(c.smallest_angle.first->point(),
                              c.smallest_angle.second->point())/FT(4);
  cluster_map.insert(Cluster_map_value_type(va,c));
}

template <typename Tr>
bool Clusters<Tr>::
get_cluster(Vertex_handle va, Vertex_handle vb, Cluster& c,
            const_iterator& it) const
{
  typedef std::pair<const_iterator, const_iterator> Range;

  Range range = cluster_map.equal_range(va);

  for(it = range.first; it != range.second; it++)
    {
      const Cluster &cl = it->second;
      if(cl.vertices.find(vb)!=cl.vertices.end()) {
        c = it->second;
        return true;
      }
    }
  return false;
}

template <typename Tr>
bool Clusters<Tr>::
get_cluster(Vertex_handle va, Vertex_handle vb, Cluster& c,
            iterator& it) 
{
  typedef std::pair<iterator, iterator> Range;

  Range range = cluster_map.equal_range(va);

  for(it = range.first; it != range.second; it++)
    {
      const Cluster &cl = it->second;
      if(cl.vertices.find(vb)!=cl.vertices.end()) {
        c = it->second;
        return true;
      }
    }
  return false;
}

template <typename Tr>
void Clusters<Tr>::
create_clusters()
{
  cluster_map.clear();
#ifndef CGAL_IT_IS_A_CONSTRAINED_TRIANGULATION_PLUS
  for(typename Tr::Finite_vertices_iterator vit = tr.finite_vertices_begin();
      vit != tr.finite_vertices_end();
      vit++)
    create_clusters_of_vertex(vit);
#else
  Unique_hash_map<Vertex_handle,bool> created(false);
  for(typename Tr::Subconstraint_iterator it = tr.subconstraints_begin();
      it != tr.subconstraints_end(); ++it) {
    Vertex_handle vh = it->first.first;
    if(!created[vh]){
      created[vh] = true;
      create_clusters_of_vertex(vh);
    }

    vh = it->first.second;
    if(!created[vh]){
      created[vh] = true;
      create_clusters_of_vertex(vh);
    }
  }
#endif
}

template <typename Tr>
void Clusters<Tr>::
create_clusters_of_vertex(const Vertex_handle v)
{
  details::Is_edge_constrained<Tr> test(tr);

  Constrained_edge_circulator begin(tr.incident_edges(v),test);

  // This circulator represents all constrained edges around the
  // vertex v. An edge [v,v'] is represented by the vertex v'.

  if(begin == 0) return; // if there is only one vertex

  Constrained_edge_circulator
    current(begin), next(begin), cluster_begin(begin);
  ++next; // next is always just after current.
  if(current == next) return;

  bool in_a_cluster = false;
  do
    {
      if(is_small_angle(target(current)->point(), v->point(),
                        target(next)->point()))
        {
          if(!in_a_cluster)
            {
              // at this point, current is the beginning of a cluster
              in_a_cluster = true;
              cluster_begin = current;
            }
        }
      else
        if(in_a_cluster)
          {
            // at this point, current is the end of a cluster and
            // cluster_begin is its beginning
            construct_cluster(v, cluster_begin, current);
            in_a_cluster = false;
          }
      ++next;
      ++current;
    } while( current!=begin );
  if(in_a_cluster)
    {
      Cluster c;
      iterator it;
      if(get_cluster(v, target(begin), c, it))
        {
          // get the cluster and erase it from the clusters map
          cluster_map.erase(it);
          construct_cluster(v, cluster_begin, begin, c);
        }
      else
        construct_cluster(v, cluster_begin, current);
    }
}

template <typename Tr>
void Clusters<Tr>::
construct_cluster(Vertex_handle v,
                  Constrained_edge_circulator begin,
                  const Constrained_edge_circulator& end,
                  Cluster c)
{
  typename Geom_traits::Compute_squared_distance_2 squared_distance =
    tr.geom_traits().compute_squared_distance_2_object();

  if(c.vertices.empty())
    {
      c.reduced = false;
      // c.rmin is not initialized because
      // reduced=false!
      c.minimum_squared_length =
        squared_distance(v->point(), target(begin)->point());
      Constrained_edge_circulator second(begin);
      ++second;
      c.smallest_angle.first = target(begin);
      c.smallest_angle.second = target(second);
    }

  bool all_edges_in_cluster=false; // tell if all incident edges are
  // in the cluster
  if(begin==end)
    all_edges_in_cluster=true;

  const Point& vp = v->point();

  FT greatest_cosine =
    squared_cosine_of_angle_times_4(c.smallest_angle.first->point(),
                                    v->point(),
                                    c.smallest_angle.second->point());

  Constrained_edge_circulator next(begin);
  ++next;
  do
    {
      c.vertices[target(begin)] = false;
      Squared_length l = squared_distance(vp,
                                        target(begin)->point());
      c.minimum_squared_length =
        (std::min)(l,c.minimum_squared_length);

      if(all_edges_in_cluster || begin!=end)
        {
          FT cosine =
            squared_cosine_of_angle_times_4(target(begin)->point(),
                                            v->point(),
                                            target(next)->point());
          if(cosine>greatest_cosine)
            {
              greatest_cosine = cosine;
              c.smallest_angle.first = target(begin);
              c.smallest_angle.second = target(next);
            }
        }
    }
  while(next++,begin++!=end);
  typedef typename Cluster_map::value_type Value_key_pair;
  cluster_map.insert(Value_key_pair(v,c));
}

template <typename Tr>
bool Clusters<Tr>::
is_small_angle(const Point& pleft,
               const Point& pmiddle,
               const Point& pright) const
{
  typename Geom_traits::Angle_2 angle = 
    tr.geom_traits().angle_2_object();
  typename Geom_traits::Orientation_2 orient =
    tr.geom_traits().orientation_2_object();

  if( angle(pleft, pmiddle, pright)==OBTUSE )
    return false;
  if( orient(pmiddle,pleft,pright)==RIGHT_TURN)
    return false;

  FT cos_alpha = squared_cosine_of_angle_times_4(pleft, pmiddle,
                                                 pright);

  if(cos_alpha > 1)
    {
      return true; //the same cluster
    }
  else
    {
      return false; //another cluster
    }
}

template <typename Tr>
typename Clusters<Tr>::FT
Clusters<Tr>::
squared_cosine_of_angle_times_4(const Point& pb, const Point& pa,
                                const Point& pc) const
{
  typename Geom_traits::Compute_squared_distance_2 squared_distance =
    tr.geom_traits().compute_squared_distance_2_object();

  const FT
    a = squared_distance(pb, pc),
    b = squared_distance(pa, pb),
    c = squared_distance(pa, pc);

  const FT num = a-(b+c);

  return (num*num)/(b*c);
}
  
} // end namespace Mesh_2

} // end namespace CGAL

#endif // CGAL_MESH_2_CLUSTERS_H