delaunay_triangulation_base_3.h

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

    for (unsigned int i=0; i<4; ++i) {
      Point_key k= f.first->vertex(i)->point();
      if (!k.is_valid()) {
	hinf=i;
	break;
      }
    }
    if (hinf==-1) {
      Point_key k= triangulation_.mirror_vertex(f.first, f.second)->point();
      if ( ! k.is_valid() ) {
	hinf=4;
      }
    }
    if (hinf != -1) {
      CGAL_KINETIC_LOG(LOG_LOTS, "hinf is " << hinf << std::endl);
      if (hinf ==4) {
	for (unsigned int i=0; i<4; ++i) {
	  Point_key k= f.first->vertex(i)->point();
	  *out= k;
	  ++out;
	}
	return;
      }
      else {
	//Facet ff(f.first, hinf);
	if (hinf%2!=0) {
	  Point_key k= triangulation_.mirror_vertex(f.first, f.second)->point();
	  *out= k;
	  ++out;
	}
	for (int i=0; i<4; ++i) {
	  // CGAL infinite cells seem to be misoriented
	  if (i==hinf) continue;
	  *out= f.first->vertex(i)->point();
	  ++out;
	}
	if (hinf%2==0) {
	  Point_key k= triangulation_.mirror_vertex(f.first, f.second)->point();
	  *out= k;
	  ++out;
	}
      }
    }
    else {
      for (unsigned int i=0; i<4; ++i) {
	Point_key k= f.first->vertex(i)->point();
	*out= k;
	++out;
      }
      Point_key k= triangulation_.mirror_vertex(f.first, f.second)->point();
      *out= k;
      ++out;
    }
  }




private:
  Certificate root_stack(const Edge &e, const typename Simulator::Time &st) const
  {
    return root_stack(*triangulation_.incident_facets(e), st);
  }

  

  void make_no_events(const Edge &e) {
    CGAL_precondition(triangulation_.has_degree_3(e));
    Facet_circulator fc= triangulation_.incident_facets(e);
    Facet_circulator fe= fc;
    do {
      if (has_event(*fc)) {
	Event_key k= triangulation_.label(*fc);
	simulator()->delete_event(k);
	triangulation_.set_label(*fc, Event_key());
      }
    } while(++fc != fe);
  }
public:
  Point_key replace_vertex(Vertex_handle vh, Point_key k) {
    Point_key ok= vh->point();
    vh->point()=k;
    vhs_[ok.to_index()]= Vertex_handle();
    vhs_[k.to_index()]= vh;
    return ok;
  }
