delaunay_triangulation_base_3.h

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

// Copyright (c) 2005  Stanford University (USA).
// 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 Lesser General Public License as
// published by the Free Software Foundation; version 2.1 of the License.
// See the file LICENSE.LGPL 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/Kinetic_data_structures/include/CGAL/Kinetic/internal/Delaunay_triangulation_base_3.h $
// $Id: Delaunay_triangulation_base_3.h 37995 2007-04-07 19:07:06Z drussel $ $Date: 2007-04-07 21:07:06 +0200 (Sat, 07 Apr 2007) $
// 
//
// Author(s)     : Daniel Russel <[e-mail]>

#ifndef CGAL_KINETIC_INTERNAL_DELAUNAY_BASE_3_H
#define CGAL_KINETIC_INTERNAL_DELAUNAY_BASE_3_H

#include <CGAL/Kinetic/basic.h>

#include <CGAL/Kinetic/internal/Triangulation_helper_3.h>
#include <CGAL/Kinetic/internal/triangulation_helpers_3.h>
#include <CGAL/Kinetic/Event_base.h>
//#include <CGAL/Kinetic/internal/Delaunay_cache_3.h>

// STL
#include <map>
#include <set>
#include <iterator>
#include <ostream>
#include <iostream>

//extern int too_late__;
//extern int filtered__;

CGAL_KINETIC_BEGIN_INTERNAL_NAMESPACE

template <class KD, class RS>
class Delaunay_event_base_3: public Event_base<KD*>
{
  typedef Event_base<KD*> P;
public:
  Delaunay_event_base_3(const RS &s,
			KD *kdel):P(kdel), s_(s){
  }
  //! Default constructor
  /*!  This really should not be used, but I need it for the
    Simulator::event() call, due to the apparent gcc compiler bug.
  */
  Delaunay_event_base_3() {

  }
  void process() {
    // for some reason VC insists that this be there
    CGAL_assertion(0 && "Process called in Delaunay_event_base_3");
  }
  const RS& root_stack() const
  {
    return s_;
  }
  std::ostream& write(std::ostream &out) const
  {
    out << "Delaunay_event";
    return out;
  }
  KD* kdel() const {
    return P::kds();
  }
protected:
  const RS s_;
};

/*template <class K, class R>
std::ostream& operator<<(std::ostream &out, const Delaunay_event_base_3<K, R> &e)
{
  e.write(out);
  return out;
  }*/


template <class KD, class RS>
class Delaunay_3_edge_flip_event: public Delaunay_event_base_3<KD, RS>
{
public:
  typedef Delaunay_event_base_3<KD, RS>   P;
  Delaunay_3_edge_flip_event(const RS &s,
			     const typename KD::Edge &e,
			     KD *kdel):P(s, kdel), e_(e) {
#ifndef NDEBUG
    o_= edge_point(e_,0);
    d_= edge_point(e_,1);
#endif
  }
  void process() {
    P::kdel()->flip(e_);
  }

  static typename KD::Point_key edge_point(const typename KD::Edge &e, int i) {
    return vertex_of_edge(e, i)->point();
  }

  std::ostream& write(std::ostream &out) const
  {
    out << "Flip ";
    internal::write_edge(e_, out);
#if 0
    out << "(" << o_ << d_<<")" << std::flush;
    CGAL_postcondition(o_== vertex_of_edge(e_,static_cast<int>(0))->point());
    CGAL_postcondition(d_== vertex_of_edge(e_,1)->point());
#endif
    return out;
  }
 

  void audit(typename KD::Event_key k) const {
    if (e_.first->edge_label(e_.second,e_.third) != k) {
      CGAL_KINETIC_ERROR("Mismatch, for label " << k << " had event " << e_.first->edge_label(e_.second,e_.third));
    }
    CGAL_assertion(e_.first->edge_label(e_.second,e_.third) == k);
  }
protected:
  const typename KD::Edge e_;
#ifndef NDEBUG
  typename KD::Point_key o_, d_;
#endif
};

/*template <class B, class C, class D>
std::ostream& operator<<(std::ostream &out, const Delaunay_3_edge_flip_event<B, C, D> &e)
{
  e.write(out);
  return out;
  }*/


