snc_point_locator.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/SNC_point_locator.h $
// $Id: SNC_point_locator.h 37431 2007-03-23 21:19:49Z hachenb $
// 
//
// Author(s)     : Miguel Granados <[e-mail]>

#ifndef CGAL_NEF_SNC_POINT_LOCATOR_H
#define CGAL_NEF_SNC_POINT_LOCATOR_H

#include <CGAL/basic.h>
#include <CGAL/Nef_3/SNC_intersection.h>
#ifdef CGAL_NEF3_POINT_LOCATOR_NAIVE
#include <CGAL/Nef_3/SNC_ray_shooter.h>
#endif
#include <CGAL/Nef_3/SNC_k3_tree_traits.h>
#include <CGAL/Nef_3/K3_tree.h>
#include <CGAL/Unique_hash_map.h>
#include <CGAL/Timer.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

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

#undef _CGAL_NEF_TRACEN
#define _CGAL_NEF_TRACEN(msg) CGAL_NEF_TRACEN( "SNC_point_locator: " << msg);

// TODO: find out the proper CGAL replacement for this macro and remove it
#define CGAL_for_each( i, C) for( i = C.begin(); i != C.end(); ++i)

// #define CGAL_NEF_TIMER(instruction) instruction
#define CGAL_NEF_TIMER(instruction)

// #define CGAL_NEF_CLOG(t) std::clog <<" "<<t<<std::endl; std::clog.flush()
#define CGAL_NEF_CLOG(t)

CGAL_BEGIN_NAMESPACE

template <typename SNC_decorator>
class SNC_point_locator
{
 public:
  class Intersection_call_back;
  typedef SNC_point_locator<SNC_decorator> Self;
  typedef typename SNC_decorator::Decorator_traits Decorator_traits;
  typedef typename SNC_decorator::SNC_structure SNC_structure;
protected:
  char version_[64];
  // time for construction, point location, ray shooting and intersection test
  mutable Timer ct_t, pl_t, rs_t, it_t; 

public: 
  typedef typename SNC_structure::Object_handle Object_handle;

  typedef typename SNC_structure::Point_3 Point_3;
  typedef typename SNC_structure::Segment_3 Segment_3;
  typedef typename SNC_structure::Ray_3 Ray_3;
  typedef typename SNC_structure::Vector_3 Vector_3;
  typedef typename SNC_structure::Aff_transformation_3 
                                  Aff_transformation_3;

  typedef typename Decorator_traits::Vertex_handle Vertex_handle;
  typedef typename Decorator_traits::Halfedge_handle Halfedge_handle;
  typedef typename Decorator_traits::Halffacet_handle Halffacet_handle;
  typedef typename Decorator_traits::Volume_handle Volume_handle;
  typedef typename Decorator_traits::Vertex_iterator Vertex_iterator;
  typedef typename Decorator_traits::Halfedge_iterator Halfedge_iterator;
  typedef typename Decorator_traits::Halffacet_iterator Halffacet_iterator;

  const char* version() { return version_; }

  virtual Object_handle locate(const Point_3& p) const = 0;

  virtual Object_handle shoot(const Ray_3& s, int mask=255) const = 0;

  virtual void intersect_with_edges( Halfedge_handle edge,
                                     const Intersection_call_back& call_back) 
    const = 0;

  virtual void intersect_with_facets( Halfedge_handle edge,
                                      const Intersection_call_back& call_back)
    const = 0;

  virtual void intersect_with_edges_and_facets( Halfedge_handle edge,
	const Intersection_call_back& call_back) const = 0;

  class Intersection_call_back 
  {
  public:
    virtual void operator()( Halfedge_handle edge, Object_handle object, 
                             const Point_3& intersection_point) const = 0;
    
    virtual ~Intersection_call_back() {}
  };

  virtual void initialize(SNC_structure* W) = 0;

  virtual Self* clone() const = 0;

  virtual void transform(const Aff_transformation_3& t) = 0;

  //virtual bool update( Unique_hash_map<Vertex_handle, bool>& V, 
  //                   Unique_hash_map<Halfedge_handle, bool>& E, 
  //                   Unique_hash_map<Halffacet_handle, bool>& F) = 0;

  virtual void add_facet(Halffacet_handle) {}

  virtual void add_edge(Halfedge_handle) {}

  virtual void add_vertex(Vertex_handle) {}

  virtual ~SNC_point_locator() {
    CGAL_NEF_CLOG("");
    CGAL_NEF_CLOG("construction_time:  "<<ct_t.time());
    CGAL_NEF_CLOG("pointlocation_time: "<<pl_t.time());
    CGAL_NEF_CLOG("rayshooting_time:   "<<rs_t.time());
    CGAL_NEF_CLOG("intersection_time:  "<<it_t.time());
    // warning: the total time showed here could be actually larger
    // that the real time used by this class, since point location
    // and intersection test use the ray shooter and so the same time 
    // could be account to more than one timer
    CGAL_NEF_CLOG("pltotal_time:       "<<
      ct_t.time()+pl_t.time()+rs_t.time()+it_t.time());
  };
};

