snc_external_structure.h

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

    SHalfedge_iterator e;
    CGAL_forall_shalfedges(e,*this->sncp()) {
      Sphere_circle c(e->circle());
      Plane_3 h = c.plane_through(e->source()->source()->point());
      CGAL_NEF_TRACEN("\n" << e->source()->twin()->source()->point() <<" - "
		      << e->source()->source()->point() <<" - "<< 
		      e->twin()->source()->twin()->source()->point() << 
		      " has plane " << h << " has circle " << e->circle() << 
		      " has signum " << sign_of(h));
      if ( sign_of(h)<0 ) continue;
      M[normalized(h)].push_back(Object_handle(e->twin())); 
      CGAL_NEF_TRACEN(" normalized as " << normalized(h));
      /*
	Unique_hash_map<SHalfedge_handle, bool> Done(false);
	SHalfedge_iterator ei;
	CGAL_forall_sedges(ei,*this->sncp()) {
	if(Done[ei]) continue;
	Sphere_circle c(ei->circle());
	Plane_3 h = c.plane_through(ei->source()->source()->point());
	CGAL_NEF_TRACEN("\n" << ei->source()->twin()->source()->point() <<" - "
	<< ei->source()->source()->point() <<" - "<< 
		      ei->twin()->source()->twin()->source()->point() << 
		      " has plane " << h << " has circle " << ei->circle() << 
		      " has signum " << sign_of(h));
      SHalfedge_handle e(ei);
      if ( sign_of(h)<0 ) {
	h = h.opposite();
	e = e->twin();
      }
      SHalfedge_around_facet_circulator sfc(e), send(sfc);
      CGAL_For_all(sfc, send) {
	M[normalized(h)].push_back(Object_handle(e->twin())); 
	Done[sfc] = true;
	Done[sfc->twin()] = true;
	CGAL_NEF_TRACEN(" normalized as " << normalized(h)); 
      }
      */
    }
    SHalfloop_iterator l;
    CGAL_forall_shalfloops(l,*this->sncp()) {
      Sphere_circle c(l->circle());
      Plane_3 h = c.plane_through(l->incident_sface()->center_vertex()->point()); 
      if ( sign_of(h)<0 ) continue;
      // CGAL_assertion( h == normalized(h));
      M[normalized(h)].push_back(Object_handle(l->twin()));
    }
    
#ifdef CGAL_NEF3_TIMER_PLANE_SWEEPS
  number_of_plane_sweeps=0;
  timer_plane_sweeps.reset();
