envelope_divide_and_conquer_3.h

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

// Copyright (c) 2005  Tel-Aviv University (Israel).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you may redistribute it under
// the terms of the Q Public License version 1.0.
// See the file LICENSE.QPL distributed with CGAL.
//
// 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: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.3-branch/Envelope_3/include/CGAL/Envelope_3/Envelope_divide_and_conquer_3.h $
// $Id: Envelope_divide_and_conquer_3.h 37935 2007-04-04 17:12:45Z efif $
//
// Author(s)     : Michal Meyerovitch     <[e-mail]>
//                 Baruch Zukerman        <[e-mail]>
//                 Efi Fogel              <[e-mail]>

#ifndef CGAL_ENVELOPE_DIVIDE_AND_CONQUER_3_H
#define CGAL_ENVELOPE_DIVIDE_AND_CONQUER_3_H

#define CGAL_ENVELOPE_SAVE_COMPARISONS
#define CGAL_ENVELOPE_USE_BFS_FACE_ORDER

#include <iostream>
#include <cassert>
#include <list>
#include <set>
#include <vector>
#include <map>
#include <time.h>

#include <CGAL/Object.h>
#include <CGAL/enum.h>
#include <CGAL/Arr_observer.h>
#include <CGAL/Envelope_3/Envelope_base.h>
#include <CGAL/Envelope_3/Env_overlay_2.h>
#include <CGAL/Envelope_3/Envelope_element_visitor_3.h>
#include <CGAL/Envelope_3/set_dividors.h>

#ifdef CGAL_ENVELOPE_USE_BFS_FACE_ORDER
#include <CGAL/Arr_face_map.h>
#include <CGAL/graph_traits_Dual_Arrangement_2.h>
#include <boost/graph/dijkstra_shortest_paths.hpp>
#endif

// this base divide & conquer algorithm splits the input into 2 groups,
// calculates the result over the 2 groups, and then merges the results like
// this:
// - overlays the maps, and set at each vertex, edge and face the data from
//   the 2 maps
// - foreach vertex, decide from the 2 envelopes data, which is the envelope
//   of all surfaces
// - foreach edge, do the same, using the Resolver class. an edge can split to
//   a constant number of parts, each with its own envelope data.
// - foreach face, do the same as for edges, using the Resolver class.
//   a face can split to a number of sub-faces that are linear with the face's
//   size, each with its own envelope data.
// - remove edges between faces with the same envelope data, which do not
//   contribute to the shape of the envelope (i.e. have the same envelope data
//   as their adjacent faces)
// - remove unneccessary vertices of two kinds:
//   a. vertices which have degree 2, the 2 incident edges can be geometrically
//      merged, and has the same envelope data as both these edges
//   b. isolated vertices which have the same envelope data as their incident
//      face

// the algorithm deals with some degenerate input including:
// 1. more than one surface on faces, edges, vertices
//    (the minimization diagram should also support this)
// 2. all degenerate cases in 2d (the minimization diagram model is
//    responsible for)

// some degenerate cases in the responsibility of the geometric traits
// 1. overlapping surfaces
// 2. a vertical surface that contributes only edge (or edges) to the envelope

CGAL_BEGIN_NAMESPACE

// The algorithm has 5 template parameters:
// 1. EnvelopeTraits_3        - the geometric traits class
// 2. MinimizationDiagram_2   - the type of the output, which is an arrangement
//                              with additional information (list of surfaces)
//                              in vertices, edges & faces
// 3. EnvelopeResolver_3      - part of the algorithm that solves the shape of
//                              the envelope between 2 surfaces over a feature
//                              of the arrangement
// 4. Overlay_2               - overlay of 2 MinimizationDiagram_2
template <class EnvelopeTraits_3, 
          class MinimizationDiagram_2,
          class EnvelopeResolver_3 =
            Envelope_element_visitor_3<EnvelopeTraits_3, MinimizationDiagram_2>,
          class Overlay_2 = Envelope_overlay_2<MinimizationDiagram_2> >
