delaunay_triangulation_2.h

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

	  CGAL_DELAUNAY_2_DEBUG(std::cout << "New cert for "
				<< TDS_helper::origin(*ec)->point() << " " 
				<< TDS_helper::destination(*ec)->point() 
				<< std::endl);
	} else {
	  CGAL_DELAUNAY_2_DEBUG(std::cout << "Not creating cert for "
				<< TDS_helper::origin(*ec)->point() << " " 
				<< TDS_helper::destination(*ec)->point() 
				<< std::endl);
	}
      }
      ++ec;
    } while (ec != del_.incident_edges(vh));
  }

  void update_neighbors(Vertex_handle vh) {
    unsigned int deg= vh->degree();
    if (deg == vh->neighbors()) return;
    if (deg ==3) {
      vh->set_neighbors(3);
      //update_vertex_to_degree_3(vh);
    } else if (vh->neighbors()==3) {
      vh->set_neighbors(deg);
      //update_vertex_to_not_degree_3(vh);
    }  else {
      vh->set_neighbors(deg);
    }
  }
  
  void decrease_neighbors(Vertex_handle vh) {
    vh->set_neighbors(vh->neighbors()-1);
    CGAL_assertion(neighbors_ok(vh));
    //if (vh->neighbors() == 3) update_vertex_to_degree_3(vh);
  }

  void increase_neighbors(Vertex_handle vh) {
    vh->set_neighbors(vh->neighbors()+1);
    CGAL_assertion(neighbors_ok(vh));
    //if (vh->neighbors() == 4) update_vertex_to_not_degree_3(vh);
  }

  void update_vertex(Vertex_handle vh) {
    CGAL_precondition(neighbors_ok(vh));
    if (vh->neighbors() == 3) {
      update_vertex_to_degree_3(vh);
    } else if (vh->neighbors() == 4) {
      update_vertex_to_not_degree_3(vh);
    }
  }

  void update_decreased_vertex(Vertex_handle vh) {
    if (vh->neighbors() == 3) {
      update_vertex_to_degree_3(vh);
    } 
  }

 void update_increased_vertex(Vertex_handle vh) {
    if (vh->neighbors() == 4) {
      update_vertex_to_not_degree_3(vh);
    } 
  }

  bool neighbors_ok(Vertex_handle vh) const {
    return vh->degree() == vh->neighbors();
    //(vh->neighbors() == vh->degree() && vh->degree() <6 
    //|| vh->degree() >=5 && vh->neighbors() >=5);
  }

  void update_edge_no_batch(const Edge &e) {
    CGAL_DELAUNAY_2_DEBUG(std::cout << "Updating edge " 
			  << TDS_helper::origin(e)->point() << " " 
			  << TDS_helper::destination(e)->point() 
			  << std::endl);
    Vertex_handle ov=TDS_helper::origin(e);
    Vertex_handle dv=TDS_helper::destination(e);
    
    CGAL_precondition(neighbors_ok(ov));
    CGAL_precondition(neighbors_ok(dv));

    if (get_undirected_edge_label(e) != Event_key()) {
      CGAL_DELAUNAY_2_DEBUG(std::cout << "Already has event " << std::endl);
      // can't do this since I create all edges around vertex of degree 4 at once
      // CGAL_assertion(0);
    } else if (ov->neighbors() ==3 
	       || dv->neighbors() ==3) {
      CGAL_DELAUNAY_2_DEBUG(std::cout << "One end has 3 " << std::endl);
    } else {
      //CGAL_DELAUNAY_2_DEBUG(std::cout << "New certificate" << std::endl);
      new_certificate(e);
    }
  }



  void update_edge(const Edge &e) {
    if (batching_) {
      //delete_certificate(e);
      batched_certs_.push_back(e);
    } else if (get_undirected_edge_label(e) == Event_key()) {
      update_edge_no_batch(e);
    }
  }

  // return true if hull
  bool points(const Edge &e, Point_key ks[4]) const {
    ks[0]= TDS_helper::origin(e)->point();
    ks[1]= TDS_helper::third_vertex(e)->point();
    ks[2]= TDS_helper::destination(e)->point();
    ks[3]= TDS_helper::mirror_vertex(e)->point();

    bool odd_parity=false;
    bool infinity=false;
    for (unsigned int i=0; i<4; ++i) {
      if (infinity) {
	ks[i-1]=ks[i];
      }
      else {
	if (!ks[i].is_valid()) {
	  infinity=true;
	  odd_parity= ((i%2)==1);
	} 
      }
    }
    if (odd_parity) {
      std::swap(ks[0], ks[1]);
    }
    return infinity;
  }


  /*bool is_hull_edge(const Edge &e) const {
    return ! TDS_helper::mirror_vertex(e)->point().is_valid()
      || ! TDS_helper::third_vertex(e)->point().is_valid()
      || ! TDS_helper::origin(e)->point().is_valid()
      || ! TDS_helper::destination(e)->point().is_valid();
  }


 
  void edge_points(const Edge &e, Point_key ks[4]) const {
    ks[0]= TDS_helper::origin(e)->point();
    ks[1]= TDS_helper::third_vertex(e)->point();
    ks[2]= TDS_helper::destination(e)->point();
    ks[3]= TDS_helper::mirror_vertex(e)->point();
  }

  // very dangerous

  void hull_points(const Edge &e,
		   Point_key ks[4]) const {
    ks[0]= TDS_helper::origin(e)->point();
    ks[1]= TDS_helper::third_vertex(e)->point();
    ks[2]= TDS_helper::destination(e)->point();
    ks[3]= TDS_helper::mirror_vertex(e)->point();

    bool odd_parity=false;
    bool infinity=false;
    for (unsigned int i=0; i<4; ++i) {
      if (infinity) {
	ks[i-1]=ks[i];
      }
      else {
	if (!ks[i].is_valid()) {
	  infinity=true;
	  odd_parity= ((i%2)==1);
	} 
      }
    }
    if (odd_parity) {
      std::swap(ks[0], ks[1]);
    }
    }*/


  void new_certificate(Edge e) {
    CGAL_precondition(get_undirected_edge_label(e) == Event_key());
    CGAL_DELAUNAY_2_DEBUG(std::cout << "\nMaking certificate for " << TDS_helper::origin(e)->point() << " " 
			  << TDS_helper::destination(e)->point() 
			  << " which would make " << TDS_helper::mirror_vertex(e)->point() << " " 
			  << TDS_helper::third_vertex(e)->point()
			  << std::endl);

    Time t; Certificate_data cd;
    Point_key ks[4];
    if (points(e,ks)) {
      if (traits_.hull_certificate_failure_time(e, ks, t, cd)) {
	Event_key k= traits_.simulator_handle()->new_event(t, Event(cd, e, this));
	set_undirected_edge_label(e, k);
      } else {
	set_undirected_edge_label(e, traits_.simulator_handle()->null_event());
      }
    } else {
      if (traits_.internal_certificate_failure_time(e, ks, t, cd)) {
	Event_key k= traits_.simulator_handle()->new_event(t, Event(cd, e, this));
	set_undirected_edge_label(e, k);
      } else {
	set_undirected_edge_label(e, traits_.simulator_handle()->null_event());
      }
    }
  }

  void delete_certificate(Edge e) {
    CGAL_DELAUNAY_2_DEBUG(std::cout << "Cleaning edge " << TDS_helper::origin(e)->point() << " " << TDS_helper::destination(e)->point() << std::endl);
    Event_key k=  get_undirected_edge_label(e);
    if (k != Event_key()) {
      traits_.simulator_handle()->delete_event(k);
      set_undirected_edge_label(e, Event_key());
    }
  }


  Edge canonicalize(Edge e) const {
    if (e.first->neighbor(e.second) < e.first) {
      return TDS_helper::mirror_edge(e);
    } else {
      return e;
    }
  }

  static Event_key get_directed_edge_label(const Edge &e) {
    return e.first->get_edge_label(e.second);
  }

  static Event_key get_undirected_edge_label(const Edge &e) {
#ifndef NDEBUG
    if (get_directed_edge_label(e) != get_directed_edge_label(TDS_helper::mirror_edge(e))) {
      std::cerr << "FAILURE Edge from " << TDS_helper::origin(e)->point() << " to " 
		<< TDS_helper::destination(e)->point() << " is screwed." << std::endl;
      std::cerr << get_directed_edge_label(e) << " "
                <<  get_directed_edge_label(TDS_helper::mirror_edge(e)) << std::endl;
      CGAL_precondition(get_directed_edge_label(e)
			== get_directed_edge_label(TDS_helper::mirror_edge(e)));
    }
#endif
    return e.first->get_edge_label(e.second);
  }

  static void set_directed_edge_label(const Edge &e,
				      Event_key l) {
    e.first->set_edge_label(e.second, l);
  }

  static void set_undirected_edge_label(const Edge &e,
					Event_key l) {
    set_directed_edge_label(e,l);
    set_directed_edge_label(TDS_helper::mirror_edge(e),l);
  }

};

