snc_constructor.h

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

    if(on_sface[0]) {
      std::cerr << "found sf0 " << std::endl;
      sv[2]->mark() = BOP(sf0->mark(), e1->mark(), inv);
      SFace_cycle_iterator sfci(sf0->sface_cycles_begin());
      CGAL_assertion_msg(sfci.is_shalfedge(), "not implemented, yet");
      SHalfedge_handle se_tmp(sfci);
      see[0] = se_tmp;
      if(see[0]->source() != sv[0]) see[0] = see[0]->snext();
      seb[1] = see[0];
    } else {
      CGAL::assign(se0, o2);
      CGAL_assertion(CGAL::assign(se0, o2));
      std::cerr << "found se0 " << se0->source()->point()
		<< "->" << se0->twin()->source()->point() << std::endl;
      std::cerr << "insert sv " << sv[2]->point() << std::endl;
      if(se0->source() != sv[0]) se0 = se0->twin();
      CGAL_assertion(se0->source() == sv[0]);
      sv[2]->mark() = BOP(se0->mark(), e1->mark(), inv);
      SHalfedge_handle se_new =
	D.new_shalfedge_pair();
      se_new->source() = sv[2];
      se_new->twin()->source() = se0->twin()->source();
      se0->twin()->source() = sv[2];

      se_new->snext() = se0->snext();
      se0->snext()->sprev() = se_new;
      se_new->twin()->sprev() = se0->twin()->sprev();
      se0->twin()->sprev()->snext() = se_new->twin();

      se0->snext() = se_new;
      se0->twin()->sprev() = se_new->twin();
      se_new->sprev() = se0;
      se_new->twin()->snext() = se0->twin();

      se_new->mark() = se_new->twin()->mark() = se0->mark();
      se_new->incident_sface() = se0->incident_sface();
      se_new->twin()->incident_sface() = se0->twin()->incident_sface();
      se_new->circle() = se0->circle();
      se_new->twin()->circle() = se0->twin()->circle();

      SM_io_parser<SM_decorator>::dump(D, std::cerr);

      see[0] = seb[1]= se0;
      ++seb[1];
    } 

    if(on_sface[1]) {
      std::cerr << "found sf1 " << std::endl;
      sv[3]->mark() = BOP(sf1->mark(), e1->mark(), inv);
      SFace_cycle_iterator sfci(sf1->sface_cycles_begin());
      CGAL_assertion_msg(sfci.is_shalfedge(), "not implemented, yet");
      SHalfedge_handle se_tmp(sfci);
      see[1] = se_tmp;
      if(see[1]->source() != sv[0]) see[1] = see[1]->snext();
      seb[0] = see[1];
    } else {
      CGAL::assign(se1, o3);
      CGAL_assertion(CGAL::assign(se1, o3));
      std::cerr << "found se1 " << se1->source()->point()
		<< "->" << se1->twin()->source()->point() << std::endl;
      std::cerr << "insert sv " << sv[3]->point() << std::endl;
      if(se1->source() != sv[0]) se1 = se1->twin();
      CGAL_assertion(se1->source() == sv[0]);
      sv[3]->mark() = BOP(se1->mark(), e1->mark(), inv);
      SHalfedge_handle se_new = 
	D.new_shalfedge_pair();
      se_new->source() = sv[3];
      se_new->twin()->source() = se1->twin()->source();
      se1->twin()->source() = sv[3];

      se_new->snext() = se1->snext();
      se1->snext()->sprev() = se_new;
      se_new->twin()->sprev() = se1->twin()->sprev();
      se1->twin()->sprev()->snext() = se_new->twin();

      se1->snext() = se_new;
      se1->twin()->sprev() = se_new->twin();
      se_new->sprev() = se1;
      se_new->twin()->snext() = se1->twin();

      se_new->mark() = se_new->twin()->mark() = se1->mark();
      se_new->incident_sface() = se1->incident_sface();
      se_new->twin()->incident_sface() = se1->twin()->incident_sface();
      se_new->circle() = se1->circle();
      se_new->twin()->circle() = se1->twin()->circle();

      SM_io_parser<SM_decorator>::dump(D, std::cerr);

      see[1] = seb[0] = se1;
      ++seb[1];
    } 

    if(E.is_isolated(e1)) {
      CGAL_assertion_msg(false, "not implemented, yet");
      return D.sphere_map();
    }

    SHalfedge_around_svertex_const_circulator 
      svc(e1->out_sedge()), send(svc);
    CGAL_assertion(svc != svc->twin()->snext());
    // TODO: one sedge
    
    int i=0;
    bool done = false;
    Oriented_side os0 =
      svc->circle().oriented_side(sv[0]->point());
    Oriented_side os1 = os0;
    while(os0 == os1 &&
	  os1 != ON_ORIENTED_BOUNDARY &&
	  ++svc != send) 
      os1 = svc->circle().oriented_side(sv[0]->point());

    if(os1 == ON_ORIENTED_BOUNDARY) {

      Sphere_segment seg(sv[2]->point(), sv[3]->point(),
			 svc->circle());
      int sv_index =
	seg.has_on(sv[0]->point()) ? 0 : 1;
      
      std::cerr << "sv " << sv[0]->point() << ", " << sv[1]->point() << std::endl;

      if(on_sface[0]) {
	SFace_cycle_iterator sfci = sf0->sface_cycles_begin();
	CGAL_assertion(sfci.is_shalfedge());
	SHalfedge_handle se_tgt(sfci);
	while(se_tgt->source() != sv[sv_index])
	  se_tgt = se_tgt->snext();
	CGAL_assertion(se_tgt->source() == sv[sv_index]);
	std::cerr << "new_shalfedge_pair" << std::endl;
	SHalfedge_handle se_new = 
	  D.new_shalfedge_pair(sv[2], se_tgt);
	std::cerr << "done" << std::endl;
	se_new->mark() = se_new->twin()->mark() = 
	  BOP(sf0->mark(), svc->mark(), inv);
	splitted[0] = true;
      }

      if(on_sface[1]) {
	SFace_cycle_iterator sfci = sf1->sface_cycles_begin();
	CGAL_assertion(sfci.is_shalfedge());
	SHalfedge_handle se_tgt(sfci);
	while(se_tgt->source() != sv[sv_index]) {
	  se_tgt = se_tgt->snext();
	  std::cerr << se_tgt->source()->point() << " " << sv[sv_index]->point() << std::endl;
	  std::cerr << &(se_tgt->source()) << " " << &sv[sv_index] << std::endl;
	}
	CGAL_assertion(se_tgt->source() == sv[sv_index]);
	std::cerr << "new_shalfedge_pair" << std::endl;
	SHalfedge_handle se_new =
	  D.new_shalfedge_pair(sv[3], se_tgt);
	std::cerr << "done" << std::endl;
	se_new->mark() = se_new->twin()->mark() = 
	  BOP(sf1->mark(), svc->mark(), inv);
	splitted[1] = true;
      }
      
      ++svc;
      os1 = svc->circle().oriented_side(sv[0]->point());
      if(os1 == ON_ORIENTED_BOUNDARY) {
	CGAL_assertion_msg(false, "not implemented, yet");
	++svc;
	scb[0] = scb[1] = sce[0] = sce[1] = svc;
	if(svc != send) {
	  os1 = svc->circle().oriented_side(sv[0]->point());
	  i = os1 == ON_POSITIVE_SIDE ? 0 : 1;
	  --sce[1-i];
	}
	done = true;
      } else {
	i = os1 == ON_POSITIVE_SIDE ? 0 : 1;
	scb[i] = svc;
	--svc;
	sce[1-i] = svc;
	--svc;
	if(svc->circle().oriented_side(sv[0]->point()) == os1) {
	  sce[i] = scb[1-i] = sce[1-i];
	  done = true;
	}
      }
    } else if(svc == send) {
      i = os1 == ON_POSITIVE_SIDE ? 0 : 1;
      CGAL_assertion_msg(false, "not implemented, yet");
      done = true;
    } else {
      CGAL_assertion(os0 != os1);
      i = os1 == ON_POSITIVE_SIDE ? 0 : 1;
      sce[1-i] = scb[i] = svc;
    }
    CGAL_assertion(scb[i] == svc);

    if(!done) {
      os0 = svc->circle().oriented_side(sv[0]->point());
      CGAL_assertion(os0 != ON_ORIENTED_BOUNDARY);
      do {
	++svc;
	os1 = svc->circle().oriented_side(sv[0]->point());
      } while(os1 == os0);

      sce[i] = scb[1-i] = svc;
      if(os1 == ON_ORIENTED_BOUNDARY) {
	CGAL_assertion_msg(false, "degenerate case not handled");
	++scb[1-i];
      }
    }

    for(; scb[0] != sce[0]; ++scb[0]) {
      for(; seb[0] != see[0]; ++see[0]) {
	Sphere_segment seg0(sv[0]->point(), sv[1]->point(), seb[0]->circle());
	Sphere_segment seg1(sv[2]->point(), sv[3]->point(), scb[0]->circle());
	Sphere_point sp = seg0.intersection(seg1);
	std::cerr << "intersections oben " << sp << std::endl;
      }
    }
    
    for(; scb[1] != sce[1]; ++scb[1]) {
      for(; seb[1] != see[1]; ++see[1]) {
	Sphere_segment seg0(sv[0]->point(), sv[1]->point(), seb[1]->circle());
	Sphere_segment seg1(sv[2]->point(), sv[3]->point(), scb[1]->circle());
	Sphere_point sp = seg0.intersection(seg1);
	std::cerr << "intersections unten " << sp << std::endl;
      }
    }

    return D.sphere_map();
  }
