k3_tree.h

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

	  if( !f_mark[f]) {
	    O.push_back(*o);
	    f_mark[f] = true;
	  }
	}
#ifdef CGAL_NEF3_TRIANGULATE_FACETS
	else if(CGAL::assign(t, *o)) {
	  Triangle_3 tr = t.get_triangle();
	  if( !t_mark[tr]) {
	    O.push_back(*o);
	    t_mark[tr] = true;
	  }
	}
#endif
#ifdef CGAL_NEF3_FACET_WITH_BOX
	else if(CGAL::assign(pf, *o)) {
	  CGAL_assertion_msg(false, "wrong type");
	}
#endif
	else
	  CGAL_assertion_msg( 0, "wrong handle");
      }
    }
    return O;
  }

  bool is_point_on_cell( const Point_3& p, const typename Objects_around_segment::Iterator& target) const {
    return is_point_on_cell( p, target.get_node(), root);
  }

  void add_facet(Halffacet_handle f) {
    root->add_facet(f,0);
  }

  void add_edge(Halfedge_handle e) {
    root->add_edge(e,0);	
  }

  void add_vertex(Vertex_handle v) {
    root->add_vertex(v,0);
  }

  class BBox_updater {  
    Bounding_box_3 b;

  public:	
    BBox_updater() {}

    void pre_visit(const Node*) {}
    void post_visit(const Node* n) {
      typename Object_list::const_iterator o;
      for( o = n->objects().begin(); 
	   o != n->objects().end(); ++o) {
        Vertex_handle v;
        if( CGAL::assign( v, *o))
          b.extend(v->point());
      }
    }

    Bounding_box_3 box() const{
      return b;
    }
    
  };

  template <typename Visitor>
  void visit_k3tree(const Node* current, Visitor& V) const {
    V.pre_visit(current);
    if(current->left() != 0) {
      visit_k3tree(current->left(), V);
      visit_k3tree(current->right(), V);
    }
    V.post_visit(current);
  }

  size_t bytes() { return root->bytes();}
  size_t leafs(int mask = 255, int lower_limit=0) { return root->leafs(mask, lower_limit);}

  void transform(const Aff_transformation_3& t) {
    // TODO: Bounding box must be updated/transformed, too
    if(root != 0)
      root->transform(t);

    BBox_updater bbup;
    visit_k3tree(root, bbup);
    bounding_box = bbup.box();
  }

  
#ifdef CODE_DOES_NOT_WORK_WITH_BOTH_KERNELS_AT_THE_SAME_TIME
template <typename T>
friend std::ostream& operator<<
(std::ostream& os, const K3_tree<T>& k3_tree) {
  os << (const Node*)k3_tree.root; // no default conversion to const Node*?
  return os;
}
#endif
std::string dump_object_list( const Object_list& O, int level = 0) {
  std::stringstream os; 
  typename Object_list::size_type v_count = 0, e_count = 0, f_count = 0;
  typename Object_list::const_iterator o;
  Vertex_handle v;
  Halfedge_handle e;
  Halffacet_handle f;
#ifdef CGAL_NEF3_TRIANGULATE_FACETS
  typename Object_list::size_type t_count = 0;
  Halffacet_triangle_handle t;
#endif
#ifdef CGAL_NEF3_FACET_WITH_BOX
  typename Object_list::size_type p_count = 0;
  Partial_facet pf;
#endif
  for( o = O.begin(); o != O.end(); ++o) {
    if( CGAL::assign( v, *o)) {
      if( level) os << v->point() << std::endl;
      ++v_count;
    }
    else if( CGAL::assign( e, *o)) {
      if( level) os << e->source()->point() << "->"
		    << e->twin()->source()->point() << std::endl;
      ++e_count;
    }
    else if( CGAL::assign( f, *o)) {
      if( level) os << "facet" << std::endl;
      ++f_count;
    }
#ifdef CGAL_NEF3_TRIANGULATE_FACETS
    else if( CGAL::assign(t, *o)) {
      if( level) os << "triangle" << std::endl;
      ++t_count;
    }	
#endif
#ifdef CGAL_NEF3_FACET_WITH_BOX
    else if( CGAL::assign(pf, *o)) {
      if( level) pf.debug();
      ++p_count;      
    }
#endif
    else 
      CGAL_assertion_msg( 0, "wrong handle");
  }
  os << v_count << "v " << e_count << "e " << f_count << "f ";
#ifdef CGAL_NEF3_TRIANGULATE_FACETS
  os  << t_count << "t ";
#endif
#ifdef CGAL_NEF3_FACET_WITH_BOX
  os  << p_count << "p ";
#endif
  return os.str();
 }

