snc_const_decorator.h

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

	CGAL_assertion_msg(0,"Damn wrong handle");
    }
    return f_visible;
  }

  Halffacet_const_handle get_visible_facet( const Halfedge_const_handle e,
                                      const Ray_3& ray) const {
   //{\Mop when one shoot a ray |ray| in order to find the facet below to
   //  an object, and an edge |e| is hit, we need to choose one of the two 
   //  facets in the adjacency list of |e| that could be 'seen'  from the
   //  piercing point of the |ray| on the local (virtual) view  of |e|
   //  \precondition |ray| target belongs to |e|. } 

    SM_const_decorator SD(&*e->source());
    if( SD.is_isolated(e))
      return Halffacet_const_handle();
    
    // We search for the plane in the adjacency list of e, which is closest
    // to the ray. The cross product of the direction of e and the orthogonal
    // vector of a plane gives us a vector vec0/vec1 on the plane of the facet 
    // and orthogonal to e, pointing inside of the facet.

    Vector_3 ev(segment(e).to_vector()), rv(ray.to_vector());
    SHalfedge_around_svertex_const_circulator sh(SD.first_out_edge(e));
    Halffacet_const_handle res = sh->facet(); 
    Vector_3 vec0(cross_product(ev,res->plane().orthogonal_vector()));
    /* // probably incorrect assertion
    CGAL_assertion_code
      (Sphere_segment _ess( SD.source(sh)->point(), 
			    SD.source(next(sh))->point(),
			    SD.circle(sh)));
    CGAL_assertion( _ess.has_on(vec0));
    */
    SHalfedge_around_svertex_const_circulator send(sh);
    CGAL_NEF_TRACEN("initial face candidate "<< res->plane()<<" with vector  "<<vec0);

    // We compare the vectors vec0/vec1 of the facets. The one that is nearest 
    // to pointing in the opposite direction of the ray, is chosen. The 
    // respective facet is the nearest to the ray.

    sh++;
    CGAL_For_all(sh,send) {
      Vector_3 vec1(cross_product(ev,sh->facet()->plane().orthogonal_vector()));
      CGAL_NEF_TRACEN("test face candidate "<< sh->facet()->plane()<<" with vector  "<<vec1);
      FT sk0(rv*vec0),  sk1(rv*vec1);
      if(sk0<=FT(0) && sk1>=FT(0))
        continue;
      if(sk0>=FT(0) && sk1<=FT(0)) {
        res = sh->facet(); 
	vec0 = vec1;
	continue;
      }

      // We have to comapare the two skalar products sk0 and sk1. Therefore 
      // we have to normalize the input vectors vec0 and vec1, which means
      // that we have to divide them by their lengths len0 and len1. 
      // To cicumvent irrational numbers, we sqaure the whole inequality.

      FT len0 = vec0.x()*vec0.x()+vec0.y()*vec0.y()+vec0.z()*vec0.z();
      FT len1 = vec1.x()*vec1.x()+vec1.y()*vec1.y()+vec1.z()*vec1.z();
      FT diff = len0*sk1*sk1 - len1*sk0*sk0;

      // if sk0<0 (and therefore sk1<0) both vectors point in a good direction.
      // Therefore we take the one pointing more in the good direction.
      // if sk0>0 (and therefore sk1>0) both vectors point in a bad direction.
      // Therefore we take the one pointing less in the bad direction.

      if((sk0>FT(0) && diff<FT(0)) || (sk0<FT(0) && diff>FT(0))) {
        res = sh->facet();
	vec0 = vec1;
      }
    }
 
    // We have to check which of the two halffacet is visible from
    // the ray. 

    if(rv*res->plane().orthogonal_vector() > FT(0))
      res = res->twin();

    CGAL_NEF_TRACEN("return "<<res->plane());
    return res; // never reached
  }
  
  Halffacet_const_handle get_visible_facet(Halffacet_const_handle f,
				      const Ray_3& ray) const 
    /*{\Mop when one shoot a ray |ray| in order to find the facet below to
      an object, and a facet |f| is hit, we need to choose the right facet
      from the halffacet pair |f| that  could be 'seen'  from the
      piercing point of the |ray| on the local (virtual) view  of |f|.
      \precondition |ray| target belongs to |f| and the intersection between
      |ray| and is not coplanar with |f|. }*/ {

    Halffacet_const_handle f_visible = f;
    if( f_visible->plane().has_on_negative_side(ray.source()))
      f_visible = f_visible->twin();
    CGAL_assertion( f_visible->plane().has_on_positive_side(ray.source()));
    return f_visible;
  }

  Halffacet_const_handle get_visible_facet( const Vertex_const_handle v, 
				      const Segment_3& ray) const 
    /*{\Mop when one shoots a ray |ray| in order to find the facet below to
      an object, and vertex |v| is hit, we need to choose one of the facets
      in the adjacency list of |v| such that it could be 'seen' from the
      piercing point of the |ray| on the sphere map on |v|.  We make it just
      locating the sphere facet |sf| pierced by |ray| and taking the adjacent 
      facet to one of the sphere segments on the boundary of |sf|.
      \precondition |ray| target is on |v| and the intersection between
      |ray| and the 2-skeleton incident to v is empty. }*/ {

    Halffacet_const_handle f_visible;
    CGAL_assertion( ray.source() != v->point());
    CGAL_assertion( ray.has_on(v->point()));
    Sphere_point sp(ray.source() - v->point());
    CGAL_NEF_TRACEN( "Locating "<<sp <<" in "<<v->point());
    CGAL_assertion(Infi_box::degree(sp.hx()) < 2 && 
		   Infi_box::degree(sp.hy()) < 2 && 
		   Infi_box::degree(sp.hz()) < 2 && 
		   Infi_box::degree(sp.hw()) == 0);
    sp = Infi_box::simplify(sp);
    CGAL_NEF_TRACEN( "Locating "<<sp <<" in "<<v->point());
    SM_point_locator L(v);
    Object_handle o = L.locate(sp);

    SFace_const_handle sf;
    CGAL_assertion(CGAL::assign(sf,o));
    CGAL::assign( sf, o);

    SFace_cycle_const_iterator fc = sf->sface_cycles_begin(),
      fce = sf->sface_cycles_end();
    if( is_empty_range( fc, fce)) {
	CGAL_NEF_TRACEN( "no adjacent facets were found.");
	f_visible =  Halffacet_const_handle();
    }
    else {
      if (fc.is_shalfedge()) {
      SHalfedge_const_handle se(fc); 
      CGAL_NEF_TRACEN( "adjacent facet found (SEdges cycle).");
	CGAL_NEF_TRACEN("se"<<PH(se));
	f_visible = se->twin()->facet();
	CGAL_NEF_TRACEN("f_visible"<<&f_visible);
      }
      else if (fc.is_shalfloop()) {
      SHalfloop_const_handle sl(fc);
      CGAL_NEF_TRACEN( "adjacent facet found (SHalfloop cycle).");
      f_visible = sl->twin()->facet();
      }
      else if(fc.is_svertex()) {
	CGAL_NEF_TRACEN( "no adjacent facets were found (but incident edge(s)).");
	f_visible = Halffacet_const_handle();
      }
      else
	CGAL_assertion_msg(0,"Damn wrong handle");
    }
    return f_visible;
  }

  /*
  Halffacet_const_handle get_visible_facet( const Halfedge_const_handle e,
                                      const Segment_3& ray) const {
    //{\Mop when one shoot a ray |ray| in order to find the facet below to
    //  an object, and an edge |e| is hit, we need to choose one of the two 
   //   facets in the adjacency list of |e| that could be 'seen'  from the
   //   piercing point of the |ray| on the local (virtual) view  of |e|
   //   \precondition |ray| target belongs to |e|. } 

    CGAL_assertion(false);

    SM_const_decorator SD;
    if( SD.is_isolated(e))
      return Halffacet_const_handle();

    Halffacet_const_handle res = facet(sh);    

    Vector_3 ed(segment(e).to_vector());
    Vector_3 ev(segment(e).to_vector()), rv(ray.to_vector());
    SHalfedge_around_svertex_const_circulator sh(SD.first_out_edge(e)), send(sh);
    Vector_3 vec0(cross_product(ev,res->plane().orthogonal_vector()));
    CGAL_NEF_TRACEN("initial face candidate "<< res->plane());
   
    sh++;
    CGAL_For_all(sh,send) {
    Vector_3 vec1(cross_product(ev,sh->plane().orthogonal_vector()));
      RT sk0(rv*vec0),  sk1(rv*vec1);
      if(sk0<0 && sk1>0)
        continue;
      if(sk0>0 && sk1<0) {
        res = facet(sh);
	continue;
      }

      RT len0 = vec0.x()*vec0.x()+vec0.y()*vec0.y()+vec0.z()*vec0.z();
      RT len1 = vec1.x()*vec1.x()+vec1.y()*vec1.y()+vec1.z()*vec1.z();
      RT sq0 = sk0 * sk0;
      RT sq1 = sk1 * sk1;     
      RT diff = len0*sq1 - len1*sq0;

      if((sk0 > 0 && diff<0) || (sk0 < 0 && diff>0))
        res = facet(sh);
    }

    return Halffacet_const_handle(); // never reached
  }
  */

  Halffacet_const_handle get_visible_facet( const Halffacet_const_handle f,
				      const Segment_3& ray) const 
    //{\Mop when one shoot a ray |ray| in order to find the facet below to
    //  an object, and a facet |f| is hit, we need to choose the right facet
    //  from the halffacet pair |f| that  could be 'seen'  from the
    //  piercing point of the |ray| on the local (virtual) view  of |f|.
    //  \precondition |ray| target belongs to |f| and the intersection between
    //  |ray| and is not coplanar with |f|. } 
    {
      Halffacet_const_handle f_visible = f;
      CGAL_assertion( !f_visible->plane().has_on(ray.source()));
      if( f_visible->plane().has_on_negative_side(ray.source()))
	f_visible = f->twin();
      CGAL_assertion( f_visible->plane().has_on_positive_side(ray.source()));
      return f_visible;
    }
};

