snc_io_parser.h

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

    }

    CGAL_NEF_TRACEN("  sface 2");
  
    SHalfedge_handle se2;
    SHalfloop_handle sl2;
    CGAL_forall_sface_cycles_of(fc,sf2) {
      if(fc.is_shalfedge()) {
	SHalfedge_handle se(fc);
	SHalfedge_around_sface_circulator ec(se),ee(se);
	CGAL_For_all(ec,ee) { 
	  CGAL_NEF_TRACEN("     " << ec->source()->point() << 
		 " | " << ec->circle().orthogonal_vector());
	  if(ml(ec, se2) == -1)
	    se2 = ec;
	}
      }
      else if(fc.is_shalfloop())
	sl2 = SHalfloop_handle(fc);
      else
	CGAL_assertion(fc.is_svertex());
    }
  
    CGAL_NEF_TRACEN("  sedge cycles existing? " << (se1 != SHalfedge_handle()) 
	   << " , " << (se2 != SHalfedge_handle()));

    if(se1 != SHalfedge_handle() && se2 == SHalfedge_handle())
      return true;
    if(se1 == SHalfedge_handle() && se2 != SHalfedge_handle())
      return false;  
  
    if(se1 == SHalfedge_handle() && se2 == SHalfedge_handle()) {
      Vector_3 vec1 = sl1->circle().orthogonal_vector();
      Vector_3 vec2 = sl2->circle().orthogonal_vector();
      CGAL_NEF_TRACEN("  sloops " << vec1 << " , " << vec2);
      if(vec1.x() != vec2.x())
	return vec1.x() < vec2.x();
      else if(vec1.y() != vec2.y())
	return vec1.y() < vec2.y();
      else if(vec1.z() != vec2.z())
	return vec1.z() < vec2.z();     
    }
    
    CGAL_assertion(se1 != SHalfedge_handle() && se2 != SHalfedge_handle());

    CGAL_NEF_TRACEN("  minimal sedge in sface 1:" << se1->source()->point() << 
	   " , " << se1->circle().orthogonal_vector());
    CGAL_NEF_TRACEN("  minimal sedge in sface 2:" << se2->source()->point() << 
	   " , " << se2->circle().orthogonal_vector());
    CGAL_NEF_TRACEN("result " << ml(se1,se2));
    switch(ml(se1, se2)) {
    case -1: return true;
    case  1: return false;
    }
    sort_sface_cycle_entries<T> SORTSFC(*this->sncp());
    return SORTSFC(*sf1->sface_cycles_begin(), *sf2->sface_cycles_begin());
  }
};

template <typename T>
class sort_volumes : public SNC_decorator<T> {
  
  typedef T SNC_structure;
  typedef CGAL::SNC_decorator<T>    Base;
  typedef typename T::Volume_handle Volume_handle;
  typedef typename T::SFace_handle  SFace_handle;
  
 public:
  sort_volumes(T& D) : Base(D) {}
  
  bool operator() (Volume_handle c1, Volume_handle c2) const {
    CGAL_NEF_TRACEN("sort volumes");
    SFace_handle sf1 = SFace_handle(c1->shells_begin()); 
    SFace_handle sf2 = SFace_handle(c2->shells_begin()); 

    sort_sfaces<T> SORT(*this->sncp());
    return SORT(sf1, sf2);
  }
};

template <typename T>
class sort_facet_cycle_entries : public T {
  
  typedef typename T::SNC_structure     SNC_structure;
  typedef typename T::SM_decorator      SM_decorator;
  typedef typename T::Object_handle     Object_handle;
  typedef typename T::SHalfedge_handle  SHalfedge_handle;
  typedef typename T::SHalfloop_handle  SHalfloop_handle;
  typedef typename T::SFace_handle      SFace_handle;
  typedef typename T::Point_3           Point_3;
  typedef typename T::Vector_3          Vector_3;
  
 public:
  sort_facet_cycle_entries(T D) : T(D) {}
  
  bool operator() (Object_handle o1, Object_handle o2) const {
    
    SHalfedge_handle se1, se2;
    SHalfloop_handle sl1, sl2;
    
    if(!CGAL::assign(se1,o1) && !CGAL::assign(sl1,o1))
      CGAL_assertion_msg(0,"wrong handle");
    
    if(!CGAL::assign(se2,o2) && !CGAL::assign(sl2,o2))
      CGAL_assertion_msg(0,"wrong handle");    
    
    if(se1 != SHalfedge_handle() && se2 != SHalfedge_handle()) {
      sort_vertices<SNC_structure> SORT(*this->sncp());
      return SORT(se1->source()->source(), se2->source()->source());
    }

    if(se1 != SHalfedge_handle())
      return true;
    if(se2 != SHalfedge_handle())
      return false;

    CGAL_assertion(sl1 != SHalfloop_handle() && 
			sl2 != SHalfloop_handle());

    SM_decorator SD(&*sl1->incident_sface()->center_vertex());
    Vector_3 vec1(sl1->circle().orthogonal_vector());
    Vector_3 vec2(sl2->circle().orthogonal_vector());
    //    CGAL_assertion(vec1 == vec2.antipode());
    if(vec1.x() != vec2.x())
      return vec1.x() < vec2.x();
    else if(vec1.y() != vec2.y())
      return vec1.y() < vec2.y();
    else
      return vec1.z() < vec2.z();          
  }
};

