%------------------------------------------------------------------------------
%KILLSTART REP
%LDEL TRACE.*?\)\;
\documentclass[a4paper]{article}
\usepackage{MyLweb,version}
\excludeversion{ignoreindiss}
\excludeversion{ignore}
\input{defs}
\input{nef_macros}

\begin{document}

\title{Point location in plane maps}
\author{Michael Seel}
\date{\today}
\maketitle 
\tableofcontents 

%KILLEND REP

\section{The Manual Page}

\input manpages/PM_point_locator.man

\section{Implementation}

This section describes two point location and ray shooting modules
working in plane maps that are decorated according to our plane map
decorator model. The first module offers naive point location without
any preprocessing of the plane map. The second module implements point
location and ray shooting for the case of iterated queries and trades
query time for a one time preprocessing phase. In the center of our
approach is a constrained triangulation of the obstacle set (the input
plane map). For the point location we offer an optimal solution with
respect to space and runtime based on a slap subdivision with
persistent dictionaries from LEDA. As CGAL demands to minimize cross
library dependency we offer a default fall-back solution based on
walks in the triangulation. For ray shooting in general obstacle sets
the theory recommends to use minimum weight triangulations. However as
minimum weight triangulations are difficult to calculate we restrict
our effort to optimize the weight of the constrained triangulation
only locally. We will present the corresponding theoretical basics in
the section where we explain the construction of the triangulation.
We start with the naive module and append the triangulation module.


\subsection{Point location and ray shooting done naively}

In this section we implement the class |PM_naive_point_locator<>| that
wraps naive point location functionality. We first sketch the class
design consisting of the local types and the interface. Then we
present the two main interface operations for point location and ray
shooting.

A point locator |PM_naive_point_locator<>| is a generalization of a
plane map decorator |PM_decorator_| where we inherit the decorator
interface.  We obtain the geometry used for point location from
|Geometry_|.
<<Naive point locator>>=
/*{\Moptions print_title=yes }*/ 
/*{\Msubst 
PM_decorator_#PMD
Geometry_#GEO
}*/
/*{\Manpage {PM_naive_point_locator}{PMD,GEO}
{Naive point location in plane maps}{PL}}*/
/*{\Mdefinition An instance |\Mvar| of data type |\Mname|
encapsulates naive point location queries within a plane map |P|.  The
two template parameters are specified via concepts. |PM_decorator_|
must be a model of the concept |PMDecorator| as described in the
appendix.  |Geometry_| 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}*/

template <typename PM_decorator_, typename Geometry_>
class PM_naive_point_locator : public PM_decorator_ {
protected:
  typedef PM_decorator_ Base;
  typedef PM_naive_point_locator<PM_decorator_,Geometry_> Self;

  const Geometry_& K;
public:
  <<PL naive local types>>
  <<PL naive helpers>>
  <<PL naive interface>>
}; // PM_naive_point_locator<PM_decorator_,Geometry_>


@ We transport types from |PM_decorator_| and |Geometry_| into the
local scope.  A CGAL |Object| is a general wrapper class for any
type. Unwrapping can be done typesave by means of C++ dynamic
casts. The |Object_handle| type is just a generalized CGAL |Object|
where we add the |NULL| equality test.
<<PL naive local types>>=
/*{\Mtypes 5}*/
typedef PM_decorator_                 Decorator;
/*{\Mtypemember equals |PM_decorator_|.}*/
typedef typename Decorator::Plane_map Plane_map;
/*{\Mtypemember the plane map type decorated by |Decorator|.}*/
typedef typename Decorator::Mark      Mark;
/*{\Mtypemember the attribute of all objects (vertices, edges, 
faces).}*/

typedef Geometry_                       Geometry;
/*{\Mtypemember equals |Geometry_|.}*/
typedef typename Geometry_::Point_2     Point;
/*{\Mtypemember the point type of the geometry kernel.\\
\require |Geometry::Point_2| equals |Plane_map::Point|.}*/
typedef typename Geometry_::Segment_2   Segment;
/*{\Mtypemember the segment type of the geometry kernel.}*/
typedef typename Geometry_::Direction_2 Direction;

/*{\Mtext Local types are handles, iterators and circulators of the 
following kind: |Vertex_const_handle|, |Vertex_const_iterator|, 
|Halfedge_const_handle|, |Halfedge_const_iterator|, |Face_const_handle|, 
|Face_const_iterator|.}*/

typedef CGAL::Object_handle Object_handle;
/*{\Mtypemember a generic handle to an object of the underlying plane
map. The kind of the object |(vertex, halfedge,face)| can be determined and
the object assigned by the three functions:\\ 
|bool assign(Vertex_const_handle& h, Object_handle o)|\\ 
|bool assign(Halfedge_const_handle& h, Object_handle o)|\\ 
|bool assign(Face_const_handle& h, Object_handle o)|\\ where each 
function returns |true| iff the assignment of |o| to |h| was valid.}*/

@ \begin{ignoreindiss}
<<PL naive local types>>=
#define CGAL_USING(t) typedef typename PM_decorator_::t t
CGAL_USING(Vertex_handle);   
CGAL_USING(Halfedge_handle);   
CGAL_USING(Face_handle);     
CGAL_USING(Vertex_const_handle); 
CGAL_USING(Halfedge_const_handle); 
CGAL_USING(Face_const_handle); 
CGAL_USING(Vertex_iterator);
CGAL_USING(Halfedge_iterator);
CGAL_USING(Face_iterator);
CGAL_USING(Vertex_const_iterator);
CGAL_USING(Halfedge_const_iterator);
CGAL_USING(Face_const_iterator);
CGAL_USING(Halfedge_around_vertex_circulator);
CGAL_USING(Halfedge_around_vertex_const_circulator);
CGAL_USING(Halfedge_around_face_circulator);
CGAL_USING(Halfedge_around_face_const_circulator);
#undef CGAL_USING

@ 
\end{ignoreindiss}%
The naive interface provides construction, access to the attributes of
type |Mark|, point location, and ray shooting. Note that the ray
shooting is parameterized by a template predicate that allows
adaptation of the termination criterion of the ray-shot.
<<PL naive interface>>=
/*{\Mcreation 3}*/

PM_naive_point_locator() : Base() {}

/*{\Moptions constref=yes}*/
PM_naive_point_locator(const Plane_map& P, const Geometry& k = Geometry()) : 
  Base(const_cast<Plane_map&>(P)), K(k) {}
/*{\Mcreate constructs a point locator working on |P|.}*/
/*{\Moptions constref=no}*/
/*{\Moperations 2.5 0.5}*/

const Mark& mark(Object_handle h) const
/*{\Mop returns the mark associated to the object |h|.}*/
{ Vertex_const_handle v; 
  Halfedge_const_handle e; 
  Face_const_handle f;
  if ( assign(v,h) ) return mark(v);
  if ( assign(e,h) ) return mark(e);
  if ( assign(f,h) ) return mark(f);
  CGAL_assertion_msg(0,
  "PM_point_locator::mark: Object_handle holds no object.");
#if !defined(__BORLANDC__)
  return mark(v); // never reached
#endif
}

<<NPL location method>>
<<NPL ray shooting method>>

@ \begin{ignoreindiss}
<<PL naive interface>>=
// 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
@ \end{ignoreindiss}
\subsubsection*{Point location} 
The following point location scheme runs in linear time without any
preprocessing.  The segment |s| is required to intersect the skeleton
of our plane map |P|. We just check if the location |p = source(s)| is
part of the embedding of any vertex or uedge of the $1$-skeleton of
|P|.
<<NPL location method>>=

Object_handle locate(const Segment& s) const
/*{\Mop returns a generic handle |h| to an object (vertex, halfedge,
face) of the underlying plane map |P| which contains the point |p =
s.source()| in its relative interior. |s.target()| must be a point
such that |s| intersects the $1$-skeleton of |P|.}*/
{ TRACEN("locate naivly "<<s);
  if (number_of_vertices() == 0) 
    CGAL_assertion_msg(0,"PM_naive_point_locator: plane map is empty.");
  Point p = K.source(s);
  Vertex_const_iterator vit;
  for(vit = vertices_begin(); vit != vertices_end(); ++vit) {
    if ( p == point(vit) ) return Object_handle(vit);
  }
  Halfedge_const_iterator eit;
  for(eit = halfedges_begin(); eit != halfedges_end(); ++(++eit)) {
    // we only have to check each second halfedge
    if ( K.contains(segment(eit),p) ) 
      return Object_handle(eit);
  }
  Vertex_const_handle v_res;
  Halfedge_const_handle e_res;
  <<NPL determine closest object intersecting s>>
  if ( e_res != Halfedge_const_handle() )
    return Object_handle((Face_const_handle)(face(e_res)));
  else
    return Object_handle((Face_const_handle)(face(v_res)));
}

