bezier_x_monotone_2.h

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

                                                                   bound2));
  }
  
  // Iddo: Implement slope comparison at a Rational point
  //       (using de Casteljau).
  
  // The slope of (X(t), Y(t)) at t0 is Y'(t)/X'(t).
  // Compute the slope of (*this).
  // Note that we take special care of the case X'(t) = 0, when the tangent
  // is vertical and its slope is +/- oo.
  CGAL_assertion (org1->has_parameter() && org2->has_parameter());
  
  Nt_traits         nt_traits;
  const Algebraic&  t1 = org1->parameter();
  Polynomial        derivX = nt_traits.derive (_curve.x_polynomial());
  Polynomial        derivY = nt_traits.derive (_curve.y_polynomial());
  Algebraic         numer1 = nt_traits.evaluate_at (derivY, t1) *
                             nt_traits.convert (_curve.x_norm());
  Algebraic         denom1 = nt_traits.evaluate_at (derivX, t1) *
                             nt_traits.convert (_curve.y_norm());
  CGAL::Sign        inf_slope1 = CGAL::ZERO;
  Algebraic         slope1;

  if (CGAL::sign (denom1) == CGAL::ZERO)
  {
    inf_slope1 = CGAL::sign (numer1);

    // If also Y'(t) = 0, we cannot perform the comparison:
    if (inf_slope1 == CGAL::ZERO)
      return (EQUAL);
  }
  else
  {
    slope1 = numer1 / denom1;
  }

  // Compute the slope of the other subcurve.
  const Algebraic&  t2 = org2->parameter();
  derivX = nt_traits.derive (cv._curve.x_polynomial());
  derivY = nt_traits.derive (cv._curve.y_polynomial());
  Algebraic         numer2 = nt_traits.evaluate_at (derivY, t2) *
                             nt_traits.convert (cv._curve.x_norm());
  Algebraic         denom2 = nt_traits.evaluate_at (derivX, t2) *
                             nt_traits.convert (cv._curve.y_norm());
  CGAL::Sign        inf_slope2 = CGAL::ZERO;
  Algebraic         slope2;

  if (CGAL::sign (denom2) == CGAL::ZERO)
  {
    inf_slope2 = CGAL::sign (numer2);

    // If also Y'(t) = 0, we cannot perform the comparison:
    if (inf_slope2 == CGAL::ZERO)
      return (EQUAL);
  }
  else
  {
    slope2 = numer2 / denom2;
  }

  // Handle the comparison when one slope (or both) is +/- oo.
  if (inf_slope1 == CGAL::POSITIVE)
    return (inf_slope2 == CGAL::POSITIVE ? EQUAL : LARGER);
  
  if (inf_slope1 == CGAL::NEGATIVE)
    return (inf_slope2 == CGAL::NEGATIVE ? EQUAL : SMALLER);

  if (inf_slope2 == CGAL::POSITIVE)
    return (SMALLER);

  if (inf_slope2 == CGAL::NEGATIVE)
    return (LARGER);

  // Compare the slopes.
  return (CGAL::compare (slope1, slope2));
}

// ---------------------------------------------------------------------------
// Get the range of t-value over which the subcurve is defined.
//
template <class RatKer, class AlgKer, class NtTrt, class BndTrt>
std::pair<typename _Bezier_x_monotone_2<RatKer, AlgKer, NtTrt,
                                        BndTrt>::Algebraic, 
          typename _Bezier_x_monotone_2<RatKer, AlgKer, NtTrt,
                                        BndTrt>::Algebraic>
_Bezier_x_monotone_2<RatKer, AlgKer, NtTrt, BndTrt>::_t_range
        (Bezier_cache& cache) const
{
  Originator_iterator  ps_org = _ps.get_originator(_curve);
  CGAL_assertion(ps_org != _ps.originators_end());
  
  Originator_iterator  pt_org = _pt.get_originator(_curve);
  CGAL_assertion(pt_org != _pt.originators_end());
  
  // Make sure that the two endpoints are exact.
  if (! ps_org->has_parameter())
    _ps.make_exact (cache);
  
  if (! pt_org->has_parameter())
    _pt.make_exact (cache);
  
  return (std::make_pair (ps_org->parameter(),
                          pt_org->parameter()));
}