template <typename T>
class sort_shell_entries : public T {
  
  typedef typename T::Object_handle Object_handle;
  typedef typename T::Shell_entry_iterator  Shell_entry_iterator;
  typedef typename T::SFace_handle  SFace_handle;
  typedef typename T::Point_3       Point_3;
  
 public:
  sort_shell_entries(T D) : T(D) {}
  
  bool operator() (Object_handle o1, Object_handle o2) const {
    SFace_handle sf1, sf2;
    CGAL::assign(sf1, o1);
    CGAL::assign(sf2, o2);
    Point_3 p1(sf1->center_vertex()->point()), p2(sf2->center_vertex()->point());
    if(p1.x() != p2.x())
      return p1.x() < p2.x();
    else if(p1.y() != p2.y())
      return p1.y() < p2.y();
    return p1.z() < p2.z();
  }
};

template<typename T>
struct find_minimal_sface_of_shell : public SNC_decorator<T> {
  
  typedef T                               SNC_structure;
  typedef CGAL::SNC_decorator<T>          Base;
  typedef typename T::Vertex_handle       Vertex_handle;
  typedef typename T::Halfedge_handle     Halfedge_handle;
  typedef typename T::Halffacet_handle    Halffacet_handle;
  typedef typename T::SFace_handle        SFace_handle;
  typedef typename T::SHalfedge_handle    SHalfedge_handle;
  typedef typename T::SHalfloop_handle    SHalfloop_handle;
  typedef CGAL::Unique_hash_map<SFace_handle,bool> SFace_visited_hash;

  SFace_visited_hash& Done;
  SFace_handle sf_min;
  sort_sfaces<T> SORT;
  
  find_minimal_sface_of_shell(T& D, SFace_visited_hash& Vi) 
    : Base(D), Done(Vi), SORT(D) {}
  
  void visit(SFace_handle h) { 
    Done[h]=true;
    if(sf_min == SFace_handle())
      sf_min = h;
    else {
      if(SORT(h,sf_min))
	sf_min = h;
    }
  }
  
  void visit(Vertex_handle ) {}
  void visit(Halfedge_handle ) {}
  void visit(Halffacet_handle ) {}
  void visit(SHalfedge_handle ) {}
  void visit(SHalfloop_handle ) {}

  SFace_handle& minimal_sface() { return sf_min; }
};

class Homogeneous_tag;
class Cartesian_tag;
template<typename Tag, typename Kernel> class Geometry_io;

template<typename ET>
class Geometry_io<Cartesian_tag, CGAL::Lazy_kernel<CGAL::Simple_cartesian<ET> > > {
 public:
  template <typename EK, typename K> static
  typename EK::Point_3 
  read_point(std::istream& in) {
    typedef Fraction_traits<typename K::FT> FracTraits;
    typename FracTraits::Type hx, hy, hz, hw;
    typename FracTraits::Numerator_type num;
    typename FracTraits::Denominator_type denom(1);
    typename FracTraits::Compose composer;
    in >> num;
    hx = composer(num, denom);
    in >> num;
    hy = composer(num, denom);
    in >> num;
    hz = composer(num, denom);
    in >> num;
    hw = composer(num, denom);
    return typename EK::Point_3(hx,hy,hz,hw);
  }

  template <typename EK, typename K> static
  typename EK::Plane_3 read_plane(std::istream& in) {
    typedef Fraction_traits<typename K::FT> FracTraits;
    typename FracTraits::Type a, b, c, d;
    typename FracTraits::Numerator_type num;
    typename FracTraits::Denominator_type denom(1);
    typename FracTraits::Compose composer;
    in >> num;
    a = composer(num, denom);
    in >> num;
    b = composer(num, denom);
    in >> num;
    c = composer(num, denom);
    in >> num;
    d = composer(num, denom);
    return typename EK::Plane_3(a,b,c,d);
  }

  template <typename R> static
  void print_point(std::ostream& out, const CGAL::Point_3<R> p) {
    typedef typename CGAL::Simple_cartesian<ET> SC;
    typedef typename SC::Point_3 Exact_point;
    typedef Geometry_io<Cartesian_tag, SC> Gio;
    
    Exact_point ep(p.x().exact(), p.y().exact(), p.z().exact());
    Gio::print_point(out, ep);
  }  

