sm_triangulator.h

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

    From[--L.end()] = Object_handle(v);
  }
  SHalfedge_const_iterator e;
  CGAL_forall_sedges(e,*E_) {
    if ( e->source() == e->twin()->source() ) {
      Seg_pair p = two_segments(E_,e);
      L.push_back(p.first); L.push_back(p.second);
      From[--L.end()] = From[--(--L.end())] = Object_handle(e);
    } else {
      L.push_back(segment(E_,e));
      From[--L.end()] = Object_handle(e);
    }
  }
  if ( E_->has_shalfloop() ) {
    Seg_pair p = two_segments(E_,E_->shalfloop());
    L.push_back(p.first); L.push_back(p.second);
    From[--L.end()] = From[--(--L.end())] = Object_handle(E_->shalfloop());
  }

  // partition segments from L to positive and negative hemisphere
  Seg_list L_pos,L_neg;
  partition_to_halfsphere(L.begin(), L.end(), L_pos, From, +1);
  partition_to_halfsphere(L.begin(), L.end(), L_neg, From, -1);

  //  typename Seg_list::iterator it;
  //    std::cerr << "L_pos" << std::endl;
  //    CGAL_forall_iterators(it,L_pos) std::cerr << *it << std::endl;
  //    std::cerr << "L_neg" << std::endl;
  //    CGAL_forall_iterators(it,L_neg) std::cerr << *it << std::endl;

  // sweep the hemispheres to create two half sphere maps
  typedef SM_subdivision<Self,Seg_iterator,Object_handle> SM_output;
  typedef typename Sphere_kernel::Positive_halfsphere_geometry PH_geometry;
  typedef CGAL::Segment_overlay_traits< 
            Seg_iterator, SM_output, PH_geometry>  PHS_traits;
  typedef CGAL::generic_sweep<PHS_traits> Positive_halfsphere_sweep;

  typedef typename Sphere_kernel::Negative_halfsphere_geometry NH_geometry;
  typedef CGAL::Segment_overlay_traits< 
            Seg_iterator, SM_output, NH_geometry> NHS_traits;
  typedef CGAL::generic_sweep<NHS_traits> Negative_halfsphere_sweep;

  SVertex_handle v_sep;
  SHalfedge_handle e_sep;
  SM_output O(*this,From); 

  typedef typename PHS_traits::INPUT Input_range;
  Positive_halfsphere_sweep SP(
    Input_range(L_pos.begin(),L_pos.end()),O,
    PH_geometry());
  SP.sweep();
  v_sep=--this->svertices_end(); e_sep=--this->shalfedges_end();

  Negative_halfsphere_sweep SM(
    Input_range(L_neg.begin(),L_neg.end()),O,
    NH_geometry());
  SM.sweep();
  ++v_sep; ++e_sep;
  // now two CCs of sphere graph are calculated
  // v_sep = first vertex of CC in negative x-sphere
  // e_sep = first edge of CC in negative x-sphere
   
  SHalfedge_iterator u;
  CGAL_forall_sedges(u,*this) {
    Sphere_segment s(u->source()->point(),u->twin()->source()->point());
    u->circle() = s.sphere_circle(); 
    u->twin()->circle() = s.sphere_circle().opposite();
  }

  Mark lower, upper;
  SM_point_locator PL(E_->sphere_map());
  PL.marks_of_halfspheres(lower,upper);
  complete_support(this->svertices_begin(), v_sep, lower);
  complete_support(v_sep, this->svertices_end(), upper);

  /*
  CGAL_forall_sedges(u,*this) {
    std::cerr << point(source(u)) << "->" << point(target(u)) << std::endl;
  }
  */

  // triangulate per hemisphere
  typedef SM_constrained_triang_traits<Self,PH_geometry>  PCT_traits;
  typedef CGAL::generic_sweep<PCT_traits> Positive_halfsphere_ct_sweep;
  typedef SM_constrained_triang_traits<Self,NH_geometry>  NCT_traits;
  typedef CGAL::generic_sweep<NCT_traits> Negative_halfsphere_ct_sweep;
  typedef std::pair<SVertex_iterator,SVertex_iterator> SVertex_pair;

  SVertex_pair vpp(this->svertices_begin(),v_sep);
  Positive_halfsphere_ct_sweep PCTS(vpp, *this,
    PH_geometry());
  PCTS.sweep();
  SVertex_pair vpn(v_sep,this->svertices_end());
  Negative_halfsphere_ct_sweep NCTS(vpn, *this,
    NH_geometry());
  NCTS.sweep();

  /*
  std::cerr << std::endl;
  CGAL_forall_sedges(u,*this) {
    std::cerr << point(source(u)) << "->" << point(target(u)) << std::endl;
  }
  */

  /* Note the we divide the world along the xy equator and 
     split the equator at y- and y+. We treat the halfcircle
     at x+ as if perturbed slightly up. This makes triangles
     that have y- or y+ as a vertex degenerate. if such triangles
     appear we repair it by flipping the edge opposite to the
     vertex y-(y+).
  */

  correct_triangle_at(this->svertices_begin());
  correct_triangle_at(--SVertex_iterator(v_sep));
  correct_triangle_at(v_sep);
  correct_triangle_at(--this->svertices_end());

  CGAL_forall_sedges(u,*this) {
    Sphere_segment s(u->source()->point(),u->twin()->source()->point());
    u->circle() = s.sphere_circle();
    u->twin()->circle() = s.sphere_circle().opposite();
  }

  // merge the hemisphere maps into one sphere map
  merge_halfsphere_maps(this->svertices_begin(),v_sep);
  this->check_integrity_and_topological_planarity(false);
}