@ If |p| is located in the interior of a face |f| we can identify |f|
only by any object in its boundary. We do this by a thought ray shoot
along |s| through all obstacles of the $1$-skeleton.  Note that we are
sure to hit at least one object of the $1$-skeleton of the plane
map. Let |q = s.target()|. We check the intersection points of the
segment |s| with all skeleton objects. As soon as we find an object we
prune |s|, store the object that gives possibly access to the face
containing |p|, and iterate.

An edge |e| that intersects |s| in its relative interior is trivial to
handle. In this case we just clip |s| by the intersection point and go
on with our search to possibly find a closer edge to |p|.

Isolated vertices are easy to treat too. They can be checked directly
as they carry a link to their incident face. Non-isolated vertices
don't carry the face information directly. There we have to find an
incident edge whose incident face intersects |s|. Intuitively when |s|
contains a non-isolated vertex |v| we use a symbolic perturbation
scheme. Let's assume that our direction is down along the negative
$y$-axis.  We just shift the whole scene infinitesimally left. Then we
only hit edges. But of course the closest edge in $A(v)$\footnote{the
adjacency list of $v$.} has to be determined not by proximity but by
the direction of edges in $A(v)$. Now we are interested in an edge |e|
of $A(v)$ that has a minimal counterclockwise angle to |s|. |e| bounds
the wedge containing |p|. See figure \figref{walkseg}.

\displayeps{walkseg}{Configurations of |s| and the $1$-skeleton of |P|}{9cm}

Assume $d$ to be a direction anchored in $v$.  The operation
|out_wedge(v,d,collinear)| returns the halfedge |e| bounding a wedge
in between two neighbored edges in the adjacency list of |v| that
contains |d|.  If |d| extends along an edge then |e| is this edge. If
|d| extends into the interior of such a wedge then |e| is the first
edge hit when |d| is rotated clockwise. As a precondition we ask |v|
not to be isolated.
<<PL naive helpers>>=
Halfedge_const_handle out_wedge(Vertex_const_handle v, 
  const Direction& d, bool& collinear) const
/*{\Xop returns a halfedge |e| bounding a wedge in between two
neighbored edges in the adjacency list of |v| which contains |d|.
If |d| extends along a edge then |e| is this edge. If |d| extends
into the interior of such a wedge then |e| is the first edge hit
when |d| is rotated clockwise. \precond |v| is not isolated.}*/
{ TRACEN("out_wedge "<<PV(v));
  assert(!is_isolated(v));
  collinear=false;
  Point p = point(v);
  Halfedge_const_handle e_res = first_out_edge(v);
  Direction d_res = direction(e_res);
  Halfedge_around_vertex_const_circulator el(e_res),ee(el);
  CGAL_For_all(el,ee) {
    if ( K.strictly_ordered_ccw(d_res, direction(el), d) )
      e_res = el; d_res = direction(e_res);
  }
  TRACEN("  determined "<<PE(e_res)<<" "<<d_res);
  if ( direction(cyclic_adj_succ(e_res)) == d ) {
    e_res = cyclic_adj_succ(e_res);
    collinear=true;
  }
  TRACEN("  wedge = "<<PE(e_res)<<" "<<collinear);
  return e_res;
}

@ We use the above operation when the vertex is not isolated.
<<NPL determine closest object intersecting s>>=
Segment ss = s; // we shorten the segment iteratively
Direction dso = K.construct_direction(K.target(s),p), d_res;
CGAL::Unique_hash_map<Halfedge_const_handle,bool> visited(false);
for(vit = vertices_begin(); vit != vertices_end(); ++vit) {
  Point p_res, vp = point(vit);
  if ( K.contains(ss,vp) ) {
    TRACEN(" location via vertex at "<<vp);
    ss = K.construct_segment(p,vp); // we shrink the segment
    if ( is_isolated(vit) ) {
      v_res = vit; e_res = Halfedge_const_handle();
    } else { // not isolated
      bool dummy;
      e_res = out_wedge(vit,dso,dummy);
      Halfedge_around_vertex_const_circulator el(e_res),ee(el);
      CGAL_For_all(el,ee) 
        visited[el] = visited[twin(el)] = true;
      /* e_res is now the counterclockwise maximal halfedge out
         of v just before s */
      if ( K.orientation(p,vp,point(target(e_res))) < 0 ) // right turn
        e_res = previous(e_res);
      // correction to make e_res visible from p
      TRACEN("  determined "<<PE(e_res));
    }
  }
}

@ Now we treat the left case of figure \figref{walkseg}. We check if
the current segment |ss| intersects a halfedge and shorten it
accordingly.  Note that the edge |e_res| is chosen such that |p| is
left of it. The |visited| flags allow us to skip all edges that have
already been examined incident to vertices on |s| which saves the
geometric computations.
<<NPL determine closest object intersecting s>>=
for (eit = halfedges_begin(); eit != halfedges_end(); ++eit) {
  if ( visited[eit] ) continue;
  Point se = point(source(eit)),
        te = point(target(eit));
  int o1 = K.orientation(ss,se);
  int o2 = K.orientation(ss,te);
  if ( o1 == -o2 && // internal intersection
       K.orientation(se,te,K.source(ss)) != 
       K.orientation(se,te,K.target(ss)) ) { 
      TRACEN(" location via halfedge "<<segment(eit));
    Point p_res = K.intersection(s,segment(eit));
    ss = K.construct_segment(p,p_res); 
    e_res = (o2 > 0 ? eit : twin(eit));
    // o2>0 => te left of s and se right of s => p left of e
    visited[eit] = visited[twin(eit)] = true;
    TRACEN("  determined "<<PE(e_res)<<" "<<mark(e_res));
    TRACEN("             "<<mark(face(e_res)));
  }
}

@ \subsubsection*{Ray shooting} 

The ray shooting has to determine a closest object |h| (vertex,
halfedge, face) of our plane map |P| which fulfills the predicate
|M|. We first locate the source point |p = s.source()| by our point
location method. If we already hit an object that fulfills |M| then we
are done. Otherwise we look for vertices or halfedges along |s|.
<<NPL ray shooting method>>=

template <typename Object_predicate>
Object_handle ray_shoot(const Segment& s, const Object_predicate& M) const
/*{\Mop returns an |Object_handle o| which can be converted to a
|Vertex_const_handle|, |Halfedge_const_handle|, |Face_const_handle|
|h| as described above.  The object predicate |M| has to have function
operators \\
|bool operator() (const Vertex_/Halfedge_/Face_const_handle&)|.\\
The object returned is intersected by the segment |s| and has
minimal distance to |s.source()| and |M(h)| holds on the converted
object. The operation returns the null handle |NULL| if the ray shoot
along |s| does not hit any object |h| of |P| with |M(h)|.}*/
{ TRACEN("naive ray_shoot "<<s);
  assert( !K.is_degenerate(s) );
  Point p = K.source(s);
  Segment ss(s);
  Direction d = K.construct_direction(K.source(s),K.target(s));
  Object_handle h = locate(s);
  Vertex_const_handle v; 
  Halfedge_const_handle e; 
  Face_const_handle f;
  if ( assign(v,h) && M(v) ||
       assign(e,h) && M(e) ||
       assign(f,h) && M(f) ) return h;
  h = Object_handle(); 
  TRACEN("not contained");
  <<NPL segment s contains vertex>>
  <<NPL segment s intersects halfedge>>
  return h;
}

@ Look at a vertex |v| at position |pv| on segment |s|. If |M(v)| we
can shorten |s| and iterate. If not |M(v)| then a ray shoot along |s|
can still hit an object at |v| where |M| is fulfilled in form of an
outgoing edge or the face in a wedge between two edges. At each vertex
we look only for the part of |s| between |pv| and |s.target()|.
<<NPL segment s contains vertex>>=
for (v = vertices_begin(); v != vertices_end(); ++v) {
  Point pv = point(v);
  if ( !K.contains(ss,pv) ) continue;
  TRACEN("candidate "<<pv);
  if ( M(v) ) {
    h = Object_handle(v);     // store vertex
    ss = K.construct_segment(p,pv); // shorten
    continue;
  }
  // now we know that v is not marked but on s
  bool collinear;
  Halfedge_const_handle e = out_wedge(v,d,collinear);
  if ( collinear ) { 
    if ( M(e) ) {
      h = Object_handle(e);
      ss = K.construct_segment(p,pv);
    }
    continue;
  }
  if ( M(face(e)) ) {
    h = Object_handle(face(e));
    ss = K.construct_segment(p,pv);
  }
} // all vertices

