snc_k3_tree_traits.h

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

  typedef typename Kernel::RT RT;
  typedef typename Kernel::Kernel_tag Kernel_tag;
  //  typedef CGAL::Bounding_box_3
  //    <typename Is_extended_kernel<Kernel>::value_type, Kernel> 
  //    Bounding_box_3;
  typedef CGAL::Bounding_box_3<Tag_true, Kernel> 
    Bounding_box_3;

  typedef typename Kernel::Intersect_3 Intersect;
  typedef CGAL::Objects_bbox<SNC_decorator> Objects_bbox;
  typedef CGAL::Side_of_plane<SNC_decorator> Side_of_plane;

  Intersect intersect_object() const {
    return Intersect();
  }

  Objects_bbox objects_bbox_object() const {
    return Objects_bbox();
  }
};

template <class SNC_decorator>
template <typename Depth>
Oriented_side 
Side_of_plane<SNC_decorator>::operator()
  (const Point_3& pop, Object_handle o, Depth depth) {
  Vertex_handle v;
  Halfedge_handle e;
  Halffacet_handle f;
#ifdef CGAL_NEF3_TRIANGULATE_FACETS
  Halffacet_triangle_handle t;
#endif
#ifdef CGAL_NEF3_FACET_WITH_BOX
  Partial_facet pf;
#endif
  if( CGAL::assign( v, o))
    return (*this)(pop, v, depth);
  else if( CGAL::assign( e, o))
    return (*this)(pop, e, depth);
  else if( CGAL::assign( f, o))
    return (*this)(pop, f, depth);
#ifdef CGAL_NEF3_FACET_WITH_BOX
  else if( CGAL::assign(pf, o))
    return (*this)(pop, pf, depth);
#endif
#ifdef CGAL_NEF3_TRIANGULATE_FACETS
  else if( CGAL::assign( t, o))
    return (*this)(pop, t, depth);
#endif
  else
    CGAL_assertion_msg( 0, "wrong handle");
  
  return Oriented_side(); // never reached
}

template <class SNC_decorator>
template <typename Depth>
Oriented_side 
Side_of_plane<SNC_decorator>::operator()
( const Point_3& pop, Vertex_handle v, Depth depth) {
  Comparison_result cr;
#ifdef CGAL_NEF_EXPLOIT_REFERENCE_COUNTING
  if(reference_counted) {
    if(!OnSideMapRC.is_defined(&(v->point().hw())))
      switch(depth%3) {
      case 0: 
        cr = CGAL::compare_x(v->point(), pop);
        OnSideMapRC[&(v->point().hw())] = cr == LARGER ? ON_POSITIVE_SIDE :
                         cr == SMALLER ? ON_NEGATIVE_SIDE : ON_ORIENTED_BOUNDARY;
        break;
      case 1:
        cr = CGAL::compare_y(v->point(), pop);
        OnSideMapRC[&(v->point().hw())] = cr == LARGER ? ON_POSITIVE_SIDE :
                         cr == SMALLER ? ON_NEGATIVE_SIDE : ON_ORIENTED_BOUNDARY;
        break;
      case 2:
        cr = CGAL::compare_z(v->point(), pop);
        OnSideMapRC[&(v->point().hw())] = cr == LARGER ? ON_POSITIVE_SIDE :
                         cr == SMALLER ? ON_NEGATIVE_SIDE : ON_ORIENTED_BOUNDARY;
        break;
      default: CGAL_assertion_msg(false, "wrong value");
      }
    return OnSideMapRC[&(v->point().hw())];
  } else {
#endif
    ComparePoints_ compare(depth%3);
    if(!OnSideMap.is_defined(v)) {
      cr = compare(v->point(), pop);
      OnSideMap[v] = cr == LARGER ? ON_POSITIVE_SIDE :
	cr == SMALLER ? ON_NEGATIVE_SIDE : ON_ORIENTED_BOUNDARY;
    }
    return OnSideMap[v];
#ifdef CGAL_NEF_EXPLOIT_REFERENCE_COUNTING
  }
#endif
  CGAL_assertion_msg(false, "should not be reached");
}

/* 
   An edge is considered intersecting a plane if its endpoints lie on the
   plane or if they lie on diferent sides.  Partial tangency is not considered
   as intersection, due the fact that a lower dimensional face (the vertex)
   should be already reported as an object intersecting the plane.
 */

