sm_overlayer.h

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

    number_of_sphere_sweeps++;
#endif
  }

  if(compute_halfsphere[cs][1]) {
    NH_geometry nhg(cs);

    /*
    SHalfedge_const_handle se;
    SHalfloop_const_handle sl;
    for(Seg_iterator it=L_neg.begin(); it!=L_neg.end();++it) {
      CGAL_NEF_TRACEN("neg " << *it);
      if(nhg.compare_xy(it->target(),it->source())<0) {
	Object_handle o = From[it]._o;
	if(CGAL::assign(se, o)) {
	  if(it->sphere_circle() == se->circle())
	    From[it] = Seg_info(se->twin(), From[it]._from);
	} else if(CGAL::assign(sl, o)) {
	  if(it->sphere_circle() == sl->circle())
	  From[it] = Seg_info(sl->twin(), From[it]._from);
	}
      }
    }
    */

    Negative_halfsphere_sweep SM(
        Input_range(L_neg.begin(),L_neg.end()),O,
        nhg);
    SM.sweep();
#ifdef CGAL_NEF3_TIMER_SPHERE_SWEEPS
    number_of_sphere_sweeps++;
#endif
  }

#ifdef CGAL_NEF3_TIMER_SPHERE_SWEEPS
  timer_sphere_sweeps.stop();
#endif

  if(compute_halfsphere[cs][0]) {
    ++v; 
    ++e;
  }
  else {
    v = this->svertices_begin(); 
    e = this->shalfedges_begin();
  }

  if(compute_halfsphere[cs][0])
    create_face_objects(this->shalfedges_begin(), e, this->svertices_begin(), v, O,
                        PH_geometry(cs));
  if(compute_halfsphere[cs][1])
    create_face_objects(e, this->shalfedges_end(), v, this->svertices_end(), O,
			NH_geometry(cs));

  CGAL_forall_sedges(e,*this) {
    e->circle() = Sphere_circle(e->source()->point(), e->twin()->source()->point()); 
    e->twin()->circle() = e->circle().opposite();
    CGAL_NEF_TRACEN(PH(e) << " with circle " << e->circle());
  }

  std::vector<Mark> mohs(4);
  SM_point_locator L0(M0);
  SM_point_locator L1(M1);
  
  L0.marks_of_halfspheres(mohs, 0, cs);
  L1.marks_of_halfspheres(mohs, 2, cs);

  CGAL_NEF_TRACEN("mohs[0]=" << mohs[0]);
  CGAL_NEF_TRACEN("mohs[1]=" << mohs[1]);
  CGAL_NEF_TRACEN("mohs[2]=" << mohs[2]);
  CGAL_NEF_TRACEN("mohs[3]=" << mohs[3]);

  CGAL_NEF_TRACEN("compute_halfsphrere\n  cs = " << cs << 
	 "\n  [cs][0] = " << compute_halfsphere[cs][0] <<
	 "\n  [cs][1] = " << compute_halfsphere[cs][1]);
  /*
  SVertex_iterator svii;
  CGAL_forall_svertices(svii, *this) {
    SVertex_const_handle vs;
    SHalfedge_const_handle es;
    if(CGAL::assign(vs, supp_object(svii,0)))
      std::cerr << svii->point() << " supported by svertex " << vs->point() << std::endl;
    else if(CGAL::assign(es, supp_object(svii,0)))
      std::cerr << svii->point() << " supported by sedge" << std::endl;
    else
      std::cerr << svii->point() << " is neither supported by svertex or sedge" << std::endl;
    if(CGAL::assign(vs, supp_object(svii,1)))
      std::cerr << svii->point() << " supported by svertex " << vs->point() << std::endl;
    else if(CGAL::assign(es, supp_object(svii,1)))
      std::cerr << svii->point() << " supported by sedge" << std::endl;
    else
      std::cerr << svii->point() << " is neither supported by svertex or sedge" << std::endl;
  }
  */
  if(compute_halfsphere[cs][0])
    complete_face_support(this->svertices_begin(), v, O, mohs, 0, 
			  compute_halfsphere[cs][1]);
  if(compute_halfsphere[cs][1])
    complete_face_support(v, this->svertices_end(), O, mohs, 1,
			  compute_halfsphere[cs][0]);

  complete_sface_marks();

  // DEBUG CODE: to do: have all svertices a halfedge below associated?
  CGAL_NEF_TRACEN("Vertex info after swep");
  CGAL_assertion_code(SVertex_iterator svi);
  CGAL_assertion_code(
    for(svi=this->svertices_begin(); svi!=this->svertices_end(); svi++) {
      CGAL_NEF_TRACEN("vertex "<<svi->point()<<" info "<<info(svi)<<
	     " marks "<<mark(svi,0)<<" "<<mark(svi,1));
    }
  )

  transfer_data(A);

  if(compute_halfsphere[cs][0] && compute_halfsphere[cs][1])
    merge_halfsphere_maps(this->svertices_begin(),v,O);
  else
    set_outer_face_mark(compute_halfsphere[cs][1], mohs);
  
  CGAL_assertion_code(this->check_integrity_and_topological_planarity());
  
  CGAL_NEF_TRACEN("subdivided");
  CGAL_assertion_code(CGAL_forall_svertices(v,*this) CGAL_NEF_TRACEN(PH(v)));
}

