sm_point_locator.h

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

    Sphere_circle c(d.circle());
    Sphere_segment s;
    bool s_init(false);
    Object_handle h = locate(p);
    SVertex_handle v; 
    SHalfedge_handle e; 
    SHalfloop_handle l; 
    SFace_handle f;
    if ( CGAL::assign(v,h) && M(v) ||
         CGAL::assign(e,h) && M(e) ||
	 CGAL::assign(l,h) && M(l) ||
         CGAL::assign(f,h) && M(f) ) return h;
    h = Object_handle(); 
    CGAL_NEF_TRACEN("not contained");
#if 0
    HASEN: s am anfang circle, ab wann segment ?
	   wo loop ?

    CGAL_forall_svertices (v,*this) {
      Point pv = v->point();
      if ( !(s_init && s.has_on(pv) ||
	    !s_init && c.has_on(pv)) ) continue;
      CGAL_NEF_TRACEN("candidate "<<pv);
      if ( M(v) ) {
        h = Object_handle(v);     // store vertex
        s = Sphere_segment(p,pv,c); // shorten
        continue;
      }
      // now we know that v is not marked but on s
      bool collinear;
      SHalfedge_handle e = out_wedge(v,d,collinear);
      if ( collinear ) { 
        if ( M(e) ) {
          h = Object_handle(e);
          s = Sphere_segment(p,pv,c);
        }
        continue;
      }
      if ( M(e->incident_sface()) ) {
        h = Object_handle(e->incident_sface());
        s = Sphere_segment(p,pv,c);
      }
    } // all vertices

    CGAL::Unique_hash_map<SHalfedge_handle,bool> visited(false);
    SHalfedge_iterator e_res;
    CGAL_forall_sedges(e,*this) {
      Sphere_segment se = segment(e);
      Sphere_point p_res;
      if ( do_intersect_internally(se,s,p_res) ) {
        // internal intersection
        CGAL_NEF_TRACEN("candidate "<<se); 
        e_res = e;
        Sphere_segment s_cand = Sphere_segment(p,p_res,c);  
        if ( s_cand.is_short() && e->circle().has_on_negative_side(p) ||
	     s_cand.is_long() && e->circle().has_on_positive_side(p) ||
	     s_cand.is_halfcircle() && 
	       strictly_ordered_ccw_at(p.antipode(),
				       direction(e),d,direction(e->twin())) )
	  e_res = e->twin();
        if ( M(e_res) ) {
          h = Object_handle(e_res); s = s_cand;
        } else if ( M(face(twin(e_res))) ) {
          h = Object_handle(face(twin(e_res))); s = s_cand;
        }
      }
    }