template <typename SNC_decorator>
class SNC_point_locator_by_spatial_subdivision : 
  public SNC_point_locator<SNC_decorator>,
  public SNC_decorator
{ 
  typedef SNC_decorator Base;
  typedef CGAL::SNC_point_locator<SNC_decorator> SNC_point_locator;
  typedef CGAL::SNC_point_locator_by_spatial_subdivision<SNC_decorator> Self;
  typedef typename SNC_decorator::SNC_structure SNC_structure;
  typedef typename SNC_decorator::Decorator_traits Decorator_traits;
  typedef typename Decorator_traits::SM_decorator SM_decorator;
  typedef CGAL::SNC_intersection<SNC_structure> SNC_intersection;
  typedef typename SNC_decorator::SM_point_locator SM_point_locator;
  typedef typename SNC_decorator::Kernel Kernel;

public:
  typedef typename CGAL::SNC_k3_tree_traits<SNC_decorator> K3_tree_traits;
  typedef typename CGAL::K3_tree<K3_tree_traits> K3_tree;
  typedef K3_tree SNC_candidate_provider;
  
  typedef typename SNC_structure::Object_handle Object_handle;
 #ifdef CGAL_NEF3_TRIANGULATE_FACETS
  typedef typename Decorator_traits::Halffacet_triangle_handle 
                                     Halffacet_triangle_handle;
#endif
#ifdef CGAL_NEF3_FACET_WITH_BOX
  typedef typename SNC_structure::Partial_facet Partial_facet;
#endif
  typedef typename SNC_structure::Point_3 Point_3;
  typedef typename SNC_structure::Segment_3 Segment_3;
  typedef typename SNC_structure::Ray_3 Ray_3;
  typedef typename SNC_structure::Vector_3 Vector_3;
  typedef typename SNC_structure::Triangle_3 Triangle_3;
  typedef typename SNC_structure::Aff_transformation_3 
                                  Aff_transformation_3;

  typedef typename SNC_structure::Infi_box Infi_box;

  typedef typename Decorator_traits::Vertex_handle Vertex_handle;
  typedef typename Decorator_traits::Halfedge_handle Halfedge_handle;
  typedef typename Decorator_traits::Halffacet_handle Halffacet_handle;
  typedef typename Decorator_traits::SHalfedge_handle SHalfedge_handle;
  typedef typename Decorator_traits::SFace_handle SFace_handle;
  typedef typename Decorator_traits::Vertex_iterator Vertex_iterator;
  typedef typename Decorator_traits::Halfedge_iterator Halfedge_iterator;
  typedef typename Decorator_traits::Halffacet_iterator Halffacet_iterator;
  typedef typename Decorator_traits::Volume_handle Volume_handle;

  typedef typename Decorator_traits::Halffacet_cycle_iterator 
                                     Halffacet_cycle_iterator;
  typedef typename Decorator_traits::SHalfedge_around_facet_circulator 
	                             SHalfedge_around_facet_circulator;

  typedef typename SNC_candidate_provider::Object_list Object_list;
  typedef typename Object_list::iterator Object_list_iterator;
  typedef typename SNC_candidate_provider::Objects_along_ray Objects_along_ray;
  typedef typename Objects_along_ray::Iterator Objects_along_ray_iterator;

public:
  SNC_point_locator_by_spatial_subdivision() : 
    initialized(false), candidate_provider(0) {}

#ifdef CGAL_NEF3_TRIANGULATE_FACETS	
  template<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;

    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) {

      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) {
            CGAL_NEF_TRACEN("  insert constraint" << sfc->source()->source()->point()
	                     << "->" << sfc->source()->twin()->source()->point());
	    ct.insert_constraint(sfc->source()->source()->point(),
	                         sfc->source()->twin()->source()->point());
          }
        }
      }
      CGAL_assertion(ct.is_valid());

      CGAL_NEF_TRACEN("number of finite triangles " << ct.number_of_faces());

      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;

      CGAL_assertion(!ct.is_infinite(first));
      traverse_triangulation(first, first->index(opposite));

      /*
      for(fi = ct.finite_faces_begin(); fi != ct.finite_faces_end(); ++fi)
        CGAL_NEF_TRACEN("  finite face " 
	  << Triangle_3(fi->vertex(0)->point(), fi->vertex(1)->point(), fi->vertex(2)->point())
	   << "was visited " << visited[fi]);
      */

      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.ccw(parent))) && !visited[f->neighbor(ct.ccw(parent))]) {
	Face_handle child(f->neighbor(ct.ccw(parent)));
	traverse_triangulation(child, child->index(f));
      } 
    } 
 
    template<typename Triangle_3>
    bool get_next_triangle(Triangle_3& tr) {
      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());
      ++fi;
      return true;
    }
  };