#endif

  template<typename Selection, typename Association>
  Sphere_map* create_edge_facet_overlay( Halfedge_const_handle e, 
					 Halffacet_const_handle f,
					 const Point_3& p,
					 const Selection& BOP, bool inv,
					 Association& ) {

#ifdef CGAL_NEF3_TIMER_SPHERE_SWEEPS
    ++number_of_edge_facet_overlays;
#endif

    CGAL_NEF_TRACEN("edge facet overlay " << p);

    Unique_hash_map<SHalfedge_handle, Mark> mark_of_right_sface;
    
    SM_decorator D(&*this->sncp()->new_vertex(p, BOP(e->mark(), f->mark())));
    SM_const_decorator E(&*e->source());

    Sphere_point ps = e->point();
    ps = normalized(ps);
    SVertex_handle v1 = D.new_svertex(ps);
    SVertex_handle v2 = D.new_svertex(ps.antipode());

    CGAL_NEF_TRACEN("new svertex 1 " << ps);
    CGAL_NEF_TRACEN("new svertex 2 " << ps.antipode());
    Halffacet_const_handle faces_p(f);
    Vector_3 vec(ps-CGAL::ORIGIN);
    if(faces_p->plane().oriented_side(p+vec) == ON_NEGATIVE_SIDE)
      faces_p = faces_p->twin();
    v1->mark() = BOP(e->mark(), faces_p->incident_volume()->mark(), inv);
    v2->mark() = BOP(e->mark(), faces_p->twin()->incident_volume()->mark(), inv);
    CGAL_NEF_TRACEN("svertex 1 " << ps << " has mark " << v1->mark());
    CGAL_NEF_TRACEN("svertex 2 " << ps.antipode() << " has mark " << v2->mark());

    if(E.is_isolated(e)) {
      CGAL_NEF_TRACEN("edge is isolated");
      Mark mf1 = BOP(e->incident_sface()->mark(), faces_p->incident_volume()->mark(), inv);
      Mark mf2 = BOP(e->incident_sface()->mark(), faces_p->twin()->incident_volume()->mark(), inv);
      Mark ml = BOP(e->incident_sface()->mark(), faces_p->mark(), inv);
      
      SFace_handle f1 = D.new_sface();
      D.link_as_isolated_vertex(v1, f1);
      f1->mark() = mf1;
      
      if(mf1 == mf2 && mf1 == ml) {
	D.link_as_isolated_vertex(v2, f1);
      }
      else {
	SHalfloop_handle l = D.new_shalfloop_pair();
	SFace_handle f2 = D.new_sface();    
	D.link_as_isolated_vertex(v2, f2);
	D.link_as_loop(l,f1);
	D.link_as_loop(l->twin(),f2);
	l->circle() = Sphere_circle(faces_p->plane()); 
	l->twin()->circle() = l->circle().opposite();
	f2->mark() = mf2;
	l->mark() = l->twin()->mark() = ml;
      }
    }
    else {
      CGAL_NEF_TRACEN("edge is not isolated");
      SVertex_handle sv;
      SHalfedge_handle se1;
      SHalfedge_handle se2;
      SFace_handle sf;
      Sphere_circle c(f->plane());
      
      SHalfedge_handle next_edge;
      SHalfedge_around_svertex_const_circulator ec(E.out_edges(e)), ee(ec);
      CGAL_For_all(ec,ee) {
	Sphere_segment seg(ec->source()->point(), 
			   ec->source()->point().antipode(), 
			   ec->circle());
	Sphere_point sp(intersection(c, seg.sphere_circle()));
	CGAL_NEF_TRACEN(seg <<" has_on " << sp);
	if(!seg.has_on(sp))
	  sp = sp.antipode();
	sv = D.new_svertex(sp);
	CGAL_NEF_TRACEN("new svertex 3 " << normalized(sp));
	sv->mark() = BOP(ec->mark(), f->mark(), inv);
	se1 = D.new_shalfedge_pair(v1, sv);
	if(next_edge == SHalfedge_handle())
	  se2 = D.new_shalfedge_pair(sv, v2); 
	else
	  se2 = D.new_shalfedge_pair(sv, next_edge, -1);
	next_edge = se2->twin();
	se1->mark() = se1->twin()->mark() = BOP(ec->mark(), faces_p->incident_volume()->mark(), inv);
	se2->mark() = se2->twin()->mark() = BOP(ec->mark(), faces_p->twin()->incident_volume()->mark(), inv);
	mark_of_right_sface[se1] = ec->incident_sface()->mark();
	se1->circle() = se2->circle() = ec->circle();
	se1->twin()->circle() = se2->twin()->circle() = se1->circle().opposite();
      }
      
      SHalfedge_around_svertex_circulator ec2(D.out_edges(v1)), ee2(ec2);
      CGAL_For_all(ec2,ee2) {
	SHalfedge_around_svertex_circulator en(ec2);
	++en;
	se1 = D.new_shalfedge_pair(ec2->twin(), en->twin(), -1, 1);
	CGAL_NEF_TRACEN("new edge pair " << ec2->twin()->source()->vector() << 
			" -> " << en->twin()->source()->vector());
	se1->circle() = Sphere_circle(faces_p->plane());
	se1->twin()->circle() = se1->circle().opposite();
        se1->mark() = se1->twin()->mark() = BOP(mark_of_right_sface[ec2], faces_p->mark(), inv);

	sf = D.new_sface();
	sf->mark() = BOP(mark_of_right_sface[ec2], faces_p->incident_volume()->mark(), inv);
	D.link_as_face_cycle(se1,sf);
	sf = D.new_sface();
	sf->mark() = BOP(mark_of_right_sface[ec2], faces_p->twin()->incident_volume()->mark(), inv);
	D.link_as_face_cycle(se1->twin(),sf);
      }   
    }
    
    return D.sphere_map();
  }

 public:
  Point_3 get_transformed_coords_of_vertex(const Point_3& p,
					   const std::list<Point_3>& segs1,
  					   const std::list<Point_3>& segs2) {
    CGAL_assertion(!segs1.empty() && !segs2.empty());
    
    int side_of_point=1;
    RT max = p.hx();
    if(CGAL_NTS abs(p.hy()) > CGAL_NTS abs(max)) {
      max = p.hy();
      side_of_point=2;
    }
    if(CGAL_NTS abs(p.hz()) > CGAL_NTS abs(max)) {
      max = p.hz();
      side_of_point=3;
    }
    if(max < RT(0)) 
      side_of_point = -side_of_point;
    
    typename std::list<Point_3>::const_iterator s1,s2,t1,t2;
    s1 = Infi_box::segment_on_side(side_of_point, segs1);
    t1 = s1;
    ++t1;
    if(t1 == segs1.end()) t1=segs1.begin();
    s2 = Infi_box::segment_on_side(side_of_point, segs2);
    t2 = s2;
    ++t2;
    if(t2 == segs2.end()) t2=segs2.begin();
    
    SNC_intersection is;
    Point_3 ip;
    bool flag=is.does_intersect_internally(Segment_3(*s1,*t1),Segment_3(*s2,*t2),ip);
    if(!flag) {
      if(*s1 == *s2) return normalized(*s1);
      else if(*s1 == *t2) return normalized(*s1);
      else if(*t1 == *s2) return normalized(*t1);
      else if(*t1 == *t2) return normalized(*t1);
    }
    return normalized(ip);
  }