template <class KD, class RS>
class Delaunay_3_facet_flip_event:  public Delaunay_event_base_3<KD, RS>
{
public:
  typedef Delaunay_event_base_3<KD, RS>   P;
  Delaunay_3_facet_flip_event(const RS &s,
			      const typename KD::Facet &e,
			      KD *kdel):  P(s, kdel),e_(e) {
#ifndef NDEBUG
    a_= vertex_of_facet(e_,0)->point();
    b_= vertex_of_facet(e_,1)->point();
    c_= vertex_of_facet(e_,2)->point();
#endif
  }
  void process() {
    P::kdel()->flip(e_);
  }
  std::ostream& write(std::ostream &out) const
  {
    out << "Flip ";
    write_facet(e_, out);
#if 0
    out << "(" << a_ << b_<<c_ << ")";
    CGAL_postcondition(a_== vertex_of_facet(e_,0)->point());
    CGAL_postcondition(b_== vertex_of_facet(e_,1)->point());
    CGAL_postcondition(c_== vertex_of_facet(e_,2)->point());
#endif
    return out;
  }
 
  void audit(typename KD::Event_key k) const {
    if (e_.first->facet_label(e_.second) != k) {
      CGAL_KINETIC_ERROR("Mismatch, for label " << k << " had event " << e_.first->facet_label(e_.second));
    }
    CGAL_assertion(e_.first->facet_label(e_.second) == k);
  }
protected:
  const typename KD::Facet e_;
#ifndef NDEBUG
  typename KD::Point_key a_, b_, c_;
#endif
};

/*template <class T, class C, class D>
std::ostream& operator<<(std::ostream &out, const Delaunay_3_facet_flip_event<T, C, D> &e)
{
  e.write(out);
  return out;
  }*/


template <class TraitsT, class Visitor>
class Delaunay_triangulation_base_3;

template <class TraitsT, class Visitor>
std::ostream &operator<<(std::ostream &out,
			 const Delaunay_triangulation_base_3<TraitsT, Visitor> &dt);

// Common base class for Delaunay and regular triangulations
template <class TraitsT, class Visitor>
class Delaunay_triangulation_base_3
{
public:
  typedef Delaunay_triangulation_base_3<TraitsT, Visitor> This;

  // KDS typedefs
  typedef typename TraitsT::Simulator Simulator;
  typedef typename TraitsT::Active_points_3_table Moving_object_table;
  typedef typename TraitsT::Kinetic_kernel Kinetic_kernel;
  typedef typename TraitsT::Instantaneous_kernel Instantaneous_kernel;
  //typedef typename Simulator::Time Time;
  typedef typename Moving_object_table::Key Point_key;
  typedef typename Moving_object_table::Data Point;
  typedef typename Simulator::Event_key Event_key;
  typedef typename TraitsT::Side_of_oriented_sphere_3::result_type Certificate;
  //typedef typename Simulator::NT NT;

  // Delaunay typedefs
  typedef Triangulation_helper_3<typename TraitsT::Triangulation> Tri;

  // Triangulation members
  typedef typename Tri::Facet Facet;
  typedef typename Tri::Edge Edge;
  typedef typename Tri::Facet_circulator Facet_circulator;
  typedef typename Tri::Cell_circulator Cell_circulator;
  typedef typename Tri::Finite_edges_iterator Finite_edges_iterator;
  typedef typename Tri::Finite_facets_iterator Finite_facets_iterator;
  typedef typename Tri::Edge_iterator Edge_iterator;
  //typedef typename Tri::Facet_iterator Facet_iterator;
  typedef typename Tri::Vertex_handle Vertex_handle;
  typedef typename Tri::Cell_handle Cell_handle;
  typedef typename Tri::All_cells_iterator All_cells_iterator;
  typedef typename Tri::All_edges_iterator All_edges_iterator;
  typedef typename Tri::All_facets_iterator All_facets_iterator;
  typedef typename Tri::Finite_vertices_iterator Finite_vertices_iterator;

  // accessory types
  typedef typename TraitsT::Facet_flip Facet_flip_event;
  typedef typename Facet_flip_event::P Event_base;
  typedef typename TraitsT::Edge_flip Edge_flip_event;
  typedef Tri Triangulation;

