CGAL users discussion list

["Sebastien Loriot (GeometryFactory)" <>]

It's a bug.

Try with the attached file instead of the one in CGAL/.

Sebastien.

On 08/05/2013 02:20 PM, rcasero wrote:
> Hi,
>
> As a follow up, I have commited the Matlab MEX function that uses CGAL's
> Alpha_Shape_3 here [1] to share
>
> https://code.google.com/p/gerardus/source/browse/trunk/matlab/CgalToolbox/CgalAlphaShape3.cpp?r=1284
>
> That code builds and works. But as I have point sets with 400,000+
> points, I tried to add the flag CGAL::Fast_location to the Delaunay
> triangulation, replacing
>
> typedefCGAL::Delaunay_triangulation_3<K,Tds>Delaunay;
>
> with
>
> typedef CGAL::Delaunay_triangulation_3<K, Tds, CGAL::Fast_location>
> Delaunay;
>
> as in example"43.5.2   Building Basic Alpha Shapes for Many Points"
>
> http://www.cgal.org/Manual/latest/doc_html/cgal_manual/Alpha_shapes_3/Chapter_main.html#Subsection_43.5.2
>
> However, this seems to give trouble with the PointWithIndex class, and I
> get a compilation error. The full error is given below, but digging a
> bit into it, it seems that the problem is with
>
> /home/rcasero/Documents/gerardus/matlab/CgalToolbox/CgalAlphaShape3.cpp:175:
>    instantiated from here
> /usr/include/c++/4.4/bits/stl_uninitialized.h:74: error: no matching
> function for call to
> ‘CGAL::Point_3<CGAL::Epick>::Point_3(std::pair<CGAL::Point_3<CGAL::Epick>,
> long unsigned int>&)’
>
> Line 175 of CgalAlphaShape3.cpp is
>
> https://code.google.com/p/gerardus/source/browse/trunk/matlab/CgalToolbox/CgalAlphaShape3.cpp?r=1284#175
>
> <CODE>
> Delaunay delaunay(x.begin(), x.end());
> </CODE>
>
> So is this a limitation of the Delaunay triangulation class when using a
> PointWithIndex, that only happens with CGAL::Fast_location? Or am I
> reading the error the wrong way?
>
> Best regards,
>
> Ramon.
>
>
> <FULL ERROR>
> [100%] Building CXX object
> matlab/CgalToolbox/CMakeFiles/cgal_alpha_shape3.dir/CgalAlphaShape3.cpp.o
> In file included from /usr/include/c++/4.4/vector:64,
>                   from
> /usr/local/include/ITK-4.3/itkMetaDataDictionary.h:22,
>                   from /usr/local/include/ITK-4.3/itkObject.h:33,
>                   from /usr/local/include/ITK-4.3/itkDataObject.h:31,
>                   from /usr/local/include/ITK-4.3/itkProcessObject.h:31,
>                   from /usr/local/include/ITK-4.3/itkImageSource.h:31,
>                   from /usr/local/include/ITK-4.3/itkImportImageFilter.h:21,
>                   from
> /home/rcasero/Documents/gerardus/matlab/CgalToolbox/../MatlabImportFilter.h:52,
>                   from
> /home/rcasero/Documents/gerardus/matlab/CgalToolbox/CgalAlphaShape3.cpp:81:
> /usr/include/c++/4.4/bits/stl_uninitialized.h: In static member function
> ‘static _ForwardIterator std::__uninitialized_copy<<anonymous>
>  >::uninitialized_copy(_InputIterator, _InputIterator, _ForwardIterator)
> [with _InputIterator =
> __gnu_cxx::__normal_iterator<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int>*, std::vector<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int>, std::allocator<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int> > > >, _ForwardIterator = CGAL::Point_3<CGAL::Epick>*,
> bool <anonymous> = false]’:
> /usr/include/c++/4.4/bits/stl_uninitialized.h:117:   instantiated from
> ‘_ForwardIterator std::uninitialized_copy(_InputIterator,
> _InputIterator, _ForwardIterator) [with _InputIterator =
> __gnu_cxx::__normal_iterator<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int>*, std::vector<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int>, std::allocator<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int> > > >, _ForwardIterator = CGAL::Point_3<CGAL::Epick>*]’
> /usr/include/c++/4.4/bits/stl_uninitialized.h:257:   instantiated from
> ‘_ForwardIterator std::__uninitialized_copy_a(_InputIterator,
> _InputIterator, _ForwardIterator, std::allocator<_Tp>&) [with
> _InputIterator =
> __gnu_cxx::__normal_iterator<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int>*, std::vector<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int>, std::allocator<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int> > > >, _ForwardIterator = CGAL::Point_3<CGAL::Epick>*, _Tp
> = CGAL::Point_3<CGAL::Epick>]’
> /usr/include/c++/4.4/bits/stl_vector.h:1024:   instantiated from ‘void
> std::vector<_Tp, _Alloc>::_M_range_initialize(_ForwardIterator,
> _ForwardIterator, std::forward_iterator_tag) [with _ForwardIterator =
> __gnu_cxx::__normal_iterator<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int>*, std::vector<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int>, std::allocator<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int> > > >, _Tp = CGAL::Point_3<CGAL::Epick>, _Alloc =
> std::allocator<CGAL::Point_3<CGAL::Epick> >]’
> /usr/include/c++/4.4/bits/stl_vector.h:1002:   instantiated from ‘void
> std::vector<_Tp, _Alloc>::_M_initialize_dispatch(_InputIterator,
> _InputIterator, std::__false_type) [with _InputIterator =
> __gnu_cxx::__normal_iterator<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int>*, std::vector<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int>, std::allocator<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int> > > >, _Tp = CGAL::Point_3<CGAL::Epick>, _Alloc =
> std::allocator<CGAL::Point_3<CGAL::Epick> >]’
> /usr/include/c++/4.4/bits/stl_vector.h:303:   instantiated from
> ‘std::vector<_Tp, _Alloc>::vector(_InputIterator, _InputIterator, const
> _Alloc&) [with _InputIterator =
> __gnu_cxx::__normal_iterator<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int>*, std::vector<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int>, std::allocator<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int> > > >, _Tp = CGAL::Point_3<CGAL::Epick>, _Alloc =
> std::allocator<CGAL::Point_3<CGAL::Epick> >]’
> /home/rcasero/Documents/gerardus/cpp/src/third-party/CGAL-4.2/include/CGAL/Triangulation_hierarchy_3.h:133:
>    instantiated from ‘ptrdiff_t
> CGAL::Triangulation_hierarchy_3<Tr>::insert(InputIterator,
> InputIterator) [with InputIterator =
> __gnu_cxx::__normal_iterator<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int>*, std::vector<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int>, std::allocator<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int> > > >, Tr = CGAL::Delaunay_triangulation_3<CGAL::Epick,
> CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick,
> CGAL::Triangulation_vertex_base_with_info_3<long unsigned int,
> CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick,
> CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick,
> CGAL::Triangulation_vertex_base_with_info_3<long unsigned int,
> CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick,
> CGAL::Triangulation_ds_vertex_base_3<void> > >,
> CGAL::Boolean_tag<false>, CGAL::Boolean_tag<false> >,
> CGAL::Alpha_shape_cell_base_3<CGAL::Epick, CGAL::Default,
> CGAL::Boolean_tag<false>, CGAL::Boolean_tag<false> > > > > >,
> CGAL::Boolean_tag<false>, CGAL::Boolean_tag<false> > >,
> CGAL::Alpha_shape_cell_base_3<CGAL::Epick, CGAL::Default,
> CGAL::Boolean_tag<false>, CGAL::Boolean_tag<false> > >, CGAL::Default>]’
> /home/rcasero/Documents/gerardus/cpp/src/third-party/CGAL-4.2/include/CGAL/Triangulation_hierarchy_3.h:98:
>    instantiated from
> ‘CGAL::Triangulation_hierarchy_3<Tr>::Triangulation_hierarchy_3(InputIterator,
> InputIterator, const typename Tr::Geom_traits&) [with InputIterator =
> __gnu_cxx::__normal_iterator<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int>*, std::vector<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int>, std::allocator<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int> > > >, Tr = CGAL::Delaunay_triangulation_3<CGAL::Epick,
> CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick,
> CGAL::Triangulation_vertex_base_with_info_3<long unsigned int,
> CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick,
> CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick,
> CGAL::Triangulation_vertex_base_with_info_3<long unsigned int,
> CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick,
> CGAL::Triangulation_ds_vertex_base_3<void> > >,
> CGAL::Boolean_tag<false>, CGAL::Boolean_tag<false> >,
> CGAL::Alpha_shape_cell_base_3<CGAL::Epick, CGAL::Default,
> CGAL::Boolean_tag<false>, CGAL::Boolean_tag<false> > > > > >,
> CGAL::Boolean_tag<false>, CGAL::Boolean_tag<false> > >,
> CGAL::Alpha_shape_cell_base_3<CGAL::Epick, CGAL::Default,
> CGAL::Boolean_tag<false>, CGAL::Boolean_tag<false> > >, CGAL::Default>]’
> /home/rcasero/Documents/gerardus/cpp/src/third-party/CGAL-4.2/include/CGAL/internal/Delaunay_triangulation_hierarchy_3.h:55:
>    instantiated from ‘CGAL::Delaunay_triangulation_3<Gt, Tds_,
> CGAL::Location_policy<CGAL::Fast>
>  >::Delaunay_triangulation_3(InputIterator, InputIterator, const Gt&)
> [with InputIterator =
> __gnu_cxx::__normal_iterator<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int>*, std::vector<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int>, std::allocator<std::pair<CGAL::Point_3<CGAL::Epick>, long
> unsigned int> > > >, Gt = CGAL::Epick, Tds_ =
> CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick,
> CGAL::Triangulation_vertex_base_with_info_3<long unsigned int,
> CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick,
> CGAL::Triangulation_ds_vertex_base_3<void> > >,
> CGAL::Boolean_tag<false>, CGAL::Boolean_tag<false> >,
> CGAL::Alpha_shape_cell_base_3<CGAL::Epick, CGAL::Default,
> CGAL::Boolean_tag<false>, CGAL::Boolean_tag<false> > >]’
> /home/rcasero/Documents/gerardus/matlab/CgalToolbox/CgalAlphaShape3.cpp:175:
>    instantiated from here
> /usr/include/c++/4.4/bits/stl_uninitialized.h:74: error: no matching
> function for call to
> ‘CGAL::Point_3<CGAL::Epick>::Point_3(std::pair<CGAL::Point_3<CGAL::Epick>,
> long unsigned int>&)’
> /home/rcasero/Documents/gerardus/cpp/src/third-party/CGAL-4.2/include/CGAL/Point_3.h:83:
> note: candidates are: CGAL::Point_3<R_>::Point_3(const typename R_::RT&,
> const typename R_::RT&, const typename R_::RT&, const typename R_::RT&)
> [with R_ = CGAL::Epick]
> /home/rcasero/Documents/gerardus/cpp/src/third-party/CGAL-4.2/include/CGAL/Point_3.h:75:
> note:                 CGAL::Point_3<R_>::Point_3(const typename
> R_::Kernel_base::Point_3&) [with R_ = CGAL::Epick]
> /home/rcasero/Documents/gerardus/cpp/src/third-party/CGAL-4.2/include/CGAL/Point_3.h:71:
> note:                 CGAL::Point_3<R_>::Point_3(const CGAL::Origin&)
> [with R_ = CGAL::Epick]
> /home/rcasero/Documents/gerardus/cpp/src/third-party/CGAL-4.2/include/CGAL/Point_3.h:69:
> note:                 CGAL::Point_3<R_>::Point_3() [with R_ = CGAL::Epick]
> /home/rcasero/Documents/gerardus/cpp/src/third-party/CGAL-4.2/include/CGAL/Point_3.h:40:
> note:                 CGAL::Point_3<CGAL::Epick>::Point_3(const
> CGAL::Point_3<CGAL::Epick>&)
> make[2]: ***
> [matlab/CgalToolbox/CMakeFiles/cgal_alpha_shape3.dir/CgalAlphaShape3.cpp.o]
> Error 1
> make[1]: *** [matlab/CgalToolbox/CMakeFiles/cgal_alpha_shape3.dir/all]
> Error 2
> make: *** [all] Error 2
> </FULL ERROR>
>
>
>
>
> On 2 August 2013 12:55, Ramón Casero Cañas <[hidden email]
> </user/SendEmail.jtp?type=node&node=4657880&i=0>> wrote:
>
>     Hi Sebastien,
>
>     Many thanks for your indications. I struggled a bit to find the
>     correct way to combine Triangulation_vertex_base_with_info_3
>     and Alpha_shape_3, but it seems to work in the end. Here's a code
>     snippet in case it's useful for somebody else (I'm programming a
>     Matlab MEX C++ function, so some types are specific to Matlab)
>
>     <CODE>
>     /* CGAL headers */
>     #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
>     #include <CGAL/Delaunay_triangulation_3.h>
>     #include <CGAL/Triangulation_vertex_base_with_info_3.h>
>     #include <CGAL/Alpha_shape_3.h>
>
>     typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
>     // vertex
>     typedef CGAL::Triangulation_vertex_base_with_info_3<mwSize, K>  Vb;
>     typedef CGAL::Alpha_shape_vertex_base_3<K, Vb>       AsVb;
>     // cell
>     typedef CGAL::Alpha_shape_cell_base_3<K>             Fb;
>     // triangulation structure: vertex and cell
>     typedef CGAL::Triangulation_data_structure_3<AsVb, Fb> Tds;
>     typedef CGAL::Delaunay_triangulation_3<K, Tds>       Delaunay;
>     typedef CGAL::Alpha_shape_3<Delaunay>                Alpha_shape_3;
>
>     typedef K::Point_3                                   Point;
>     typedef std::pair<Point, mwIndex>                    PointWithIndex;
>     typedef Alpha_shape_3::Alpha_iterator                Alpha_iterator;
>
>     typedef Alpha_shape_3::Facet                         Facet;
>
>        // read points from function
>        std::vector<PointWithIndex> x(nrowsX);
>        for (mwIndex i = 0; i < nrowsX; ++i) {
>          x[i] = std::make_pair(
>        some_function_to_read_a_point(i),
>       i+1 // because this will be a row index in Matlab, 1, ..., Nrows
>       );
>        }
>
>        // Delaunay triangulation
>        //
>     
> http://www.cgal.org/Manual/latest/doc_html/cgal_manual/Triangulation_3/Chapter_main.html#Subsection_39.5.3
>        Delaunay delaunay(x.begin(), x.end());
>        CGAL_assertion(delaunay.number_of_vertices() == nrowsX);
>
>        // compute alpha shape
>        Alpha_shape_3 as(delaunay);
>     </CODE>
>
>     Best regards,
>
>     Ramón.
>
>
>
>     On 2 August 2013 06:29, Sebastien Loriot (GeometryFactory) [via
>     cgal-discuss] <[hidden email]
>     </user/SendEmail.jtp?type=node&node=4657880&i=1>> wrote:
>
>         Use one of the following examples [1] to build the triangulation
>         of your
>         points and set their indices. Then build the alpha shape from the
>         triangulation.
>
>         then you can do:
>
>         typedef CGAL::Delaunay_triangulation_3<.....> DT3;
>
>         ....
>
>         for (std::list<Facet>::iterator it = facets.begin(); it !=
>         facets.end();
>            ++it)
>         {
>             it->first->vertex( DT3::vertex_triple_index(it->second,0)
>         )->info();
>             it->first->vertex( DT3::vertex_triple_index(it->second,1)
>         )->info();
>             it->first->vertex( DT3::vertex_triple_index(it->second,2)
>         )->info();
>         }
>
>         Sebastien.
>
>         [1]
>         
> http://www.cgal.org/Manual/latest/doc_html/cgal_manual/Triangulation_3/Chapter_main.html#Subsection_39.5.3
>
>         On 08/01/2013 08:44 PM, rcasero wrote:
>
>          > Dear all,
>          >
>          > I'm computing alpha shapes on a set of points
>          >
>          >    // read points from function
>          >    std::vector<Point>  x =
>          >      some_function_to_read_points();
>          >
>          >    // compute alpha shape
>          >    Alpha_shape_3 as(x.begin(), x.end());
>          >
>          > based on the example from the "3D Alpha shapes" documentation:
>          >
>          >
>         
> http://www.cgal.org/Manual/latest/doc_html/cgal_manual/Alpha_shapes_3/Chapter_main.html#Subsection_43.5.1
>          >
>
>          > What I would like to do now is extract the triangulation of a
>         given alpha
>          > shape in index form, e.g.
>          >
>          > tri = [3 2 8
>          >          2 1 3
>          >         ...]
>          >
>          > where the indices are the indices of the points in the x
>         vector above.
>          >
>          > This is similar to a previous question I got help from
>         Sebastien Loriot to
>          >
>          >
>         
> http://cgal-discuss.949826.n4.nabble.com/Getting-index-of-facet-closest-to-a-point-tp4657635.html
>          >
>
>          > but I don't manage to find the right pointers in this case so
>         that I can do
>          > something like
>          >
>          >     std::distance(x.begin(), some_iterator_here);
>          >
>          > to get the index.
>          >
>          > I can get a list of the facets
>          >
>          >    // get alpha-shape surface
>          >    std::list<Facet>        facets;
>          >    as.get_alpha_shape_facets(std::back_inserter(facets),
>          >                              Alpha_shape_3::REGULAR);
>          >
>          > and I can iterate said facets
>          >
>          >    for (std::list<Facet>::iterator it = facets.begin(); it !=
>         facets.end();
>          > ++it) {
>          >
>          >    }
>          >
>          >  From the documentation and an old question replied to by
>         sloriot
>          >
>          >
>         
> http://stackoverflow.com/questions/7938311/cgal-help-getting-triangles-coordinates-from-delaunay-triangulation
>          >
>
>          > I know I can the coordinates of the facet vertices as
>          >
>          >    for (std::list<Facet>::iterator it = facets.begin(); it !=
>         facets.end();
>          > ++it) {
>          >
>          >      it->first->vertex((it->second+1)%4)->point();
>          >      it->first->vertex((it->second+2)%4)->point();
>          >      it->first->vertex((it->second+3)%4)->point();
>          >
>          >    }
>          >
>          > but the it->...->point(); above are not pointing to elements
>         in x, neither
>          > it->first->vertex((it->second+1)%4).
>          >
>          > Many thanks for any help or, ahem, pointers.
>          >
>          > Best regards,
>          >
>          > Ramon.
>          >
>          >
>          >
>          > --
>          > View this message in context:
>         
> http://cgal-discuss.949826.n4.nabble.com/Extracting-triangulation-from-alpha-shape-tp4657870.html
>          > Sent from the cgal-discuss mailing list archive at Nabble.com.
>          >
>
>
>         --
>         You are currently subscribed to cgal-discuss.
>         To unsubscribe or access the archives, go to
>         https://sympa.inria.fr/sympa/info/cgal-discuss
>
>
>
>
>         
> ------------------------------------------------------------------------
>         If you reply to this email, your message will be added to the
>         discussion below:
>         
> http://cgal-discuss.949826.n4.nabble.com/Extracting-triangulation-from-alpha-shape-tp4657870p4657872.html
>
>         To unsubscribe from Extracting triangulation from alpha shape,
>         click here.
>         NAML
>         
> <http://cgal-discuss.949826.n4.nabble.com/template/NamlServlet.jtp?macro=macro_viewer&id=instant_html%21nabble%3Aemail.naml&base=nabble.naml.namespaces.BasicNamespace-nabble.view.web.template.NabbleNamespace-nabble.view.web.template.NodeNamespace&breadcrumbs=notify_subscribers%21nabble%3Aemail.naml-instant_emails%21nabble%3Aemail.naml-send_instant_email%21nabble%3Aemail.naml>
>
>
>
>
>
>     --
>     Dr. Ramón Casero Cañas
>
>     Oxford e-Research Centre (OeRC)
>     University of Oxford
>     7 Keble Rd
>     Oxford OX1 3QG
>
>     tlf <a href="tel:%2B44%20%280%29%201865%20610739"
>     value="+441865610739" target="_blank">+44 (0) 1865 610739
>     web http://www.cs.ox.ac.uk/people/Ramon.CaseroCanas
>     photos http://www.flickr.com/photos/rcasero/
>
>
>
>
> --
> Dr. Ramón Casero Cañas
>
> Oxford e-Research Centre (OeRC)
> University of Oxford
> 7 Keble Rd
> Oxford OX1 3QG
>
> tlf     +44 (0) 1865 610739
> web http://www.cs.ox.ac.uk/people/Ramon.CaseroCanas
> photos http://www.flickr.com/photos/rcasero/
>
> ------------------------------------------------------------------------
> View this message in context: Re: Extracting triangulation from alpha
> shape
> <http://cgal-discuss.949826.n4.nabble.com/Extracting-triangulation-from-alpha-shape-tp4657870p4657880.html>
> Sent from the cgal-discuss mailing list archive
> <http://cgal-discuss.949826.n4.nabble.com/> at Nabble.com.