protected:
  bool has_event(const Edge &e) const
  {
    return triangulation_.label(e) != Event_key();
  }

  bool has_event(const Facet &e) const
  {
    
    return triangulation_.label(e) != Event_key();
  }


  void create_edge_flips(Vertex_handle v) {
    CGAL_precondition(!is_degree_4(v));
    std::vector<Cell_handle> ics;
    triangulation().incident_cells(v, std::back_inserter(ics));
    for (unsigned int i=0; i< ics.size(); ++i) {
      int j=-1;// disable warning
      bool ret=ics[i]->has_vertex(v, j); // initializes j
      CGAL_assertion(j != -1);
      CGAL_assertion(ret);
      for (int k=0; k<4 ; ++k) {
	if (k==j) continue;
	Edge e(ics[i], j, k);
	if (is_degree_3(e) && !has_event(e) && !has_degree_4_vertex(e)) {
	  // rather than make_certificate due to ordering dependencies
	  make_edge_flip(e);
	}
      }
    }
  }

  void suppress_edge_flips(Vertex_handle v) {
    CGAL_precondition(is_degree_4(v));
    std::vector<Cell_handle> ics;
    triangulation().incident_cells(v, std::back_inserter(ics));
    for (unsigned int i=0; i< ics.size(); ++i) {
      int j=-1; // keep some dumb compiler happy
      bool ret=ics[i]->has_vertex(v, j);
      CGAL_assertion(ret);
      for (int k=0; k<4 ; ++k) {
	if (k==j) continue;
	Edge e(ics[i], j, k);
	if (has_event(e)) {
	  simulator()->delete_event(triangulation_.label(e));
	  triangulation().set_label(e, Event_key());
	}
      }
    }
  }

  void make_edge_flip(Edge &edge) {
    CGAL_KINETIC_LOG(LOG_LOTS, "Making edge flip ");
    CGAL_KINETIC_LOG_WRITE(LOG_LOTS, triangulation_.write_labeled_edge(edge, LOG_STREAM ));
    CGAL_KINETIC_LOG(LOG_LOTS,std::endl);
    CGAL_assertion(triangulation_.has_degree_3(edge));
    typename Simulator::Event_key k= typename Simulator::Event_key();
    Facet_circulator fc= triangulation_.incident_facets(edge), fe=fc;
    do {
      if (has_event(*fc)) {
	CGAL_assertion( k == Event_key() );
	k=triangulation_.label(*fc);
	triangulation_.set_label(*fc, typename Simulator::Event_key());
	Event_key kk=change_to_edge_flip(edge, k);
	triangulation_.set_label(edge, kk);
	simulator()->audit_event(kk);
	return;
      }
    } while (++fc != fe);


    CGAL_KINETIC_LOG(LOG_LOTS, "Making up edge event.\n");
    make_certificate(edge);
  }

  void make_not_edge_flip(Edge &edge, Cell_handle h) {
    if (true) {
      CGAL_KINETIC_LOG(LOG_LOTS, "Making edge ");
      CGAL_KINETIC_LOG_WRITE(LOG_LOTS, triangulation_.write_labeled_edge(edge, LOG_STREAM ));
      CGAL_KINETIC_LOG(LOG_LOTS, " not an edge flip.\n");
    }
    CGAL_assertion(is_degree_3(edge) || print());
    typename Simulator::Event_key index = triangulation_.label(edge);
    CGAL_precondition( index != Event_key() );

    triangulation_.set_label(edge, typename Simulator::Event_key());

    Cell_circulator fc=triangulation_.incident_cells(edge), pfc=fc, ef;
    ++fc;
    ef=fc;
    do {
      if( (h != fc) && ( h != pfc) ) {
	for (unsigned int i=0; i< 4; ++i) {
	  if (pfc->neighbor(i) == fc){
	    Facet f(pfc, i);
	    CGAL_precondition(!has_event(f));
	    triangulation_.set_label(f, change_to_facet_flip(f, index));
	    simulator()->audit_event(triangulation_.label(f));
	  }
	}
      }
      ++pfc; ++fc;
    } while (ef != fc);
  }

  template <class C>
  std::back_insert_iterator<C> back_inserter(C &c) const
  {
    return std::back_insert_iterator<C>(c);
  }

  bool contains(typename std::vector<Facet>::iterator beg,
		typename std::vector<Facet>::iterator end, Facet f) {
    Facet of = triangulation_.opposite(f);
    for (; beg != end; ++beg) {
      if (f.first == beg->first) {
	if (f.second == beg->second) return true;
      }
      else if (of.first == beg->first) {
	if (of.second == beg->second) return true;
      }
    }
    return false;
  }

  bool labeling_is_valid() const
  {
    for (All_facets_iterator eit= triangulation_.all_facets_begin(); eit != triangulation_.all_facets_end();
	 ++eit) {
      Facet f= *eit;
      Facet of= triangulation_.opposite(f);
      if (triangulation_.label(f) != triangulation_.label(of)) {
	triangulation_.write_labeled_facet(f, std::cerr);
	std::cerr << " does not match ("<< triangulation_.label(f) << ", "
		  << triangulation_.label(of) << ")\n";
      }
    }
    for (All_edges_iterator eit= triangulation_.all_edges_begin(); eit != triangulation_.all_edges_end();
	 ++eit) {
      Edge e=*eit;
      triangulation_.label(e);
    }
    return true;
  }



  bool audit_structure() const
  {
    if (triangulation_.dimension() != 3) return true;

    std::set<Point_key> pks;
    for (Finite_vertices_iterator eit = triangulation_.finite_vertices_begin();
	  eit != triangulation_.finite_vertices_end(); ++eit) {
      CGAL_assertion_code(Point_key k= eit->point());
      CGAL_assertion(pks.find(k) == pks.end());
      CGAL_assertion_code(pks.insert(k));
    }
    
    if (!has_certificates()) {
      for (All_edges_iterator eit = triangulation_.all_edges_begin();
	   eit != triangulation_.all_edges_end(); ++eit) {
	CGAL_assertion(!has_event(*eit));
      }
      for (All_facets_iterator eit = triangulation_.all_facets_begin();
	   eit != triangulation_.all_facets_end(); ++eit) {
	CGAL_assertion(!has_event(*eit));
      }
    } else {
      CGAL_KINETIC_LOG(LOG_SOME, "Auditing structure" << std::endl);
      //print();
      
      for (All_edges_iterator eit = triangulation_.all_edges_begin();
	   eit != triangulation_.all_edges_end(); ++eit) {
	bool isd3= is_degree_3(*eit);
	bool hd4= has_degree_4_vertex(*eit);
	if (!isd3 || hd4) {
	  if (has_event(*eit)) {
	    std::cerr << "Edge should not have certificate ";
	    triangulation_.write_labeled_edge(*eit, std::cerr);
	    std::cerr << std::endl;
	    simulator()->audit_event(triangulation_.label(*eit));
	    CGAL_assertion(0);
	  }
	} else if ( isd3) {
	  if (!has_event(*eit)) {
	    std::cerr << "Edge should have certificate ";
	    triangulation_.write_labeled_edge(*eit, std::cerr);
	    std::cerr << std::endl;
	    CGAL_assertion(0);
	  } else {
	    simulator()->audit_event(triangulation_.label(*eit));
	  }
	}
      }

      for (All_facets_iterator eit = triangulation_.all_facets_begin();
	   eit != triangulation_.all_facets_end(); ++eit) {
	bool hsd3= has_degree_3_edge(*eit);
	bool hd4= has_degree_4_vertex(*eit);
	if (hsd3 || hd4) {
	  if (has_event(*eit)) {
	    std::cerr << "Facet should not have certificate ";
	    triangulation_.write_labeled_facet(*eit, std::cerr);
	    std::cerr << std::endl;
	    simulator()->audit_event(triangulation_.label(*eit));
	    CGAL_assertion(0);
	  }
	} else {
	  if (!has_event(*eit)) {
	    std::cerr << "Facet should have certificate ";
	    triangulation_.write_labeled_facet(*eit, std::cerr);
	    std::cerr << std::endl;
	    CGAL_assertion(0);
	  } else {
	    simulator()->audit_event(triangulation_.label(*eit));
	  }
	}
      }
    }
    return true;
  }

  Event_key change_to_edge_flip(const Edge &e, Event_key k) {
    if (k== simulator()->null_event()) return k;
    Certificate s= extract_root_stack(k);
    return simulator()->set_event(k, Edge_flip_event(s, e, tr_.wrapper_handle()));
  }

  Event_key change_to_facet_flip(const Facet &f, Event_key k) {
    if (k== simulator()->null_event()) return k;
    Certificate s= extract_root_stack(k);
    return simulator()->set_event(k, Facet_flip_event(s, f, tr_.wrapper_handle()));
  }

  /*Moving_object_table* moving_object_table() {
    return mpt_.pointer();
    }*/