template <typename Map>
template <typename Association>
void SM_overlayer<Map>::
transfer_data(Association& A) {

  SVertex_iterator sv;
  SHalfedge_handle se;
  SVertex_const_handle sv0,sv1;
  SHalfedge_const_handle se0, se1;
  SHalfloop_const_handle sl0, sl1;

  CGAL_forall_svertices(sv, *this) {
    //    std::cerr << "svertex " << sv->point() << std::endl;
    Object_handle o0 = supp_object(sv,0), o1 = supp_object(sv,1);
    if(o0 == NULL) {
      if(CGAL::assign(sv1, o1))
	A.handle_support(sv, sv1);
      else
	continue;
    } else if(CGAL::assign(se0, o0)) {
      if(o1 == NULL) 
	continue;
      else if(assign(se1, o1))
	A.handle_support(sv, se0, se1);
      else if(CGAL::assign(sv1, o1))
	A.handle_support(sv, se0, sv1);
      else if(CGAL::assign(sl1, o1))
	A.handle_support(sv, se0, sl1);
      else
	CGAL_assertion_msg(false, "wrong handle");
    } else if(CGAL::assign(sv0, o0)) {
      if(o1 == NULL)
	A.handle_support(sv, sv0);
      else if(CGAL::assign(se1, o1))
	A.handle_support(sv, sv0, se1);
      else if(CGAL::assign(sv1, o1))
	A.handle_support(sv, sv0, sv1);
      else if(CGAL::assign(sl1, o1))
	A.handle_support(sv, sv0, sl1);
      else
	CGAL_assertion_msg(false, "wrong handle");
    } else if(CGAL::assign(sl0, o0)) {
      if(o1 == NULL)
	continue;
      else if(CGAL::assign(sv1, o1))
	A.handle_support(sv, sl0, sv1);
      else if(CGAL::assign(se1, o1))
	A.handle_support(sv, sl0, se1);
      else if(CGAL::assign(sl1, o1))
	A.handle_support(sv, sl0, sl1);
    } else 
      CGAL_assertion_msg(false, "wrong handle");
  }

  CGAL_forall_sedges(se, *this) {
    CGAL_assertion(is_forward(se));
    Object_handle o0 = supp_object(se,0), o1 = supp_object(se,1);   
    if(o0 == NULL) {
      if(assign(se1, o1))
	A.handle_support(se, se1);
      else if(assign(sl1, o1))
	A.handle_support(se, sl1);
      else
	continue; // CGAL_assertion_msg(false, "wrong handle");
    } else if(assign(se0, o0)) {
      if(o1 == NULL)
	A.handle_support(se, se0);
      else if(assign(se1, o1))
	A.handle_support(se, se0, se1);
      else if(assign(sl1, o1))
	A.handle_support(se, se0, sl1);
      else
	CGAL_assertion_msg(false, "wrong handle");    
    } else if(assign(sl0, o0)) {
      if(o1 == NULL)
	A.handle_support(se, sl0);
      else if(assign(se1, o1))
	A.handle_support(se, sl0, se1);
      else if(assign(sl1, o1))
	A.handle_support(se, sl0, sl1);
      else
	CGAL_assertion_msg(false, "wrong handle");
    } else
      CGAL_assertion_msg(false, "wrong handle");	      
  }
}