#endif

  virtual void initialize(SNC_structure* W) {
#ifdef CGAL_NEF_LIST_OF_TRIANGLES
    set_snc(*W);
    candidate_provider = new SNC_candidate_provider(W);
#else // CGAL_NEF_LIST_OF_TRIANGLES
    CGAL_NEF_TIMER(ct_t.start());
    strcpy( this->version_, "Point Locator by Spatial Subdivision (tm)");
#ifdef CGAL_NEF3_TRIANGULATE_FACETS
    CGAL_NEF_CLOG(version()<<" (with triangulated facets)");
#else
    CGAL_NEF_CLOG(version());
#endif
    CGAL_assertion( W != NULL);
//    (Base) *this = SNC_decorator(*W);
	set_snc(*W);
    Object_list objects;
    Vertex_iterator v;
    Halfedge_iterator e;
    Halffacet_iterator f;
    CGAL_forall_vertices( v, *this->sncp())
      objects.push_back(Object_handle(Vertex_handle(v)));
    typename Object_list::difference_type v_end = objects.size();
    CGAL_forall_edges( e, *this->sncp())
      objects.push_back(Object_handle(Halfedge_handle(e)));
    CGAL_forall_facets( f, *this->sncp()) {
#ifdef CGAL_NEF3_TRIANGULATE_FACETS

#ifndef CGAL_NEF3_TRIANGULATION_MINIMUM
#define CGAL_NEF3_TRIANGULATION_MINIMUM 25
#endif

      Halffacet_cycle_iterator fci = f->facet_cycles_begin();
      CGAL_assertion(fci.is_shalfedge());
      SHalfedge_around_facet_circulator safc(fci), send(safc);
      int length = 0;
      int stop = CGAL_NEF3_TRIANGULATION_MINIMUM;
      while(++length < stop && ++safc != send);
      if(length >= stop) {

	CGAL_NEF_TRACEN("triangulate facet " << f->plane());
      
	typedef typename CGAL::Triangulation_euclidean_traits_xy_3<Kernel>       XY;
	typedef typename CGAL::Triangulation_euclidean_traits_yz_3<Kernel>       YZ;
	typedef typename CGAL::Triangulation_euclidean_traits_xz_3<Kernel>       XZ;
	
	Triangle_3 tr;
	
	Vector_3 orth = f->plane().orthogonal_vector();
	int c = CGAL::abs(orth[0]) > CGAL::abs(orth[1]) ? 0 : 1;
	c = CGAL::abs(orth[2]) > CGAL::abs(orth[c]) ? 2 : c;
	
	if(c == 0) {
	  Triangulation_handler<YZ> th(f);
	  while(th.get_next_triangle(tr)) {
	    Halffacet_triangle_handle th( f, tr);
	    objects.push_back(Object_handle(th));
	    CGAL_NEF_TRACEN("add triangle " << tr);
	  }
	} else if(c == 1) {
	  Triangulation_handler<XZ> th(f);
	  while(th.get_next_triangle(tr)) {
	    Halffacet_triangle_handle th( f, tr);
	    objects.push_back(Object_handle(th));
	    CGAL_NEF_TRACEN("add triangle " << tr);
	  }
	} else if(c == 2) {
	  Triangulation_handler<XY> th(f);
	  while(th.get_next_triangle(tr)) {
	    Halffacet_triangle_handle th( f, tr);
	    objects.push_back(Object_handle(th));
	    CGAL_NEF_TRACEN("add triangle " << tr);
	  }
	} else
	  CGAL_assertion_msg(false, "wrong value");
      } else {
        CGAL_NEF_TRACEN("add facet " << f->plane());
        objects.push_back(Object_handle(Halffacet_handle(f)));
      }
#elif defined CGAL_NEF3_FACET_WITH_BOX
#ifndef CGAL_NEF3_PARTITION_MINIMUM
#define CGAL_NEF3_PARTITION_MINIMUM 6
#endif
      Halffacet_cycle_iterator fci = f->facet_cycles_begin();
      CGAL_assertion(fci.is_shalfedge());
      SHalfedge_around_facet_circulator safc(fci), send(safc);
      int length = 0;
      int stop = CGAL_NEF3_PARTITION_MINIMUM;
      while(++length < stop && ++safc != send);
      if(length >= stop) {
	CGAL_NEF_TRACEN("use Partial facets ");
	Partial_facet pf(f);
	objects.push_back(Object_handle(pf));
      } else {
        objects.push_back(Object_handle(Halffacet_handle(f)));
      }
#else
      objects.push_back(Object_handle(Halffacet_handle(f)));
#endif
    }
    Object_list_iterator oli=objects.begin()+v_end;
    if(initialized)
      delete candidate_provider;
    candidate_provider = new SNC_candidate_provider(objects,oli);
    // CGAL_NEF_TRACEN(*candidate_provider);
    CGAL_NEF_TIMER(ct_t.stop());
#endif // CGAL_NEF_LIST_OF_TRIANGLES
    initialized = true;
  }

  virtual Self* clone() const { 
    return new Self; 
  }