class Envelope_divide_and_conquer_3
{
public:
  typedef EnvelopeTraits_3                              Traits;
  typedef typename Traits::Surface_3                    Surface_3;
  typedef typename Traits::Xy_monotone_surface_3        Xy_monotone_surface_3;

  typedef MinimizationDiagram_2                         Minimization_diagram_2;

  typedef typename Traits::Point_2                      Point_2;
  typedef typename Traits::X_monotone_curve_2           X_monotone_curve_2;
  typedef typename Traits::Curve_2                      Curve_2;

  typedef EnvelopeResolver_3                            Envelope_resolver;

  typedef Envelope_divide_and_conquer_3<Traits, Minimization_diagram_2,
                                        EnvelopeResolver_3,
                                        Overlay_2>      Self;

protected:

  typedef typename Minimization_diagram_2::Halfedge_const_iterator
    Halfedge_const_iterator;
  typedef typename Minimization_diagram_2::Halfedge_handle
    Halfedge_handle;
  typedef typename Minimization_diagram_2::Halfedge_iterator
    Halfedge_iterator;
  typedef typename Minimization_diagram_2::Face_handle
    Face_handle;
  typedef typename Minimization_diagram_2::Edge_iterator
    Edge_iterator;
  typedef typename Minimization_diagram_2::Face_iterator
    Face_iterator;
  typedef typename Minimization_diagram_2::Vertex_handle
    Vertex_handle;
  typedef typename Minimization_diagram_2::Vertex_iterator
    Vertex_iterator;
  typedef typename Minimization_diagram_2::Hole_iterator
    Hole_iterator;
  typedef typename Minimization_diagram_2::Ccb_halfedge_circulator
    Ccb_halfedge_circulator;
  typedef typename Minimization_diagram_2::Halfedge_around_vertex_circulator

    Halfedge_around_vertex_circulator;

  typedef Arr_observer<Minimization_diagram_2>          Md_observer;
  typedef typename Minimization_diagram_2::Dcel::Dcel_data_iterator
    Envelope_data_iterator;

#ifdef CGAL_ENVELOPE_USE_BFS_FACE_ORDER
  typedef CGAL::Dual<Minimization_diagram_2>
    Dual_Minimization_diagram_2;
#endif

public:
  // c'tor
  Envelope_divide_and_conquer_3(Envelope_type type = LOWER)
  {    
    // Allocate the traits.
    traits = new Traits;                                                     

    own_traits = true;

    // Allocate the Envelope resolver with our traits
    resolver = new Envelope_resolver(traits, type);

    m_is_lower = ((type == LOWER) ? true : false);
  }

  Envelope_divide_and_conquer_3(Traits* tr, Envelope_type type = LOWER)
  {
    // Set the traits.
    traits = tr;
    own_traits = false;

    // Allocate the Envelope resolver with our traits
    resolver = new Envelope_resolver(traits, type);

    m_is_lower = ((type == LOWER) ? true : false);
  }
  
  // virtual destructor.
  virtual ~Envelope_divide_and_conquer_3()
  {
    // Free the traits object, if necessary.
    if (own_traits)
      delete traits;

    // Free the resolver
    delete resolver;
  }

  // compute the envelope of surfaces in 3D, using the default arbitrary
  // dividor
  template <class SurfaceIterator>
  void construct_lu_envelope(SurfaceIterator begin, SurfaceIterator end,
                             Minimization_diagram_2 &result)
  {
    Arbitrary_dividor dividor;
    construct_lu_envelope(begin, end, result, dividor);
  }
    

  // compute the envelope of surfaces in 3D using the given set dividor
  template <class SurfaceIterator, class SetDividor>
  void construct_lu_envelope(SurfaceIterator begin, SurfaceIterator end, 
                             Minimization_diagram_2 &result,
                             SetDividor& dividor)
  {
    if (begin == end)
    {
      return; // result is empty
    }
    
    // make the general surfaces xy-monotone
    std::list<Xy_monotone_surface_3> xy_monotones;
    for (; begin != end; ++begin)
      traits->make_xy_monotone_3_object()(*begin, m_is_lower,
                                          std::back_inserter(xy_monotones));

    // recursively construct the envelope of the xy-monotone parts
    construct_lu_envelope_xy_monotones(xy_monotones.begin(), 
                                       xy_monotones.end(), result, dividor);

    CGAL_assertion(is_envelope_valid(result));     
  }

