delaunay_triangulation_2.h

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

	  ++fc;
	} while (fc != fe);
      }
    }
    // remove from triangulation
    del_.geom_traits().set_time(traits_.rational_current_time());
    del_.remove(vh);
    //new_edges_.clear();
    if (del_.dimension()==2 && has_certificates_) {
      std::vector<Face_handle> faces;

      del_.get_conflicts(k,std::back_inserter(faces));
      for (unsigned int i=0; i< faces.size(); ++i) {
	for (unsigned int j=0; j<3; ++j) {
	  update_neighbors(faces[i]->vertex(j));
	}
      }
      for (unsigned int i=0; i< faces.size(); ++i) {
	for (unsigned int j=0; j<3; ++j) {
	  update_vertex(faces[i]->vertex(j));
	}
      }
      for (unsigned int i=0; i< faces.size(); ++i) {
	for (unsigned int j=0; j<3; ++j) {
	  Edge e(faces[i],j);
	  //Event_key k= get_undirected_edge_label(e);
	  // a bit redundant for certificates which don't fail
	  update_edge(e);
	  //new_edges_.insert(e);
	}
	watcher_.create_face(faces[i]);
      }
      
    }
    watcher_.post_remove_vertex(k);
  }

  //! The assertion will catch that the object is in the same sorted order
  void set(Point_key k) {
    //std::cout << "Object changed " << k << std::endl;

    //new_edges_.clear();
    traits_.point_changed(k);
    if (del_.dimension() != 2) {
      CGAL_KINETIC_LOG(CGAL::Kinetic::LOG_SOME,"Triangulation is still 1D.\n");
      return;
    }

    Vertex_handle vh=vertex_handle(k);
    if (vh == Vertex_handle()) {
      CGAL_KINETIC_LOG(CGAL::Kinetic::LOG_SOME, "Point " << k << " is not in triangulation on set."<< std::endl);
      return;
    }

    if (has_certificates_) {
      Face_handle f= vh->face(), fe= vh->face();
      int i= f->index(vh);
      do {
	//int i= fc->index(vh);
	Edge e0(f, i);
	delete_certificate(e0);
	update_edge(e0);
	Edge e1=Edge(f, (i+1)%3);
	delete_certificate(e1);
	update_edge(e1);
	f= f->neighbor((i+1)%3);
	i= f->index(vh);
      } while (f != fe);
    }

   
    watcher_.change_vertex(vh);

    /*if (has_certificates_) {
      Face_handle f= vh->face(), fe= vh->face();
      int i= f->index(vh);
      do {
	//int i= fc->index(vh);
	Edge e0(f, i);
	update_edge(e0);
	Edge e1=Edge(f, (i+1)%3);
	update_edge(e1);
	f= f->neighbor((i+1)%3);
	i= f->index(vh);
      } while (f != fe);
      }*/
    //write(std::cout);
  }


  void insert(Point_key k) {
    // evil hack
    CGAL_precondition(k.to_index() >= vhs_.size() || vertex_handle(k) == Vertex_handle());
    CGAL_DELAUNAY_2_DEBUG(std::cout << "Inserting " << k << std::endl);
    bool was_2d= (del_.dimension()==2);

    del_.geom_traits().set_time(traits_.rational_current_time());
    if (was_2d && has_certificates_) {
      //std::cout << "removing extra certificates.\n";
      std::vector<Face_handle> faces;
      del_.get_conflicts(k, std::back_inserter(faces));
      for (unsigned int i=0; i< faces.size(); ++i) {
	Face_handle f= faces[i];
	for (unsigned int j=0; j<3; ++j) {
	  Edge e(f, j);
	  Event_key k= get_undirected_edge_label(e);
	  delete_certificate(e);
	}
	watcher_.destroy_face(faces[i]);
      }

      if (faces.empty()) {
	CGAL_KINETIC_LOG(CGAL::Kinetic::LOG_SOME, "DELAUNAY vertex not successfully inserted " << k << std::endl);
	return;
      }
    }
    watcher_.pre_insert_vertex(k);
    Vertex_handle vh= del_.insert(k);
    set_vertex_handle(k, vh);
    
    CGAL_assertion(vertex_handle(k) != Vertex_handle());
    
    CGAL_DELAUNAY_2_DEBUG(std::cout << "Vertex " << vertex_handle(k)->point() 
			  << " has " << vertex_handle(k)->neighbors() << std::endl);
   

    // now have to update
    if (has_certificates_) {
   

      if (!was_2d && del_.dimension()==2) {
	vh->set_neighbors(vh->degree());
	has_certificates_=false;
	set_has_certificates(true, 0); 
      } else if (del_.dimension() == 2) {
	update_neighbors(vh);
	update_vertex(vh);
	// update vertices
	{
	  Face_handle f= vh->face(), fe= vh->face();
	  int i= f->index(vh);
	  do {
	    update_neighbors(f->vertex((i+1)%3));
	    f= f->neighbor((i+1)%3);
	    i= f->index(vh);
	  } while (f != fe);
	} 

	{
	  Face_handle f= vh->face(), fe= vh->face();
	  int i= f->index(vh);
	  do {
	    update_vertex(f->vertex((i+1)%3));
	    f= f->neighbor((i+1)%3);
	    i= f->index(vh);
	  } while (f != fe);
	} 
	// update edges 
	{
	  Face_handle f= vh->face(), fe= vh->face();
	  int i= f->index(vh);
	  do {
	    //int i= fc->index(vh);
	    Edge e0(f, i);
	    update_edge(e0);
	    Edge e1=Edge(f, (i+1)%3);
	    update_edge(e1);
	    f= f->neighbor((i+1)%3);
	    i= f->index(vh);
	  } while (f != fe);
	}
      }
    } else {
      vertex_handle(k)->set_neighbors(vh->degree());
    }
    watcher_.post_insert_vertex(vh);
    //write(std::cout);
    //if (del_.dimension()==2) audit();
  }



  Comparison_result compare_concurrent(Event_key a, Event_key b) const {
    Edge ea= traits_.simulator_handle()->template event<Event>(a).edge();
    Edge eb= traits_.simulator_handle()->template event<Event>(b).edge();
    return traits_.compare_concurrent(a, ea, b, eb);
  }




  Edge flip(const Edge &e, Certificate_data cert) {
    ++num_events_;
    CGAL_precondition(!batching_);
    CGAL_KINETIC_LOG(CGAL::Kinetic::LOG_SOME, "\n\n\n\n\n\nDELAUNAY Flipping edge "
		     << TDS_helper::origin(e)->point()
		     << TDS_helper::destination(e)->point() 
		     << " to get " << TDS_helper::third_vertex(e)->point()
		     << ", " << TDS_helper::mirror_vertex(e)->point()<< std::endl);
    //CGAL_KINETIC_LOG(CGAL::Kinetic::LOG_NONE, TDS_helper::destination(e)->point() << std::endl);
    //CGAL_KINETIC_LOG(CGAL::Kinetic::LOG_SOME, " at "  << traits_.simulator()->current_time() << std::endl);

    

    Face_handle face= e.first;
    int index= e.second;
    Face_handle mirror_face = face->neighbor(index);
    int mirror_index =face->mirror_index(index);
    Edge em(mirror_face,mirror_index);
    CGAL_precondition(mirror_face->neighbor(mirror_index) == face);

    Face_handle bef;
    int bei;

    if (!traits_.is_exact()
	&& del_.is_edge(TDS_helper::third_vertex(e), TDS_helper::mirror_vertex(e),
			bef, bei)) {
      // we have a numeric error, lets try to rebuild the neighboring certificates
      CGAL_KINETIC_LOG(CGAL::Kinetic::LOG_SOME,
		       "DELAUNAY ERROR not flipping unflippable edge" << std::endl);
      CGAL_KINETIC_LOG(CGAL::Kinetic::LOG_SOME, 
		       *traits_.simulator_handle() << std::endl;);
      //make this better
      //double ub=to_interval(traits_.simulator_handle()->next_event_time()).second;
      Edge bad_edge(bef, bei);
      Event_key bek= get_undirected_edge_label(bad_edge);
      
      if (bek == traits_.simulator_handle()->null_event()) {
	CGAL_KINETIC_LOG(CGAL::Kinetic::LOG_SOME,
			 "Dropping the event." << std::endl);
	set_undirected_edge_label(e, Event_key());
	return e;
      } else {

	double ub = CGAL::to_interval(traits_.simulator_handle()->event_time(bek)).second;
	
	ub= (std::max)(ub+.0000001, 
		       nextafter(ub, (std::numeric_limits<double>::max)()));
	Time t(ub);
	CGAL_precondition(CGAL::compare(t, traits_.simulator_handle()->next_event_time()) == CGAL::LARGER);
	Event_key k =traits_.simulator_handle()->new_event(t, Event(cert, e, this));
	set_undirected_edge_label(e, k);
	return e;
      }
    }
    CGAL_precondition(!del_.is_edge(TDS_helper::third_vertex(e), TDS_helper::mirror_vertex(e),
				    bef, bei));

    set_directed_edge_label(e, Event_key());
    set_directed_edge_label(em, Event_key());

    delete_certificate(Edge(face, (index+1)%3));
    delete_certificate(Edge(face, (index+2)%3));
    delete_certificate(Edge(mirror_face, (mirror_index+1)%3));
    delete_certificate(Edge(mirror_face, (mirror_index+2)%3));

   
    watcher_.pre_flip(e);
    del_.tds().flip(face,index);
   
    // we also know that CGAL preserves the edge index of the flipped edge
    mirror_index = mirror_face->index(face);
    index= face->index(mirror_face);

    Edge flipped_edge(face,index);
    Edge mirror_flipped_edge(face->neighbor(index), mirror_index);
    CGAL_assertion(mirror_flipped_edge == TDS_helper::mirror_edge(flipped_edge));
    //CGAL_postcondition(del_.is_face(face));

    CGAL_assertion(mirror_index == face->mirror_index(index));
    CGAL_assertion(mirror_face == face->neighbor(index));

    watcher_.post_flip(flipped_edge);

    decrease_neighbors(flipped_edge.first->vertex(flipped_edge.second));
    increase_neighbors(flipped_edge.first->vertex((flipped_edge.second+1)%3));
    increase_neighbors(flipped_edge.first->vertex((flipped_edge.second+2)%3));
    decrease_neighbors(mirror_flipped_edge.first->vertex(mirror_flipped_edge.second));

    // make sure it is not created by a 3-4 transition
    set_directed_edge_label(flipped_edge, traits_.simulator_handle()->null_event());
    set_directed_edge_label(mirror_flipped_edge, traits_.simulator_handle()->null_event());

    update_decreased_vertex(flipped_edge.first->vertex(flipped_edge.second));
    update_increased_vertex(flipped_edge.first->vertex((flipped_edge.second+1)%3));
    update_increased_vertex(flipped_edge.first->vertex((flipped_edge.second+2)%3));
    update_decreased_vertex(mirror_flipped_edge.first->vertex(mirror_flipped_edge.second));

    update_edge(Edge(flipped_edge.first, (flipped_edge.second+1)%3));
    update_edge(Edge(flipped_edge.first, (flipped_edge.second+2)%3));
    update_edge(Edge(mirror_flipped_edge.first, (mirror_flipped_edge.second+1)%3));
    update_edge(Edge(mirror_flipped_edge.first, (mirror_flipped_edge.second+2)%3));

    {

      Time t; Certificate_data cd;
      if (traits_.certificate_failure_time(flipped_edge,cert, t, cd)){
	Event_key k =traits_.simulator_handle()->new_event(t,
							   Event(cd, flipped_edge, this));
	set_directed_edge_label(flipped_edge, k);
	set_directed_edge_label(mirror_flipped_edge, k);
      } else {
	set_directed_edge_label(flipped_edge, traits_.simulator_handle()->null_event());
	set_directed_edge_label(mirror_flipped_edge, traits_.simulator_handle()->null_event());
      }
    }
    //write(std::cout);
    //new_edges_.clear();
    //new_edges_.insert(flipped_edge);

    CGAL_KINETIC_LOG(CGAL::Kinetic::LOG_SOME, "Created " << TDS_helper::origin(flipped_edge)->point());
    CGAL_KINETIC_LOG(CGAL::Kinetic::LOG_SOME, TDS_helper::destination(flipped_edge)->point() << std::endl);

    return flipped_edge;
  }

  Visitor &visitor() {
    return watcher_;
  }

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


  bool has_event(const Edge &e) const {
    return get_directed_edge_label(e) != Event_key();
  }

  bool has_finite_event(const Edge &e) const {
    return has_event(e) && get_directed_edge_label(e) != traits_.simulator_handle()->null_event();
  }

