pm_point_locator.h

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

    Halfedge_const_iterator e_res;
    for(e = this->halfedges_begin(); e != this->halfedges_end(); ++(++e)) {
      Segment es = segment(e);
      int o1 = K.orientation(ss,K.source(es));
      int o2 = K.orientation(ss,K.target(es));
      if ( o1 == -o2 && o1 != 0 &&
           K.orientation(es, K.source(ss)) ==
           - K.orientation(es, K.target(ss)) ) {
        // internal intersection
        CGAL_NEF_TRACEN("candidate "<<es);
        Point p_res = K.intersection(s,es);
        e_res = (o2 > 0 ? e : twin(e));
        // o2 > 0 => te left of s and se right of s => p left of e
        if ( M(e_res) ) {
          h = Object_handle(e_res);
          ss = K.construct_segment(p,p_res);  
        } else if ( M(face(twin(e_res))) ) {
          h = Object_handle(face(twin(e_res)));
          ss = K.construct_segment(p,p_res);  
        }
      }
    }

    return h;
  }


  // C++ is really friendly:
  #define USECMARK(t) const Mark& mark(t h) const { return Base::mark(h); }
  #define USEMARK(t)  Mark& mark(t h) const { return Base::mark(h); }
  USEMARK(Vertex_handle)
  USEMARK(Halfedge_handle)
  USEMARK(Face_handle)
  USECMARK(Vertex_const_handle)
  USECMARK(Halfedge_const_handle)
  USECMARK(Face_const_handle)
  #undef USEMARK
  #undef USECMARK
  /*{\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.}*/
}; // PM_naive_point_locator<PM_decorator_,Geometry_>


/*{\Moptions print_title=yes }*/ 
/*{\Msubst 
PM_decorator_#PMD
Geometry_#GEO
}*/
/*{\Manpage {PM_point_locator}{PMD,GEO} 
{Point location in plane maps via LMWT}{PL}}*/ 
/*{\Mdefinition An instance |\Mvar| of data type |\Mname|
encapsulates point location queries within a plane map |P|. The two
template parameters are specified via concepts. |PMD| must be a model
of the concept |PMDecorator| as described in the appendix.  |GEO| must
be a model of the concept |AffineGeometryTraits_2| as described in the
appendix. For a specification of plane maps see also the concept of
|PMConstDecorator|.}*/

/*{\Mgeneralization PMD^#PM_naive_point_locator<PMD,GEO>}*/