  Delaunay_triangulation_base_3(TraitsT tr, Visitor v): tr_(tr),
							triangulation_(tr.instantaneous_kernel_object()),
							soc_(tr.side_of_oriented_sphere_3_object()),
							o3_(tr.orientation_3_object()),
							v_(v) {
    if (0) print();
    has_certificates_=false;
  }

  const TraitsT& simulation_traits_object() const {return tr_;}


  const Visitor& visitor() const
  {
    return v_;
  }
  Visitor& visitor() {
    return v_;
  }

  //! Just write the objects in order;
  template <class Stream>
  void write(Stream &out) const
  {
    //out << triangulation_;
    for (All_cells_iterator cit= triangulation_.all_cells_begin(); cit != triangulation_.all_cells_end();
	 ++cit) {
      Cell_handle h= cit;
      for (unsigned int i=0; i<4; ++i) {
	out << h->vertex(i)->point() << " ";
      }
      out << "--";
      for (unsigned int i=0; i<4; ++i) {
	Facet f(h, i);
	out << triangulation_.label(f) << " ";
      }
      out << std::endl;
    }
    for (All_edges_iterator eit = triangulation_.all_edges_begin(); eit != triangulation_.all_edges_end(); ++eit) {
      triangulation_.write_labeled_edge(*eit, out);
      if (is_degree_3(*eit)) {
	//out << " " << triangulation_.label(*eit);
      }
      else {
	if (has_event(*eit)) {
	  out << " ?? " ;           //<< triangulation_.label(*eit);
	}
	//CGAL_assertion(triangulation_.label(*eit)== Event_key::null());
      }
      out << std::endl;
    }
    for (All_facets_iterator eit = triangulation_.all_facets_begin(); eit != triangulation_.all_facets_end(); ++eit) {
      triangulation_.write_labeled_facet(*eit, out);
      if (!has_degree_3_edge(*eit)) {
	//out << " " << triangulation_.label(*eit);
      }
      else {
	//CGAL_assertion(triangulation_.label(*eit)== Event_key::null());
	if (has_event(*eit)) {
	  out << " ?? ";            //<< triangulation_.label(*eit);
	}
      }
      out << std::endl;
    }
  }
  bool is_degree_3(const Edge &e) const {
    return triangulation_.has_degree_3(e);
  }

  bool is_degree_4(Vertex_handle vh) const {
    return triangulation_.degree(vh) == 4;
  }

  bool has_degree_3_edge(const Facet& f) const {
    return triangulation_.has_degree_3(f);
  }

  bool has_degree_4_vertex(const Facet& f) const {
    return is_degree_4(triangulation_.vertex(f,0))
      || is_degree_4(triangulation_.vertex(f,1))
      || is_degree_4(triangulation_.vertex(f,2));
  }

  bool has_degree_4_vertex(const Edge& f) const {
    return is_degree_4(triangulation_.vertex(f,0))
      || is_degree_4(triangulation_.vertex(f,1));
  }

  bool print() const
  {
    write(std::cout);
    return true;
  }

  /*!
    Delete all incident face and edge events (have to push edge events on to facets).
    Insert point.
    Set edges to never fails. Set outside facets.
  */
  /*Vertex_handle push_vertex(Point_key k, Cell_handle c) {
    clean_cell(c);
    v_.pre_insert_vertex(k, c);
    // into max dim simplex?
    Vertex_handle vh=triangulation_.tds().insert_in_cell( c);
    vh->set_point(k);
    set_vertex_handle(k, vh);
    std::vector<Cell_handle> ics;
    triangulation_.incident_cells(vh, std::back_insert_iterator<std::vector<Cell_handle> >(ics));
    CGAL_postcondition(ics.size() == 4);
    for (unsigned int j=0; j< ics.size(); ++j) {
      Cell_handle cc= ics[j];
      handle_new_cell(cc);
    }
    v_.post_insert_vertex(vh);
    return vh;
    }*/

