bezier_bounding_rational_traits.h

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

  }


  // Returns true if endpoints coincide, and adds the Bound_pair
  // to the list if the endpoint hasn't already been inserted there.
  // Precondition !spansOverlap
  bool _endpoints_coincide(const Control_point_vec& cv1, const NT& l1, const NT& r1,
			        const Control_point_vec& cv2, const NT& l2, const NT& r2,
              std::set<NT>& endpoint_intersections,
              Bound_pair_lst& intersection_pairs) const
  {
    const Point_2& s1 = cv1.front();
    const Point_2& t1 = cv1.back();
    const Point_2& s2 = cv2.front();
    const Point_2& t2 = cv2.back();

    // What happens when an endpoint of one is an inner point (e.g., 1/3) of another?
    // We will not be able to discover this endcase and therefore get to can_refine()==false.
    int endpoint_intersection_num = 0; // For sanity check and return value.

    NT x, y, t_param1, t_param2;
    if (s1==s2) // Using operator== of kernel Point_2
    {
      ++endpoint_intersection_num;
      x = s1.x();
      y = s1.y();
      t_param1 = l1;
      t_param2 = l2;
    }
    if (s1==t2)
    {
      ++endpoint_intersection_num;
      x = s1.x();
      y = s1.y();
      t_param1 = l1;
      t_param2 = r2;
    }
    if (t1==s2)
    {
      ++endpoint_intersection_num;
      x = t1.x();
      y = t1.y();
      t_param1 = r1;
      t_param2 = l2;
    }
    if (t1==t2)
    {
      ++endpoint_intersection_num;
      x = t1.x();
      y = t1.y();
      t_param1 = r1;
      t_param2 = r2;
    }

    CGAL_assertion(endpoint_intersection_num <= 1);
    if (endpoint_intersection_num == 0)
      return false;

    // Construct a degenerate rational Bound_pair
    if (endpoint_intersections.find(t_param1) == endpoint_intersections.end())
    {
      // Found a new intersection point that wasn't inserted yet.
      endpoint_intersections.insert(t_param1); 
      intersection_pairs.push_back(Bound_pair(
                                    Bez_point_bound(cv1, t_param1, t_param1, Bez_point_bound::RATIONAL_PT, true),
                                    Bez_point_bound(cv2, t_param2, t_param2, Bez_point_bound::RATIONAL_PT, true),
                                    Bez_point_bbox(x,x,y,y)
                                  )
                                );
    }
    return true;
  }


  // Auxilary recursive function for intersecting two curves.
  void _CrvCrvInterAux
      (const Control_point_vec& cv1, const NT& l1, const NT& r1,
       const Control_point_vec& cv2, const NT& l2, const NT& r2,
       bool should_check_span,
       std::set<NT>& endpoint_intersections,
       Bound_pair_lst& intersection_pairs)
  {
    //iddo debug
    //static unsigned int dbg = 0;
    //std::cout << ++dbg << ", " << recursion_level << "; ";

    // Check if we got to subdivision termination criteria. 
    bool can_refine1 = can_refine(cv1, l1, r1);
    bool can_refine2 = can_refine(cv2, l2, r2);;
    if (!can_refine1 || !can_refine2) 
    {
      //Failed in the bounded approach - stop the subdivision.
      //Construct output with can_refine=false...  We don't really need a
      // meaningful bbox since we can't refine, so we take cv1 bbox (which
      // certainly contains the point).
      Bez_point_bbox bbox1;

      cp_bbox(cv1, bbox1);
      intersection_pairs.push_back
        (Bound_pair (Bez_point_bound (cv1, l1, r1,
                                      Bez_point_bound::INTERSECTION_PT, false),
                     Bez_point_bound (cv2, l2, r2,
                                      Bez_point_bound::INTERSECTION_PT, false),
                     bbox1));

      return;
    }

    Bez_point_bbox bbox1, bbox2;

    //_construct_bboxes (cv1, cv2,
    //                   bbox1, bbox2);
    cp_bbox(cv1, bbox1);
    cp_bbox(cv2, bbox2);

    // The bounding boxes do not overlap - return.
    if (!bbox1.Overlaps(bbox2))
    {
      //iddo debug
      //std::cout << "boxes do not overlap - curve do not intersect"  << std::endl;
      return;
    }
    else // Debug print
    {
      //std::cout << "boxes overlap - continuing..." << std::endl;
    }

    Vector_2 dir1, dir2, leftMost1, rightMost1, leftMost2, rightMost2;

    bool cv1SpanOK = true;
	  bool cv2SpanOK = true;
	  bool spansOverlap = false;

    if (should_check_span)
    {
      cv1SpanOK = _CrvTanAngularSpan(cv1, dir1, leftMost1, rightMost1);
	    cv2SpanOK = _CrvTanAngularSpan(cv2, dir2, leftMost2, rightMost2);
      if (cv1SpanOK && cv2SpanOK)
      {
  	    spansOverlap = _AngularSpansOverlap(leftMost1, rightMost1, leftMost2, rightMost2);
      }
      else
      {
        spansOverlap = true;
      }
    }

    if (!spansOverlap)
    {
      // Debug print.
      //std::cout << "Spans do NOT Overlap." << std::endl;

      // Debug print
      /*
      std::cout << "CV1:\n";
      Control_point_vec::const_iterator debug_iter = cv1.begin();
      while (debug_iter != cv1.end())
      {
        std::cout << "[" << *debug_iter << "]";
        ++debug_iter;
      }
      std::cout << std::endl;
      std::cout << "CV2:\n";
      debug_iter = cv2.begin();
      while (debug_iter != cv2.end())
      {
        std::cout << "[" << *debug_iter << "]";
        ++debug_iter;
      }
      std::cout << std::endl;
      */

      // Checking for endpoint intersections.
      bool endpoint_coincidence = _endpoints_coincide(cv1, l1, r1,
			        cv2, l2, r2,
              endpoint_intersections,
              intersection_pairs);

      if (endpoint_coincidence)
      {
        // Debug print
        //std::cout << "endpoint coincidence intersection." << std::endl; 
        //std::cout << "curves intersect! intersectionIntervals.size() == " << intersection_pairs.size() << std::endl;

        //Since spans do not overlap and there is an endpoint coincidence
        //there is no need for further checking (the endpoint is the intersection).
        return;
      }

      const Point_2& s1 = cv1.front();
      const Point_2& t1 = cv1.back();
      const Point_2& s2 = cv2.front();
      const Point_2& t2 = cv2.back();

      Line_2 skew1a, skew1b;
      _cvSkewBbox(cv1, skew1a, skew1b);

      Line_2 skew2a, skew2b;
      _cvSkewBbox(cv2, skew2a, skew2b);

      // Check for opposite orientations between: 
      // (skew2a, s1)(skew2a, t1)
      // (Skew2b, s1)(Skew2b, t1)
      // and vice versa
      //RI: TODO not use < 0 (count on numeric value) on Oriented_side
      Oriented_side or_2a_s1 = skew2a.oriented_side(s1);
      Oriented_side or_2a_t1 = skew2a.oriented_side(t1);

      Oriented_side or_2b_s1 = skew2b.oriented_side(s1);
      Oriented_side or_2b_t1 = skew2b.oriented_side(t1);

      bool s1_t1_are_opposite = ((or_2a_s1*or_2a_t1 < 0) &&
								      (or_2b_s1*or_2b_t1 < 0));


      Oriented_side or_1a_s2 = skew1a.oriented_side(s2);
      Oriented_side or_1a_t2 = skew1a.oriented_side(t2);

      Oriented_side or_1b_s2 = skew1b.oriented_side(s2);
      Oriented_side or_1b_t2 = skew1b.oriented_side(t2);

      bool s2_t2_are_opposite = ((or_1a_s2*or_1a_t2 < 0) &&
								      (or_1b_s2*or_1b_t2 < 0));

      if (s1_t1_are_opposite && s2_t2_are_opposite) 
      {
        // Construct the bbox from intersection of skew lines.
        Bez_point_bbox inter_bbox;
        Control_point_vec aux_vec;
        Object  res;
        Point_2 p;

        res = intersection (skew1a, skew2a);
        if (!assign(p, res)) 
        {
          CGAL_assertion(false);
        }
        aux_vec.push_back(p);

        res = intersection (skew1a, skew2b);
        if (!assign(p, res)) 
        {
          CGAL_assertion(false);
        }
        aux_vec.push_back(p);

        res = intersection (skew1b, skew2a);
        if (!assign(p, res)) 
        {
          CGAL_assertion(false);
        }
        aux_vec.push_back(p);

        res = intersection (skew1b, skew2b);
        if (!assign(p, res)) 
        {
          CGAL_assertion(false);
        }
        aux_vec.push_back(p);

        cp_bbox(aux_vec, inter_bbox);

        intersection_pairs.push_back(Bound_pair(
                                        Bez_point_bound(cv1, l1, r1, Bez_point_bound::INTERSECTION_PT, true),
                                        Bez_point_bound(cv2, l2, r2, Bez_point_bound::INTERSECTION_PT, true),
                                        inter_bbox
                                      )
                                    );

        // Debug print
        //std::cout << "curves intersect! intersectionIntervals.size() == " << intersection_pairs.size() << std::endl;
        // gets(aux);

	      return;		
      }
    }

    // Subdivide the two curves and recurse. 
    Control_point_vec cv1a, cv1b; 
    
    bisect_control_polygon_2 (cv1.begin(), cv1.end(),
                              std::back_inserter(cv1a),
                              std::front_inserter(cv1b));
    //DeCasteljau(cv1, 0.5, cv1a, cv1b);
    NT t_mid1 = NT(0.5) * (l1 + r1);

    Control_point_vec cv2a, cv2b; 
    
    bisect_control_polygon_2 (cv2.begin(), cv2.end(),
                              std::back_inserter(cv2a),
                              std::front_inserter(cv2b));
    //DeCasteljau(cv2, 0.5, cv2a, cv2b);
    NT t_mid2 = NT(0.5) * (l2 + r2);

    // Recursion:
    _CrvCrvInterAux (cv1a, l1, t_mid1, cv2a, l2, t_mid2, 
                     spansOverlap, endpoint_intersections,
                     intersection_pairs);
    _CrvCrvInterAux (cv1a, l1, t_mid1, cv2b, t_mid2, r2, 
                     spansOverlap, endpoint_intersections,
                     intersection_pairs);
    _CrvCrvInterAux (cv1b, t_mid1, r1, cv2a, l2, t_mid2, 
                     spansOverlap, endpoint_intersections,
                     intersection_pairs);
    _CrvCrvInterAux (cv1b, t_mid1, r1, cv2b, t_mid2, r2,
                     spansOverlap, endpoint_intersections,
                     intersection_pairs);
  }

};

