k3_tree.h

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

// Copyright (c) 1997-2000  Max-Planck-Institute Saarbruecken (Germany).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you may redistribute it under
// the terms of the Q Public License version 1.0.
// See the file LICENSE.QPL 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/Nef_3/include/CGAL/Nef_3/K3_tree.h $
// $Id: K3_tree.h 38349 2007-04-19 17:32:41Z hachenb $
// 
//
// Author(s)     : Miguel Granados <[e-mail]>
//                 Peter Hachenberger <[e-mail]>

#ifndef CGAL_NEF_K3_TREE_H
#define CGAL_NEF_K3_TREE_H

#include <CGAL/basic.h>
#include <CGAL/Unique_hash_map.h>
#include <CGAL/Nef_3/quotient_coordinates_to_homogeneous_point.h>
#include <CGAL/Lazy_kernel.h>
#include <CGAL/Cartesian.h>

#ifdef CGAL_NEF3_TRIANGULATE_FACETS
#include <CGAL/Constrained_triangulation_2.h>
#include <CGAL/Triangulation_data_structure_2.h>
#include <CGAL/Triangulation_euclidean_traits_xy_3.h>
#include <CGAL/Triangulation_euclidean_traits_yz_3.h>
#include <CGAL/Triangulation_euclidean_traits_xz_3.h>
#include <CGAL/Constrained_triangulation_face_base_2.h>
#endif

#include <deque>
#include <sstream>
#include <string>
#include <map>

#undef CGAL_NEF_DEBUG
#define CGAL_NEF_DEBUG 503
#include <CGAL/Nef_2/debug.h>

CGAL_BEGIN_NAMESPACE

template <typename Triangle_3>
void sort_triangle_by_lexicographically_smaller_vertex
(const Triangle_3& t, int& a, int& b, int& c) {
  typedef typename Triangle_3::R Kernel;
  a = 0;
  for( int i = 1; i < 3; ++i) {
    if( compare_xyz<Kernel>( t[a], t[i]) == SMALLER)
      a = i;
  }
  b = (a + 1) % 3;
  c = (b + 1) % 3;
  if( compare_xyz<Kernel>( t[b], t[c]) == LARGER)
    std::swap( b, c);
  return;
}

template <typename Triangle_3>
struct Compare_triangle_3 {
  typedef typename Triangle_3::R Kernel;
  bool operator()( const Triangle_3& t1, const Triangle_3& t2) const {
    int v1[3], v2[3];
    sort_triangle_by_lexicographically_smaller_vertex
      ( t1, v1[0], v1[1], v1[2]);
    sort_triangle_by_lexicographically_smaller_vertex
      ( t2, v2[0], v2[1], v2[2]);
    for( int i = 0; i < 3; ++i) {
      Comparison_result c = compare_xyz<Kernel>( t1[v1[i]], t2[v2[i]]);
      if( c == SMALLER)
	return true;
      else if( c == LARGER)
	return false;
    }
    return false; // the two triangles are equivalent
  }
};

template <class Traits>
class K3_tree
{

template <typename Kernel, typename Object, 
          typename Vertex, typename Coordinate>
class Smaller_than
{
public:
  Smaller_than(Coordinate c) : coord(c) {
    CGAL_assertion( c >= 0 && c <=2);
  }
  bool operator()( const Vertex& v1, const Vertex& v2) {
    switch(coord) {
    case 0: return CGAL::compare_x(v1->point(), v2->point()) == SMALLER;
    case 1: return CGAL::compare_y(v1->point(), v2->point()) == SMALLER;
    case 2: return CGAL::compare_z(v1->point(), v2->point()) == SMALLER;
    default: CGAL_assertion(false);
    }
    return false;
  }

  bool operator()( const Object& o1, const Object& o2) {
    Vertex v1,v2;
    CGAL::assign(v1,o1);
    CGAL::assign(v2,o2);
    switch(coord) {
    case 0: return CGAL::compare_x(v1->point(), v2->point()) == SMALLER;
    case 1: return CGAL::compare_y(v1->point(), v2->point()) == SMALLER;
    case 2: return CGAL::compare_z(v1->point(), v2->point()) == SMALLER;
    default: CGAL_assertion(false);
    }
    return false;
  }
private:
  Coordinate coord;
};


template <typename Object, typename Vertex, 
          typename Coordinate, typename EK>
  class Smaller_than<CGAL::Lazy_kernel<EK>, Object, Vertex, Coordinate>
{
public:
  Smaller_than(Coordinate c) : coord(c) {
    CGAL_assertion( c >= 0 && c <=2);
  }
  bool operator()( const Vertex& v1, const Vertex& v2) {
    switch(coord) {
    case 0: return CGAL::to_interval(v1->point().x()).second <
			   CGAL::to_interval(v2->point().x()).first;
    case 1: return CGAL::to_interval(v1->point().y()).second <
			   CGAL::to_interval(v2->point().y()).first;
    case 2: return CGAL::to_interval(v1->point().z()).second <
			   CGAL::to_interval(v2->point().z()).first;
    default: CGAL_assertion(false);
    }
    return false;
  }