@ All halfedges can shorten |ss| when intersecting the segment in its
interior and when the halfedge fulfills |M|.
<<NPL segment s intersects halfedge>>=
Halfedge_const_iterator e_res;
for(e = halfedges_begin(); e != 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
    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);  
    }
  }
}

@ \subsubsection*{Correctness and Runtimes}

We summarize the facts in a lemma.

\begin{lemma}
If the segment |s| intersects the skeleton of |P| then
|point_location(s)| returns a generic handle |h| to the object
(vertex, halfedge, face) whose embedding in the plane contains the
source of |s|. The running time is linear in the size of the plane map
|P|.

If the segment |s| intersects the skeleton of |P| then
|ray_shoot(s,M)| determines a generic handle |h| which is |NULL| if no
object of |P| that intersects |s| fulfills |M(h)|\footnote {Note that
the predicate |M| is defined on vertex, halfedge, and face handles for
performance reasons. However we sloppily write |M(h)| on a generic handle
wrapping one of the above.}  and which can be converted to a
corresponding vertex, halfedge, or face otherwise. In the latter case
the object fulfills |M(h)| and additionally has minimal distance to
|source(s)|.  The running time is linear in the size of the plane map.
\end{lemma}

\begin{proof}
For the point location case note that we iterate over all objects of
the skeleton of |P| twice. If |p = source(s)| is contained in a
skeleton object then |h| is determined already in $\langle$\textit{NPL
locate}$\rangle$ in one of the two iterations. If |p| is contained in the
plane in between the $1$-skeleton then the closest skeleton object
along |s| defines the face. If the closest object is a vertex we
determine this vertex in the vertex iteration of $\langle$\textit{NPL
determine closest object intersecting s}$\rangle$.  For non-isolated
vertices we also iterate over all edges in the adjacency list of the
vertices (in |out_wedge()|) to determine a halfedge visible from |p|
along |s|. If the closest object is a halfedge intersecting |s| in its
iterior the halfedge iteration of $\langle$\textit{NPL determine
closest object intersecting s}$\rangle$ does the job.  As all
geometric tests take constant time this gives us the linear runtime
bound.  Note that starting from |s| we only shorten |s|. Due to our
precondition |s| indeed is shortened at least once. This ensures that
we obtain an object which allows us to obtain the face containing |p|.

For the ray shooting note that if |p| already is contained in an object
marked correctly then the point location already determines this
object. If not then a similar iteration over all objects now only
considering the predicate |M| does the job. If no object is determined
then |h| is the |NULL| handle. The runtime stays linear.
\end{proof}

\begin{ignoreindiss}
<<PL naive helpers>>=
Segment segment(Halfedge_const_handle e) const
{ return K.construct_segment(point(source(e)), point(target(e))); }

Direction direction(Halfedge_const_handle e) const
{ return K.construct_direction(point(source(e)),point(target(e))); }

<<PL naive interface>>=  
/*{\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.}*/
@ \end{ignoreindiss}
\subsection{Point location and ray shooting based on constrained 
triangulations}

Our second point location structure is based on a convex subdivision
of the plane map by means of a constrained triangulation\footnote{The
constrained triangulation of plane map is a triangulation of the
vertices such that all edges are part of the triangulation}. We create
the point location structure starting from a plane map. The edges and
isolated vertices of the plane map are our obstacle set. We reference
the plane map again via a decorator.  We present the class layout, and
the implementation of the two main interface operations |locate()| and
|ray_shoot()|.  The construction of the constrained triangulation and
the integration of persistent dictionaries as a location structure is
shown in the following sections.

The main class is derived from the naive point location class to
inherit its primitive methods. The class |PM_point_locator| decorates
the input plane map |P|. We store the additional triangulation via an
additional decorator |CT| in the local scope. Note that we talk about
two plane map structures here. The plane map decorated by the |*this|
object which we also call the input structure and the plane map
storing the constrained triangulation of the input structure. We have
to maintain links such that each object in |CT| knows the object of
|*this| that supports\footnote{an object $a$ supports $b$ if
|embedding(a)| contains |embedding(b)|.} it.  We use an optional point
location structure based on persistent dictionaries that is
configuration dependent. All properties of that optional module are
explained in section \ref{point location with persistent
dictionaries}.
<<Triangulation point locator>>= 
/*{\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;
  <<TPL persistent module members>>
public:
  <<TPL local types>>
protected:
  <<TPL protected members>>
public:
  <<TPL interface>>
}; // PM_point_locator<PM_decorator_,Geometry_>


@ \begin{ignoreindiss}
<<TPL local types>>=
#define CGAL_USING(t) typedef typename Base::t t
CGAL_USING(Decorator);
CGAL_USING(Plane_map);
CGAL_USING(Mark);
CGAL_USING(Geometry);
CGAL_USING(Point);
CGAL_USING(Segment);
CGAL_USING(Direction);
CGAL_USING(Object_handle);
CGAL_USING(Vertex_handle);   
CGAL_USING(Halfedge_handle);   
CGAL_USING(Face_handle);     
CGAL_USING(Vertex_const_handle); 
CGAL_USING(Halfedge_const_handle); 
CGAL_USING(Face_const_handle); 
CGAL_USING(Vertex_iterator);
CGAL_USING(Halfedge_iterator);
CGAL_USING(Face_iterator);
CGAL_USING(Vertex_const_iterator);
CGAL_USING(Halfedge_const_iterator);
CGAL_USING(Face_const_iterator);
CGAL_USING(Halfedge_around_vertex_circulator);
CGAL_USING(Halfedge_around_vertex_const_circulator);
CGAL_USING(Halfedge_around_face_circulator);
CGAL_USING(Halfedge_around_face_const_circulator);
#undef CGAL_USING

@ \end{ignoreindiss} 
Our constrained triangulation |CT| decorates the $1$-skeleton of a
plane map consisting just of vertices and edges. We save the face
objects there. Each vertex |v| of |CT| is linked to an object of type
|VF_pair| which stores the input vertex supporting |v| and potentially
an incident face of the input plane map if |v| is isolated (before the
triangulation process). Each halfedge |e| of |CT| is linked to an
object of type |EF_pair| which stores the input halfedge supporting
|e| and the face incident to |e| in the input structure.
<<TPL local types>>=
/*{\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;

@ We associate the above pair structures to the skeleton objects via
the available |GenPtr| slot. We have one additional pointer in each
object that we can use for temporary information. We use a special
class |geninfo| for the expansion and the access. See the manual
page of |geninfo| in the appendix for its usage.
<<TPL protected members>>=
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; }


@ The public interface provides construction, access to the underlying
triangulation, point location, and ray shooting. Note that the ray
shooting is parameterized by a template predicate that allows
adaptation of the termination criterion of the ray shooting walk.
<<TPL interface>>= 
/*{\Mcreation 3}*/

PM_point_locator() { 
  <<TPL persistent module default init>> 
}

/*{\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}*/

<<TPL location method>>
<<TPL ray shooting method>>

Object_handle walk_in_triangulation(const Point& p) const;


@ \subsubsection*{Construction and Destruction}

We use a constrained triangulation of the input plane map to do ray
shooting walks. We leave the input structure |P| unchanged. Therefore
we clone the structure into one which we extend to a constrained
triangulation. The cloning can be enriched by actions on the cloned
objects via a data accessor. See the description of the operation
|clone_skeleton| in the manual page of |PMDecorator| for the concept
of this data accessor.  Our data accessor |CT_link_to_original| does
the linkage of input objects to cloned objects. Additionally it takes
care that the cloned objects know the faces incident to their
supporting input objects. For the used data association method see the
manual page of |geninfo|.
<<TPL local types>>=
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); 
    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));
    TRACEN("linking to org "<<PE(hn));
  }
};

@ Note that we have two decorators: the |*this| object decorates the
plane map |P|, the internal decorator |CT| works on a cloned plane map
on the heap.  We extend the cloned representation to a constrained
triangulation. Afterwards we locally optimize the weight of the
constrained triangulation to obtain a better expected performance for
the ray shooting walk. See the implementation of
|minimize_weight_CT()| for more information on this.  Depending on the
presence of |LEDA| we also create an additional point location
structure based on the slap method using LEDA's persistent
dictionaries.  See $\langle$\textit{TPL persistent module pm
init}$\rangle$ in section \ref{point location with persistent
dictionaries} for further information.
<<TPL construction>>=
template <typename PMD, typename GEO>
PM_point_locator<PMD,GEO>::
PM_point_locator(const Plane_map& P, const Geometry& k) :
  Base(P,k), CT(*(new Plane_map),k)