//TODO - These classes can both go outside to another file.
//if we decide to take some of the functionality to a template algo file.
template <class _NT> 
struct _Bez_point_bound {
  enum Point_type {RATIONAL_PT, VERTICAL_TANGENCY_PT, INTERSECTION_PT, UNDEFINED};

  typedef _NT                    NT;
  typedef typename Cartesian<NT>::Point_2 Point_2;
  typedef std::deque<Point_2>                   Control_point_vec;

  Control_point_vec bounding_polyline;

  NT t_min; // iddo: in the future, maybe a better rep for these numbers
  NT t_max; // possibly a mantissa+exp representation.
  Point_type point_type;
  bool can_refine;

  _Bez_point_bound() : bounding_polyline(),
    t_min(), t_max(), point_type(UNDEFINED), can_refine(false)
  {}

  _Bez_point_bound(const _Bez_point_bound& o) : bounding_polyline(o.bounding_polyline),
    t_min(o.t_min), t_max(o.t_max), point_type(o.point_type), can_refine(o.can_refine)
  {}

  _Bez_point_bound(const Control_point_vec& bp,
    const NT& min_t, const NT& max_t, 
    Point_type pt, bool cr) : bounding_polyline(bp),
    t_min(min_t), t_max(max_t), point_type(pt), can_refine(cr)
  {}

};