template <typename Decorator_>
template <typename Iterator, typename T>
void SM_triangulator<Decorator_>::
partition_to_halfsphere(Iterator start, Iterator beyond, Seg_list& L, 
  CGAL::Unique_hash_map<Iterator,T>& M, int pos) const
{ CGAL_NEF_TRACEN("partition_to_halfsphere ");
  CGAL_assertion(pos!=0);
  bool add_cross = true;
  Sphere_segment s1,s2;
  Sphere_circle xycircle(0,0,pos);
  while ( start != beyond ) { 
    if(start->source().hz() * pos > 0 || start->target().hz() * pos > 0)
      add_cross = false;
    int i = start->intersection(xycircle,s1,s2);
    if (i>1) { L.push_back(s2); M[--L.end()] = M[start]; }
    if (i>0) { L.push_back(s1); M[--L.end()] = M[start]; }
    ++start;
  }
  // now all segments are split into halfspheres
  // we still have to:
  // - split segments containing our special poles y^-, y^+
  // - split halfcircles
  // - add four equator segments 
  Sphere_point S(0,-1,0),N(0,1,0);
  Sphere_circle yzcircle(1,0,0);
  typename Seg_list::iterator it, itl;

  bool part_in_hemisphere(false);
  CGAL_forall_iterators(it,L) { CGAL_NEF_TRACEN("  "<<*it);
    if ( equal_as_sets(it->sphere_circle(),xycircle) ) {
      CGAL_NEF_TRACEN("  splitting xy seg "<<*it);
      int n1 =  it->intersection(yzcircle,s1,s2);
      if (n1 > 1 && !s2.is_degenerate()) 
      { M[ L.insert(it,s2) ] = M[it]; }
      if (n1 > 0 && !s1.is_degenerate()) 
      { M[ L.insert(it,s1) ] = M[it]; }
      int n2 =  it->intersection(yzcircle.opposite(),s1,s2);
      if (n2 > 1 && !s2.is_degenerate()) 
      { M[ L.insert(it,s2) ] = M[it]; }
      if (n2 > 0 && !s1.is_degenerate()) 
      { M[ L.insert(it,s1) ] = M[it]; }
      itl = it; --it; L.erase(itl); M[itl] = T();
      // at least one item was appended
    } else {
      part_in_hemisphere = true;
    }
  }
  CGAL_forall_iterators(it,L) {
    if ( it->is_halfcircle() ) {
      CGAL_NEF_TRACEN("  splitting halfcircle "<<*it);
      Sphere_segment s1,s2;
      it->split_halfcircle(s1,s2);
      *it = s2; 
      M[ L.insert(it,s1) ] = M[it];
    }
  }
  // append 4 xy-equator segments:
  Sphere_segment sp(S,N,xycircle);
  Sphere_segment sm(S,N,xycircle.opposite());
  Sphere_segment s[4];
  sp.split_halfcircle(s[0],s[1]);
  sm.split_halfcircle(s[2],s[3]);
  L.insert(L.end(),s,s+4);
  /* if no segment is covering the interior of the hemisphere
     we have to add a trivial segment to allow for a correct
     triangulation */
  if ( !part_in_hemisphere || add_cross) {
    Sphere_point p(0,0,pos);
    Sphere_circle c(1,0,0);
    L.push_back(Sphere_segment(p,p,c));
  }
}

template <typename Decorator_>
void SM_triangulator<Decorator_>::
merge_nodes(SHalfedge_handle e1, SHalfedge_handle e2)
{
  SVertex_handle v1 = e1->source(), v2 = e2->twin()->source();
  CGAL_NEF_TRACEN("merge_nodes "<<PH(v1)<<PH(v2));
  CGAL_assertion(v1->point()==v2->point());
  SHalfedge_handle ep1 = e1->sprev(), en2 = e2->snext();
  SHalfedge_around_svertex_circulator eav(out_edges(v2)),ee(eav);
  CGAL_For_all(eav,ee) { set_source(eav,v1); }
  link_as_prev_next_pair(e2,e1);  
  link_as_prev_next_pair(ep1,en2); 
  assert_equal_marks(v1,v2);
  discard_info(v2);
  delete_vertex_only(v2);
}