template <class SNC_decorator>
template <typename Depth>
Oriented_side 
Side_of_plane<SNC_decorator>::operator()
( const Point_3& pop, Halfedge_handle e, Depth depth) {
  Vertex_handle v = e->source();
  Vertex_handle vt = e->twin()->source();
  /*
  Comparison_result cr;
  if(!OnSideMap.is_defined(v))
    switch(depth%3) {
    case 0: 
      cr = CGAL::compare_x(v->point(), pop);
      OnSideMap[v] = cr == LARGER ? ON_POSITIVE_SIDE :
                       cr == SMALLER ? ON_NEGATIVE_SIDE : ON_ORIENTED_BOUNDARY;
      break;
    case 1:
      cr = CGAL::compare_y(v->point(), pop);
      OnSideMap[v] = cr == LARGER ? ON_POSITIVE_SIDE :
                       cr == SMALLER ? ON_NEGATIVE_SIDE : ON_ORIENTED_BOUNDARY;
      break;
    case 2:
      cr = CGAL::compare_z(v->point(), pop);
      OnSideMap[v] = cr == LARGER ? ON_POSITIVE_SIDE :
                       cr == SMALLER ? ON_NEGATIVE_SIDE : ON_ORIENTED_BOUNDARY;
      break;
    default: CGAL_assertion_msg(false, "wrong value");
    }
  if(!OnSideMap.is_defined(vt))
    switch(depth%3) {
    case 0: 
      cr = CGAL::compare_x(vt->point(), pop);
      OnSideMap[vt] = cr == LARGER ? ON_POSITIVE_SIDE :
                       cr == SMALLER ? ON_NEGATIVE_SIDE : ON_ORIENTED_BOUNDARY;
      break;
    case 1:
      cr = CGAL::compare_y(vt->point(), pop);
      OnSideMap[vt] = cr == LARGER ? ON_POSITIVE_SIDE :
                       cr == SMALLER ? ON_NEGATIVE_SIDE : ON_ORIENTED_BOUNDARY;
      break;
    case 2:
      cr = CGAL::compare_z(vt->point(), pop);
      OnSideMap[vt] = cr == LARGER ? ON_POSITIVE_SIDE :
                       cr == SMALLER ? ON_NEGATIVE_SIDE : ON_ORIENTED_BOUNDARY;
      break;
    default: CGAL_assertion_msg(false, "wrong value");
    }
  Oriented_side src_side = OnSideMap[v];
  Oriented_side tgt_side = OnSideMap[vt];
*/
  Oriented_side src_side = (*this) (pop, v, depth);
  Oriented_side tgt_side = (*this) (pop, vt, depth);
  if( src_side == tgt_side)
    return src_side;
  if( src_side == ON_ORIENTED_BOUNDARY)
    return tgt_side;
  if( tgt_side == ON_ORIENTED_BOUNDARY)
    return src_side;
  return ON_ORIENTED_BOUNDARY;
}

#ifdef CGAL_NEF3_TRIANGULATE_FACETS
template <typename SNC_decorator>
template <typename Depth>
Oriented_side
Side_of_plane<SNC_decorator>::operator()
( const Point_3& pop, Halffacet_triangle_handle t, Depth depth) {
  bool on_positive_side = false, on_negative_side = false;
  Triangle_3 tr(t.get_triangle());
  for( int i = 0; i < 3; ++i) {
    Oriented_side side = ON_ORIENTED_BOUNDARY;
    Comparison_result cr;
    if(
#ifdef CGAL_NEF_EXPLOIT_REFERENCE_COUNTING
       !reference_counted || 
#endif
       !OnSideMapRC.is_defined(&(tr[i].hw()))) {
      switch(depth%3) {
      case 0: 
        cr = CGAL::compare_x(tr[i], pop);
        side = cr == LARGER ? ON_POSITIVE_SIDE :
                 cr == SMALLER ? ON_NEGATIVE_SIDE : ON_ORIENTED_BOUNDARY;
        break;
      case 1:
        cr = CGAL::compare_y(tr[i], pop);
        side = cr == LARGER ? ON_POSITIVE_SIDE :
                 cr == SMALLER ? ON_NEGATIVE_SIDE : ON_ORIENTED_BOUNDARY;
        break;
      case 2:
        cr = CGAL::compare_z(tr[i], pop);
        side = cr == LARGER ? ON_POSITIVE_SIDE :
                 cr == SMALLER ? ON_NEGATIVE_SIDE : ON_ORIENTED_BOUNDARY;
        break;
      default: CGAL_assertion_msg(false, "wrong value");
      }
#ifdef CGAL_NEF_EXPLOIT_REFERENCE_COUNTING
      if(reference_counted) 
        OnSideMapRC[&(tr[i].hw())] = side;
    } else if(reference_counted)
      side = OnSideMapRC[&(tr[i].hw())];
#endif
    if( side == ON_POSITIVE_SIDE)
      on_positive_side = true;
    else if( side == ON_NEGATIVE_SIDE)
      on_negative_side = true;
  }
  if( on_positive_side && on_negative_side)
    return ON_ORIENTED_BOUNDARY;
  if( !on_positive_side && !on_negative_side)
    return ON_ORIENTED_BOUNDARY;
  if( on_positive_side) {
    CGAL_assertion( !on_negative_side);
    return ON_POSITIVE_SIDE;
  }
  CGAL_assertion( on_negative_side);
  return ON_NEGATIVE_SIDE;
}
#endif