template <typename EW>
template <typename Visitor>
void SNC_const_decorator<EW>::
visit_shell_objects(SFace_const_handle f, Visitor& V) const
{ 
  std::list<SFace_const_handle> SFaceCandidates;
  std::list<Halffacet_const_handle> FacetCandidates;
  CGAL::Unique_hash_map<SFace_const_handle,bool> DoneSF(false);
  CGAL::Unique_hash_map<Vertex_const_handle,bool> DoneV(false);
  CGAL::Unique_hash_map<SVertex_const_handle,bool> DoneSV(false);
  CGAL::Unique_hash_map<Halffacet_const_handle,bool> DoneF(false);
  SFaceCandidates.push_back(f);  DoneSF[f] = true;
  while ( true ) {
    if ( SFaceCandidates.empty() && FacetCandidates.empty() ) break;
    if ( !FacetCandidates.empty() ) {
      Halffacet_const_handle f = *FacetCandidates.begin();
      FacetCandidates.pop_front();
      V.visit(f); // report facet
      Halffacet_cycle_const_iterator fc;
      CGAL_forall_facet_cycles_of(fc,f) {
        if (fc.is_shalfedge() ) {
	  SHalfedge_const_handle e(fc);
	  SHalfedge_const_handle she;
          SHalfedge_around_facet_const_circulator ec(e),ee(e);
          CGAL_For_all(ec,ee) { she = ec->twin();
            if ( DoneSF[she->incident_sface()] ) continue;
            SFaceCandidates.push_back(she->incident_sface());
            DoneSF[she->incident_sface()] = true;
          }
        } else if (fc.is_shalfloop() ) { 
	  SHalfloop_const_handle l(fc);
	  SHalfloop_const_handle ll = l->twin();
          if ( DoneSF[ll->incident_sface()] ) continue;
          SFaceCandidates.push_back(ll->incident_sface());
          DoneSF[ll->incident_sface()] = true;
        } else CGAL_assertion_msg(0,"Damn wrong handle.");
      }
    }
    if ( !SFaceCandidates.empty() ) {
      SFace_const_handle sf = *SFaceCandidates.begin();
      SFaceCandidates.pop_front();
      V.visit(sf);
      if ( !DoneV[sf->center_vertex()] )
        V.visit(sf->center_vertex()); // report vertex
      DoneV[sf->center_vertex()] = true;
      //      SVertex_const_handle sv;
      SM_const_decorator SD(&*sf->center_vertex());
      /*      
      CGAL_forall_svertices(sv,SD){
	if(SD.is_isolated(sv) && !DoneSV[sv])
	  V.visit(sv);
      }
      */
      SFace_cycle_const_iterator fc;
      CGAL_forall_sface_cycles_of(fc,sf) {
        if (fc.is_shalfedge() ) {
	  SHalfedge_const_handle e(fc);
	  SHalfedge_around_sface_const_circulator ec(e),ee(e);
          CGAL_For_all(ec,ee) { 
	    V.visit(SHalfedge_const_handle(ec));
            SVertex_const_handle vv = ec->twin()->source();
            if ( !SD.is_isolated(vv) && !DoneSV[vv] ) {
              V.visit(vv); // report edge
              DoneSV[vv] = DoneSV[vv->twin()] = true;
            }
            Halffacet_const_handle f = ec->twin()->facet();
            if ( DoneF[f] ) continue;
            FacetCandidates.push_back(f); DoneF[f] = true;
          }
        } else if (fc.is_svertex() ) {
	  SVertex_const_handle v(fc);
          if ( DoneSV[v] ) continue; 
          V.visit(v); // report edge
	  V.visit(v->twin());
          DoneSV[v] = DoneSV[v->twin()] = true;
	  CGAL_assertion(SD.is_isolated(v));
	  SFaceCandidates.push_back(v->twin()->incident_sface());
	  DoneSF[v->twin()->incident_sface()]=true;
          // note that v is isolated, thus twin(v) is isolated too
	  //	  SM_const_decorator SD;
	  //	  SFace_const_handle fo;
	  //	  fo = v->twin()->incident_sface();
	  /*
	  if(SD.is_isolated(v)) 
	    fo = v->source()->sfaces_begin();
	  else
	    fo = v->twin()->incident_sface();
	  */
        } else if (fc.is_shalfloop() ) {
	  SHalfloop_const_handle l(fc);
	  V.visit(l);
          Halffacet_const_handle f = l->twin()->facet();
          if ( DoneF[f] ) continue;
          FacetCandidates.push_back(f);  DoneF[f] = true;
        } else CGAL_assertion_msg(0,"Damn wrong handle.");
      }
    }
  }
}


CGAL_END_NAMESPACE
#endif //CGAL_SNC_CONST_DECORATOR_H