template <typename Decorator_>
void SM_triangulator<Decorator_>::
merge_halfsphere_maps(SVertex_handle v1, SVertex_handle v2)
{ CGAL_NEF_TRACEN("merging halfspheres "<<PH(v1)<<PH(v2));
  CGAL_assertion(v1->point()==v2->point());
  std::list<SHalfedge_pair> L_equator;
  SHalfedge_around_sface_circulator 
    ep(last_out_edge(v1)), en(first_out_edge(v2)->twin());
  do { 
   L_equator.push_back(SHalfedge_pair(ep,en));
   merge_nodes(ep,en); ++ep; --en; 
  } while ( ep->source() != v1 );
  
  typename std::list<SHalfedge_pair>::iterator it;
  CGAL_forall_iterators(it,L_equator) { 
    SHalfedge_handle e1 = it->first, e2 = it->second;
    SHalfedge_handle e1t = e1->twin(), e2t = e2->twin();
    CGAL_NEF_TRACEV(PH(e1));CGAL_NEF_TRACEV(PH(e2));
    SHalfedge_handle e2tp = e2t->sprev();
    SHalfedge_handle e2tn = e2t->snext();
    link_as_prev_next_pair(e2tp,e1);
    link_as_prev_next_pair(e1,e2tn);
    SFace_handle f = e2t->incident_sface();
    if ( is_sm_boundary_object(e2t) )
    { undo_sm_boundary_object(e2t,f); store_sm_boundary_object(e1,f); }
    set_face(e1,f);
    if ( e2 == first_out_edge(e2->source()) )
      set_first_out_edge(e2->source(),e1t);
    discard_info(e2);
    delete_edge_pair_only(e2);
  }
}

template <typename Decorator_>
void SM_triangulator<Decorator_>::
complete_support(SVertex_iterator v_start, SVertex_iterator v_end,
		 Mark mohs) const
{ CGAL_NEF_TRACEN("complete_support");
  Mark m_buffer(mohs); 
  for (SVertex_iterator v = v_start; v != v_end; ++v) { 
    CGAL_NEF_TRACEN(" vertex = "<<PH(v));
    SHalfedge_handle e_below = halfedge_below(v);
    if ( v != v_start ) {
      if ( e_below != SHalfedge_handle() ) {
	m_buffer = incident_mark(e_below); 
      } else { // e_below does not exist
	/* this is only the case for a vertex v on the final equatorial
	   halfcircle; there we take the mark from an inedge edge into v */
	//	CGAL_assertion( point(v).z() == 0 && 
	//	  ( pos > 0 ? (point(v).x() >= 0) : (point(v).x()<=0)) );
	m_buffer = incident_mark(first_out_edge(v)->sprev());
      }
    }
    CGAL_NEF_TRACEN(" face mark below "<<m_buffer);

    Object_handle o = support(v);
    SVertex_const_handle vs;
    SHalfedge_const_handle es;
    SHalfloop_const_handle ls;
    if ( o == NULL ) { v->mark() = m_buffer; }
    else if ( CGAL::assign(vs,o) ) { v->mark() = vs->mark(); }
    else if ( CGAL::assign(es,o) ) {
      if ( es->source()->point() == v->point() ) 
	{ v->mark() = es->source()->mark(); }
      else if ( es->twin()->source()->point() == v->point() ) 
	{ v->mark() = es->twin()->source()->mark(); }
      else { v->mark() = es->mark(); }
    }
    else if ( CGAL::assign(ls,o) ) { v->mark() = ls->mark(); }
    else CGAL_assertion_msg(0,"damn wrong support.");
    CGAL_NEF_TRACEN(" face mark at "<<v->mark());

    if ( is_isolated(v) ) continue;

    SHalfedge_around_svertex_circulator e(first_out_edge(v)), hend(e);
    CGAL_For_all(e,hend) {
      CGAL_NEF_TRACEN("  edge "<<PH(e));
      if ( !is_forward(e) ) break;
      if ( support(e) != NULL ) {
        SHalfedge_const_handle ei;
        if ( CGAL::assign(ei,support(e)) ) { 
          if ( ei->circle() != e->circle() ) { ei = ei->twin(); }
          CGAL_assertion( ei->circle() == e->circle() ); 
          CGAL_NEF_TRACEN("  supporting edge "<<PH(ei));
          incident_mark(e->twin()) = ei->twin()->incident_sface()->mark();
          e->mark() = ei->mark();
          incident_mark(e) = m_buffer = ei->incident_sface()->mark(); 
        }
        SHalfloop_const_handle li;
        if ( CGAL::assign(li,support(e)) ) { 
          if ( li->circle() != e->circle() ) { li = li->twin(); }
          CGAL_assertion( li->circle() == e->circle() ); 
          CGAL_NEF_TRACEN("  supporting loop "<<PH(li));
          incident_mark(e->twin()) = li->twin()->incident_sface()->mark();
          e->mark() = li->mark();
          incident_mark(e) = m_buffer = li->incident_sface()->mark();
        }
      } else { CGAL_NEF_TRACEN("  support from face below ");
        incident_mark(e->twin()) = e->mark() = 
        incident_mark(e) = m_buffer;
      }
      CGAL_NEF_TRACEN("  new face mark "<<m_buffer);

    } // CGAL_For_all(e,hend)

    CGAL_NEF_TRACEN(" mark of "<<PH(v));
  }

}




CGAL_END_NAMESPACE
#undef CGAL_USING
#endif //CGAL_SM_TRIANGULATOR_H