/* 
   As for the edges, if a facet is tanget to the plane it is not considered as
   a interesection since lower dimensional faces, like the edges and vertices
   where the tangency occurrs, should be reported as the objects intersecting 
   the plane.
   So, an intersection is reported if all vertices of the facet lie on plane,
   for which it is only necessary to check three vertices, or if the facet 
   has vertices on both sides of the plane, so the intersection is known
   as far as two vertices located on different sides of the plane.
*/

#ifdef CGAL_NEF3_FACET_WITH_BOX
template <class SNC_decorator>
template <typename Depth>
Oriented_side
Side_of_plane<SNC_decorator>::operator()
  (const Point_3& pop, Partial_facet& pf, Depth depth) {
  CGAL_assertion_msg(false, "not implemented yet");
  
  return ON_ORIENTED_BOUNDARY;
}
#endif

template <class SNC_decorator>
template <typename Depth>
Oriented_side 
Side_of_plane<SNC_decorator>::operator()
  (const Point_3& pop, Halffacet_handle f, Depth depth) {
    CGAL_assertion( std::distance( f->facet_cycles_begin(), f->facet_cycles_end()) > 0);
    /*
#ifdef CGAL_NEF3_FACET_WITH_BOX
    switch(depth%3) {
    case 0:
      if(f->b.min_coord(0) > pop.x())
	return ON_POSITIVE_SIDE;
      if(f->b.max_coord(0) < pop.x())
	return ON_NEGATIVE_SIDE;
      break;
    case 1:
      if(f->b.min_coord(1) > pop.y())
	return ON_POSITIVE_SIDE;
      if(f->b.max_coord(1) < pop.y())
	return ON_NEGATIVE_SIDE;
      break;
    case 2:
      if(f->b.min_coord(2) > pop.z())
	return ON_POSITIVE_SIDE;
      if(f->b.max_coord(2) < pop.z())
	return ON_NEGATIVE_SIDE;
      break;
    default: CGAL_assertion_msg(false, "wrong value");
    }    
    return ON_ORIENTED_BOUNDARY;
#else  
    */
  Halffacet_cycle_iterator fc(f->facet_cycles_begin());
  SHalfedge_handle e;
  CGAL_assertion(fc.is_shalfedge());
  e = SHalfedge_handle(fc);
  SHalfedge_around_facet_circulator sc(e), send(sc);
  //CGAL_assertion( iterator_distance( sc, send) >= 3); // TODO: facet with 2 vertices was found, is it possible?

  Oriented_side facet_side;
  Vertex_handle v;
  do {
    v = sc->source()->center_vertex();
    facet_side = (*this) (pop, v, depth);
    ++sc;
  }
  while( facet_side == ON_ORIENTED_BOUNDARY && sc != send);
  if( facet_side == ON_ORIENTED_BOUNDARY)
    return ON_ORIENTED_BOUNDARY;
  CGAL_assertion( facet_side != ON_ORIENTED_BOUNDARY);
  Oriented_side point_side;
  while( sc != send) {
    v = sc->source()->center_vertex();
    point_side = (*this) (pop, v, depth);
    ++sc;
    if( point_side == ON_ORIENTED_BOUNDARY)
      continue;
    if( point_side != facet_side)
      return ON_ORIENTED_BOUNDARY;
  }
  return facet_side;
  //#endif
}

CGAL_END_NAMESPACE

#endif // CGAL_NEF_SNC_K3_TREE_TRAITS_H