template <typename Map>
void SM_overlayer<Map>::
set_outer_face_mark(int offset, const std::vector<Mark>& mohs) {

  SFace_handle sf = this->new_sface();
  assoc_info(sf);
  mark(sf, 0) = mohs[offset];
  mark(sf, 1) = mohs[offset+2];

  SHalfedge_iterator e;
  CGAL_forall_shalfedges(e, *this) {
    if ( e->incident_sface() != SFace_handle() ) continue;
    link_as_face_cycle(e,sf); 
  }
  
  SVertex_handle v;
  CGAL_forall_svertices(v, *this) {
    if(!is_isolated(v) || v->incident_sface() != SFace_handle()) continue;
    link_as_isolated_vertex(v,sf);
  }
}

template <typename Map>
template <typename Iterator, typename T>
void SM_overlayer<Map>::
partition_to_halfsphere(Iterator start, Iterator beyond, Seg_list& L, 
			CGAL::Unique_hash_map<Iterator,T>& M, 
			Sphere_circle xycircle, Sphere_circle yzcircle, 
			bool include_equator) const
{ CGAL_NEF_TRACEN("partition_to_halfsphere ");
//  CGAL_assertion(pos!=0);
  Sphere_segment s1,s2;
  //  Sphere_circle xycircle(0,0,pos);
  if(include_equator) {
    while ( start != beyond) {
      int i = start->intersection(xycircle,s1,s2);
      CGAL_NEF_TRACEN("segment " << start->source() << " " << start->target());
      if (i>1) { 
	L.push_back(s2); M[--L.end()] = M[start];    
	CGAL_NEF_TRACEN(">1 " << s2.source() << " " << s2.target()); 
      }
      if (i>0) { 
	L.push_back(s1); M[--L.end()] = M[start];    
	CGAL_NEF_TRACEN(">0 " << s1.source() << " " << s1.target()); 
      }
      ++start;
    }
  }
  else {
    while(start != beyond) {
      L.push_back(*start);
      M[--L.end()] = M[start];
      ++start;
    }
  }

  // now all segments are split into hemispheres
  // we still have to:
  // - split segments containing our special poles y^-, y^+
  // - split halfcircles
  // - add four equator segments 

  //  Sphere_circle yzcircle(1,0,0);
  typename Seg_list::iterator it, itl;
  
  CGAL_forall_iterators(it,L) { CGAL_NEF_TRACEN("  "<<*it);
    if ( equal_as_sets(it->sphere_circle(),xycircle) ) {
      CGAL_NEF_TRACEN("  splitting xy seg "<<*it);
      bool added=false;
      int n1 =  it->intersection(yzcircle,s1,s2);
      if (n1 > 1 && !s2.is_degenerate()) { 
	M[ L.insert(it,s2) ] = M[it]; 
	added=true; 
	CGAL_NEF_TRACEN(">1 " << s2.source() << " " << s2.target()); 
      }
      if (n1 > 0 && !s1.is_degenerate()) { 
	M[ L.insert(it,s1) ] = M[it]; 
	added = true; 
	CGAL_NEF_TRACEN(">1 " << s1.source() << " " << s1.target()); 
      }
      int n2 =  it->intersection(yzcircle.opposite(),s1,s2);
      if (n2 > 1 && !s2.is_degenerate()) { 
	M[ L.insert(it,s2) ] = M[it]; 
	added=true; 
	CGAL_NEF_TRACEN(">1 " << s2.source() << " " << s2.target()); 
      }
      if (n2 > 0 && !s1.is_degenerate()) { 
	M[ L.insert(it,s1) ] = M[it]; 
	added=true; 
	CGAL_NEF_TRACEN(">1 " << s1.source() << " " << s1.target()); 
      }
      if(added) {
	itl = it; --it; M[itl] = T(); L.erase(itl);
      }
      // at least one item was appended
    }
  }

  CGAL_forall_iterators(it,L) {
    if ( it->is_halfcircle() ) {
      CGAL_NEF_TRACEN("  splitting halfcircle "<<*it);
      Sphere_segment s1,s2;
      it->split_halfcircle(s1,s2);
      CGAL_NEF_TRACEN("    into " << s1);
      CGAL_NEF_TRACEN("    into " << s2);
      *it = s2; 
      M[ L.insert(it,s1) ] = M[it];
    }
  }

  if(include_equator) {
    // append 4 xy-equator segments:
    Sphere_point S(0,-1,0),N(0,1,0);
    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);
  }
}