public:
  void clean_cell(Cell_handle h) {
    CGAL_precondition(has_certificates_);
    for (unsigned int i=0; i< 4; ++i) {
      for (unsigned int j=0; j<i; ++j) {
	Edge e(h, i, j);
	if (has_event(e)) {
	  make_not_edge_flip(e, h);
	  //delete_event(triangulation_.label(e));
	  //triangulation_.set_label(e, Event_key());
	}
      }
    }
    for (unsigned int i=0; i< 4; ++i) {
      Facet f(h, i);
      if (has_event(f)) {
	simulator()->delete_event(triangulation_.label(f));
	triangulation_.set_label(f, Event_key());
      }
    }
    v_.destroy_cell(h);
  }

  void handle_new_cell(Cell_handle h) {
    CGAL_precondition(has_certificates_);
    for (unsigned int i=0; i<4; ++i) {
      for (unsigned int j=0; j< i; ++j) {
	Edge e(h, i, j);
	if (is_degree_3(e) && !has_event(e) && !has_degree_4_vertex(e)) {
	  make_edge_flip(e);
	}
      }
      Facet f(h, i);
      if (!has_event(f) && !has_degree_3_edge(f)) {
	make_certificate(f);
      }
     
    }
    for (unsigned int i=0; i<4; ++i) {
      Vertex_handle vh= h->vertex(i);
      if (is_degree_4(vh)) {
	suppress_edge_flips(vh);
      } else {
	create_edge_flips(vh);
      }
    }

    v_.create_cell(h);
  }