  template <typename R> static
  void print_plane(std::ostream& out, const CGAL::Plane_3<R> p) {
    typedef typename CGAL::Simple_cartesian<ET> SC;
    typedef typename SC::Plane_3 Exact_plane;
    typedef Geometry_io<Cartesian_tag, SC> Gio;
    
    Exact_plane ep(p.a().exact(), p.b().exact(), 
		   p.c().exact(), p.d().exact());
    Gio::print_plane(out, ep);
  }  
};

template<typename Kernel>
class Geometry_io<Cartesian_tag, Kernel> {
 public:
  template <typename EK, typename K> static
  typename EK::Point_3 read_point(std::istream& in) {
    typedef Fraction_traits<typename K::FT> FracTraits;
    typename FracTraits::Type hx, hy, hz, hw;
    typename FracTraits::Numerator_type num;
    typename FracTraits::Denominator_type denom(1);
    typename FracTraits::Compose composer;
    in >> num;
    hx = composer(num, denom);
    in >> num;
    hy = composer(num, denom);
    in >> num;
    hz = composer(num, denom);
    in >> num;
    hw = composer(num, denom);
    return typename EK::Point_3(hx,hy,hz,hw);
  }

  template <typename EK, typename K> static
  typename EK::Plane_3 read_plane(std::istream& in) {
    typedef Fraction_traits<typename K::FT> FracTraits;
    typename FracTraits::Type a, b, c, d;
    typename FracTraits::Numerator_type num;
    typename FracTraits::Denominator_type denom(1);
    typename FracTraits::Compose composer;
    in >> num;
    a = composer(num, denom);
    in >> num;
    b = composer(num, denom);
    in >> num;
    c = composer(num, denom);
    in >> num;
    d = composer(num, denom);
    return typename EK::Plane_3(a,b,c,d);
  }

  template <typename R> static
  void print_point(std::ostream& out, const CGAL::Point_3<R> p) {
    typedef Fraction_traits<typename R::FT> FracTraits;
    typedef std::vector<typename FracTraits::Numerator_type> NV;
    typedef typename NV::iterator NV_iter;

    typename FracTraits::Numerator_type num;
    typename FracTraits::Denominator_type denom;
    typename FracTraits::Decompose decomposer;
    NV vec;

    decomposer(p.hx(),num,denom);
    vec.push_back(num);
    vec.push_back(denom);
    vec.push_back(denom);
    vec.push_back(denom);
    decomposer(p.hy(),num,denom);
    vec[0]*=denom;
    vec[1]*=num;
    vec[2]*=denom;
    vec[3]*=denom;
    decomposer(p.hz(),num,denom);
    vec[0]*=denom;
    vec[1]*=denom;
    vec[2]*=num;
    vec[3]*=denom;
    Normalizing<Homogeneous_tag>::
      normalized(vec.begin(),vec.end());
    out << vec[0] << " " << vec[1] << " "
	<< vec[2] << " " << vec[3];
  }

  template <typename R> static
  void print_plane(std::ostream& out, const CGAL::Plane_3<R> p) {
    typedef Fraction_traits<typename R::FT> FracTraits;
    typedef std::vector<typename FracTraits::Numerator_type> NV;
    typedef typename NV::iterator NV_iter;

    typename FracTraits::Numerator_type num;
    typename FracTraits::Denominator_type denom;
    typename FracTraits::Decompose decomposer;
    NV vec;

    decomposer(p.a(),num,denom);
    vec.push_back(num);
    vec.push_back(denom);
    vec.push_back(denom);
    vec.push_back(denom);
    decomposer(p.b(),num,denom);
    vec[0]*=denom;
    vec[1]*=num;
    vec[2]*=denom;
    vec[3]*=denom;
    decomposer(p.c(),num,denom);
    vec[0]*=denom;
    vec[1]*=denom;
    vec[2]*=num;
    vec[3]*=denom;
    decomposer(p.d(),num,denom);
    vec[0]*=denom;
    vec[1]*=denom;
    vec[2]*=denom;
    vec[3]*=num;
    Normalizing<Homogeneous_tag>::
      normalized(vec.begin(),vec.end());

    out << vec[0] << " " << vec[1] << " "
	<< vec[2] << " " << vec[3];
  }
};

template<typename Kernel>
class Geometry_io<Homogeneous_tag, Kernel> {
 public:
  template <typename EK, typename K> static
  typename EK::Point_3 read_point(std::istream& in) {
    typename K::RT hx, hy, hz, hw;
    in >> hx >> hy >> hz >> hw;
    return typename EK::Point_3(hx, hy, hz, hw);
  }

  template <typename EK, typename K> static
  typename EK::Plane_3 read_plane(std::istream& in) {
    typename K::RT a, b, c, d;
    in >> a >> b >> c >> d;
    return typename EK::Plane_3(a, b, c, d);
  }

  template <typename R> static
  void print_point(std::ostream& out, const CGAL::Point_3<R>& p) {
    out << p;
  }

  template <typename R> static
  void print_plane(std::ostream& out, const CGAL::Plane_3<R>& p) {