trapezoidal_decomposition_2.h

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

	    if (!lb_circ.is_right_rotation())
	      end_point.set_lb(0);
	  }
    }
    else
    {
      end_point.set_rt(0);
      end_point.set_lb(0);
    }
  }
  
  /*update
    tr.bottom()
    vertical_ray_shoot downward from tr
    tr.top()
    vertical_ray_shoot upward from tr
  */
  reference
  insert_curve_at_point_using_geometry(const X_curve &     cv,
                                       const Point&        p,
                                       pointer &           tr,
                                       const Locate_type &
                                       CGAL_precondition_code(lt))
  {
    CGAL_assertion(traits);
    CGAL_precondition(lt == POINT);
    
    if (traits->compare_cw_around_point_2_object()
        (cv, CGAL_CURVE_IS_TO_RIGHT(cv,p),
         tr->top(), CGAL_CURVE_IS_TO_RIGHT(tr->top(),p), p) == SMALLER)
      tr->set_top(cv);
    if (traits->compare_cw_around_point_2_object()
        (cv, CGAL_CURVE_IS_TO_RIGHT(cv,p), tr->bottom(),
         CGAL_CURVE_IS_TO_RIGHT(tr->bottom(),p), p, false) == SMALLER)
      tr->set_bottom(cv);
    return *tr;
  }
  
  reference
  insert_curve_at_point_using_data_structure(const X_curve & cv,
                                             const Point & p,
                                             pointer & tr,
                                             const Locate_type &
                                             CGAL_precondition_code(lt))
  {
    CGAL_precondition(lt==TRAPEZOID || lt==UNBOUNDED_TRAPEZOID);
    
    Data_structure *tt=tr->get_node();
    
#ifdef CGAL_TD_DEBUG
    
    CGAL_assertion(tr->get_node());
    
#endif
    
    return *split_trapezoid_by_point(*tt,p,cv,cv);
    //return *tr;
  }
        
  void replace_curve_at_point_using_geometry(const X_curve& old_cv,
                                             const X_curve& new_cv,
                                             reference sep)
  {
    
#ifdef CGAL_TD_DEBUG
    
    CGAL_precondition(traits->is_degenerate_point(sep));
    
#endif
    
    if (!sep.is_top_unbounded() && traits->equal_2_object()(sep.top(), old_cv))
      sep.set_top(new_cv);
    if (!sep.is_bottom_unbounded() &&
        traits->equal_2_object()(sep.bottom(), old_cv)) sep.set_bottom(new_cv);
  }
  
  /* update geometric boundary(top and bottom) for trapezoids
     traveled along an iterator till end reached
     precondition:
     end==0 or end is on the path of the iterator
     postcondition:
     end is pointer to the last trapezoid encountered,if any
  */
  void replace_curve_using_geometry(In_face_iterator & it,
                                    const X_curve & old_cv,
                                    const X_curve & new_cv,pointer & end)
  {
    pointer last=0;
    while (it.operator->()!=end)
    {
      if (!it->is_top_unbounded() &&
          traits->equal_2_object()(it->top(), old_cv))
        it->set_top(new_cv);
      if (!it->is_bottom_unbounded() &&
          traits->equal_2_object()(it->bottom(),old_cv))
        it->set_bottom(new_cv);
      last=it.operator->();
      ++it;
    }
    end=last;
  }
  
  /*
    description:
    advances input Data structure using data structure,input point p and 
    possibly X_curve cv till
    p is found(if cv hadn't been given)
    cv is found(if cv was given)
    or
    leaf node reached
    postcondition:
    output is the closest active trapezoid to p/cv
    remark:
    use this function with care!
  */
  /*static */
  Locate_type search_using_data_structure(Data_structure& curr,const_Traits_ptr traits,
                                          const Point& p,const X_curve* cv,
                                          Comparison_result up = EQUAL) const
  {
    const Point* pp;
    const X_curve* pc;
    
#ifdef CGAL_TD_DEBUG
    pointer old = NULL;
#endif
    
    while(true)
    {
#ifdef CGAL_TD_DEBUG
	// unbounded loop
      CGAL_assertion(curr.operator->() != old);
      old = curr.operator->();
#endif

      if (traits->is_degenerate_point(*curr))
        // point node conditional (separation)
	  {
	    // extract point from trapezoid
	    pp = &curr->left();
	    if (CGAL_POINT_IS_LEFT_LOW(p, *pp))
        {
          curr = curr.left();
          continue;
        }
	    else if (CGAL_POINT_IS_LEFT_LOW(*pp, p))
        {
          curr = curr.right();
          continue;
        }
	    else if (traits->equal_2_object()(*pp, p))
        {
          if (!cv)
		  {
		    if ( up == EQUAL ) {                // point found!
		      if (curr->is_active()) return POINT;
		      curr = curr.left();
		    }
		    else if ( up == LARGER ) {          // vertical ray shut up
		      curr = curr.right();                      
		    }
		    else /*if ( up == SMALLER ) */ {
		      curr = curr.left();               // vertical ray shut down
		    }
		    continue;
		  }
          else
		  {

#ifndef CGAL_TD_DEBUG
		    CGAL_warning(traits->equal_2_object()
                         (traits->construct_min_vertex_2_object()(*cv), p) ||
                         traits->equal_2_object()
                         (traits->construct_max_vertex_2_object()(*cv), p));
#else
		    CGAL_assertion(traits->equal_2_object()
                           (traits->construct_min_vertex_2_object()(*cv), p) ||
                           traits->equal_2_object()
                           (traits->construct_max_vertex_2_object()(*cv), p));
#endif
		    curr = traits->equal_2_object()
		      (traits->construct_min_vertex_2_object()(*cv), p) ?
		      curr.right() : curr.left();
		    // (Oren 14/4/02) ??
                    
		    continue;
		  }
        }
	    else
        {
                
#ifndef CGAL_TD_DEBUG
          CGAL_warning(CGAL_POINT_IS_LEFT_LOW(p,*pp) ||
                       CGAL_POINT_IS_LEFT_LOW(*pp,p) ||
                       traits->equal_2_object()(*pp,p));
#else
		CGAL_assertion(CGAL_POINT_IS_LEFT_LOW(p,*pp) ||
                       CGAL_POINT_IS_LEFT_LOW(*pp,p) ||
                       traits->equal_2_object()(*pp,p));
#endif

		return Locate_type();
        }
	  }
      if (traits->is_degenerate_curve(*curr))
	  {
	    // CURVE SEPRATION
	    pc = &curr->top();
	    Comparison_result cres = traits->compare_y_at_x_2_object()(p, *pc);
	    if (cres == SMALLER)
        {
          curr = curr.left();
          continue;
        }
	    else if (cres == LARGER)
        {
          curr = curr.right();
          continue;
        }
	    else
        {  
          // p on CURVE  
#ifndef CGAL_TD_DEBUG      
          CGAL_warning((cres == EQUAL) &&
               !(traits->compare_x_2_object()(p,
                  traits->construct_max_vertex_2_object()(*pc)) == LARGER) &&
                       !(traits->compare_x_2_object()(p,
                  traits->construct_min_vertex_2_object()(*pc)) == SMALLER));
#else
                
          CGAL_postcondition((cres == EQUAL) &&
                             !(traits->compare_x_2_object()(p,
                  traits->construct_max_vertex_2_object()(*pc)) == LARGER) &&
                             !(traits->compare_x_2_object()(p,
                  traits->construct_min_vertex_2_object()(*pc)) == SMALLER));
#endif
          if (!cv)
		  {
		    // For a vertical curve, we always visit it after visiting
		    // one of its endpoints.
		    if ((up == EQUAL) || traits->is_vertical(*curr)) {
		      //std::cout << "EQUAL or VERTICAL" << std::endl;
		      if (curr->is_active()) return CURVE;
		      curr = curr.left();
		    }
		    else if (up == LARGER) {
		      curr = curr.right();
		    }
		    else /* if (up==SMALLER) */ {
		      curr = curr.left();
		    }
		    continue;
		  }
          else
		  {
                    
#ifndef CGAL_TD_DEBUG          
		    CGAL_warning(traits->equal_2_object()
                         (traits->construct_min_vertex_2_object()(*cv),
                          traits->construct_min_vertex_2_object()(*pc)) ||
                         traits->equal_2_object()
                         (traits->construct_max_vertex_2_object()(*cv),
                          traits->construct_max_vertex_2_object()(*pc)));
#else
		    if (!(traits->equal_2_object()
                  (traits->construct_min_vertex_2_object()(*cv),
                   traits->construct_min_vertex_2_object()(*pc))||
                  traits->equal_2_object()
                  (traits->construct_max_vertex_2_object()(*cv),
                   traits->construct_max_vertex_2_object()(*pc))))
            {
              std::cerr << "\npc " << *pc;
              std::cerr << "\ncv " << *cv << std::endl;
              CGAL_assertion(traits->equal_2_object()
                             (traits->construct_min_vertex_2_object()(*cv),
                              traits->construct_min_vertex_2_object()(*pc)) ||
                             traits->equal_2_object()
                             (traits->construct_max_vertex_2_object()(*cv),
                              traits->construct_max_vertex_2_object()(*pc)));
            }
#endif

		    Comparison_result res =
		      traits->equal_2_object()
		      (traits->construct_min_vertex_2_object()(*cv),
		       traits->construct_min_vertex_2_object()(*pc)) ?
		      traits->compare_cw_around_point_2_object()
		      (*pc, CGAL_CURVE_IS_TO_RIGHT(*pc,p),
		       *cv, CGAL_CURVE_IS_TO_RIGHT(*cv,p), p) :
		      traits->compare_cw_around_point_2_object()
		      (*cv, CGAL_CURVE_IS_TO_RIGHT(*cv,p),
		       *pc, CGAL_CURVE_IS_TO_RIGHT(*pc,p), p ,false);
                    
		    switch(res)
            {
		      case LARGER:
               curr = curr.right();
               break;
             case SMALLER:
              curr = curr.left();
              break;
             case EQUAL:
              switch(up)
			  {
               case LARGER:
			    curr = curr.right();
			    break;
               case SMALLER:
			    curr = curr.left();
			    break;
               case EQUAL:
			    if (curr->is_active()) return CURVE;
			    curr = curr.left();
			    break;
                            
#ifdef CGAL_TD_DEBUG
               default:
			    CGAL_assertion(up==LARGER||up==SMALLER||up==EQUAL);
			    return Locate_type();
#endif
			  }
              break;

#ifdef CGAL_TD_DEBUG
             default:
              CGAL_assertion(res == LARGER || res == SMALLER || res == EQUAL);
              return Locate_type();
#endif

            }
		  }
        } 
	  }
      else
	  {
	    // !is_degenerate()
	    if (curr->is_active())
	      return curr->is_unbounded() ? UNBOUNDED_TRAPEZOID : TRAPEZOID;
	    curr = curr.left();
	    continue;
	  }
    }
  }

  Data_structure container2data_structure(map_nodes& ar, int left,
                                          int right) const
  {

#ifndef CGAL_TD_DEBUG
    
    CGAL_warning(traits);
    
#else
    
    CGAL_assertion(traits);
    
#endif
    
    if (right>left)
    {
      int d=(int)CGAL_CLIB_STD::floor((double(right+left))/2);
      // Replacing operator [] of map with find to please MSVC 7
      Point p = (ar.find(d)->second)->right();
      //Point p=ar[d]->right();
      Data_structure curr=
        Data_structure(X_trapezoid(&p,&p,0,0),
                       container2data_structure(ar,left,d),
                       container2data_structure(ar,d+1,right));
      curr.left()->set_node((Data_structure*)&curr.left());
      curr.right()->set_node((Data_structure*)&curr.right());
      curr->set_node(&curr);// fake temporary node
      curr->remove(); // mark as deleted
      curr->set_node(0);

      return curr;
    }
    else
      // Replacing operator [] of map with find to please MSVC 7
      return ar.find(left)->second;
    //return ar[left];
  }
  /*==============================================
    Trapezoidal_decomposition_2 public member functions
    ==============================================*/