bool update( Unique_hash_map<Vertex_handle, bool>& V, 
             Unique_hash_map<Halfedge_handle, bool>& E, 
             Unique_hash_map<Halffacet_handle, bool>& F) {
  return update( root, V, E, F);
}

bool update( Node* node,
             Unique_hash_map<Vertex_handle, bool>& V, 
             Unique_hash_map<Halfedge_handle, bool>& E, 
             Unique_hash_map<Halffacet_handle, bool>& F) {
  CGAL_assertion( node != 0);
  if( node->is_leaf()) {
    bool updated = false;
    Object_list* O = &node->object_list;
    typename Object_list::iterator onext, o = O->begin();
    while( o != O->end()) {
      onext = o;
      onext++;
      Vertex_handle v;
      Halfedge_handle e;
      Halffacet_handle f;
      if( CGAL::assign( v, *o)) {
        if( !V[v]) {
          O->erase(o);
          updated = true;
        }
      }
      else if( CGAL::assign( e, *o)) {
        if( !E[e]) {
          O->erase(o);
          updated = true;
        }         
      }
      else if( CGAL::assign( f, *o)) {
        if( !F[f]) {
          O->erase(o);
          updated = true;
        }
      }
      else CGAL_assertion_msg( 0, "wrong handle");
      o = onext;
    }
    return updated;
  }
  // TODO: protect the code below from optimizations!
  bool left_updated = update( node->left_node, V, E, F);
  CGAL_NEF_TRACEN("k3_tree::update(): left node updated? "<<left_updated);
  bool right_updated = update( node->right_node, V, E, F);
  CGAL_NEF_TRACEN("k3_tree::update(): right node updated? "<<right_updated);
  return (left_updated || right_updated);
}

~K3_tree() {
  CGAL_NEF_TRACEN("~K3_tree: deleting root...");
  delete root;
}

private:
  
template <typename Depth>
Node* build_kdtree(Object_list& O, Object_iterator v_end, 
	           Depth depth, Node* parent=0, int non_efective_splits=0) {
  CGAL_precondition( depth >= 0);
  CGAL_NEF_TRACEN( "build_kdtree: "<<O.size()<<" objects, "<<"depth "<<depth);
  CGAL_NEF_TRACEN( "build_kdtree: "<<dump_object_list(O,1));
  if( !can_set_be_divided(O.begin(), v_end, depth)) {
    CGAL_NEF_TRACEN("build_kdtree: set cannot be divided");
    return new Node( parent, 0, 0, Plane_3(), O);
  }
  Object_iterator median; 
  Plane_3 partition_plane = construct_splitting_plane(O.begin(), v_end, median, depth);
  CGAL_NEF_TRACEN("build_kdtree: plane: "<<partition_plane<< " " << partition_plane.point());

  Object_list O1, O2;
  Vertex_handle vm,vx;
  CGAL::assign(vm,*median);
  Smaller_ smaller(depth%3);
  for(Object_iterator oi=O.begin();oi!=median;++oi) {
    O1.push_back(*oi);
    CGAL::assign(vx,*oi);
    if(!smaller(vx, vm))
      O2.push_back(*oi);
  }

  O1.push_back(*median);
  O2.push_back(*median);
  
  for(Object_iterator oi=median+1;oi!=v_end;++oi) {
    O2.push_back(*oi);
    CGAL::assign(vx,*oi);
    if(!smaller(vm, vx))
      O1.push_back(*oi);
  }

#ifdef CGAL_NEF_EXPLOIT_REFERENCE_COUNTING
  Side_of_plane sop(reference_counted);
#else
  Side_of_plane sop;
#endif
  typename Object_list::size_type v_end1 = O1.size();
  typename Object_list::size_type v_end2 = O2.size();
  bool splitted = classify_objects( v_end, O.end(), partition_plane, sop,
                                    std::back_inserter(O1), 
                                    std::back_inserter(O2), depth);

  bool non_efective_split = false;
  if( !splitted) {
    CGAL_NEF_TRACEN("build_kdtree: splitting plane not found");
    //    if(depth > max_depth)
    return new Node( parent, 0, 0, Plane_3(), O);
    non_efective_split = true;
  } else {
    CGAL_NEF_TRACEN("Sizes " << O1.size() << ", " << O2.size() << ", " << O.size());
    CGAL_assertion( O1.size() <= O.size() && O2.size() <= O.size());
    CGAL_assertion( O1.size() + O2.size() >= O.size());
    non_efective_split = ((O1.size() == O.size()) || (O2.size() == O.size()));
  }
  if( non_efective_split)
    non_efective_splits++;
  else
    non_efective_splits = 0;
  if( non_efective_splits > 2) {
    CGAL_NEF_TRACEN("build_kdtree: non efective splits reached maximum");
    return new Node( parent, 0, 0, Plane_3(), O);
  }
  Node *node = new Node( parent, 0, 0, partition_plane, Object_list());
  node->left_node = build_kdtree( O1, O1.begin()+v_end1, depth + 1, node, non_efective_splits);
  node->right_node = build_kdtree( O2, O2.begin()+v_end2, depth + 1, node, non_efective_splits);
  return node;
}