#endif
    CGAL_assertion_msg(0,"not yet correct");
    return h;
  }

  Object_handle ray_shoot(const Sphere_point& p, 
			  const Sphere_circle& c,
			  Sphere_point& ip,
			  bool start_inclusive = false) { 
    Sphere_segment seg(p, p.antipode(), c);
    return ray_shoot(seg, ip, start_inclusive);
  }

  Object_handle ray_shoot(const Sphere_segment& d, 
			  Sphere_point& ip,
			  bool start_inclusive = false,
			  bool beyond_end = true,
			  bool end_inclusive = false) { 

    // TODO: end_inclusive=true does not work properly for sedges and sloops

    CGAL_NEF_TRACEN("ray shoot");
    Sphere_circle c(d.sphere_circle());
    Sphere_point p(d.source());
    Sphere_segment s;
    bool s_init(false);
    
    if(!beyond_end) {
      s = d;
      s_init = true;
    }

    Object_handle h = Object_handle();

    if(s_init) {
      CGAL_NEF_TRACEN(" at begin " << s_init << ":" << s);
    } else {
      CGAL_NEF_TRACEN(" at begin " << s_init << ":" << c);
    }

    SVertex_iterator vi;
    CGAL_forall_svertices (vi,*this) {
      Sphere_point pv = vi->point();
      if ((s_init && !s.has_on(pv)) ||
	  (!s_init && !c.has_on(pv))) continue;
      CGAL_NEF_TRACEN("candidate "<<pv);
      if ((start_inclusive || p != pv) && 
	  (end_inclusive || !s_init || s.target() != pv)) {
        h = Object_handle(vi);     // store vertex
        s = Sphere_segment(p,pv,c); // shorten
	ip = pv;
	s_init = true;
      }
    }

    // TODO: edges on the ray.

    SHalfedge_iterator ei;
    CGAL_forall_sedges(ei,*this) {
      Sphere_segment se = segment(ei);
      CGAL_NEF_TRACEN("ray_shoot " << s_init);
      if(s_init)
	CGAL_NEF_TRACEN("  " << s.source() << "->" << s.target() << " | " << s.sphere_circle() << " is long " << s.is_long());
      CGAL_NEF_TRACEN("  " << se.source() << "->" << se.target() << " | " << se.sphere_circle() << " is long " << se.is_long());

      // TODO: start-end point of s on se or c

      Sphere_point p_res;
      if(se.source() == se.target()) {

	if(s_init) {
	  if(s.is_long()) {
	    Sphere_segment first_half(p,p.antipode(),c);
	    Sphere_segment second_part(p.antipode(), s.target(), c);
	    if(!do_intersect_internally(ei->circle(), first_half, p_res)) {
	      do_intersect_internally(ei->circle(), second_part, p_res);
	    }
	  } else {
	    if(!do_intersect_internally(ei->circle(), s, p_res)) continue;
	  }
	} else {
	  Sphere_segment first_half(p,p.antipode(),c);
	  Sphere_segment second_part(p.antipode(),p,c);
	  if(!do_intersect_internally(ei->circle(), first_half, p_res)) {
	    do_intersect_internally(ei->circle(), second_part, p_res);
	  }
	}

      } else {
	if(s_init) {	  
	  if(s.is_long() && se.is_long()) {
	    Sphere_segment first_half(p,p.antipode(),c);
	    Sphere_segment second_part(p.antipode(), s.target(), c);
	    if(!do_intersect_internally(se, first_half, p_res) &&
	       !do_intersect_internally(se, second_part, p_res)) {
	      if(se.has_on(p.antipode()))
		p_res = p.antipode();
	      else
		continue;
	    }	    
	  } else {
	    if(!do_intersect_internally(se, s, p_res)) continue;
	  }  
	} else {
	  if(se.is_long()) {
	    Sphere_segment first_half(p,p.antipode(),c);
	    Sphere_segment second_half(p.antipode(),p,c); 
	    if(!do_intersect_internally(se, first_half, p_res)) {
	      if(!do_intersect_internally(se, second_half, p_res)) {
		if(start_inclusive)
		  p_res = p;
		else
		  p_res = p.antipode();
	      }
	    }
	  } else {
	    if(!do_intersect_internally(c, se, p_res)) continue;
	  }
	}
      }
      
      CGAL_NEF_TRACEN("candidate "<<se); 
      if (start_inclusive || p != p_res) {
	h = Object_handle(ei); 
	s = Sphere_segment(p,p_res,c);
	ip = p_res;
	s_init = true;
      }
    }

    // TODO: start-end point of s on cl

    if(this->has_shalfloop()) {
      Sphere_circle cl(this->shalfloop()->circle());
      if(!s_init || s.is_long()) {
	if(cl.has_on(p)) {
	  ip = p.antipode();
	  return Object_handle(SHalfloop_handle(this->shalfloop()));
	} else 	  
	  s = Sphere_segment(p,p.antipode(),c);
      }
      Sphere_point p_res;
      CGAL_NEF_TRACEN("do intersect " << cl << ", " << s);
      if(!do_intersect_internally(cl,s,p_res))
	return h;
      /*
      if(p_res == p.antipode()) // does this happen ? test has_on for p/p.antipode ?
	p_res = p;
      */
      CGAL_NEF_TRACEN("found intersection point " << p_res);
      CGAL_assertion_code(Sphere_segment testseg(p,p_res,c));
      CGAL_assertion(!testseg.is_long());
      if (start_inclusive || p != p_res) {
	ip = p_res;
	return Object_handle(SHalfloop_handle(this->shalfloop()));
      }
    }

    return h;
  }

  void marks_of_halfspheres(std::vector<Mark>& mohs, int offset, 
			    int axis=2);
  void marks_of_halfspheres(Mark& unten, Mark& oben, int axis=2);

  /*{\Mimplementation Naive query operations are realized by checking
  the intersection points of the $1$-skeleton of the plane map |P| with
  the query segments $s$. This method takes time linear in the size $n$
  of the underlying plane map without any preprocessing.}*/
}; // SM_point_locator<SM_decorator>