template <class Sim, class Del, class W, class T>
std::ostream &operator<<(std::ostream &out, const Delaunay_triangulation_2<Sim, Del, W, T> &kd)
{
  kd.write(out);
  return out;
}

template <class Sim, class Del, class W, class T>
void Delaunay_triangulation_2<Sim, Del, W, T>::audit() const
  {
    if (!has_certificates_) return;
    CGAL_KINETIC_LOG(CGAL::Kinetic::LOG_SOME, "Auditing delaunay" << std::endl);
    
    if (del_.number_of_vertices() < 50) {
      CGAL_KINETIC_LOG(CGAL::Kinetic::LOG_LOTS, *this);
    }
    if (del_.dimension() != 2) return;
    Basic_Delaunay sdel(traits_.instantaneous_kernel_object());
    sdel.geom_traits().set_time(traits_.simulator_handle()->audit_time());
    for (typename Triangulation::Finite_vertices_iterator cit= del_.finite_vertices_begin();
	 cit != del_.finite_vertices_end(); ++cit){
      sdel.insert(cit->point());
    }
    /*    sdel.insert(traits_.active_points_2_table_handle()->keys_begin(),
	  traits_.active_points_2_table_handle()->keys_end());*/

    //CGAL_KINETIC_LOG(CGAL::Kinetic::LOG_LOTS, sdel << std::endl);

    if (del_.dimension() != sdel.dimension()) {
      CGAL_KINETIC_LOG(CGAL::Kinetic::LOG_NONE, "AUDIT FAILURE Dimensions don't match in audit" << std::endl);
      return;
    }
    CGAL_exactness_assertion(del_.dimension() == sdel.dimension());
   
   for (typename Triangulation::All_vertices_iterator vit = del_.all_vertices_begin();
	vit != del_.all_vertices_end(); ++vit) {
     if (vit->point() != Point_key()) {
       
       if (!neighbors_ok(vit)) {
	 CGAL_KINETIC_LOG(CGAL::Kinetic::LOG_NONE, "AUDIT FAILURE stored degree is " << vit->neighbors() 
			  << " and actual is " << vit->degree() << " for " << vit->point() << std::endl);
	 CGAL_exactness_assertion(neighbors_ok(vit));
       }
     }
   }
   for (typename Basic_Delaunay::All_vertices_iterator vit = sdel.all_vertices_begin();
	 vit != sdel.all_vertices_end(); ++vit) {
      bool found=false;
    
      //Object_key k= vit->point();
      for (typename Triangulation::All_vertices_iterator vit2= del_.all_vertices_begin();
	   vit2 != del_.all_vertices_end(); ++vit2) {
	//Object_key k2= vit2->point();
	if (vit->point() == vit2->point()) {
	  found=true;
	  //int d= vit->degree();
	  //int d2= vit2->degree();
	  if (vit->degree() != vit2->degree()) {
	    CGAL_KINETIC_LOG(CGAL::Kinetic::LOG_NONE, "AUDIT FAILURE Degrees don't match in: " 
			     << vit->point() << std::endl);
	  }
	  CGAL_exactness_assertion(vit->degree() == vit2->degree());


	}
      }
      if (!found) {
	CGAL_KINETIC_LOG(CGAL::Kinetic::LOG_NONE, "AUDIT FAILURE Matching vertex not found: " 
			 << vit->point() << std::endl);
      }
      CGAL_exactness_assertion(found);
     
    }

  

    typename Simulation_traits::Instantaneous_kernel ik= traits_.instantaneous_kernel_object();
    ik.set_time(traits_.simulator_handle()->audit_time());

    typename Simulation_traits::Instantaneous_kernel::Side_of_oriented_circle_2 ic2
      = ik.side_of_oriented_circle_2_object();
    for (typename Triangulation::Finite_edges_iterator fit = del_.finite_edges_begin(); 
	 fit != del_.finite_edges_end(); ++fit){
      Point_key k0= fit->first->vertex((fit->second+1)%3)->point();
      Point_key k2= fit->first->vertex((fit->second+2)%3)->point();
      Point_key k3= TDS_helper::mirror_vertex(*fit)->point();
      Point_key k1= TDS_helper::third_vertex(*fit)->point();
      if (k1== Point_key() || k3== Point_key()) continue;
      typename Triangulation::Geom_traits::Current_coordinates cc= del_.geom_traits().current_coordinates_object();
      /*typedef typename Triangulation::Geom_traits::Current_coordinates::result_type P2;
      P2 p0= cc(k0);
      P2 p1= cc(k1);
      P2 p2= cc(k2);
      P2 p3= cc(k3);*/
      if (ic2(k0, k1, k2, k3) != CGAL::ON_POSITIVE_SIDE) {
	CGAL_KINETIC_LOG(CGAL::Kinetic::LOG_NONE, "AUDIT FAILURE Failed certificate: " << k0 << " " << k1 << " " 
			 << k2 << " " << k3 << std::endl);
	CGAL_KINETIC_LOG(CGAL::Kinetic::LOG_NONE, "AUDIT FAILURE Points are: " << cc(k0) << ": " << cc(k1) << ": " << cc(k2) 
			 << ": " << cc(k3) << std::endl);
      }
      CGAL_exactness_assertion(ic2(k0, k1, k2, k3) == CGAL::ON_POSITIVE_SIDE);
      
    }

    for (typename Triangulation::Edge_iterator fit = del_.edges_begin(); fit != del_.edges_end(); ++fit){
      if (get_directed_edge_label(*fit) != 
	  get_directed_edge_label(TDS_helper::mirror_edge(*fit))) {
	CGAL_KINETIC_LOG(CGAL::Kinetic::LOG_NONE, "AUDIT FAILURE mismatched labels on " 
			 << TDS_helper::origin(*fit)->point() << " " 
			 << TDS_helper::destination(*fit)->point() 
			 << " front has " << get_directed_edge_label(*fit)
			 << " and back has " << get_directed_edge_label(TDS_helper::mirror_edge(*fit))<< std::endl);
      }
      if (TDS_helper::origin(*fit)->degree()==3 || TDS_helper::destination(*fit)->degree()==3) {
	if (get_undirected_edge_label(*fit) != Event_key()) {
	  CGAL_KINETIC_LOG(CGAL::Kinetic::LOG_NONE, "AUDIT FAILURE certificate on degree 3 edge: " 
			   << TDS_helper::origin(*fit)->point()
			   << " " <<  TDS_helper::destination(*fit)->point() 
			   << TDS_helper::origin(*fit)->degree() <<  " "
			   << TDS_helper::destination(*fit)->degree() << std::endl);
	} 
	CGAL_exactness_assertion(get_undirected_edge_label(*fit) == Event_key());
      } else {
	if (get_undirected_edge_label(*fit) == Event_key()) {
	  CGAL_KINETIC_LOG(CGAL::Kinetic::LOG_NONE, "AUDIT FAILURE no certificate on edge: " 
			   << TDS_helper::origin(*fit)->point()
			   << " " <<  TDS_helper::destination(*fit)->point() << std::endl);
	  CGAL_KINETIC_LOG(CGAL::Kinetic::LOG_NONE, "AUDIT FAILURE degrees are: " 
			   << TDS_helper::origin(*fit)->degree()
			   << " " <<  TDS_helper::destination(*fit)->degree() << std::endl);
	  
	} 
	CGAL_exactness_assertion(get_undirected_edge_label(*fit) != Event_key());
      }
    }

  }

CGAL_KINETIC_END_NAMESPACE
#endif