template <class NT> 
struct _Bez_point_bbox {
  NT min_x, max_x, min_y, max_y;

  _Bez_point_bbox() : min_x(), max_x(), min_y(), max_y() {}
  _Bez_point_bbox(const _Bez_point_bbox& other) : 
    min_x(other.min_x), max_x(other.max_x), 
    min_y(other.min_y), max_y(other.max_y) 
  {}

  _Bez_point_bbox(const NT& _min_x, const NT& _max_x,
    const NT& _min_y, const NT& _max_y) : 
  min_x(_min_x), max_x(_max_x), min_y(_min_y), max_y(_max_y) {}

  void AddS(const _Bez_point_bbox& other) {
    min_x = std::min(min_x, other.min_x);
    min_y = std::min(min_y, other.min_y);
    max_x = std::max(max_x, other.max_x);
    max_y = std::max(max_y, other.max_y);
  }

  _Bez_point_bbox Add(const _Bez_point_bbox& other) const {
    _Bez_point_bbox res(other);
    res.AddS(*this);
    return res;
  }

  bool Overlaps(const _Bez_point_bbox& other) const
  {
    if (max_x < other.min_x ||
        max_y < other.min_y ||
        other.max_x < min_x ||
        other.max_y < min_y)
        return false;
    return true;
  }

  bool Overlaps_x(const _Bez_point_bbox& other) const
  {
    if (max_x < other.min_x ||
        other.max_x < min_x)
        return false;
    return true;
  }

};

CGAL_END_NAMESPACE

#endif //CGAL_BEZIER_BOUNDING_RATIONAL_TRAITS_H