template <typename Depth>
bool can_set_be_divided(Object_iterator start, Object_iterator end, Depth depth) {
  if(depth >= max_depth)
    return false;
  if(std::distance(start,end)<2)
    return false;
  return true;
}

template <typename OutputIterator, typename Depth>
bool classify_objects(Object_iterator start, Object_iterator end, 
                      Plane_3 partition_plane, Side_of_plane& sop,
                      OutputIterator o1, OutputIterator o2, Depth depth) {
  typename Object_list::difference_type on_oriented_boundary = 0;
  typename Object_list::const_iterator o;
  
  Point_3 point_on_plane(partition_plane.point());

  for( o = start; o != end; ++o) {
#ifdef CGAL_NEF3_FACET_WITH_BOX
    Partial_facet pf;
    if(CGAL::assign(pf, *o)) {
      Partial_facet pfn,pfp;
      if(pf.divide(partition_plane, pfn, pfp)) {
	*o1 = Object_handle(pfn);
	++o1;
	*o2 = Object_handle(pfp);
	++o2;
	continue;
      }
    }
#endif
    Oriented_side side = sop( point_on_plane, *o, depth);
    if( side == ON_NEGATIVE_SIDE || side == ON_ORIENTED_BOUNDARY) {
      *o1 = *o;
      ++o1;
    }
    if( side == ON_POSITIVE_SIDE || side == ON_ORIENTED_BOUNDARY) {
      *o2 = *o;
      ++o2;
    }
    if( side == ON_ORIENTED_BOUNDARY)
      ++on_oriented_boundary;
  }
  return (on_oriented_boundary != std::distance(start,end));
}


template <typename Depth>
Plane_3 construct_splitting_plane(Object_iterator start, Object_iterator end,
                                  Object_iterator& median, Depth depth) {
  CGAL_precondition( depth >= 0);
  typename Object_list::difference_type n=std::distance(start,end);
  CGAL_assertion(n>1);

  std::nth_element(start, start+n/2, end,
  	           Smaller_(depth%3));

  Vertex_handle v;
  median = start+n/2;
  CGAL::assign(v,*median);
  switch( depth % 3) {
  case 0: return Plane_3( v->point(), Vector_3( 1, 0, 0)); break;
  case 1: return Plane_3( v->point(), Vector_3( 0, 1, 0)); break;
  case 2: return Plane_3( v->point(), Vector_3( 0, 0, 1)); break;
  }

  CGAL_assertion_msg( 0, "never reached");
  return Plane_3();
}

const Node *locate_cell_containing( const Point_3& p, const Node* node) const {
  CGAL_precondition( node != 0);
  if( node->is_leaf())
    return node;
  else {
    Oriented_side side = node->plane().oriented_side(p);
    if( side == ON_NEGATIVE_SIDE || side == ON_ORIENTED_BOUNDARY)
      return locate_cell_containing( p, node->left());
    else { // side == ON_POSITIVE_SIDE 
      CGAL_assertion( side == ON_POSITIVE_SIDE);
      return locate_cell_containing( p, node->right());
    }
  }
}

const Object_list& locate( const Point_3& p, const Node* node) const {
  CGAL_precondition( node != 0);
  return locate_cell_containing( p, node)->objects();
}

bool is_point_on_cell( const Point_3& p, const Node* target, const Node* current) const {
  CGAL_precondition( target != 0 && current != 0);
  if( current->is_leaf())
    return (current == target); 
  Oriented_side side = current->plane().oriented_side(p);
  if( side == ON_NEGATIVE_SIDE)
    return is_point_on_cell( p, target, current->left());
  else if( side == ON_POSITIVE_SIDE)
    return is_point_on_cell( p, target, current->right());
  CGAL_assertion( side == ON_ORIENTED_BOUNDARY);
  return (is_point_on_cell( p, target, current->left()) ||
          is_point_on_cell( p, target, current->right()));
}

};

CGAL_END_NAMESPACE

#endif // CGAL_NEF_K3_TREE_H