template <typename D>
void SM_point_locator<D>::
marks_of_halfspheres(std::vector<Mark>& mohs, int offset, int axis) {
  Mark lower, upper;
  marks_of_halfspheres(lower, upper, axis);
  mohs[offset] = lower;
  mohs[offset+1] = upper;
}

template <typename D>
void SM_point_locator<D>::
marks_of_halfspheres(Mark& lower, Mark& upper, int axis) {

  CGAL_NEF_TRACEN("marks_of_halfspheres ");

  Sphere_point y_minus;
  if(axis!=1) 
    y_minus = Sphere_point(0,-1,0);
  else
    y_minus = Sphere_point(0,0,1);
  Object_handle h = locate(y_minus);
  SFace_handle f;
  if ( CGAL::assign(f,h) ) { 
    CGAL_NEF_TRACEN("on face " << mark(f));
    lower = upper = mark(f);
    return;
  }

  SHalfedge_handle e;
  if ( CGAL::assign(e,h) ) { 
    CGAL_assertion(e->circle().has_on(y_minus));
    Sphere_point op(CGAL::ORIGIN+e->circle().orthogonal_vector());
    CGAL_NEF_TRACEN("on edge "<<op);
    if (axis==0 && ((op.z() < 0) || ((op.z() == 0) && (op.x() < 0)))) e = e->twin();
    if (axis==1 && ((op.x() > 0) || ((op.x() == 0) && (op.y() < 0)))) e = e->twin();
    if (axis==2 && ((op.x() > 0) || ((op.x() == 0) && (op.z() < 0)))) e = e->twin();
    upper = e->incident_sface()->mark();
    lower = e->twin()->incident_sface()->mark();
    return;
  }

  SHalfloop_handle l;
  if ( CGAL::assign(l,h) ) {
    CGAL_assertion(l->circle().has_on(y_minus));
    Sphere_point op(CGAL::ORIGIN+l->circle().orthogonal_vector());
    CGAL_NEF_TRACEN("on loop "<<op);
    if (axis==0 && ((op.z() < 0) || ((op.z() == 0) && (op.x() < 0)))) l = l->twin();
    if (axis==1 && ((op.x() > 0) || ((op.x() == 0) && (op.y() < 0)))) l = l->twin();
    if (axis==2 && ((op.x() > 0) || ((op.x() == 0) && (op.z() < 0)))) l = l->twin();
    upper = l->incident_sface()->mark();
    lower = l->twin()->incident_sface()->mark();
    return;
  }

  Sphere_circle c;
  switch(axis) {
  case 0: c = Sphere_circle(1,0,0); break;
  case 1: c = Sphere_circle(0,1,0); break;
  case 2: c = Sphere_circle(0,0,1); break;
  }
  Sphere_direction right(c),left(c.opposite());
  bool collinear(false);
  SVertex_handle v;
  if ( CGAL::assign(v,h) ) {
    CGAL_assertion(v->point()==y_minus);
    if(is_isolated(v))
      upper = lower = mark(v->incident_sface());
    else {
      e = out_wedge(v,left,collinear); 
      if ( collinear ) upper = e->twin()->incident_sface()->mark();
      else upper = e->incident_sface()->mark();
      e = out_wedge(v,right,collinear); 
      if ( collinear ) lower = e->twin()->incident_sface()->mark();
      else lower = e->incident_sface()->mark();
    }
  }
}

CGAL_END_NAMESPACE
#endif // CGAL_SM_POINT_LOCATOR_H