// ---------------------------------------------------------------------------
// Compare the relative y-position of two x-monotone subcurve to the right
// (or to the left) of their intersection point, whose multiplicity is
// greater than 1.
//
template <class RatKer, class AlgKer, class NtTrt, class BndTrt>
Comparison_result
_Bezier_x_monotone_2<RatKer, AlgKer, NtTrt, BndTrt>::_compare_to_side
        (const Self& cv,
         const Point_2& p,
         bool to_right,
         Bezier_cache& cache) const
{
  // Get the intersection points of the two curves from the cache. Note that
  // we make sure that the ID of this->_curve is smaller than of cv's curve ID.
  const bool             no_swap_curves = (_curve.id() < cv._curve.id());
  bool                   do_ovlp;
  const Intersect_list&  inter_list =
    (no_swap_curves ?
     (cache.get_intersections (_curve.id(),
                               _curve.x_polynomial(), _curve.x_norm(),
                               _curve.y_polynomial(), _curve.y_norm(),
                               cv._curve.id(),
                               cv._curve.x_polynomial(), cv._curve.x_norm(),
                               cv._curve.y_polynomial(), cv._curve.y_norm(),
                               do_ovlp)) :
     (cache.get_intersections (cv._curve.id(),
                               cv._curve.x_polynomial(), cv._curve.x_norm(),
                               cv._curve.y_polynomial(), cv._curve.y_norm(),
                               _curve.id(),
                               _curve.x_polynomial(), _curve.x_norm(),
                               _curve.y_polynomial(), _curve.y_norm(),
                               do_ovlp)));

  if (do_ovlp)
    return (EQUAL);
  
  // Get the parameter value for the point p.
  Originator_iterator          org = p.get_originator(_curve);
  
  CGAL_assertion (org != p.originators_end());
  CGAL_assertion (org->has_parameter());
  
  const Algebraic&             t0 = org->parameter();

  // Get the parameter range of the curve.
  const std::pair<Algebraic,
                  Algebraic>&  range = _t_range (cache);
  const Algebraic&             t_src = range.first;
  const Algebraic&             t_trg = range.second;

  // Find the next intersection point that lies to the right of p.
  Intersect_iter               iit;
  Algebraic                    next_t;
  Comparison_result            res = CGAL::EQUAL;
  bool                         found = false;
  
  for (iit = inter_list.begin(); iit != inter_list.end(); ++iit)
  {
    // Check if the current point lies to the right (left) of p. We do so by
    // considering its originating parameter value s (or t, if we swapped
    // the curves).
    const Algebraic&     t = (no_swap_curves ? (iit->s) : iit->t);
    
    res = CGAL::compare (t, t0);
    if ((to_right && ((_inc_to_right && res == LARGER) ||
                      (! _inc_to_right && res == SMALLER))) ||
        (! to_right && ((_inc_to_right && res == SMALLER) ||
                        (! _inc_to_right && res == LARGER))))
    {
      if (! found)
      {
        next_t = t;
        found = true;
      }
      else
      {
        // If we have already located an intersection point to the right
        // (left) of p, choose the leftmost (rightmost) of the two points.
        res = CGAL::compare (t, next_t);
        if ((to_right && ((_inc_to_right && res == SMALLER) ||
                          (! _inc_to_right && res == LARGER))) ||
            (! to_right && ((_inc_to_right && res == LARGER) ||
                            (! _inc_to_right && res == SMALLER))))
        {
          next_t = t;
        }
      }
    }
  }
  
  // If the next intersection point occurs before the right (left) endpoint
  // of the subcurve, keep it. Otherwise, take the parameter value at
  // the endpoint.
  if (found)
  {
    if (to_right == _dir_right)
      res = CGAL::compare (t_trg, next_t);
    else
      res = CGAL::compare (t_src, next_t);
  }
  
  if (! found ||
      (to_right && ((_inc_to_right && res == SMALLER) ||
                    (! _inc_to_right && res == LARGER))) ||
      (! to_right && ((_inc_to_right && res == LARGER) ||
                      (! _inc_to_right && res == SMALLER))))
  {
    next_t = ((to_right == _dir_right) ? t_trg : t_src);
  }
  
  // Find a rational value between t0 and t_next. Using this value, we
  // a point with rational coordinates on our subcurve. We also locate a point
  // on the other curve with the same x-coordinates.
  Nt_traits           nt_traits;
  const Rational&     mid_t = nt_traits.rational_in_interval (t0, next_t);
  const Rat_point_2&  q1 = _curve (mid_t);
  const Algebraic&    y2 = cv._get_y (q1.x(), cache);

  // We now just have to compare the y-coordinates of the two points we have
  // computed.
  return (CGAL::compare (nt_traits.convert (q1.y()), y2));
}