{ TRACEN("PM_point_locator construction");
  CT.clone_skeleton(P,CT_link_to_original(CT,*this));
  triangulate_CT();
  minimize_weight_CT();
  <<TPL persistent module pm init>>
}

@ Destruction mirrors the allocations within |CT| above.  We discard
the extended information slots in all objects (vertices and
halfedges). Then we discard all objects.  Finally we free the memory
for the |CT| plane map on the heap. Optionally we clean up the memory
of the persistent module.
<<TPL destruction>>=
template <typename PMD, typename GEO>
PM_point_locator<PMD,GEO>::
~PM_point_locator()
{ TRACEN("clear_static_point_locator");
  Vertex_iterator vit, vend = CT.vertices_end();
  for (vit = CT.vertices_begin(); vit != vend; ++vit) {
    geninfo<VF_pair>::clear(CT.info(vit));
  }
  Halfedge_iterator eit, eend = CT.halfedges_end();
  for (eit = CT.halfedges_begin(); eit != eend; ++eit) {
    geninfo<EF_pair>::clear(CT.info(eit));
  }
  CT.clear();
  delete &(CT.plane_map());
  <<TPL persistent module destruction>>
}

@ \subsubsection*{Point location} 

The following chunk interfaces the basic point location layer
depending on the presence of LEDA. |LOCATE_IN_TRIANGULATION| is either
|walk_in_triangulation| or a call to the persistent point location
structure in $\langle$\textit{TPL persistent module members}$\rangle$
of section \ref{point location with persistent dictionaries}.
<<TPL location method>>=

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_ = 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));
  }
  assert(0); return h; // compiler warning
}

@ The following function |walk_in_triangulation(q)| provides a
fall-back solution to point location. It returns an |Object_handle f|
which is either a vertex |v| with |point(v) == q| or a halfedge |e|
where |e| contains |q| in its relative interior, or one of the two
triangles incident to |e| (or |twin(e)| respectively) contains
|q|. This semantics is equivalent to the function from the generic
point location framework programmed by Sven Thiel \cite{thiel-mt}
which allows us to obtain a halfedge or vertex that contains |q| or
that is vertically below |q|.

We walk from the embedding |p = point(v)| of the first vertex $v =
|vertices_begin()|$ of our structure to the point |q| that we expect to
lie in the interior of the triangulation. We examine the skeleton of
the triangulation. We can thereby hit vertices and edges, the latter
in two ways. We walk along an edge |e| when |segment(e)| and |s = pq|
overlap, or we cross |e| in its relative interior. Our walk is thus an
iteration storing a vertex or an edge and a flag how we traverse
them. We store the oriented halfedge |e| either in direction of |s| or
such that |q| is left of |e|.
<<TPL interface>>=
enum object_kind { VERTEX, EDGE_CROSSING, EDGE_COLLINEAR };
<<TPL walk in triangulation>>=
template <typename PMD, typename GEO>
typename PM_point_locator<PMD,GEO>::Object_handle  
PM_point_locator<PMD,GEO>::walk_in_triangulation(const Point& q) const
{ 
  TRACEN("walk in triangulation "<<q);
  Vertex_const_handle v = CT.vertices_begin();
  Halfedge_const_handle e;
  Point p = CT.point(v);
  if ( p == q ) return Object_handle(v);
  Segment s = K.construct_segment(p,q);
  Direction dir = K.construct_direction(p,q);
  object_kind current = VERTEX;
  while (true) switch ( current ) {
    case VERTEX:
      <<TPL segment walk via a vertex>>
      break;
    case EDGE_CROSSING:
      <<TPL segment walk via a crossing edge>>
      break;
    case EDGE_COLLINEAR:
      <<TPL segment walk via a collinear edge>>
      break;
  }
#if !defined(__BORLANDC__)
  return Object_handle(); // never reached warning acceptable
#endif
}

@ We walk along |s| and hit a vertex |v|. If we have reached |q| we
are done.  Otherwise we walk the adjacency list to find either the
wedge between two edges via which we leave the vertex or we find a
halfedge collinear to |s|. The adjacency list walk is encapsulated in
|out_wedge|.  We obtain a halfedge |e_out| out of |v|, either
collinear or bounding the wedge right of |s|. If |e_out| is collinear
to |s| it is the next object.  Otherwise we take the edge |next(e_out)|
closing the wedge to a triangle which we certainly hit next on our
walk along |s|.
<<TPL segment walk via a vertex>>=
{ 
  TRACEN("vertex "<<CT.point(v));
  if ( CT.point(v) == q ) 
    return Object_handle(v); // stop walking at q
  bool collinear;
  Halfedge_const_handle e_out = CT.out_wedge(v,dir,collinear);
  if (collinear) // ray shoot via e_out
  { e = e_out; current = EDGE_COLLINEAR; }
  else  // ray shoot in wedge left of e_out
  { e = CT.twin(CT.next(e_out)); current = EDGE_CROSSING; }
}

@ When we cross an edge |e| from right to left it could be that |q| is
on |e|, or contained in the incident triangle right of it.  If neither
is true we continue through the incident triangle left of |e|. There
we can hit one of two other edges or the vertex opposite to |e|. Note
however that this vertex can lie beyond |q| on the directed line
through |p| and |q|.
<<TPL segment walk via a crossing edge>>=
{ TRACEN("crossing edge "<<CT.segment(e));
  if ( !(K.orientation(CT.segment(e),q) > 0) ) // q not left of e
    return Object_handle(e);
  Vertex_const_handle v_cand = CT.target(CT.next(e));
  int orientation_ = K.orientation(p,q,CT.point(v_cand));
  switch( orientation_ ) {
    case 0:  // collinear 
      if ( K.strictly_ordered_along_line(p,q,CT.point(v_cand)) ) 
        return Object_handle(e);
      v = v_cand; current = VERTEX; break;
    case +1: // left_turn
      e = twin(next(e)); current = EDGE_CROSSING; break;
    case -1: 
      e = twin(previous(e)); current = EDGE_CROSSING; break;
  }
}

@ Hitting a collinear edge is easy to treat. If it does not contain
|q| its target is the next object.
<<TPL segment walk via a collinear edge>>=
{ TRACEN("collinear edge "<<CT.segment(e));
  if ( K.strictly_ordered_along_line(
         CT.point(CT.source(e)),q,CT.point(CT.target(e))) ) 
    return Object_handle(e);
  v = CT.target(e); current = VERTEX;
}

@ \subsubsection*{Correctness and Runtime}

We summarize the execution properties of the above code in a lemma.
\begin{lemma}
Point location in a plane map |P| of size $O(n)$ based on the
persistent module determines the object (vertex, halfedge, face)
containing |p| in expected time $O( \log n )$. The fall back solution
|walk_in_triangulation(p)| determines the object (vertex, halfedge,
face) containing |p| in time $O(n)$.
\end{lemma}

\begin{proof} We omit the details of the persistent point location
package.  The elaborate treatment can be found in \cite{thiel-mt}.
The constrained triangulation |CT| of the input plane map |P| is
asymptotically not larger than the plane map $O(n)$.  We sketch the
intuition about the persistent point location module (PPL). The PPL
partitions the plane into at most $2n+1$ slaps defined by vertical
lines through the $n$ vertices of |CT|. The stripes between the lines
and the lines are considered as slaps.  Each slap is subdivided into
convex patches by objects of |CT|. The patches are stored in a
(persistent) dictionary via their bounding skeleton objects.  A
location within a slap can be determined within logarithmic query
time.  Now a point location can be done in time $O(2 \log n)$ first by
a binary search for the correct slap and second by a query of the
corresponding dictionary. The construction of the PPL can be done in
$O( n \log n )$ by the plane sweep framework described in
\cite{thiel-mt} and is used in section \ref{point location with
persistent dictionaries}. The partially persistent dictionary
implementation of S. Thiel stays within the linear space bound of |P|
and |CT|.

For the |walk_in_triangulation(q)| query operation note that the
initialization is well defined. We start with the first vertex of the
triangulation. Then each iteration reaches either |q| or one of the
states |VERTEX|, |EDGE_CROSSING|, |EDGE_COLLINEAR|. Convince yourself
by revisiting the code chunks that each of the possible three
configurations covers all geometric possibilities that the segment
walk along |s| can encounter.
\end{proof}


@ \subsubsection*{Ray shooting} 