  bool operator()( const Object& o1, const Object& o2) {
    Vertex v1,v2;
    CGAL::assign(v1,o1);
    CGAL::assign(v2,o2);
    switch(coord) {
    case 0: return CGAL::to_interval(v1->point().x()).second <
			   CGAL::to_interval(v2->point().x()).first;
    case 1: return CGAL::to_interval(v1->point().y()).second <
			   CGAL::to_interval(v2->point().y()).first;
    case 2: return CGAL::to_interval(v1->point().z()).second <
			   CGAL::to_interval(v2->point().z()).first;
    default: CGAL_assertion(false);
    }
    return false;
  }
private:
  Coordinate coord;
};

#ifdef CGAL_NEF3_TRIANGULATE_FACETS
  template<typename SNC_structure, typename Kernel>
  class Triangulation_handler {

    typedef typename CGAL::Triangulation_vertex_base_2<Kernel>               Vb;
    typedef typename CGAL::Constrained_triangulation_face_base_2<Kernel>     Fb;
    typedef typename CGAL::Triangulation_data_structure_2<Vb,Fb>             TDS;
    typedef typename CGAL::No_intersection_tag                               Itag;
    typedef typename CGAL::Constrained_triangulation_2<Kernel,TDS,Itag>      CT;

    typedef typename CT::Face_handle           Face_handle;
    typedef typename CT::Finite_faces_iterator Finite_face_iterator;
    typedef typename CT::Edge                  Edge;

    typedef typename SNC_structure::Halffacet_cycle_iterator 
                                    Halffacet_cycle_iterator;
    typedef typename SNC_structure::SHalfedge_around_facet_circulator
                                    SHalfedge_around_facet_circulator;

    CT ct;
    CGAL::Unique_hash_map<Face_handle, bool> visited;
    Finite_face_iterator fi;

  public:
    template<typename Halffacet_handle>
    Triangulation_handler(Halffacet_handle f) : visited(false) {

      typedef typename SNC_structure::Halffacet_cycle_iterator 
                                      Halffacet_cycle_iterator;
      typedef typename SNC_structure::SHalfedge_around_facet_circulator
                                      SHalfedge_around_facet_circulator;

      Halffacet_cycle_iterator fci;
      for(fci=f->facet_cycles_begin(); fci!=f->facet_cycles_end(); ++fci) {
	if(fci.is_shalfedge()) {
          SHalfedge_around_facet_circulator sfc(fci), send(sfc);
	  CGAL_For_all(sfc,send) {
	    ct.insert_constraint(sfc->source()->source()->point(),
	                         sfc->source()->twin()->source()->point());
          }
        }
      }
      CGAL_assertion(ct.is_valid());

      typename CT::Face_handle infinite = ct.infinite_face();
      typename CT::Vertex_handle ctv = infinite->vertex(1);
      if(ct.is_infinite(ctv)) ctv = infinite->vertex(2);
      CGAL_assertion(!ct.is_infinite(ctv));

      typename CT::Face_handle opposite;
      typename CT::Face_circulator vc(ctv,infinite);
      do { opposite = vc++;
      } while(!ct.is_constrained(CT::Edge(vc,vc->index(opposite))));
      typename CT::Face_handle first = vc;

      traverse_triangulation(first, first->index(opposite));
      
      fi = ct.finite_faces_begin();
    }

    void traverse_triangulation(Face_handle f, int parent) {
      visited[f] = true;
      if(!ct.is_constrained(Edge(f,ct.cw(parent))) && !visited[f->neighbor(ct.cw(parent))]) {
	Face_handle child(f->neighbor(ct.cw(parent)));
	traverse_triangulation(child, child->index(f));
      } 
      if(!ct.is_constrained(Edge(f,ct.cw(parent))) && !visited[f->neighbor(ct.cw(parent))]) {
	Face_handle child(f->neighbor(ct.cw(parent)));
	traverse_triangulation(child, child->index(f));
      } 
    } 
 