// ---------------------------------------------------------------------------
// Approximate the intersection points between the two given curves.
//
template <class RatKer, class AlgKer, class NtTrt, class BndTrt>
bool _Bezier_x_monotone_2<RatKer, AlgKer, NtTrt,
                          BndTrt>::_approximate_intersection_points
        (const Curve_2& B1,
         const Curve_2& B2,
         std::list<Point_2>& inter_pts) const
{
  inter_pts.clear();

  // Make local copies of the control polygons.
  typename Bounding_traits::Control_point_vec  cp1, cp2;
  
  std::copy (B1.control_points_begin(), B1.control_points_end(),
             std::back_inserter (cp1));
  std::copy (B2.control_points_begin(), B2.control_points_end(),
             std::back_inserter(cp2));

  // Use the bounding traits to isolate the intersection points.
  Bounding_traits                              bound_tr;
  typename Bounding_traits::Bound_pair_lst     inter_pairs;

  bound_tr.CrvCrvInter (cp1, cp2, 
                        inter_pairs);
  
  // Construct the approximated points.
  typename Bounding_traits::Bound_pair_lst::const_iterator iter;
  
  for (iter = inter_pairs.begin(); iter != inter_pairs.end(); ++iter)
  {
    const Bez_point_bound&  bound1 = iter->bound1; 
    const Bez_point_bound&  bound2 = iter->bound2; 
    const Bez_point_bbox&   bbox = iter->bbox; 
    
    // In case it is impossible to further refine the point, stop here.
    if (! bound1.can_refine || ! bound2.can_refine)
      return (false);
    
    // Create the approximated intersection point.
    Point_2                 pt;
    
    if (bound1.point_type == Bounding_traits::Bez_point_bound::RATIONAL_PT &&
        bound2.point_type == Bounding_traits::Bez_point_bound::RATIONAL_PT)
    {
      CGAL_assertion (CGAL::compare (bound1.t_min, bound1.t_max) == EQUAL); 
      CGAL_assertion (CGAL::compare (bound2.t_min, bound2.t_max) == EQUAL); 
      Rational  t1 = bound1.t_min;
      Rational  t2 = bound2.t_min;

      pt = Point_2 (B1, t1);
      pt.add_originator (Originator (B2, t2));
    }
    else
    {
      pt.add_originator (Originator (B1, bound1));
      pt.add_originator (Originator (B2, bound2));
    }
    pt.set_bbox (bbox);
    
    inter_pts.push_back (pt);
  }

  // The approximation process ended OK.
  return (true);
}

// ---------------------------------------------------------------------------
// Compute the intersections with the given subcurve.
//
template <class RatKer, class AlgKer, class NtTrt, class BndTrt>
bool _Bezier_x_monotone_2<RatKer, AlgKer, NtTrt, BndTrt>::_intersect
    (const Self& cv,
     Intersection_map& inter_map,
     Bezier_cache& cache,
     std::vector<Intersection_point_2>& ipts,
     Self& /* ovlp_cv */) const
{
  CGAL_precondition (_curve.id() < cv._curve.id());

  ipts.clear();

  // Construct the pair of curve IDs and look for it in the intersection map.
  Curve_pair                 curve_pair (_curve.id(), cv._curve.id());
  Intersection_map_iterator  map_iter = inter_map.find (curve_pair);
  std::list<Point_2>         inter_pts;
  bool                       app_ok = true;

  if (map_iter != inter_map.end())
  {
    // Get the intersection points between the two supporting curves as stored
    // in the map.
    inter_pts = map_iter->second;
  }
  else
  {
    // Approximate the intersection points and store them in the map.
    app_ok = _approximate_intersection_points (_curve,
                                               cv._curve,
                                               inter_pts);

    if (app_ok)
      inter_map[curve_pair] = inter_pts;
  }

  // Try to approximate the intersection points.
  bool                in_range1, in_range2;
  bool                correct_res;

  if (app_ok)
  {
    // If the approximation went OK, then we know that we just have simple
    // intersection points (with multiplicity 1). We go over the points
    // and report the ones lying in the parameter ranges of both curves.
    typename std::list<Point_2>::iterator  pit;
    
    for (pit = inter_pts.begin(); pit != inter_pts.end(); ++pit)
    {
      // Check if the point is in the range of this curve - first using
      // its parameter bounds, and if we fail we perform an exact check.
      in_range1 = _is_in_range (*pit, correct_res);
      
      if (! correct_res)
      {
        if (! pit->is_exact())
          pit->make_exact (cache);
        
        Originator_iterator  p_org = pit->get_originator (_curve);
        CGAL_assertion (p_org != pit->originators_end());
        
        in_range1 = _is_in_range (p_org->parameter(), cache);
      }
      
      if (! in_range1)
        continue;

      // Check if the point is in the range of the other curve - first using
      // its parameter bounds, and if we fail we perform an exact check.
      in_range2 = cv._is_in_range (*pit, correct_res);
      
      if (! correct_res)
      {
        if (! pit->is_exact())
          pit->make_exact (cache);
        
        Originator_iterator  p_org = pit->get_originator (cv._curve);
        CGAL_assertion (p_org != pit->originators_end());
        
        in_range2 = cv._is_in_range (p_o