template <typename PM_decorator_, typename Geometry_>
class PM_point_locator : public 
  PM_naive_point_locator<PM_decorator_,Geometry_> {
protected:
  typedef PM_naive_point_locator<PM_decorator_,Geometry_> Base;
  typedef PM_point_locator<PM_decorator_,Geometry_> Self;
  Base CT;
  #ifdef CGAL_USING_PPL
  typedef PM_persistent_PL_traits<Base>  PMPPLT; 
  typedef PointLocator<PMPPLT>           PMPP_locator;
  PMPP_locator* pPPL;
  #define LOCATE_IN_TRIANGULATION pPPL->locate_down
  #else
  #define LOCATE_IN_TRIANGULATION walk_in_triangulation
  #endif

public:

  typedef typename Base::Decorator Decorator;
  typedef typename Base::Plane_map Plane_map;
  typedef typename Base::Mark Mark;
  typedef typename Base::Geometry Geometry;
  typedef typename Base::Point Point;
  typedef typename Base::Segment Segment;
  typedef typename Base::Direction Direction;
  typedef typename Base::Object_handle Object_handle;
  typedef typename Base::Vertex_handle Vertex_handle;   
  typedef typename Base::Halfedge_handle Halfedge_handle;   
  typedef typename Base::Face_handle Face_handle;     
  typedef typename Base::Vertex_const_handle Vertex_const_handle; 
  typedef typename Base::Halfedge_const_handle Halfedge_const_handle; 
  typedef typename Base::Face_const_handle Face_const_handle; 
  typedef typename Base::Vertex_iterator Vertex_iterator;
  typedef typename Base::Halfedge_iterator Halfedge_iterator;
  typedef typename Base::Face_iterator Face_iterator;
  typedef typename Base::Vertex_const_iterator Vertex_const_iterator;
  typedef typename Base::Halfedge_const_iterator Halfedge_const_iterator;
  typedef typename Base::Face_const_iterator Face_const_iterator;
  typedef typename Base::Halfedge_around_vertex_circulator Halfedge_around_vertex_circulator;
  typedef typename Base::Halfedge_around_vertex_const_circulator Halfedge_around_vertex_const_circulator;
  typedef typename Base::Halfedge_around_face_circulator Halfedge_around_face_circulator;
  typedef typename Base::Halfedge_around_face_const_circulator Halfedge_around_face_const_circulator;

#ifndef CGAL_CFG_USING_BASE_MEMBER_BUG_3
  using Base::K;
  using Base::number_of_vertices;
  using Base::faces_begin;
#endif

  /*{\Mtypes 2}*/
  /*{\Mtext All local types of |PM_naive_point_locator| are inherited.}*/
  typedef std::pair<Vertex_const_handle,Face_const_handle>   VF_pair;
  typedef std::pair<Halfedge_const_handle,Face_const_handle> EF_pair;

  struct CT_link_to_original : Decorator { // CT decorator
    const Decorator& Po;
    CT_link_to_original(const Decorator& P, const Decorator& Poi) 
      : Decorator(P), Po(Poi) {}
    void operator()(Vertex_handle vn, Vertex_const_handle vo) const
    { Face_const_handle f;
      if ( Po.is_isolated(vo) ) f = Po.face(vo);
      geninfo<VF_pair>::create(info(vn));
      geninfo<VF_pair>::access(info(vn)) = VF_pair(vo,f); 
      CGAL_NEF_TRACEN("linking to org "<<PV(vn));
    }
    void operator()(Halfedge_handle hn, Halfedge_const_handle ho) const
    { geninfo<EF_pair>::create(info(hn));
      geninfo<EF_pair>::access(info(hn)) = EF_pair(ho,Po.face(ho));
      CGAL_NEF_TRACEN("linking to org "<<PE(hn));
    }
  };

protected:
  Vertex_const_handle input_vertex(Vertex_const_handle v) const
  { return geninfo<VF_pair>::const_access(CT.info(v)).first; }

  Halfedge_const_handle input_halfedge(Halfedge_const_handle e) const
  { return geninfo<EF_pair>::const_access(CT.info(e)).first; }

  Face_const_handle input_face(Halfedge_const_handle e) const
  { return geninfo<EF_pair>::const_access(CT.info(e)).second; }


  Object_handle input_object(Vertex_const_handle v) const
  { return Object_handle(input_vertex(v)); }

  Object_handle input_object(Halfedge_const_handle e) const 
  { Halfedge_const_handle e_org = input_halfedge(e);
    if ( e_org != Halfedge_const_handle() )
      return Object_handle( e_org );
    // now e_org is not existing
    return Object_handle( input_face(e) );
  }

  /*{\Mimplementation 
  The efficiency of this point location module is mostly based on
  heuristics. Therefore worst case bounds are not very expressive. The
  query operations take up to linear time for subsequent query
  operations though they are better in practise. They trigger a one-time
  initialization which needs worst case $O(n^2)$ time though runtime
  tests often show subquadratic results. The necessary space for the
  query structure is subsumed in the storage space $O(n)$ of the input
  plane map. The query times are configuration dependent. If LEDA is
  present then point location is done via the slap method based on
  persistent dictionaries.  Then $T_{pl}(n) = O( \log(n) )$. If CGAL is
  not configured to use LEDA then point location is done via a segment
  walk in the underlying convex subdivision of $P$. In this case
  $T_{pl}(n)$ is the number of triangles crossed by a walk from the
  boundary of the structure to the query point. The time for the ray
  shooting operation $T_{rs}(n)$ is the time for the point location
  $T_{pl}(n)$ plus the time for the walk in the triangulation that is
  superimposed to the plane map. Let's consider the plane map edges as
  obstacles and the additional triangulation edges as non-obstacle
  edges. Let's call the sum of the lengths of all edges of the
  triangulation its weight. If the calculated triangulation
  approximates\footnote{The calculation of general
  minimum-weight-triangulations is conjectured to be NP-complete and
  locally-minimum-weight-triangulations that we use are considered good
  approximations.} the minimum weight triangulation of the obstacle set
  then the stepping quotient\footnote {The number of non-obstacle edges
  crossed until an obstacle edge is hit.} for a random direction of the
  ray shot is expected to be $O( \sqrt{n} )$.}*/


  struct CT_new_edge : Decorator { 
    const Decorator& _DP;
    CT_new_edge(const Decorator& CT, const Decorator& DP) : 
      Decorator(CT), _DP(DP) {}
    void operator()(Halfedge_handle& e) const
    { Halfedge_handle e_from = previous(e);
      Face_const_handle f;
      if ( is_closed_at_source(e) ) // source(e) was isolated before
        f = geninfo<VF_pair>::access(info(source(e))).second;
      else
        f = geninfo<EF_pair>::access(info(e_from)).second;
      mark(e) = _DP.mark(f);
      geninfo<EF_pair>::create(info(e));
      geninfo<EF_pair>::create(info(twin(e)));

      geninfo<EF_pair>::access(info(e)).first = 
      geninfo<EF_pair>::access(info(twin(e))).first = 
        Halfedge_const_handle();

      geninfo<EF_pair>::access(info(e)).second = 
      geninfo<EF_pair>::access(info(twin(e))).second = f;
      CGAL_NEF_TRACEN("CT_new_edge "<<PE(e));
    }
  };

  void triangulate_CT() const
  {
    CGAL_NEF_TRACEN("triangulate_CT");
    typedef CGAL::Constrained_triang_traits<
      Decorator,Geometry,CT_new_edge> NCTT;
    typedef CGAL::generic_sweep<NCTT> Constrained_triang_sweep;
    CT_new_edge NE(CT,*this);
    Constrained_triang_sweep T(NE,CT.plane_map(),this->K); T.sweep();
  }

  void minimize_weight_CT() const
  { CGAL_NEF_TRACEN("minimize_weight_CT");
    if ( this->number_of_vertices() < 2 ) return;
    std::list<Halfedge_handle> S;
    /* We maintain a stack |S| of edges containing diagonals 
       which might have to be flipped. */
    int flip_count = 0;
    Halfedge_iterator e;
    for (e = CT.halfedges_begin(); e != CT.halfedges_end(); ++(++e)) {
      Halfedge_const_handle e_org = input_halfedge(e);
      if ( e_org != Halfedge_const_handle() ) 
        continue;
      S.push_back(e);
    }

    while ( !S.empty() ) {
      Halfedge_handle e = S.front(); S.pop_front();
      Halfedge_handle r = twin(e);
      Halfedge_const_handle e_org = input_halfedge(e);
      if ( e_org != Halfedge_const_handle() ) 
        continue;
      Halfedge_handle e1 = next(r);
      Halfedge_handle e3 = next(e);
      // e1,e3: edges of quadrilateral with diagonal e

      Point a = point(source(e1));
      Point b = point(target(e1));
      Point c = point(source(e3));
      Point d = point(target(e3));

      if (! (this->K.orientation(b,d,a) > 0 && // left_turn
             this->K.orientation(b,d,c) < 0) ) // right_turn
        continue;

      if ( this->K.first_pair_closer_than_second(b,d,a,c) ) { // flip
        CGAL_NEF_TRACEN("flipping diagonal of quadilateral"<<a<<b<<c<<d);
        Halfedge_handle e2 = next(e1);
        Halfedge_handle e4 = next(e3);
        S.push_back(e1); 
        S.push_back(e2); 
        S.push_back(e3); 
        S.push_back(e4); 
        flip_diagonal(e);  
        flip_count++;
      }


    }
    CGAL_NEF_TRACEN("  flipped "<<flip_count);
  }

public:
  /*{\Mcreation 3}*/

  PM_point_locator() { 
    #ifdef CGAL_USING_PPL
    pPPL = 0;
    #endif
 
  }

  /*{\Moptions constref=yes}*/
  PM_point_locator(const Plane_map& P, const Geometry& k = Geometry());
  /*{\Mcreate constructs a point locator working on |P|.}*/

  ~PM_point_locator();

  /*{\Moperations 2.5 0.5}*/

  const Decorator& triangulation() const { return CT; }
  /*{\Mop access to the constrained triangulation structure that
  is superimposed to |P|.}*/
  /*{\Moptions constref=no}*/

  
  Object_handle locate(const Point& p) const
  /*{\Mop returns a generic handle |h| to an object (vertex, halfedge,
  face) of |P| which contains the point |p| in its relative
  interior.}*/
  { 
    Object_handle h = LOCATE_IN_TRIANGULATION(p);
    Vertex_const_handle v_triang;
    if ( assign(v_triang,h) ) {
      return input_object(v_triang);
    } 
    Halfedge_const_handle e_triang;
    if ( assign(e_triang,h) ) {
      Halfedge_const_handle e = input_halfedge(e_triang);
      if ( e == Halfedge_const_handle() ) // inserted during triangulation
        return Object_handle(input_face(e_triang)); 
      int orientation_ = this->K.orientation(segment(e),p);
      if ( orientation_ == 0 ) return Object_handle(e);
      if ( orientation_ < 0 )  return Object_handle(face(twin(e)));
      if ( orientation_ > 0 )  return Object_handle(face(e));
    } 
    CGAL_assertion(!check_tag(typename Is_extended_kernel<Geometry>::value_type()));
    return Face_const_handle(faces_begin());
    //    assert(0); return h; // compiler warning
  }

  bool ray_shoot_from_outer_facet(Segment& , object_kind& , 
				  Vertex_const_handle &, 
				  Halfedge_const_handle& ,
				  const Tag_true& ) const {
    return false;
  }