The ray shooting along a segment |s| works as follows. We first locate
|p=s.source()| in our underlying triangulation. We might hit a vertex,
an edge or a face. In the triangulation we do a segment walk starting
in |p| along |s| similar to the one in |walk_in_triangulation()|. Only
that we have the additional termination criterion of the object
predicate |M|. We examine the skeleton of the triangulation. We can
thereby hit vertices and edges, the latter in two kinds. We walk along
an edge |e| when |segment(e)| and |s| overlap, or we cross |e| in its
relative interior. Our walk is an iteration storing a vertex or an
edge and a flag how we traverse them. We store the oriented halfedge
either in direction of |s| or such that |s.target()| is left of |e|.

First we have to get to a starting point by point location.  Then we
iterate until we hit an object of our input structure as determined
by the object predicate |M|.
<<TPL ray shooting method>>=

template <typename Object_predicate>
Object_handle ray_shoot(const Segment& s, const Object_predicate& M) const
/*{\Mop returns an |Object_handle o| which can be converted to a
|Vertex_const_handle|, |Halfedge_const_handle|, |Face_const_handle|
|h| as described above.  The object predicate |M| has to have 
function operators\\
|bool operator() (const Vertex_/ Halfedge_/Face_const_handle&) const|.\\
The object returned is intersected by the segment |s| and has minimal
distance to |s.source()| and |M(h)| holds on the converted object. The
operation returns the null handle |NULL| if the ray shoot along |s|
does not hit any object |h| of |P| with |M(h)|.}*/
{ TRACEN("ray_shoot "<<s);
  CGAL_assertion( !K.is_degenerate(s) );
  Point p = K.source(s), q = K.target(s);
  Direction d = K.construct_direction(p,q); 
  Vertex_const_handle v;
  Halfedge_const_handle e;
  object_kind current;
  Object_handle h = LOCATE_IN_TRIANGULATION(p);
  <<TPL ray shoot initialization>>
  while (true) switch ( current ) {
    case VERTEX:
      <<TPL ray shoot hits a vertex>>
      break;
    case EDGE_CROSSING:
      <<TPL ray shoot hits a crossing edge>>
      break;
    case EDGE_COLLINEAR:
      <<TPL ray shoot hits a collinear edge>>
      break;
  } 
  // assert(0); return h; // compiler warning
}

@ If we hit a vertex we are done.
<<TPL ray shoot initialization>>=
if ( assign(v,h) ) {
  TRACEN("located vertex "<<PV(v));
  current = VERTEX;
}  

@ If the result is a halfedge |e|, it can contain |p| or |e| (or its
twin) can bound the triangle that contains |p| in its interior.
<<TPL ray shoot initialization>>=
if ( assign(e,h) ) {
  TRACEN("located edge "<<PE(e));
  int orientation_ = K.orientation( segment(e), p);
  if ( orientation_ == 0 ) { // p on segment
    <<edge e contains the starting point p>>
  } else { // p not on segment, thus in triangle
    if ( orientation_ < 0  ) e = CT.twin(e);
    // now p left of e
    <<edge e bounds the triangle containing p>>
  }
}

@ In this case we have to handle the cases that |s| overlaps
|segment(e)| or that we shoot into a face adjacent to |e|.
<<edge e contains the starting point p>>=
TRACEN("on edge "<<PE(e));
if ( d == CT.direction(e) ) 
{ current = EDGE_COLLINEAR; }
else if ( d == CT.direction(CT.twin(e)) ) 
{ e = CT.twin(e); current = EDGE_COLLINEAR; }
else { // crossing
  current = EDGE_CROSSING;
  if ( !(K.orientation(CT.segment(e),q)>0) ) // not left_turn
    e = CT.twin(e); 
}

@ If the triangle (underlying face) which contains |p| fulfills |M| we
return it. Otherwise we have to find the object (vertex or edge) hit
by the ray shoot from |p| along |s|. Note that it might be that |s| is
totally contained in the triangle in which case we return the |NULL|
handle.
<<edge e bounds the triangle containing p>>=
TRACEN("in face at "<<PE(e));
if ( M(input_face(e)) ) // face mark
  return Object_handle(input_face(e));

Point p1 = CT.point(CT.source(e)), 
      p2 = CT.point(CT.target(e)), 
      p3 = CT.point(CT.target(next(e)));
int or1 = K.orientation(p,q,p1);
int or2 = K.orientation(p,q,p2);
int or3 = K.orientation(p,q,p3);
if ( or1 == 0 && !K.left_turn(p1,p2,q) )
{ v = CT.source(e); current = VERTEX; }
else if ( or2 == 0 && !K.left_turn(p2,p3,q) )
{ v = CT.target(e); current = VERTEX; }
else if ( or3 == 0 && !K.left_turn(p3,p1,q) )
{ v = CT.target(CT.next(e)); current = VERTEX; }
else if ( or2 > 0 && or1 < 0 && !K.left_turn(p1,p2,q) )
{ e = CT.twin(e); current = EDGE_CROSSING; }
else if ( or3 > 0 && or2 < 0 && !K.left_turn(p2,p3,q) )
{ e = CT.twin(CT.next(e)); current = EDGE_CROSSING; }
else if ( or1 > 0 && or3 < 0 && !K.left_turn(p3,p1,q) )
{ e = CT.twin(CT.previous(e)); current = EDGE_CROSSING; }
else return Object_handle();

@ Now we are on our walk along |s|. If we hit a vertex |v| with |M(v)|
we have found the right object. Otherwise it could be that we have
reached |q|. If not we walk the adjacency list to find either the
wedge between two edges via which we leave the vertex or we find a
halfedge collinear to |s|. The adjacency list walk is encapsulated in
|out_wedge()|. We obtain a halfedge |e_out| out of |v|, either
collinear or bounding the wedge right of |s|. Note that if the wedge
represents a face |f| with |M(f)| we are done. Otherwise we take the
edge |next(e_out)| closing the wedge to a triangle which we certainly
hit next on our walk along |s|.
<<TPL ray shoot hits a vertex>>=
{ TRACEN("vertex "<<CT.point(v));
  Vertex_const_handle v_org = input_vertex(v);
  if ( M(v_org) ) return Object_handle(v_org);
  if ( CT.point(v) == q ) return Object_handle();
  // stop walking at q, or determine next object on s:
  bool collinear;
  Halfedge_const_handle e_out = CT.out_wedge(v,d,collinear);
  if (collinear) // ray shoot via e_out
  { e = e_out; current = EDGE_COLLINEAR; }
  else { // ray shoot in wedge left of e_out
    if ( M(input_face(e_out)) )
      return Object_handle(input_face(e_out));
    e = CT.twin(CT.next(e_out)); current = EDGE_CROSSING;
  }
}

@ When we cross an edge |e| we can find either the edge or its
incident face fulfilling |M|. Note however that edges that were
created during the triangulation process are supported by the input
face which they split. Thus we just ignore them in this case. If
neither the edge nor the face fulfill |M| we continue through the
triangle that |e| bounds. There we can hit one of two other edges or
the vertex opposite to |e|.
<<TPL ray shoot hits a crossing edge>>=
{ TRACEN("crossing edge "<<segment(e));
  if ( K.orientation(CT.segment(e),q) == 0 ) 
    return Object_handle();
  Halfedge_const_handle e_org = input_halfedge(e);
  if ( e_org != Halfedge_const_handle() ) { // not a CT edge
    if ( M(e_org) ) return Object_handle(e_org);
    if ( M(face(e_org)) ) return Object_handle(face(e_org));
  }
  Vertex_const_handle v_cand = CT.target(CT.next(e));
  TRACEN("v_cand "<<PV(v_cand));
  int orientation_ = K.orientation(p,q,CT.point(v_cand));
  switch( orientation_ ) {
    case 0: 
      v = v_cand; current = VERTEX; break;
    case +1: 
      e = CT.twin(CT.next(e)); current = EDGE_CROSSING; break;
    case -1: 
      e = CT.twin(CT.previous(e)); current = EDGE_CROSSING; break;
  }
}

@ If we hit a collinear edge we check if the edge was inserted 
during the triangulation process. If it is a triangulation edge
we examine the underlying face, otherwise we check the edge
itself. The next object is its target vertex.
<<TPL ray shoot hits a collinear edge>>=
{ TRACEN("collinear edge "<<CT.segment(e));
  Halfedge_const_handle e_org = input_halfedge(e);
  if ( e_org == Halfedge_const_handle() ) { // a CT edge
    if ( M(input_face(e)) ) 
      return Object_handle(input_face(e));
  } else { // e_org is not a CT edge
    if ( M(e_org) )
      return Object_handle(e_org);
  }
  if ( K.strictly_ordered_along_line(
         CT.point(CT.source(e)),q,CT.point(CT.target(e))) ) 
    return Object_handle();
  v = CT.target(e); current = VERTEX;
}