#endif

    typename Map_planes::iterator it;
    CGAL_forall_iterators(it,M) { 
      //    progress2++;
      //      CGAL_NEF_TRACEN("  plane "<<it->first<<"             "<<(it->first).point());
      FM_decorator D(*this->sncp());
      D.create_facet_objects(it->first,it->second.begin(),it->second.end());
    }
    //       CGAL_NEF_SETDTHREAD(1);
  }

  void create_volumes() {
  /*{\Mop collects all shells incident to a volume and creates the
  volumes.  \precond |categorize_facet_cycles_and_creating_facets()| was
  called before.}*/

#ifdef CGAL_NEF3_TIMER_POINT_LOCATION
    number_of_ray_shooting_queries=0;
    timer_ray_shooting.reset();
#endif 

    //    CGAL_NEF_SETDTHREAD(37*43*503*509);

    CGAL_NEF_TRACEN(">>>>>create_volumes");
    Sface_shell_hash     ShellSf(0);
    Face_shell_hash      ShellF(0);
    SFace_visited_hash Done(false);
    Shell_explorer V(*this,ShellSf,ShellF,Done);
    std::vector<SFace_handle> MinimalSFace;
    std::vector<SFace_handle> EntrySFace;
    std::vector<bool> Closed;
    
    SFace_iterator f;
    // First, we classify all the Shere Faces per Shell.  For each Shell we
    //     determine its minimum lexicographyly vertex and we check wheter the
    //     Shell encloses a region (closed surface) or not. 
    CGAL_forall_sfaces(f,*this->sncp()) {
      //    progress++;
      CGAL_NEF_TRACEN("sface in " << ShellSf[f]);
      if ( Done[f] ) 
	continue;
      V.minimal_sface() = f;
      visit_shell_objects(f,V);
      
      CGAL_NEF_TRACEN("minimal vertex " << V.minimal_sface()->center_vertex()->point());

      MinimalSFace.push_back(V.minimal_sface());
      EntrySFace.push_back(f);
      V.increment_shell_number();
      CGAL_NEF_TRACEN("sface out " << ShellSf[f]);
    }
    
    for(unsigned int i=0; i<EntrySFace.size(); ++i)
      Closed.push_back(false);
    
    Halffacet_iterator hf;
    CGAL_forall_facets(hf,*this) 
      if(ShellF[hf] != ShellF[hf->twin()]) {
	Closed[ShellF[hf]] = true;
	Closed[ShellF[hf->twin()]] = true;
      }
    
    CGAL_assertion( pl != NULL);

#ifdef CGAL_NEF3_TIMER_INITIALIZE_KDTREE
    CGAL::Timer timer_initialize_kdtree;
    timer_initialize_kdtree.start();
#endif
    pl->initialize(this->sncp()); // construct the point locator 
#ifdef CGAL_NEF3_TIMER_INITIALIZE_KDTREE
    timer_initialize_kdtree.stop();
    if(cgal_nef3_timer_on)
      std::cout << "Runtime_initialize_kdtree: " 
		<< timer_initialize_kdtree.time() << std::endl;
#endif

    //   then, we determine the Shells which correspond to Volumes via a ray
    //   shootting in the direction (-1,0,0) over the Sphere_map of the minimal 
    //   vertex.  The Shell corresponds to a Volume if the object hit belongs 
    //   to another Shell. 
    
    this->sncp()->new_volume(); // outermost volume (nirvana)
    if(MinimalSFace.size() == 0) return;
    Vertex_handle v_min = MinimalSFace[0]->center_vertex();
    for( unsigned int i = 0; i < MinimalSFace.size(); ++i) {
      //    progress2++;
      Vertex_handle v = MinimalSFace[i]->center_vertex();
      if(CGAL::lexicographically_xyz_smaller(v->point(),v_min->point()))
	v_min=v;
      CGAL_NEF_TRACEN( "Shell #" << i << " minimal vertex: " << v->point());
      SM_point_locator D((Sphere_map*) &*v);
      Object_handle o = D.locate(Sphere_point(-1,0,0));
      SFace_const_handle sfc;
      if( !CGAL::assign(sfc, o) || ShellSf[sfc] != i) {
	CGAL_NEF_TRACEN("the shell encloses a volume");
	CGAL_NEF_TRACEN("sface hit? "<<CGAL::assign(sfc,o));
	if( CGAL::assign(sfc, o)) CGAL_NEF_TRACEN("sface on another shell? "<<ShellSf[sfc]);
	SFace_handle f = MinimalSFace[i];
	CGAL_assertion( ShellSf[f] == i );
	if( Closed[i] ) {
	  CGAL_NEF_TRACEN("Shell #" << i << " is closed");
	  SM_decorator SD(&*v);
	  Volume_handle c = this->sncp()->new_volume();
	  c->mark() = f->mark(); // TODO test if line is redundant
	  link_as_inner_shell(f, c );
	  CGAL_NEF_TRACE( "Shell #" << i <<" linked as inner shell");
	  CGAL_NEF_TRACEN( "(sface" << (CGAL::assign(sfc,o)?"":" not") << " hit case)");
	}
      }
    }
    
    // finaly, we go through all the Shells which do not correspond to a Volume 
    //     and we assign them to its enclosing Volume determined via a facet below
    //     check. 
    
    CGAL_forall_sfaces(f,*this->sncp()) {
      //    progress3++;
      if ( f->volume() != Volume_handle() ) 
	continue;
      CGAL_NEF_TRACEN( "Outer shell #" << ShellSf[f] << " volume?");
      Volume_handle c = determine_volume( MinimalSFace[ShellSf[f]], 
					  MinimalSFace, ShellSf );
      c->mark() = f->mark();
      link_as_outer_shell( f, c );
    } 
  }

  Halffacet_handle get_facet_below( Vertex_handle vi, 
				    const std::vector< SFace_handle>& MinimalSFace, 
				    const Sface_shell_hash&  Shell) const {
    // {\Mop determines the facet below a vertex |vi| via ray shooting. }
    
    Halffacet_handle f_below;
    Point_3 p = vi->point();
    if(!Infi_box::is_standard(p))
      return Halffacet_handle();
    
    Ray_3 ray = Ray_3(p, Direction_3(-1,0,0));
#ifdef CGAL_NEF3_TIMER_POINT_LOCATION
    number_of_ray_shooting_queries++;
    timer_ray_shooting.start();
#endif 
    Object_handle o = pl->shoot(ray);
#ifdef CGAL_NEF3_TIMER_POINT_LOCATION
    timer_ray_shooting.stop();
#endif 
    // The ray here has an special property since it is shooted from the lowest
    // vertex in a shell, so it would be expected that the ray goes along the
    // interior of a volume before it hits a 2-skeleton element.
    // Unfortunatelly, it seems to be possible that several shells are incident
    // to this lowest vertex, and in consequence, the ray could also go along
    // an edge or a facet belonging to a different shell.
    // This fact invalidates the precondition of the get_visible_facet method,
    // (the ray must pierce the local view of the hit object in a sface).
    // This situation has not been analyzed and has to be verified. Granados.
    Vertex_handle v;
    Halfedge_handle e;
    Halffacet_handle f;
    CGAL_NEF_TRACEN("get_facet_below");
    if( CGAL::assign(v, o)) {
      CGAL_NEF_TRACEN("facet below from from vertex...");
      f_below = get_visible_facet(v, ray);
      if( f_below == Halffacet_handle()) {
	CGAL_assertion(v->sfaces_begin() == v->sfaces_last());
	f_below = get_facet_below(MinimalSFace[Shell[v->sfaces_begin()]]->center_vertex(), 
				  MinimalSFace,Shell);
      }
    }
    else if( CGAL::assign(e, o)) {
      CGAL_NEF_TRACEN("facet below from from edge...");
      f_below = get_visible_facet(e, ray);
      if( f_below == Halffacet_handle()) {
	CGAL_assertion(e->source()->sfaces_begin() == e->source()->sfaces_last());
	f_below = get_facet_below(MinimalSFace[Shell[e->source()->sfaces_begin()]]->center_vertex(), 
				  MinimalSFace, Shell);
      }
    }
    else if( CGAL::assign(f, o)) {
      CGAL_NEF_TRACEN("facet below from from facet...");
      f_below = get_visible_facet(f, ray);
      CGAL_assertion( f_below != Halffacet_handle());
    }
    else { CGAL_NEF_TRACEN("no facet below found..."); }
    return f_below;
  }  

  Volume_handle determine_volume( SFace_handle sf, 
                const std::vector< SFace_handle>& MinimalSFace, 
				  const Sface_shell_hash&  Shell ) const {
    //{\Mop determines the volume |C| that a shell |S| pointed by |sf| 
    //  belongs to.  \precondition |S| separates the volume |C| from an enclosed
    //  volume.}
    
    CGAL_NEF_TRACEN("determine volume");
    Vertex_handle v_min = MinimalSFace[Shell[sf]]->center_vertex();
    
    Halffacet_handle f_below = get_facet_below(v_min, MinimalSFace, Shell);
    if ( f_below == Halffacet_handle())
      return SNC_decorator(*this).volumes_begin();
    Volume_handle c = f_below->incident_volume();
    if( c != Volume_handle()) {
      CGAL_NEF_TRACE( "Volume " << &*c << " hit ");
      CGAL_NEF_TRACEN("(Shell #" << Shell[adjacent_sface(f_below)] << ")");
      return c;
    }
    SFace_handle sf_below = adjacent_sface(f_below);
    CGAL_NEF_TRACE( "Shell not assigned to a volume hit ");
    CGAL_NEF_TRACEN( "(Inner shell #" << Shell[sf_below] << ")");
    c = determine_volume( sf_below, MinimalSFace, Shell);
    link_as_inner_shell( sf_below, c);
    return c;
  }
   
  void build_external_structure() {

#ifdef CGAL_NEF3_TIMER_EXTERNAL_STRUCTURE
    CGAL::Timer timer_external_structure;
    timer_external_structure.start();
#endif

#ifdef CGAL_NEF3_TIMER_PLUECKER
    CGAL::Timer timer_pluecker;
    timer_pluecker.start();
#endif
    //    SNC_io_parser<SNC_structure> O0(std::cerr,*this->sncp());
    //    O0.print();
    pair_up_halfedges();
#ifdef CGAL_NEF3_TIMER_PLUECKER
    timer_pluecker.stop();
     if(cgal_nef3_timer_on)
      std::cout << "Runtime_pluecker: " 
		<< timer_pluecker.time() << std::endl;
#endif
    link_shalfedges_to_facet_cycles();
    //    SNC_io_parser<SNC_structure> O0(std::cerr,*this->sncp());
    //    O0.print();
    categorize_facet_cycles_and_create_facets();
    create_volumes();

#ifdef CGAL_NEF3_TIMER_EXTERNAL_STRUCTURE
    timer_external_structure.stop();
    if(cgal_nef3_timer_on)
      std::cout << "Runtime_external_structure: " 
		<< timer_external_structure.time() << std::endl;
#endif
  }

  template<typename Association>
  void build_after_binary_operation(Association) {

#ifdef CGAL_NEF3_TIMER_SIMPLIFICATION
    CGAL::Timer timer_simplification;
    timer_simplification.start();
#endif
    SNC_simplify simp(*this->sncp());
    simp.vertex_simplification(SNC_simplify::NO_SNC);
#ifdef CGAL_NEF3_TIMER_SIMPLIFICATION
    timer_simplification.stop();
    if(cgal_nef3_timer_on)
      std::cout << "Runtime_simplification: " 
		<< timer_simplification.time() << std::endl;
#endif
    
    CGAL_NEF_TRACEN("\nnumber of vertices (so far...) = "
		    << this->sncp()->number_of_vertices());
    
    CGAL_NEF_TRACEN("=> resultant vertices (after simplification): ");
    
    build_external_structure();
  }

  void clear_external_structure() {
    this->sncp()->clear_snc_boundary();

    while(this->sncp()->volumes_begin()!= this->sncp()->volumes_end())
      this->sncp()->delete_volume(this->sncp()->volumes_begin());

    while(this->sncp()->halffacets_begin()!= this->sncp()->halffacets_end())
      this->sncp()->delete_halffacet_pair(this->sncp()->halffacets_begin());

    SHalfedge_iterator se;
    CGAL_forall_shalfedges(se,*this)
      se->facet() = Halffacet_handle();

    SFace_iterator sf;
    CGAL_forall_sfaces(sf,*this)
      sf->volume() = Volume_handle();
  }
};