// Copyright (c) 1998, 2001, 2003  INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
// You can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation,
// either version 3 of the License, or (at your option) any later version.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the 
software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL$
// $Id$
//
// Author(s)     : Olivier Devillers 
<>
//                 Sylvain Pion

#ifndef CGAL_TRIANGULATION_HIERARCHY_3_H
#define CGAL_TRIANGULATION_HIERARCHY_3_H

#include <CGAL/basic.h>
#include <CGAL/triangulation_assertions.h>
#include <CGAL/Triangulation_hierarchy_vertex_base_3.h>
#include <CGAL/Location_policy.h>

#include <boost/random/linear_congruential.hpp>
#include <boost/random/geometric_distribution.hpp>
#include <boost/random/variate_generator.hpp>

#ifndef CGAL_TRIANGULATION_3_DONT_INSERT_RANGE_OF_POINTS_WITH_INFO
#include <CGAL/Spatial_sort_traits_adapter_3.h>
#include <CGAL/internal/info_check.h>

#include <boost/tuple/tuple.hpp>
#include <boost/iterator/zip_iterator.hpp>
#include <boost/mpl/and.hpp>
#endif //CGAL_TRIANGULATION_3_DONT_INSERT_RANGE_OF_POINTS_WITH_INFO

namespace CGAL {

// This class is deprecated, but must be kept for backward compatibility.
//
// It would be better to move its content to the Delaunay_triangulation_3
// specializations for Fast_location and make Triangulation_hierarchy_3 the
// empty nutshell instead.
//
// Then, later, maybe merge the Compact/Fast codes in a cleaner factorized 
way.

template < class Tr >
class Triangulation_hierarchy_3
  : public Tr
{
  // parameterization of the hierarchy
  // maximal number of points is 30^5 = 24 millions !
  enum { ratio = 30 };
  enum { minsize = 20};
  enum { maxlevel = 5};

public:
  typedef Tr                                   Tr_Base;
  typedef Fast_location                        Location_policy;
  typedef typename Tr_Base::Geom_traits        Geom_traits;
  typedef typename Tr_Base::Point              Point;
  typedef typename Tr_Base::size_type          size_type;
  typedef typename Tr_Base::Vertex_handle      Vertex_handle;
  typedef typename Tr_Base::Cell_handle        Cell_handle;
  typedef typename Tr_Base::Vertex_iterator    Vertex_iterator;
  typedef typename Tr_Base::Vertex             Vertex;
  typedef typename Tr_Base::Locate_type        Locate_type;
  typedef typename Tr_Base::Finite_vertices_iterator  
Finite_vertices_iterator;
  typedef typename Tr_Base::Finite_cells_iterator     Finite_cells_iterator;
  typedef typename Tr_Base::Finite_facets_iterator    Finite_facets_iterator;
  typedef typename Tr_Base::Finite_edges_iterator     Finite_edges_iterator;

  using Tr_Base::number_of_vertices;
  using Tr_Base::geom_traits;

private:

  // here is the stack of triangulations which form the hierarchy
  Tr_Base*       hierarchy[maxlevel];
  boost::rand48  random;

  void set_up_down(Vertex_handle up, Vertex_handle down)
  {
    up->set_down(down);
    down->set_up(up);
  }

public:

  Triangulation_hierarchy_3(const Geom_traits& traits = Geom_traits());

  Triangulation_hierarchy_3(const Triangulation_hierarchy_3& tr);

  template < typename InputIterator >
  Triangulation_hierarchy_3(InputIterator first, InputIterator last,
                            const Geom_traits& traits = Geom_traits())
    : Tr_Base(traits)
  {
      hierarchy[0] = this;
      for(int i=1; i<maxlevel; ++i)
          hierarchy[i] = new Tr_Base(traits);
      insert(first, last);
  }

  Triangulation_hierarchy_3 & operator=(const Triangulation_hierarchy_3& tr)
  {
    Triangulation_hierarchy_3 tmp(tr);
    swap(tmp);
    return *this;
  }

  ~Triangulation_hierarchy_3();

  void swap(Triangulation_hierarchy_3 &tr);

  void clear();

  // CHECKING
  bool is_valid(bool verbose = false, int level = 0) const;

  // INSERT REMOVE
  Vertex_handle insert(const Point &p, Vertex_handle hint)
  {
    return insert(p, hint == Vertex_handle() ? this->infinite_cell() : 
hint->cell());
  }

  Vertex_handle insert(const Point &p, Cell_handle start = Cell_handle ());

  Vertex_handle insert(const Point &p, Locate_type lt, Cell_handle loc,
                       int li, int lj);

#ifndef CGAL_TRIANGULATION_3_DONT_INSERT_RANGE_OF_POINTS_WITH_INFO
  template < class InputIterator >
  std::ptrdiff_t
  insert( InputIterator first, InputIterator last,
          typename boost::enable_if<
            boost::is_convertible<
                typename std::iterator_traits<InputIterator>::value_type,
                Point
            >
          >::type* = NULL
  )
#else
  template < class InputIterator >
  std::ptrdiff_t
  insert( InputIterator first, InputIterator last)
#endif //CGAL_TRIANGULATION_3_DONT_INSERT_RANGE_OF_POINTS_WITH_INFO
  {
    size_type n = number_of_vertices();

      std::vector<Point> points (first, last);
      spatial_sort (points.begin(), points.end(), geom_traits());

      // hints[i] is the vertex of the previously inserted point in level i.
      // Thanks to spatial sort, they are better hints than what the hierarchy
      // would give us.
      Vertex_handle hints[maxlevel];
      for (typename std::vector<Point>::const_iterator p = points.begin(), 
end = points.end();
              p != end; ++p)
      {
          int vertex_level = random_level();

          Vertex_handle v = hints[0] = hierarchy[0]->insert (*p, hints[0]);
          Vertex_handle prev = v;

          for (int level = 1; level <= vertex_level; ++level) {
              v = hints[level] = hierarchy[level]->insert (*p, hints[level]);
              set_up_down(v, prev);
              prev = v;
          }
      }
      return number_of_vertices() - n;
  }

#ifndef CGAL_TRIANGULATION_3_DONT_INSERT_RANGE_OF_POINTS_WITH_INFO
private:
  //top stands for tuple-or-pair
  template <class Info>
  const Point& top_get_first(const std::pair<Point,Info>& pair) const { 
return pair.first; }
  template <class Info>
  const Info& top_get_second(const std::pair<Point,Info>& pair) const { 
return pair.second; }
  template <class Info>
  const Point& top_get_first(const boost::tuple<Point,Info>& tuple) const { 
return boost::get<0>(tuple); }
  template <class Info>
  const Info& top_get_second(const boost::tuple<Point,Info>& tuple) const { 
return boost::get<1>(tuple); }

  template <class Tuple_or_pair,class InputIterator>
  std::ptrdiff_t insert_with_info(InputIterator first,InputIterator last)
  {
    size_type n = number_of_vertices();
    std::vector<std::ptrdiff_t> indices;
    std::vector<Point> points;
    std::vector<typename Vertex::Info> infos;
    std::ptrdiff_t index=0;
    for (InputIterator it=first;it!=last;++it){
      Tuple_or_pair value=*it;
      points.push_back( top_get_first(value)  );
      infos.push_back ( top_get_second(value) );
      indices.push_back(index++);
    }

    typedef Spatial_sort_traits_adapter_3<Geom_traits,Point*> Search_traits;

    
spatial_sort(indices.begin(),indices.end(),Search_traits(&(points[0]),geom_traits()));


    // hints[i] is the vertex of the previously inserted point in level i.
    // Thanks to spatial sort, they are better hints than what the hierarchy
    // would give us.
    Vertex_handle hints[maxlevel];
    for (typename std::vector<std::ptrdiff_t>::const_iterator
      it = indices.begin(), end = indices.end();
      it != end; ++it)
    {
        int vertex_level = random_level();

        Vertex_handle v = hints[0] = hierarchy[0]->insert (points[*it], 
hints[0]);
        v->info()=infos[*it];
        Vertex_handle prev = v;

        for (int level = 1; level <= vertex_level; ++level) {
            v = hints[level] = hierarchy[level]->insert (points[*it], 
hints[level]);
            set_up_down(v, prev);
            prev = v;
        }
    }
    return number_of_vertices() - n;
  }

public:

  template < class InputIterator >
  std::ptrdiff_t
  insert( InputIterator first,
          InputIterator last,
          typename boost::enable_if<
            boost::is_convertible<
              typename std::iterator_traits<InputIterator>::value_type,
              std::pair<Point,typename internal::Info_check<Vertex>::type>
            > >::type* =NULL
  )
  {
    return insert_with_info< std::pair<Point,typename 
internal::Info_check<Vertex>::type> >(first,last);
  }

  template <class  InputIterator_1,class InputIterator_2>
  std::ptrdiff_t
  insert( boost::zip_iterator< boost::tuple<InputIterator_1,InputIterator_2> 
> first,
          boost::zip_iterator< boost::tuple<InputIterator_1,InputIterator_2> 
> last,
          typename boost::enable_if<
            boost::mpl::and_<
              boost::is_convertible< typename 
std::iterator_traits<InputIterator_1>::value_type, Point >,
              boost::is_convertible< typename 
std::iterator_traits<InputIterator_2>::value_type, typename 
internal::Info_check<Vertex>::type >
            >
          >::type* =NULL
  )
  {
    return insert_with_info< boost::tuple<Point,typename 
internal::Info_check<Vertex>::type> >(first,last);
  }
#endif //CGAL_TRIANGULATION_3_DONT_INSERT_RANGE_OF_POINTS_WITH_INFO

  void remove(Vertex_handle v);

  template < typename InputIterator >
  size_type remove(InputIterator first, InputIterator beyond)
  {
    size_type n = number_of_vertices();
    while (first != beyond) {
      remove(*first);
      ++first;
    }
    return n - number_of_vertices();
  }

  template < typename InputIterator >
  size_type remove_cluster(InputIterator first, InputIterator beyond)
  {
    CGAL_triangulation_precondition(!this->does_repeat_in_range(first, 
beyond));
    CGAL_triangulation_precondition(!this->infinite_vertex_in_range(first, 
beyond));
    size_type n = this->number_of_vertices();
    std::vector<Vertex_handle> vo(first, beyond), vc;
    int l=0;
    while(1) {
      size_type n = vo.size();
      if(n == 0) break;
      for(size_type i=0; i<n; i++) {
        if(vo[i]->up() != Vertex_handle()) vc.push_back(vo[i]->up());
      }
      hierarchy[l++]->remove_cluster(vo.begin(), vo.end());
      std::swap(vo,vc);
      vc.clear();
    }
    return n - this->number_of_vertices();
  }

#ifndef CGAL_NO_DEPRECATED_CODE
  CGAL_DEPRECATED Vertex_handle move_point(Vertex_handle v, const Point & p);
#endif

  Vertex_handle move_if_no_collision(Vertex_handle v, const Point &p);
  Vertex_handle move(Vertex_handle v, const Point &p);

public: // some internal methods

  // INSERT REMOVE DISPLACEMENT
  // GIVING NEW FACES

  template <class OutputItCells>
  Vertex_handle insert_and_give_new_cells(const Point  &p, 
                                          OutputItCells fit,
                                          Cell_handle start = Cell_handle() );
                
  template <class OutputItCells>
  Vertex_handle insert_and_give_new_cells(const Point& p,
                                          OutputItCells /* fit */,
                                          Vertex_handle hint)
  {
    return insert_and_give_new_cells(p, hint == Vertex_handle() ? 
                                     this->infinite_cell() : hint->cell());   
                  
  }

  template <class OutputItCells>
  Vertex_handle insert_and_give_new_cells(const Point& p,
                                          Locate_type lt,
                                          Cell_handle c, int li, int lj, 
                                          OutputItCells fit);

  template <class OutputItCells>
  void remove_and_give_new_cells(Vertex_handle v, 
                                 OutputItCells fit);

  template <class OutputItCells>
  Vertex_handle move_if_no_collision_and_give_new_cells(Vertex_handle v, 
                                                        const Point &p, 
OutputItCells fit);
        
public: 


  //LOCATE
  Cell_handle locate(const Point& p, Locate_type& lt, int& li, int& lj,
                     Vertex_handle hint) const
  {
    return locate(p, lt, li, lj, hint == Vertex_handle() ? 
this->infinite_cell() : hint->cell());
  }

  Cell_handle locate(const Point& p, Vertex_handle hint) const
  {
    return locate(p, hint == Vertex_handle() ? this->infinite_cell() : 
hint->cell());
  }

  Cell_handle locate(const Point& p, Locate_type& lt, int& li, int& lj,
                     Cell_handle start = Cell_handle ()) const;

  Cell_handle locate(const Point& p, Cell_handle start = Cell_handle ()) 
const;

  Vertex_handle
  nearest_vertex(const Point& p, Cell_handle start = Cell_handle()) const;

protected:

  struct locs {
      Cell_handle pos;
      int li, lj;
      Locate_type lt;
  };

  void locate(const Point& p, Locate_type& lt, int& li, int& lj,
              locs pos[maxlevel], Cell_handle start = Cell_handle ()) const;

  int random_level();
};


template <class Tr >
Triangulation_hierarchy_3<Tr>::
Triangulation_hierarchy_3(const Geom_traits& traits)
  : Tr_Base(traits)
{
  hierarchy[0] = this;
  for(int i=1;i<maxlevel;++i)
    hierarchy[i] = new Tr_Base(traits);
}

// copy constructor duplicates vertices and cells
template <class Tr>
Triangulation_hierarchy_3<Tr>::
Triangulation_hierarchy_3(const Triangulation_hierarchy_3<Tr> &tr)
    : Tr_Base(tr)
{
  hierarchy[0] = this;
  for(int i=1; i<maxlevel; ++i)
    hierarchy[i] = new Tr_Base(*tr.hierarchy[i]);

  // up and down have been copied in straightforward way
  // compute a map at lower level

  std::map< Vertex_handle, Vertex_handle > V;

  for( Finite_vertices_iterator it = hierarchy[0]->finite_vertices_begin(),
       end = hierarchy[0]->finite_vertices_end(); it != end; ++it)
    if (it->up() != Vertex_handle())
      V[ it->up()->down() ] = it;

  for(int j=1; j<maxlevel; ++j) {
    for( Finite_vertices_iterator it = hierarchy[j]->finite_vertices_begin(),
         end = hierarchy[j]->finite_vertices_end(); it != end; ++it) {
        // current it->down() pointer goes in original instead in copied 
triangulation
        set_up_down(it, V[it->down()]);
        // make map for next level
        if (it->up() != Vertex_handle())
            V[ it->up()->down() ] = it;
    }
  }
}

template <class Tr>
void
Triangulation_hierarchy_3<Tr>::
swap(Triangulation_hierarchy_3<Tr> &tr)
{
  Tr_Base::swap(tr);
  for(int i=1; i<maxlevel; ++i)
      std::swap(hierarchy[i], tr.hierarchy[i]);
}

template <class Tr>
Triangulation_hierarchy_3<Tr>::
~Triangulation_hierarchy_3()
{
  clear();
  for(int i=1; i<maxlevel; ++i) {
    delete hierarchy[i];
  }
}

template <class Tr>
void
Triangulation_hierarchy_3<Tr>::
clear()
{
  for(int i=0;i<maxlevel;++i)
    hierarchy[i]->clear();
}

template <class Tr>
bool
Triangulation_hierarchy_3<Tr>::
is_valid(bool verbose, int level) const
{
  bool result = true;

  // verify correctness of triangulation at all levels
  for(int i=0; i<maxlevel; ++i)
        result = result && hierarchy[i]->is_valid(verbose, level);

  // verify that lower level has no down pointers
  for( Finite_vertices_iterator it = hierarchy[0]->finite_vertices_begin(),
       end = hierarchy[0]->finite_vertices_end(); it != end; ++it)
    result = result && (it->down() == Vertex_handle());

  // verify that other levels has down pointer and reciprocal link is fine
  for(int j=1; j<maxlevel; ++j)
    for( Finite_vertices_iterator it = hierarchy[j]->finite_vertices_begin(),
         end = hierarchy[j]->finite_vertices_end(); it != end; ++it)
      result = result && &*(it) == &*(it->down()->up());

  // verify that other levels has down pointer and reciprocal link is fine
  for(int k=0; k<maxlevel-1; ++k)
    for( Finite_vertices_iterator it = hierarchy[k]->finite_vertices_begin(),
         end = hierarchy[k]->finite_vertices_end(); it != end; ++it)
      result = result && ( it->up() == Vertex_handle() ||
                &*it == &*(it->up())->down() );

  return result;
}

template <class Tr>
typename Triangulation_hierarchy_3<Tr>::Vertex_handle
Triangulation_hierarchy_3<Tr>::
insert(const Point &p, Cell_handle start)
{
  int vertex_level = random_level();
  Locate_type lt;
  int i, j;
  // locate using hierarchy
  locs positions[maxlevel];
  locate(p, lt, i, j, positions, start);
  // insert at level 0
  Vertex_handle vertex = hierarchy[0]->insert(p,
                                              positions[0].lt,
                                              positions[0].pos,
                                              positions[0].li,
                                              positions[0].lj);
  Vertex_handle previous = vertex;
  Vertex_handle first = vertex;

  int level = 1;
  while (level <= vertex_level ){
      if (positions[level].pos == Cell_handle())
          vertex = hierarchy[level]->insert(p);
      else
          vertex = hierarchy[level]->insert(p,
                                            positions[level].lt,
                                            positions[level].pos,
                                            positions[level].li,
                                            positions[level].lj);
    set_up_down(vertex, previous);
    previous=vertex;
    level++;
  }
  return first;
}

template <class Tr>
template <class OutputItCells>
typename Triangulation_hierarchy_3<Tr>::Vertex_handle
Triangulation_hierarchy_3<Tr>::
insert_and_give_new_cells(const Point &p, OutputItCells fit, Cell_handle 
start)
{
  int vertex_level = random_level();
  Locate_type lt;
  int i, j;
  // locate using hierarchy
  locs positions[maxlevel];
  locate(p, lt, i, j, positions, start);
  // insert at level 0
  Vertex_handle vertex = hierarchy[0]->insert_and_give_new_cells(p,
                                                                 
positions[0].lt,
                                                                 
positions[0].pos,
                                                                 
positions[0].li,
                                                                 
positions[0].lj,fit);
  Vertex_handle previous = vertex;
  Vertex_handle first = vertex;

  int level = 1;
  while (level <= vertex_level ){
    if (positions[level].pos == Cell_handle())
      vertex = hierarchy[level]->insert(p);
    else
      vertex = hierarchy[level]->insert(p,
                                        positions[level].lt,
                                        positions[level].pos,
                                        positions[level].li,
                                        positions[level].lj);
    set_up_down(vertex, previous);
    previous=vertex;
    level++;
  }
  return first;
}

template <class Tr>
typename Triangulation_hierarchy_3<Tr>::Vertex_handle
Triangulation_hierarchy_3<Tr>::
insert(const Point &p, Locate_type lt, Cell_handle loc, int li, int lj)
{
  int vertex_level = random_level();
  // insert at level 0
  Vertex_handle vertex = hierarchy[0]->insert(p,lt,loc,li,lj);
  Vertex_handle previous = vertex;
  Vertex_handle first = vertex;

  if (vertex_level > 0) {
    Locate_type lt;
    int i, j;
    // locate using hierarchy
    locs positions[maxlevel];
    locate(p, lt, i, j, positions, vertex->cell());

    int level = 1;
    while (level <= vertex_level ){
      if (positions[level].pos == Cell_handle())
        vertex = hierarchy[level]->insert(p);
      else
        vertex = hierarchy[level]->insert(p,
            positions[level].lt,
            positions[level].pos,
            positions[level].li,
            positions[level].lj);
      set_up_down(vertex, previous);
      previous=vertex;
      level++;
    }
  }
  return first;
}

template <class Tr>
template <class OutputItCells>
typename Triangulation_hierarchy_3<Tr>::Vertex_handle
Triangulation_hierarchy_3<Tr>::
insert_and_give_new_cells(const Point &p, Locate_type lt, Cell_handle loc, 
  int li, int lj, OutputItCells fit)
{
  int vertex_level = random_level();
  // insert at level 0
  Vertex_handle vertex = 
    hierarchy[0]->insert_and_give_new_cells(p,lt,loc,li,lj,fit);
  Vertex_handle previous = vertex;
  Vertex_handle first = vertex;

  if (vertex_level > 0) {
    Locate_type lt;
    int i, j;
    // locate using hierarchy
    locs positions[maxlevel];
    locate(p, lt, i, j, positions, vertex->cell());

    int level = 1;
    while (level <= vertex_level ){
      if (positions[level].pos == Cell_handle())
        vertex = hierarchy[level]->insert(p);
      else
        vertex = hierarchy[level]->insert(p,
                                          positions[level].lt,
                                          positions[level].pos,
                                          positions[level].li,
                                          positions[level].lj);
      set_up_down(vertex, previous);
      previous=vertex;
      level++;
    }
  }
  return first;
}

template <class Tr>
void
Triangulation_hierarchy_3<Tr>::
remove(Vertex_handle v)
{
  CGAL_triangulation_precondition(v != Vertex_handle());
  for (int l = 0; l < maxlevel; ++l) {
    Vertex_handle u = v->up();
    hierarchy[l]->remove(v);
    if (u == Vertex_handle())
        break;
    v = u;
  }
}

template <class Tr>
template <class OutputItCells>
void
Triangulation_hierarchy_3<Tr>::
remove_and_give_new_cells(Vertex_handle v, OutputItCells fit)
{
  CGAL_triangulation_precondition(v != Vertex_handle());
  CGAL_triangulation_precondition(!is_infinite(v));
  for (int l = 0; l < maxlevel; ++l) {
    Vertex_handle u = v->up();
    if(l) hierarchy[l]->remove(v);
    else hierarchy[l]->remove_and_give_new_cells(v, fit);
    if (u == Vertex_handle())
        break;
    v = u;
  }
}

#ifndef CGAL_NO_DEPRECATED_CODE
template < class Tr >
typename Triangulation_hierarchy_3<Tr>::Vertex_handle
Triangulation_hierarchy_3<Tr>::
move_point(Vertex_handle v, const Point & p)
{
  CGAL_triangulation_precondition(v != Vertex_handle());
  Vertex_handle old, ret;

  for (std::size_t l = 0; l < maxlevel; ++l) {
    Vertex_handle u = v->up();
    Vertex_handle w = hierarchy[l]->move_point(v, p);
    if (l == 0) {
        ret = w;
    }
    else {
        set_up_down(w, old);
    }
    if (u == Vertex_handle())
        break;
    old = w;
    v = u;
  }

  return ret;
}
#endif

template <class Tr>
typename Triangulation_hierarchy_3<Tr>::Vertex_handle
Triangulation_hierarchy_3<Tr>::
move_if_no_collision(Vertex_handle v, const Point & p)
{
  CGAL_triangulation_precondition(!this->is_infinite(v));       
  if(v->point() == p) return v;
  Vertex_handle ans;
  for (int l = 0; l < maxlevel; ++l) {
    Vertex_handle u = v->up();
    if(l) hierarchy[l]->move_if_no_collision(v, p);
    else ans = hierarchy[l]->move_if_no_collision(v, p);
    if(ans != v) return ans;
    if (u == Vertex_handle())
      break;
    v = u;
  }
  return ans;
}

template <class Tr>
typename Triangulation_hierarchy_3<Tr>::Vertex_handle
Triangulation_hierarchy_3<Tr>::
move(Vertex_handle v, const Point & p)
{
  CGAL_triangulation_precondition(!this->is_infinite(v));
  if(v->point() == p) return v;
  Vertex_handle w = move_if_no_collision(v,p);
  if(w != v) {
    remove(v);
    return w;
  }
  return v;
}

template <class Tr>
template <class OutputItCells>
typename Triangulation_hierarchy_3<Tr>::Vertex_handle
Triangulation_hierarchy_3<Tr>::
move_if_no_collision_and_give_new_cells(
  Vertex_handle v, const Point & p, OutputItCells fit)
{
  CGAL_triangulation_precondition(!is_infinite(v));     
  if(v->point() == p) return v;
  Vertex_handle ans;
  for (int l = 0; l < maxlevel; ++l) {
    Vertex_handle u = v->up();
    if(l) hierarchy[l]->move_if_no_collision(v, p);
    else ans = 
           hierarchy[l]->move_if_no_collision_and_give_new_cells(v, p, fit);
    if(ans != v) return ans;
    if (u == Vertex_handle())
      break;
    v = u;
  }
  return ans;
}

template <class Tr>
inline
typename Triangulation_hierarchy_3<Tr>::Cell_handle
Triangulation_hierarchy_3<Tr>::
locate(const Point& p, Locate_type& lt, int& li, int& lj, Cell_handle start) 
const
{
  if (start != Cell_handle ())
    return Tr_Base::locate (p, lt, li, lj, start);
  locs positions[maxlevel];
  locate(p, lt, li, lj, positions);
  return positions[0].pos;
}

template <class Tr>
inline
typename Triangulation_hierarchy_3<Tr>::Cell_handle
Triangulation_hierarchy_3<Tr>::
locate(const Point& p, Cell_handle start) const
{
  if (start != Cell_handle ())
    return Tr_Base::locate (p, start);
  Locate_type lt;
  int li, lj;
  return locate(p, lt, li, lj);
}

template <class Tr>
void
Triangulation_hierarchy_3<Tr>::
locate(const Point& p, Locate_type& lt, int& li, int& lj,
       locs pos[maxlevel], Cell_handle start) const
{
  int level = maxlevel;

  // find the highest level with enough vertices
  while (hierarchy[--level]->number_of_vertices() < (size_type) minsize) {
    if ( ! level)
        break;  // do not go below 0
  }

  for (int i=level+1; i<maxlevel; ++i)
      pos[i].pos = Cell_handle();

  Cell_handle position = Cell_handle();
  while(level > 0) {
    // locate at that level from "position"
    // result is stored in "position" for the next level
    pos[level].pos = position = hierarchy[level]->locate(p,
                                                         pos[level].lt,
                                                         pos[level].li,
                                                         pos[level].lj,
                                                         position);
    // find the nearest vertex.
    Vertex_handle nearest = hierarchy[level]->nearest_vertex_in_cell(p, 
position);

    // go at the same vertex on level below
    nearest = nearest->down();
    position = nearest->cell();                // incident cell
    --level;
  }

  if (start != Cell_handle())
    position = start;

  pos[0].pos = hierarchy[0]->locate(p, lt, li, lj, position); // at level 0
  pos[0].lt = lt;
  pos[0].li = li;
  pos[0].lj = lj;
}

template <class Tr>
typename Triangulation_hierarchy_3<Tr>::Vertex_handle
Triangulation_hierarchy_3<Tr>::
nearest_vertex(const Point& p, Cell_handle start) const
{
    return Tr_Base::nearest_vertex(p, start != Cell_handle() ? start : 
locate(p));
}

template <class Tr>
int
Triangulation_hierarchy_3<Tr>::
random_level()
{
  boost::geometric_distribution<> proba(1.0/ratio);
  boost::variate_generator<boost::rand48&, boost::geometric_distribution<> > 
die(random, proba);

  return (std::min)(die(), (int)maxlevel)-1;
}

} //namespace CGAL

#endif // CGAL_TRIANGULATION_HIERARCHY_3_H