@ \subsubsection*{Correctness and Runtimes}

The following lemma summarizes the properties of the ray-shooting
operation.
\begin{lemma}
If the segment |s| is non-degenerate, then |ray_shoot(s,M)| determines
a generic handle |h| which is |NULL| if no object of |P| that
intersects |s| fulfills |M(h)|\footnote {Note that the predicate |M| is
defined on vertex, halfedge, and face handles for performance
reasons. However we sloppily write |M(h)| on a generic handle wrapping one
of the above.}  and which can be converted to a corresponding vertex,
halfedge, or face otherwise. In the latter case the object fulfills
|M(h)| and additionally has minimal distance to |source(s)|.
The query time is that of the point location plus the time for the 
walk in |CT|. 
\end{lemma}
\begin{proof}
The starting point is always the point location of |p = source(s)|.
If the object |u| that contains |p| does not fulfill |M(u)|, then a
walk in the triangulation is started in $\langle$\textit{TPL ray
shoot}$\rangle$. This iteration is determined by one of the three
geometric configurations $|VERTEX|, |EDGE_CROSSING|,
|EDGE_COLLINEAR|$. We only have to ensure that the iteration terminates
correctly. In each case we either look for an object |u| hit by a ray
along |s| for which |M(u)| holds, or we check if we have to end our
walk as we have reached |q = target(s)|. 
\end{proof}
Note that we cannot bound the stepping for the ray shooting walk
besides the trivial linear bound (the size of the plane map and the
triangulation). However we will cite some theory and heuristic tests
below that provide the hope that the walk is much cheaper.

Finally there are two operations for the objects in the constrained
triangulation that return the corresponding objects from the input
plane map that support them.
<<TPL protected members>>=
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} )$.}*/


@ \subsection{The constrained triangulation}

\repdiss{We calculate the constrained triangulation of the cloned
skeleton by our constrained triangulation sweep as presented in
chapter \ref{constrained triangulations}}{}. The class
|Constrained_triang_traits<>| obtains three traits models |PMDEC|,
|GEOM|, and |NEWEDGE| on instantiation. The two former parameters are
the local types |PM_decorator| and |Geometry|. The latter is the data
accessor type implemented below.

All edges of the cloned skeleton carry links to the edges of |P|. We
mark additional triangulation edges by a default link to
|Halfedge_const_handle()| that can be checked after the triangulation
process to separate input halfedges (obstacles) from triangulation
halfedges. We forward input face links from existing edges of the
triangulation to newly created edges. Thereby each newly created edge
knows the face that it splits into triangles. This information
association and transfer is done by means of |CT_new_edge|.

The following function operator |operator()(Halfedge_handle&)|
attributes the newly created triangulation edge |e| in |CT| by an
additional |EF_pair| in a similar way as the cloning process did it
for the input edges. As |e| does not have a supporting halfedge in |P|
it obtains the default handle. |e| obtains the mark information from a
face that supports it with respect to the input structure.  The face
usually comes from an adjacent edge of |CT|. If there is none, then we
obtain is from the source which was isolated in the input structure
before |e| was created.
<<TPL protected members>>=
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;
    TRACEN("CT_new_edge "<<PE(e));
  }
};

@ Now completing |CT| into a full triangulation is a trivial
instantiation of the sweep class. We plug the types into
|Constrained_triang_traits|, and that type into our generic sweep
framework. Then we create the sweep object and trigger the sweep. Note
also how we transport references to the plane map structure and the
geometric kernel on construction.
<<TPL protected members>>=
void triangulate_CT() const
{
  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(),K); T.sweep();
}

@ \subsection{Minimizing the triangulation weight}

For a plane map or a triangulation let its weight be the sum of the
length of all edges.  After the calculation of the constrained
triangulation we minimize locally the weight of the triangulation by a
sequence of flipping operations starting from all non-constrained
edges. We first motivate the following procedure by some theory and
some experiments.
\begin{itemize}
\item Ray shooting by line walks in planar triangulations are
worst-case expensive. Agarwal et al \cite{agarwal95} show that such a
walk can have linear costs in the size of the triangulation though the
shoot does not hit any obstacle edge.  Aronov and Fortune show in
\cite{aronovfortune97} that average case stepping cost for random
lines can be bounded by $O( \sqrt{n} )$ for plane maps of size $n$ if
we use triangulations of minimal weight. Unfortunately such
triangulations are hard to construct in general. For general points
sets the problem is conjectured to be NP-hard. In simple polygons
there's a $O(n^3)$ solution by dynamic programming which is hardly
practical. Thus currently only approximation or even heuristics offer
practical solutions.

\item Note that the delaunay triangulation of a point set of size $n$
in general can be a bad approximation of the minimal weight
triangulation \cite{kirkpatrick80} but it has been shown to be a
$O(\log n)$ approximation \cite{lingas86} if the points are 
universally distributed. However we have no such result for the
constrained case.

\item Heuristic results on the approximation of minimum weight
triangulation are presented by Dickerson et al. in
\cite{dickerson95}. They present results that propose to use local
minimal weight triangulations as an approximations of the general
unsolved problem.

\item The running time of the following flipping algorithm is worst
case quadratic though runtime tests seem to show a subquadratic
runtime for random inputs\footnote{see the runtime results for the
flipping number in \cite{dickerson95}}. The quadratic bound for the
worst case runtime is based on a result of Hurtado et
al. \cite{hurtado96} where they show that the graph of triangulations
of a simple polygon is connected and has the complexity of
$O(n^2)$. Note that the termination is guaranteed due to the fact that
for each quadrangle of points we can only choose the minimal diagonal
once, thus we never cycle.
\end{itemize}
<<TPL protected members>>=
void minimize_weight_CT() const
{ TRACEN("minimize_weight_CT");
  if ( 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;
  <<insert all non constrained edges in S>>
  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;
    <<continue if quadrilateral with diagonal e is non-convex>>
    <<flip if there's a lighter edge>>
  }
  TRACEN("  flipped "<<flip_count);
}

@ We insert all non-constrained edges into our stack |S|. An edge is
not constrained if it has no link to a corresponding supporting input
edge in |P|.
<<insert all non constrained edges in S>>=
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);
}

@ A non-constrained edge can only be flipped if it is the diagonal of
a convex quadrilateral. We can continue if it is not.
<<continue if quadrilateral with diagonal e is non-convex>>=
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 (! (K.orientation(b,d,a) > 0 && // left_turn
       K.orientation(b,d,c) < 0) ) // right_turn
  continue;

@ We check if for our non-constrained edge |e| the complementary
diagonal has less weight. If yes we flip |e|.
<<flip if there's a lighter edge>>=
if ( K.first_pair_closer_than_second(b,d,a,c) ) { // flip
  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++;
}


@ \begin{ignoreindiss}
\subsection{The file wrapper}
<<PM_point_locator.h>>=
<<CGAL Header1>>
// file          : include/CGAL/Nef_2/PM_point_locator.h
<<CGAL Header2>>
#ifndef CGAL_PM_POINT_LOCATOR_H
#define CGAL_PM_POINT_LOCATOR_H

#include <CGAL/basic.h>
#include <CGAL/Unique_hash_map.h>
#include <CGAL/Nef_2/Constrained_triang_traits.h>
#include <CGAL/Nef_2/Object_handle.h>
#undef _DEBUG
#define _DEBUG 17
#include <CGAL/Nef_2/debug.h>
#include <CGAL/Nef_2/geninfo.h>
<<Conditional persistent point location inclusion>>

CGAL_BEGIN_NAMESPACE

<<Naive point locator>>
<<Triangulation point locator>>
<<TPL construction>>
<<TPL destruction>>
<<TPL walk in triangulation>>

CGAL_END_NAMESPACE

#endif // CGAL_PM_POINT_LOCATOR_H



@ \end{ignoreindiss}
\subsection{Point location with persistent dictionaries}
\label{point location with persistent dictionaries}

We add a dynamically allocated point location structure of type
|PointLocator<T>| working on the constrained triangulation of the
class |PM_point_locator|. 
\begin{ignoreindiss} The used persistent dictionaries are part
of LEDA starting version 4.2.
<<Conditional persistent point location inclusion>>=
#ifdef CGAL_USE_LEDA
# if CGAL_LEDA_VERSION < 500
#  include <LEDA/basic.h>
# else
#  include <LEDA/system/basic.h>
# endif
 