protected:
  void handle_changed_cell(Cell_handle) {
    
  }

  Certificate root_stack(const Facet &f,
			 const typename Simulator::Time &st) const
  {
    std::vector<Point_key> ids;
    point_keys(f, back_inserter(ids));
#ifndef NDEBUG
    std::vector<Point_key> mids;
    Facet of= triangulation_.opposite(f);
    point_keys(of, back_inserter(mids));
#endif
#ifndef NDEBUG
    CGAL_KINETIC_LOG(LOG_LOTS, "Creating root_stack for points ");
    for (typename std::vector<Point_key>::const_iterator cit= ids.begin(); cit != ids.end(); ++cit) {
      CGAL_KINETIC_LOG(LOG_LOTS, *cit);
    }
    CGAL_KINETIC_LOG(LOG_LOTS, std::endl);
#endif
    bool is_const=true;
    for (unsigned int i=0; i<ids.size(); ++i) {
      if (!point(ids[i]).is_constant()) {
	is_const=false;
	break;
      }
    }
    if (is_const) {
      //typename Kinetic_kernel::Certificate_function cf(1.0);
      return Certificate(); //simulator()->root_stack_object(cf);
    }
    else if (ids.size()==4) {
      /*if (point(ids[0]).is_constant()){
      // hack
      std::swap(ids[0], ids[3]);
      std::swap(ids[1], ids[2]);
      }*/
      return o3_(point(ids[0]),
		 point(ids[1]),
		 point(ids[2]),
		 point(ids[3]),
		 st,
		 simulator()->end_time());
    }
    else {
      CGAL_assertion(ids.size()==5);
      /*if (point(ids[0]).is_constant()){
      // hack for linear case
      std::swap(ids[0], ids[4]);
      std::swap(ids[1], ids[2]);
      }*/
      return soc_(point(ids[0]),
		  point(ids[1]),
		  point(ids[2]),
		  point(ids[3]),
		  point(ids[4]),
		  st,
		  simulator()->end_time());
    }
    // Some compilers give warnings without this
    CGAL_postcondition(0);
    return Certificate();
    //return simulator()->root_stack_object(typename TraitsT::Simulator::Function_kernel::Function(0));
  }

  
  TraitsT tr_;
  Triangulation triangulation_;
  std::vector<Vertex_handle> vhs_;
  typename TraitsT::Side_of_oriented_sphere_3 soc_;
  typename TraitsT::Orientation_3 o3_;
  bool has_certificates_;
  Visitor v_;
};

template <class TraitsT, class Visitor>
inline std::ostream &operator<<(std::ostream &out,
				const Delaunay_triangulation_base_3<TraitsT, Visitor> &dt)
{
  dt.write(out);
  return out;
}


CGAL_KINETIC_END_INTERNAL_NAMESPACE
#endif