public:
  
  Trapezoidal_decomposition_2(bool rebuild=true) :
    depth_threshold(CGAL_TD_DEFAULT_DEPTH_THRESHOLD),
    size_threshold(CGAL_TD_DEFAULT_SIZE_THRESHOLD) 
  {init();set_needs_update(rebuild);}
  Trapezoidal_decomposition_2(const double& depth_th,const double& size_th) : 
    depth_threshold(depth_th),size_threshold(size_th) 
  {init();set_needs_update(rebuild);}
  Trapezoidal_decomposition_2(const_Self_ref td) :
	needs_update_(td.needs_update_),
	number_of_curves_(td.number_of_curves_),    
	traits(td.traits),
	last_cv(NULL), prev_cv(NULL), 
	depth_threshold(td.depth_threshold),
	size_threshold(td.size_threshold)
  {
    hash_map_tr_ptr htr;
    /*! \todo allocate hash_map size according to content.
     * \todo change vector<> to in_place_list and pointer hash to trapezoidal
     * hash..
     */
    vector_container vtr;
    int sz;
    Td_active_trapezoid pr;
    sz=X_trapezoid_filter(vtr, &td.get_data_structure());
    //! \todo Reduce the 3 iterations to 1 (or 2) iterator.
    // First iteration: filter out the active trapezoids.
    typename vector_container::const_iterator it;
    for (it=vtr.begin(); it!=vtr.end(); ++it) {
      Data_structure* ds_copy=new Data_structure(*it);
      const X_trapezoid* cur=&*it;
      X_trapezoid* tr_copy=&*(*ds_copy);
      tr_copy->set_node(ds_copy);
      CGAL_assertion(&*(*tr_copy->get_node())==tr_copy);
      ds_copy->set_depth(cur->get_node()->depth());
      // We cheat a little with the depth.
      htr.insert(typename hash_map_tr_ptr::value_type(cur, tr_copy));
      // Second iteration: generate new copies of trapezoids and nodes.
    }
    for (it=vtr.begin(); it!=vtr.end(); ++it) {
      const X_trapezoid* cur=&*it;
      X_trapezoid* tr_copy=htr.find(cur)->second;
      const Data_structure *child;