# if __LEDA__ > 410 && __LEDA__ < 441
#  define CGAL_USING_PPL
#  include <CGAL/Nef_2/PM_persistent_PL.h>
# endif
#endif
@ \end{ignoreindiss}%
The point location framework is added to the point location class only
if the necessary modules are present. The traits class
|PM_persistent_PL_traits| is defined below. The point location
framework |PointLocator| is described in Sven Thiel's masters thesis
\cite{thiel-mt}. 
<<TPL persistent module members>>=
#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

@ The first parameter of |PMPP_locator| is the graph structure worked
on, the second is a traits class object. We forward a reference to the
geometric kernel |K| that we already obtained for the |*this| object
of type |PM_point_locator|.
<<TPL persistent module pm init>>=
#ifdef CGAL_USING_PPL
pPPL = new PMPP_locator(CT,PMPPLT(K));
#endif

<<TPL persistent module default init>>=
#ifdef CGAL_USING_PPL
pPPL = 0;
#endif

<<TPL persistent module destruction>>=
#ifdef CGAL_USING_PPL
delete pPPL; pPPL=0; 
#endif
@ \begin{ignoreindiss}
<<Triangulation point locator>>=
#ifdef CGAL_USING_PPL
static const char* pointlocationversion = "point location via pers dicts";
#else
static const char* pointlocationversion = "point location via seg walks";
#endif

@ \end{ignoreindiss}%
The rest of this section describes a traits model for the persistent
point location framework as developed by Sven Thiel \cite{thiel-mt}.
We adapt his framework for point location via persistent dictionaries
within the constrained triangulation. Thereby we obtain a better
runtime bound than the fall-back method |walk_in_triangulation()|.

The class |PM_persistent_PL_traits| transports a bunch of types and
methods into the point location data structure |PointLocator| of
S. Thiel.  Note that the class references a plane map and a geometric
kernel, that we both forward on construction. \label{persistent point
location traits}
<<class PM_persistent_PL_traits>>=
template <typename PMPL>
struct PM_persistent_PL_traits 
{
  <<PP types in local scope>>
  <<PP iterators and conversion>>
  <<PP edge categorization>>
  <<PP x-structure interface>>
  <<PP curve interface>>
  <<PP generic location and postprocessing>>
  <<PP construction>>
};

@ \begin{ignoreindiss}
<<PM_persistent_PL.h>>=
<<CGAL Header1>>
// file          : include/CGAL/Nef_2/PM_persistent_PL.h
<<CGAL Header2>>
#ifndef CGAL_PM_PERSISTENT_PL_H
#define CGAL_PM_PM_PERSISTENT_PL_H
#include <CGAL/Nef_2/gen_point_location.h>

<<class PM_persistent_PL_traits>>

#endif // CGAL_PM_PM_PERSISTENT_PL_H

@ \end{ignoreindiss}%
The type mapping is trivial. For the persistent point location we need
the interface types |Graph|, |Node|, |Edge|. We additionally introduce
the geometry.
<<PP types in local scope>>=
typedef PMPL  Graph;
typedef typename PMPL::Vertex_const_handle   Node;
typedef typename PMPL::Halfedge_const_handle Edge;
typedef typename PMPL::Face_const_handle     Face;
typedef typename PMPL::Object_handle         Object_handle;

typedef typename PMPL::Geometry  Geometry;
typedef typename PMPL::Point     Point;
typedef typename PMPL::Segment   Segment;
const Geometry* pK;

@ We store a reference to the geometric kernel in the traits class.
<<PP construction>>=
PM_persistent_PL_traits() : pK(0) {}
PM_persistent_PL_traits(const Geometry& k) : pK(&k) {}
virtual ~PM_persistent_PL_traits() {}
virtual void sweep_begin(const Graph&) {}
virtual void sweep_moveto(const XCoord&) {}
virtual void sweep_end() {}
virtual void clear() {}

<<PP iterators and conversion>>=
typedef typename PMPL::Vertex_const_iterator NodeIterator;
NodeIterator Nodes_begin(const Graph& G) const { return G.vertices_begin(); }
NodeIterator Nodes_end(const Graph& G) const { return G.vertices_end(); }
Node toNode(const NodeIterator& nit) const { return nit; }

typedef typename PMPL::Halfedge_around_vertex_const_circulator HAVC;
struct IncEdgeIterator {
  HAVC _start, _curr;
  bool met;
  IncEdgeIterator() {}
  IncEdgeIterator(HAVC c) : 
    _start(c), _curr(c), met(false) {}
  IncEdgeIterator& operator++()
  { if (_curr==_start)
      if (!met)  { met=true; ++_curr; }
      else       { _curr=HAVC(); }
    else ++_curr;
    return *this; 
  }
  bool operator==(const IncEdgeIterator& it2) const
  { return _curr==it2._curr; }
  bool operator!=(const IncEdgeIterator& it2) const
  { return !(*this==it2); }
};
Edge toEdge(const IncEdgeIterator& eit) const { return eit._curr; }

IncEdgeIterator IncEdges_begin(const Graph& G, const Node& n) 
{ return IncEdgeIterator(HAVC(G.first_out_edge(n))); }
IncEdgeIterator IncEdges_end(const Graph& G, const Node& n)   
{ return IncEdgeIterator(); }

<<PP edge categorization>>=
enum EdgeCategory 
{ StartingNonVertical, StartingVertical, EndingNonVertical, EndingVertical };

Node opposite(const Graph& G, const Edge& e, const Node& u)
{ if ( G.source(e) == u ) return G.target(e);
  else                    return G.source(e); }

EdgeCategory ClassifyEdge(const Graph& G, const Edge& e, const Node& u)
{
  Point p_u = G.point(u);
  Point p_v = G.point(opposite(G,e,u));

  int cmpX = pK->compare_x(p_u, p_v);
  if ( cmpX < 0 ) return StartingNonVertical;
  if ( cmpX > 0 ) return EndingNonVertical;

  int cmpY = pK->compare_y(p_u, p_v); 
  assert(cmpY != 0);
  if ( cmpY < 0 ) return StartingVertical;
  return EndingVertical;
}    

@ The generic persistent point location framework maintains slaps of
the graph skeleton ordered by the x-coordinates of the embedding of
the nodes. We make the |point| type representing its x-coordinate and
use |compare_x| from the kernel to sort them.
<<PP x-structure interface>>=
typedef Point XCoord;
const XCoord getXCoord(const Point& p) const 
{ return p; }
const XCoord getXCoord(const Graph& G, const Node& n) const 
{ return G.point(n); }

class PredLessThanX {
  const Geometry* pK;
public:
  PredLessThanX() : pK(0) {}
  PredLessThanX(const Geometry* pKi) : pK(pKi) {}
  PredLessThanX(const PredLessThanX& P) : pK(P.pK) 
  { TRACEN("copy PredLessThanX"); }
  int operator() (const XCoord& x1, const XCoord& x2) const
  { return pK->compare_x(x1,x2) < 0; }
};
PredLessThanX getLessThanX() const { return PredLessThanX(pK); }

@ In our case curves are segments. We obtain them from the geometric
kernel.
<<PP curve interface>>=
// Curve connected functionality:
typedef Segment  Curve;

Curve makeCurve(const Point& p) const 
{ return pK->construct_segment(p,p); }
Curve makeCurve(const Graph& G, const Node& n) const
{ return makeCurve(G.point(n)); }
Curve makeCurve(const Graph& G, const Edge& e) const
{ Point ps = G.point(G.source(e)), pt = G.point(G.target(e));
  Curve res(G.point(G.source(e)),G.point(G.target(e)));
  if ( pK->compare_xy(ps,pt) < 0 ) res = pK->construct_segment(ps,pt);
  else                             res = pK->construct_segment(pt,ps);
  return res; 
}

struct PredCompareCurves {
 const Geometry* pK;
 PredCompareCurves() : pK(0) {}
 PredCompareCurves(const Geometry* pKi) : pK(pKi) {}
 PredCompareCurves(const PredCompareCurves& P) : pK(P.pK) {}

 int cmppntseg(const Point& p, const Curve& s) const
 { 
   if ( pK->compare_x(pK->source(s),pK->target(s)) != 0 ) // !vertical
     return pK->orientation(pK->source(s),pK->target(s), p); 
   if ( pK->compare_y(p,pK->source(s)) <= 0 ) return -1;
   if ( pK->compare_y(p,pK->target(s)) >= 0 ) return +1;
   return 0;
 }