  // compute the envelope of xy-monotone surfaces in 3D,
  // using the default arbitrary dividor
  template <class SurfaceIterator>
  void construct_envelope_xy_monotone(SurfaceIterator begin,
                                      SurfaceIterator end,
                                      Minimization_diagram_2 &result)
  {
    Arbitrary_dividor dividor;
    construct_envelope_xy_monotone(begin, end, result, dividor);
  }

  // compute the envelope of xy-monotone surfaces in 3D using the given
  // set dividor
  template <class SurfaceIterator, class SetDividor>
  void construct_envelope_xy_monotone(SurfaceIterator begin,
                                      SurfaceIterator end,
                                      Minimization_diagram_2 &result,
                                      SetDividor& dividor)
  {
    if (begin == end)
      return; // result is empty

    // recursively construct the envelope of the xy-monotone parts
    construct_lu_envelope_xy_monotones(begin, end, result, dividor);
    CGAL_assertion(is_envelope_valid(result));
  }

  /*! Access the traits object (const version). */
  const Traits* get_traits() const
  {
    return traits;
  }

  /*! Access the traits object (non-const version). */
  Traits* get_traits() 
  {
    return traits;
  }

  void reset()
  {
    resolver->reset();
  }

protected:

  // compute the envelope of xy-monotone surfaces in 3D 
  template <class SurfaceIterator, class SetDividor>
  void construct_lu_envelope_xy_monotones(SurfaceIterator begin,
                                          SurfaceIterator end,
                                          Minimization_diagram_2 &result,
                                          SetDividor& dividor)
  {
    if (begin == end)
    {
      return; // result is empty
    }

    SurfaceIterator first = begin++;

    if (begin == end)
    {
      // only one surface is in the collection. insert it the result
      Xy_monotone_surface_3& surf = *first;
      
      deal_with_one_surface(surf, result);
      return;      
    }
   
    // divide the surfaces into 2 groups (insert surface to each group
    // alternately)
    std::list<Xy_monotone_surface_3> group1, group2;
    dividor(first, end,
            std::back_inserter(group1), std::back_inserter(group2));
    
    // recursively calculate the LU_envelope of the 2 groups
    Minimization_diagram_2 result1(traits), result2(traits);
    construct_lu_envelope_xy_monotones(group1.begin(), group1.end(),
                                       result1, dividor);
    construct_lu_envelope_xy_monotones(group2.begin(), group2.end(),
                                       result2, dividor);
        
    // merge the results:
    merge_envelopes(result1, result2, result);

    result1.clear();
    result2.clear();
    
    CGAL_assertion(is_envelope_valid(result));   
  }

  void deal_with_one_surface(Xy_monotone_surface_3& surf,
                             Minimization_diagram_2& result)
  {
    typedef std::list<Object>                            Boundary_list;
    typedef std::pair<X_monotone_curve_2, Oriented_side> Boundary_xcurve;
    typedef Boundary_list::iterator                      Boundary_iterator;

    Boundary_list     boundary;
    traits->
      construct_projected_boundary_2_object()(surf,
                                              std::back_inserter(boundary));
    
    if(boundary.empty())
    {
      //one infinite surface
      result.unbounded_face()->set_data(surf);
      return;
    }

    for (Boundary_iterator boundary_it = boundary.begin();
        boundary_it != boundary.end();
        ++boundary_it)
    {
      const Object& obj = *boundary_it;
      Boundary_xcurve boundary_cv;
      if(assign(boundary_cv, obj))
      {
        Oriented_side side = boundary_cv.second;
        Halfedge_handle he =
          insert_non_intersecting_curve(result, boundary_cv.first);
        
        if(side == ON_ORIENTED_BOUNDARY)
        {
          // vertical xy-surface
          he->face()->set_no_data();
          he->twin()->face()->set_no_data();
          
          continue;
        }

        if(he->face() != he->twin()->face())