template <typename Map>
template <typename Below_accessor, typename Halfsphere_geometry>
void SM_overlayer<Map>::
create_face_objects(SHalfedge_iterator e_start, SHalfedge_iterator e_end,
  SVertex_iterator v_start, SVertex_iterator v_end,
  const Below_accessor& D, 
  const Halfsphere_geometry& SG)
{
  CGAL_NEF_TRACEN("create_face_objects()");
  if(e_start != e_end) {
    CGAL::Unique_hash_map<SHalfedge_handle,int> SFaceCycle(-1);
    std::vector<SHalfedge_handle>  MinimalSHalfedge;
    SHalfedge_around_sface_circulator hfc(last_out_edge(v_start)),hend(hfc);
    CGAL_NEF_TRACEN("equator cycle "<<PH(hfc));
    CGAL_For_all(hfc,hend) SFaceCycle[hfc]=0; // outer face cycle = 0
    MinimalSHalfedge.push_back(first_out_edge(v_start)->twin());
    int i=1; 
    for (SHalfedge_iterator e = e_start; e != e_end; ++e) {
      if ( SFaceCycle[e] >= 0 ) continue; // already assigned
      SHalfedge_around_sface_circulator hfc(e),hend(hfc);
      SHalfedge_handle e_min = e;
      CGAL_NEF_TRACEN(""); 
      CGAL_NEF_TRACEN("  face cycle numbering "<<i);
      CGAL_For_all(hfc,hend) {
	SFaceCycle[hfc]=i; // assign face cycle number
	if (hfc->twin()->source() == e_min->twin()->source()) {
	  Sphere_point p1 = e->source()->point(), 
	    p2 = e->twin()->source()->point(), 
	    p3 = e->snext()->twin()->source()->point();	  
	  if ( SG.orientation(p1,p2,p3) <= 0 )
	    e_min = hfc;
	} else if ( SG.compare_xy(hfc->twin()->source()->point(), e_min->twin()->source()->point()) < 0 )
	  e_min = hfc;    
	CGAL_NEF_TRACEN(PH(hfc));
      } CGAL_NEF_TRACEN("");
      MinimalSHalfedge.push_back(e_min);
      ++i;
    }
    
    for (int j=1; j<i; ++j) {
      SHalfedge_handle e = MinimalSHalfedge[j];
      CGAL_NEF_TRACEN("  face cycle "<<j<<" minimal halfedge "<<PH(e));
      Sphere_point p1 = e->source()->point(), 
	p2 = e->twin()->source()->point(), 
	p3 = e->snext()->twin()->source()->point();
      if ( SG.or