protected:
  Traits traits_;
  Visitor watcher_;
  Triangulation del_;


  Simulator_listener siml_;
  Moving_point_table_listener motl_; 
  std::vector<Vertex_handle> vhs_;

  bool has_certificates_;
  bool batching_;
  std::vector<Edge> batched_certs_;


  mutable unsigned int num_events_;
  //mutable unsigned int num_single_certificates_;

  const typename Traits::Point_2& point(Point_key k) const
  {
    return traits_.point(k);
  }


  Vertex_handle vertex_handle(Point_key k) const {
    //if (k.to_index() >= vhs_.size()) return Vertex_handle();
    CGAL_precondition(k.to_index() < vhs_.size());
    return vhs_[k.to_index()];
  }
  void set_vertex_handle(Point_key k, Vertex_handle vh) {
    vhs_.resize(std::max BOOST_PREVENT_MACRO_SUBSTITUTION(static_cast<unsigned int>(k.to_index()+1),
			 static_cast<unsigned int>(vhs_.size())));
    vhs_[k.to_index()]=vh;
  }

  void update_vertex_to_degree_3(Vertex_handle vh) {
    CGAL_DELAUNAY_2_DEBUG(std::cout << "Degree 3 for " 
			  << vh->point() << std::endl);
    typename Triangulation::Edge_circulator ec= del_.incident_edges(vh);
    do {
      delete_certificate(*ec);
      ++ec;
    } while (ec != del_.incident_edges(vh));
  }

  void update_vertex_to_not_degree_3(Vertex_handle vh) {
    CGAL_DELAUNAY_2_DEBUG(std::cout << "Degree 4 for " 
			      << vh->point() << std::endl);
    typename Triangulation::Edge_circulator ec= del_.incident_edges(vh);
    do {
      if (get_undirected_edge_label(*ec) == Event_key()) {
	// check other vertex, it it is not changed and not 3, build
	Vertex_handle ov= ec->first->vertex((ec->second+1)%3);
	if (ov== vh) {
	  ov= ec->first->vertex((ec->second+2)%3);
	}
	
	if (ov->neighbors() != 3){
	  new_certificate(*ec);