    template<typename Triangle_3>
    bool get_next_triangle(Triangle_3& tr) {
      if(fi == ct.finite_faces_end()) return false;
      ++fi;
      while(fi != ct.finite_faces_end() && visited[fi] == false) ++fi;
      if(fi == ct.finite_faces_end()) return false;
      tr = Triangle_3(fi->vertex(0)->point(), fi->vertex(1)->point(), fi->vertex(2)->point());
      return true;
    }
  };
#endif

public:
  friend class Objects_along_ray;
  friend class Objects_around_segment;

public:
  
typedef typename Traits::SNC_decorator SNC_decorator;
typedef typename Traits::Infimaximal_box Infimaximal_box;
typedef typename Traits::Vertex_handle Vertex_handle;
typedef typename Traits::Halfedge_handle Halfedge_handle;
typedef typename Traits::Halffacet_handle Halffacet_handle;
#ifdef CGAL_NEF3_TRIANGULATE_FACETS
typedef typename Traits::Halffacet_triangle_handle Halffacet_triangle_handle;
#endif
#ifdef CGAL_NEF3_FACET_WITH_BOX
typedef typename Traits::Partial_facet Partial_facet;
#endif
typedef typename Traits::Object_handle Object_handle;
typedef std::vector<Object_handle> Object_list;
typedef typename Object_list::const_iterator Object_const_iterator;
typedef typename Object_list::iterator Object_iterator;
typedef typename Object_list::size_type size_type;

typedef typename Traits::Point_3 Point_3;
typedef typename Traits::Segment_3 Segment_3;
typedef typename Traits::Ray_3 Ray_3;
typedef typename Traits::Vector_3 Vector_3;
typedef typename Traits::Plane_3 Plane_3;
typedef typename Traits::Triangle_3 Triangle_3;
typedef typename Traits::Aff_transformation_3 Aff_transformation_3;

typedef typename Traits::Bounding_box_3 Bounding_box_3;
typedef typename Traits::Side_of_plane Side_of_plane;
typedef typename Traits::Objects_bbox Objects_bbox;

typedef typename Traits::Kernel Kernel;
typedef typename Kernel::RT RT;
typedef typename Kernel::FT FT;

typedef Smaller_than<
  Kernel,
  Object_handle,
  Vertex_handle, 
  int> Smaller_;

class Node {
  friend class K3_tree<Traits>;
public:
  Node( Node* p, Node* l, Node* r, Plane_3 pl, const Object_list& O) : 
    parent_node(p), left_node(l), right_node(r), splitting_plane(pl), 
	object_list(O) {
    if(l == 0)
      point_on_plane = Point_3();
    else
      point_on_plane = pl.point();
  }
  bool is_leaf() const {
    CGAL_assertion( (left_node != 0 && right_node != 0) ||
                    (left_node == 0 && right_node == 0));
    return (left_node == 0 && right_node == 0);
  }
  const Node* parent() const { return parent_node; }
  const Node* left() const { return left_node; }
  const Node* right() const { return right_node; }
  const Plane_3& plane() const { return splitting_plane; }
  const Object_list& objects() const { return object_list; }

  void transform(const Aff_transformation_3& t) {
    if(left_node != 0) {
	CGAL_assertion(right_node != 0);
	left_node->transform(t);
 	right_node->transform(t);
  	splitting_plane = splitting_plane.transform(t);
#ifdef CGAL_NEF3_TRIANGULATE_FACETS 
    } else {
      Halffacet_triangle_handle tri;
      typename Object_list::iterator o;
      for(o = object_list.begin(); o != object_list.end(); ++o)
        if(assign(tri,*o)) {
          tri.transform(t);
	  *o = Object_handle(tri);
        }
#endif // CGAL_NEF3_TRIANGULATE_FACETS 
    }
  }

  std::size_t bytes() {
    // bytes used for the Kd-tree
    std::size_t s = sizeof(Node);
    if(left_node != 0)
      s += left_node->bytes();
    if(right_node != 0)
      s += right_node->bytes();
    typename Object_list::iterator o;
    for(o = object_list.begin(); o != object_list.end(); ++o)
      s += sizeof(*o);
    return s;
  }

  std::size_t leafs(int mask = 255, int lower_limit = 0) {
    std::size_t s = 0;
    Halffacet_handle f;    
    Halfedge_handle e;
    Vertex_handle v;
    typename Object_list::iterator o;
    if(mask == 0) 
      s = 1;
    else {
      for(o = object_list.begin(); o != object_list.end(); ++o) {
	if((mask & 1) && assign(v,*o))
	  ++s;
	else if((mask&2) && assign(e,*o))
	  ++s;
	else if(((mask&4) || (mask&8)) && assign(f,*o)) {
	  if(mask&4)
	    ++s;
	  else {
	    int length = 0;	
	    typename Traits::SHalfedge_around_facet_circulator safc(f->facet_cycles_begin()), 
	      send(safc);
	    while(++length < lower_limit && ++safc != send);
	    if(length >= lower_limit)
	      ++s;	
	  }
	}
      }
    }
    
    if(left_node != 0)
      s += left_node->leafs(mask, lower_limit);
    if(right_node != 0)
      s += right_node->leafs(mask, lower_limit);
    return s;
  }

  template<typename Depth>
  void add_facet(Halffacet_handle f, Depth depth) {
    if(left_node == 0) {
      object_list.push_back(Object_handle(f));
      return;
    }

    Side_of_plane sop;
    Oriented_side side = sop(splitting_plane.point(), f, depth);
    if( side == ON_NEGATIVE_SIDE || side == ON_ORIENTED_BOUNDARY)
      left_node->add_facet(f, depth+1);
    if( side == ON_POSITIVE_SIDE || side == ON_ORIENTED_BOUNDARY)
      right_node->add_facet(f, depth+1);