 int operator()(const Curve& s1, const Curve& s2) const
 { 
   Point a = pK->source(s1); 
   Point b = pK->target(s1);
   Point c = pK->source(s2);
   Point d = pK->target(s2);
   if ( a==b ) 
     if ( c==d ) return pK->compare_y(a,c);
     else        return  cmppntseg(a, s2);
   if ( c==d )   return -cmppntseg(c, s1);
   // now both are non-trivial:
   int cmpX = pK->compare_x(a, c);
   if ( cmpX < 0 ) 
     return - pK->orientation(a,b,c);
   if ( cmpX > 0 ) 
     return   pK->orientation(c,d,a);

   int cmpY = pK->compare_y(a, c);
   if ( cmpY < 0 ) return -1;
   if ( cmpY > 0 )  return +1;
    
   // cmpX == cmpY == 0 => a == c
   return pK->orientation(c,d,b);
 }
};

PredCompareCurves getCompareCurves() const
{ return PredCompareCurves(pK); }

@ The generic location type |GenericLocation| is defined in the
generic point location framework. The framework allows us to introduce
a postprocessing phase after the actual location phase.  From the
location phase we obtain two locations in two slaps, which are
neigbored in the x-structure. The second location |L_plus| is non-nil
if the first slap has width zero and is defined by a node at an
x-coordinate |x|. Then |L_plus| returns an edge |e| just below $x +
\epsilon$. This allows us to extract the face that we search from |e|.
<<PP generic location and postprocessing>>=
typedef GenericLocation<Node, Edge> Location;
typedef Object_handle QueryResult;

virtual Object_handle 
PostProcess(const Location& L, const Location& L_plus, 
  const Point& p) const 
{ /* we only get an L_plus (non-nil) if L is ABOVE a vertex
     in which case we want to extract the face from the edge
     below (p+epsilon) available via L_plus. */
  if (!L_plus.is_nil()) { assert(L_plus.is_edge());
    return Object_handle(Edge(L_plus));
  } else { 
    if ( L.is_edge() ) {
      return Object_handle(Edge(L));
    }
    if ( L.is_node() ) {
      Node v(L); assert( v != Node() );
      return Object_handle(v);
    }
    return Object_handle();
  }
}

@ \begin{ignoreindiss}
\section{A Test Case}
<<PM_point_locator-test.C>>=
#include <CGAL/basic.h>
#include <CGAL/leda_integer.h>
#include <CGAL/Homogeneous.h>
#include "Affine_geometry.h"
#include <CGAL/Nef_2/HalfedgeDS_default.h>
#include <CGAL/test_macros.h>
#include <CGAL/Nef_2/HDS_items.h>
#include <CGAL/Nef_2/PM_const_decorator.h>
#include <CGAL/Nef_2/PM_decorator.h>
#include <CGAL/Nef_2/PM_io_parser.h>
#include <CGAL/Nef_2/PM_visualizor.h>
#include <CGAL/Nef_2/PM_overlayer.h>
#include <CGAL/Nef_2/PM_point_locator.h>


// KERNEL:
typedef CGAL::Homogeneous<leda_integer>   Hom_kernel;
typedef CGAL::Affine_geometry<Hom_kernel> Aff_kernel;
typedef Aff_kernel::Segment_2 Segment;

// HALFEDGE DATA STRUCTURE:
struct HDS_traits { 
  typedef Aff_kernel::Point_2 Point; 
  typedef bool Mark;
};
typedef  CGAL::HalfedgeDS_default<HDS_traits,HDS_items> Plane_map;
typedef  CGAL::PM_const_decorator<Plane_map> CDecorator;
typedef  CGAL::PM_decorator<Plane_map>       Decorator;
typedef  CGAL::PM_overlayer<Decorator,Aff_kernel> PM_aff_overlayer;

typedef CDecorator::Halfedge_const_handle  Halfedge_const_handle;
typedef CDecorator::Vertex_const_handle    Vertex_const_handle;
typedef CDecorator::Face_const_handle      Face_const_handle;
typedef Decorator::Halfedge_handle         Halfedge_handle;
typedef Decorator::Vertex_handle           Vertex_handle;

// INPUT:
typedef  std::list<Segment>::const_iterator Iterator;

struct Object_DA {
  const Decorator& D;
  Object_DA(const Decorator& Di) : D(Di) {} 
  void supporting_segment(Halfedge_handle e, Iterator it) const
  { D.mark(e) = true; }
  void trivial_segment(Vertex_handle v, Iterator it) const
  { D.mark(v) = true; }
  void starting_segment(Vertex_handle v, Iterator it) const
  { D.mark(v) = true; }
  void passing_segment(Vertex_handle v, Iterator it) const
  { D.mark(v) = true; }
  void ending_segment(Vertex_handle v, Iterator it) const
  { D.mark(v) = true; }
};

// POINT LOCATION:

typedef CGAL::PM_naive_point_locator<Decorator,Aff_kernel> PM_naive_PL;
typedef CGAL::PM_point_locator<Decorator,Aff_kernel>       PM_triang_PL;
typedef PM_naive_PL::Object_handle Object_handle;

void print(Object_handle h)
{  
  Vertex_const_handle v;
  Halfedge_const_handle e;
  Face_const_handle f;
  if (assign(v,h)) cout << " ohandle = vertex " << PV(v) << endl;
  if (assign(e,h)) cout << " ohandle = halfedge " << PE(e) << endl;
  if (assign(f,h)) cout << " ohandle = face " << PE(f->halfedge()) << endl;
}

// struct INSET {  bool operator()(bool b) const { return b; } };
struct INSET {
  const CDecorator& D;
  INSET(const CDecorator& Di) : D(Di) {}
  bool operator()(Vertex_const_handle v) const { return D.mark(v); }
  bool operator()(Halfedge_const_handle e) const { return D.mark(e); }
  bool operator()(Face_const_handle f) const { return D.mark(f); }
};

int main(int argc, char* argv[])
{
  SETDTHREAD(17);
  CGAL::set_pretty_mode(cerr);
  CGAL::Window_stream W;
  W.init(-50,50,-50,1);
  W.display();
  W.message("Insert segments to create a map.");
  Plane_map H;
  PM_aff_overlayer PMOV(H,Aff_kernel());
  CGAL::PM_visualizor<PM_aff_overlayer,Aff_kernel> V(W,PMOV);
  std::list<Segment> L;
  Segment s;
  if ( argc == 2 ) {
    std::ifstream log(argv[1]);
    while ( log >> s ) { L.push_back(s); W << s; }
  }
  while ( W >> s ) L.push_back(s);
  std::string fname(argv[0]);
  fname += ".log";
  std::ofstream log(fname.c_str());
  for (Iterator sit = L.begin(); sit != L.end(); ++sit) 
    log << *sit << " ";
  log.close();
  
  Object_DA ODA(PMOV);
  PMOV.create(L.begin(),L.end(),ODA);
  V.draw_skeleton();
  W.message("Insert one segment for point location");
  PM_naive_PL  PL1(H);
  PM_triang_PL PL2(H);
  Object_handle h;
  W >> s;
  CGAL::PM_io_parser<PM_aff_overlayer>::dump(PMOV);
  h = PL1.locate(s);
  print(h);
  h = PL2.locate(s.source());
  print(h);
  h = PL1.ray_shoot(s,INSET(PL1));
  h = PL2.ray_shoot(s,INSET(PL2));
  W.read_mouse();

  return 0;
}


@ \begin{ignore}
<<CGAL Header1>>=
// ============================================================================
//
// Copyright (c) 1997-2000 The CGAL Consortium
//
// This software and related documentation is part of an INTERNAL release
// of the Computational Geometry Algorithms Library (CGAL). It is not
// intended for general use.
//
// ----------------------------------------------------------------------------
//
// release       : $CGAL_Revision$
// release_date  : $CGAL_Date$
//
<<CGAL Header2>>=
// package       : Nef_2 
// chapter       : Nef Polyhedra
//
// source        : nef_2d/PM_point_locator_2.lw
// revision      : $Id$
// revision_date : $Date$
//
// author(s)     : Michael Seel <seel@mpi-sb.mpg.de>
// maintainer    : Michael Seel <seel@mpi-sb.mpg.de>
// coordinator   : Michael Seel <seel@mpi-sb.mpg.de>
//
// implementation: Point location module
// ============================================================================
@ \end{ignore}
\end{ignoreindiss}
%KILLSTART REP
\bibliographystyle{alpha}
\bibliography{comp_geo,geo_mod,diss}

\newpage
\section{Appendix}
\input manpages/PMConstDecorator.man
\input manpages/PMDecorator.man
\input manpages/AffineGeometryTraits_2.man
\input manpages/geninfo.man

\end{document}
%KILLEND