  Cell_handle pop_vertex(Vertex_handle v) {
    v_.pre_remove_vertex(v);
    CGAL_precondition(is_degree_4(v));
    std::vector<Cell_handle> ics;
    triangulation_.incident_cells(v, std::back_insert_iterator<std::vector<Cell_handle> >(ics));
    Point_key k= v->point();
    for (unsigned int j=0; j< ics.size(); ++j) {
      clean_cell(ics[j]);
    }
   
    set_vertex_handle(v->point(), NULL);

    Cell_handle h= triangulation_.tds().remove_from_maximal_dimension_simplex(v);
    handle_new_cell(h);
    v_.post_remove_vertex(k);

    return h;
  }

  void delete_vertex(Point_key) {
    CGAL_assertion(0);
  }


  Vertex_handle change_vertex(Point_key k) {
    if (!has_certificates()) return NULL;
    //triangulation_.geom_traits().set_time(simulator()->current_time_as_nt());
    std::vector<Cell_handle> incident_cells;
    triangulation_.incident_cells(vertex_handle(k), back_inserter(incident_cells));
    for (typename std::vector<Cell_handle>::iterator it= incident_cells.begin(); 
	 it != incident_cells.end();
	 ++it) {
      clean_cell(*it);
    }

    for (typename std::vector<Cell_handle>::iterator it= incident_cells.begin();
	 it != incident_cells.end(); ++it) {
      handle_new_cell(*it);
    }
    v_.change_vertex(vertex_handle(k));
    return vertex_handle(k);
  }

  /*typename Triangulation::Vertex_handle new_vertex_regular(Point_key k, Cell_handle h=NULL) {
    CGAL_precondition(!has_certificates_);
    typename Simulator::NT nt= simulator()->next_time_representable_as_nt();
    CGAL_precondition(simulator()->current_time() == nt);
    //simulator()->set_current_time(nt);
    triangulation_.geom_traits().set_time(nt);
    typename Triangulation::Vertex_handle vh= triangulation_.insert(k, h);
    v_.create_vertex(vh);
    return vh;
    }*/

  typename Triangulation::Vertex_handle insert(Point_key k) {
    /* NT nt= simulator()->next_time_representable_as_nt();
    simulator()->set_current_time(nt);
    triangulation_.geom_traits().set_time(nt);*/
    set_instantaneous_time();
    //typename Simulator::NT nt= simulator()->next_time_representable_as_nt();
    //simulator()->set_current_time(nt);
    //std::cout << "Locating at time " << triangulation_.geom_traits().time() << std::endl;
    Cell_handle h= triangulation_.locate(k);
    /*if (h != Cell_handle() && h->vertex(0) != Vertex_handle()
	       && h->vertex(1) != Vertex_handle()
	       && h->vertex(2) != Vertex_handle()
	       && h->vertex(3) != Vertex_handle()) {
      std::cout << "Kinetic located in " << h->vertex(0)->point() << " "
		<< h->vertex(1)->point() << " "
		<< h->vertex(2)->point() << " "
		<< h->vertex(3)->point() << std::endl;
    } else {
      std::cout << "Kinetic located outside hull" << std::endl;
      }*/
    return insert(k,h);
  }

  //!
  /*!
    Some old certificate edges will be lost, have to find all conflicts first.

    An important problem case is the edges which were previously degree 3 but are no longer
    after insertion. They can now be degree 4 and will have their certificate deleted properly.
    However, this leaves a face with no certificate (the one or more not destroyed by the new vertex).

    This occurs when a degree 3 edges has a facet added but none destroyed--i.e. a boundary edge
    with degree 3.
  */
  typename Triangulation::Vertex_handle insert(Point_key k, Cell_handle h) {
    //CGAL_precondition(h != NULL);

    std::vector<Facet> bfacets;
    std::vector<Cell_handle> cells;
    std::vector<Facet> ifacets;
    //typename Simulator::NT nt= simulator()->next_time_representable_as_nt();
    //CGAL_precondition(simulator()->current_time() == nt);
    //triangulation_.geom_traits().set_time(nt);
    set_instantaneous_time();
    CGAL_precondition(triangulation_.geom_traits().time() == simulator()->current_time());
    Vertex_handle vh;
    v_.pre_insert_vertex(k, h);
    if (triangulation_.dimension() == 3) {
      triangulation_.find_conflicts(k, h, back_inserter(bfacets), 
				    back_inserter(cells),back_inserter(ifacets));
      if (has_certificates_) {
	for (unsigned int i=0; i < cells.size(); ++i) {