bezier_point_2.h

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

  // Try to handle the comparison using the x-range of the bounding boxes.
  bool    can_refine1 = true;
  bool    can_refine2 = true;

  do
  {
    // Compare the x-ranges of the bounding boxes.
    if (CGAL::compare (_bbox.max_x, pt._bbox.min_x) == SMALLER)
      return (SMALLER);
    
    if (CGAL::compare (_bbox.min_x, pt._bbox.max_x) == LARGER)
      return (LARGER);

    // Check if only one of the points is exactly represented.
    if (is_exact())
    {
      // Compare the exact x-coordinate to pt's bounding box.
      if (CGAL::compare (*p_alg_x, Algebraic (pt._bbox.min_x)) == SMALLER)
        return (SMALLER);
      if (CGAL::compare (*p_alg_x, Algebraic (pt._bbox.max_x)) == LARGER)
        return (LARGER);
    }

    if (pt.is_exact())
    {
      // Compare the bounding box to pt's exact x-coordinate.
      if (CGAL::compare (Algebraic (_bbox.max_x), *(pt.p_alg_x)) == SMALLER)
        return (SMALLER);
      if (CGAL::compare (Algebraic (_bbox.min_x), *(pt.p_alg_x)) == LARGER)
        return (LARGER);
    }

    // Try to refine the representation of the points.
    // If we cannot refine any more, compute the exact representation of the
    // point.
    if (! is_exact())
    {
      if (can_refine1)
        can_refine1 = this->_refine();
      
      if (! can_refine1)
        this->_make_exact (cache);
    }
    else
    {
      can_refine1 = false;
    }

    if (! pt.is_exact())
    {
      if (can_refine2)
        can_refine2 = pt._refine();
      
      if (! can_refine2)
        pt._make_exact (cache);
    }
    else
    {
      can_refine2 = false;
    }

  } while (can_refine1 || can_refine2);

  // If we reached here, we have an exact representation of both points, and
  // we simply have to compare their x-coordinates.
  return (CGAL::compare (*p_alg_x, *(pt.p_alg_x)));
}

// ---------------------------------------------------------------------------
// Compare the two points xy-lexicographically.
//
template <class RatKer, class AlgKer, class NtTrt, class BndTrt>
Comparison_result
_Bezier_point_2_rep<RatKer, AlgKer, NtTrt, BndTrt>::compare_xy
         (Self& pt,
          Bezier_cache& cache)
{
  // First compare the x-coordinates.
  const Comparison_result   res_x = this->compare_x (pt,
                                                     cache);

  if (res_x != EQUAL)
    return (res_x);

  // Handle rational points separately.
  if (is_rational() && pt.is_rational())
  {
    return (CGAL::compare (*p_rat_y, *(pt.p_rat_y)));
  }

  // If we reached here, the two point should have already been computed in
  // an exact manner (otherwise we could not have determined the equality of
  // their x-coordinates).
  CGAL_assertion (is_exact() && pt.is_exact());
  return (CGAL::compare (*p_alg_y, *(pt.p_alg_y)));
}

// ---------------------------------------------------------------------------
// Determine the vertical position of the point with respect to a given curve.
//
template <class RatKer, class AlgKer, class NtTrt, class BndTrt>
Comparison_result
_Bezier_point_2_rep<RatKer, AlgKer, NtTrt, BndTrt>::vertical_position
        (const Control_point_vec& cp,
         const BoundNT& t_min,
         const BoundNT& t_max)
{
  Bounding_traits       bound_tr;
  std::list<Subcurve>   subcurves;

  subcurves.push_back(Subcurve(cp,0,1));

  // Iddo: maybe make here a comparison with an incrementor (e.g., loop 4 
  //       times first..)
  bool                  can_refine_pt = true;
  bool                  can_refine_cv = true;
  Bez_point_bbox        scv_bbox;

  while (can_refine_pt || can_refine_cv)
  {
    // Iddo: Implement Algebraic comparisons with scv_bbox if is_exact
    //       (or use the current bbox of the exact rep as is done now).

    // Go over the list of subcurves and consider only those lying in the
    // given [t_min, t_max] bound.
    typename std::list<Subcurve>::iterator  iter = subcurves.begin();
    bool                                    was_split = false;
    bool                                    is_fully_in_t_range;

    while (iter != subcurves.end())
    {
      if (CGAL::compare (iter->r, t_min) == SMALLER ||
          CGAL::compare (iter->l, t_max) == LARGER)
      {
        // Subcurve out of bounds - erase it and continue to next subcurve.
        subcurves.erase(iter++);
        continue;
      }

      // Construct the bounding box of the subcurve and compare it to
      // the bounding box of the point.
      // Iddo: In the future it might be not a bez_point_bbox
      //       but a bbox<Rational>.
      bound_tr.cp_bbox (iter->cp, scv_bbox);

      if (! _bbox.Overlaps_x(scv_bbox))
      {
        // Subcurve out of x bounds - erase it and continue to next subcurve.
        subcurves.erase(iter++);
        continue;
      }

      is_fully_in_t_range = (CGAL::compare (iter->l, t_min) != SMALLER) &&
                            (CGAL::compare (iter->r, t_max) != LARGER);

      if (_bbox.Overlaps(scv_bbox) || ! is_fully_in_t_range)
      {
        // Iddo: This is a special case of subdividing the curve
        //       not as part of an Originator. Think again if the can_refine
        //       and de Casteljau should not be different here!!
        can_refine_cv = bound_tr.can_refine (iter->cp, iter->l, iter->r);

        if (! can_refine_pt && ! can_refine_cv)
          // It is not possible to refine the point or the subcurve anymore:
          return (EQUAL);

        if (! can_refine_cv)
        {
          ++iter;
          continue;
        }

        // Subdivide the current subcurve and replace iter with the two
        // resulting subcurves.
        Control_point_vec  left_cp, right_cp;

        bisect_control_polygon_2 (iter->cp.begin(), iter->cp.end(),
                                  std::back_inserter(left_cp),
                                  std::front_inserter(right_cp));

        //bound_tr.DeCasteljau (iter->cp, 0.5, left, right);

        // Insert two new subcurves before iter and remove iter.
        BoundNT   t_mid = BoundNT(0.5) * (iter->l + iter->r);
        
        subcurves.insert (iter, Subcurve (left_cp, iter->l, t_mid));
        subcurves.insert (iter, Subcurve (right_cp, t_mid, iter->r));
        subcurves.erase(iter++);

        was_split = true;
        continue;
      }
      else
      {
        // We found a subcurve whose x-range contains our point, but whose
        // bounding box is disjoint from the bounding box of the point.
        // We can therefore compare the y-positions of the bounding boxes.
        CGAL_assertion (! _bbox.Overlaps (scv_bbox) &&
                        _bbox.Overlaps_x (scv_bbox) &&
                        is_fully_in_t_range);

        return (CGAL::compare (_bbox.max_y, scv_bbox.max_y));
      }

      // If we got here without entering one of the clauses above,
      // then iter has not been incremented yet.
      ++iter; 
    }

    // If we reached here without splitting a subcurve, then we have a
    // subcurve whose bbox (given by scv_bbox) is totally below or totally
    // above the bounding box of the point.
    if (! was_split)
      break;

    // Try to refine the bounding box of the point.
    if (! is_exact())
      can_refine_pt = _refine();

    if (! can_refine_pt && ! can_refine_cv)
      // It is not possible to refine the point or the subcurve anymore:
      return (EQUAL);
  }

  // If we reached here, we cannot refine any more:
  return (EQUAL);
}

// ---------------------------------------------------------------------------
// Refine the bounds of the point.
//
template <class RatKer, class AlgKer, class NtTrt, class BndTrt>
bool _Bezier_point_2_rep<RatKer, AlgKer, NtTrt, BndTrt>::_refine ()
{
  // Get the first originator and consider the point type.
  Bounding_traits  bound_tr;
  Originator&      orig1 = *(_origs.begin());

  if (orig1.point_bound().point_type == Bez_point_bound::VERTICAL_TANGENCY_PT)
  {
    CGAL_assertion(_origs.size() == 1);

    // Refine the vertical tangency point.
    std::pair<Bez_point_bound, Bez_point_bbox>  refined_tang_pt;

    bound_tr.refine_tangency_point (orig1.point_bound().bounding_polyline,
                                    orig1.point_bound().t_min, 
                                    orig1.point_bound().t_max,
                                    refined_tang_pt);

    if (! refined_tang_pt.first.can_refine)
      // Indicate that it was not possible to refine the point.
      return (false);

    // Update the originator and the bounding box of the point.
    orig1.update_point_bound (refined_tang_pt.first);
    _bbox = refined_tang_pt.second;
    return (true);
  }
  
  if (orig1.point_bound().point_type == Bez_point_bound::INTERSECTION_PT)
  {
    CGAL_assertion(_origs.size() == 2);

    // Obtain the other curve that originates the intersection point and use
    // it to refine its reprsentation.
    Orig_iter    org_it = _origs.begin();
    ++org_it;
    Originator&  orig2 = *org_it;

    typename Bounding_traits::Bound_pair   intersect_pt (orig1.point_bound(),
                                                         orig2.point_bound(),
                                                         _bbox);
    typename Bounding_traits::Bound_pair   refined_pt;

    bound_tr.refine_intersection_point (intersect_pt, refined_pt);

    if (! refined_pt.bound1.can_refine || ! refined_pt.bound2.can_refine)
      // Indicate that it was not possible to refine the point.
      return (false);

    // Update the originators and the bounding box of the point.
    orig1.update_point_bound (refined_pt.bound1);
    orig2.update_point_bound (refined_pt.bound2);
    _bbox = refined_pt.bbox;
    return (true);
  }

  // We should never reach here:
  //CGAL_assertion (false);
  return (false);
}

// ---------------------------------------------------------------------------
// Compute the point in an exact manner.
//
template <class RatKer, class AlgKer, class NtTrt, class BndTrt>
void _Bezier_point_2_rep<RatKer, AlgKer, NtTrt, BndTrt>::_make_exact
        (Bezier_cache& cache)
{                
  if (is_exact())
    return;

  /* Ron: For informational purposes ...
  std::cout << "***** MAKE EXACT (" << _origs.size() << " ORIGINATORS) *****"
            << std::endl;
  */

  // Check if the point is a vertical tangency point of the originator.
  if (_origs.size() == 1)
  {
    Orig_iter   org_it = _origs.begin();

    CGAL_assertion (org_it->point_bound().point_type ==
                    Bez_point_bound::VERTICAL_TANGENCY_PT);

    // Compute (using the cache) the vertical tangency parameters of
    // the current curve.
    const typename Bezier_cache::Vertical_tangency_list&          vt_list =
      cache.get_vertical_tangencies (org_it->curve().id(),
                                     org_it->curve().x_polynomial(),
                                     org_it->curve().x_norm());
    typename Bezier_cache::Vertical_tangency_iter                 vt_it;

    // Look for a parameter within the range of the bounding interval.
    const Algebraic  t_min (org_it->point_bound().t_min);
    const Algebraic  t_max (org_it->point_bound().t_max);

    for (vt_it = vt_list.begin(); vt_it != vt_list.end(); ++vt_it)
    {
      if (CGAL::compare (*vt_it, t_min) != SMALLER &&
          CGAL::compare (*vt_it, t_max) != LARGER)
      {
        // Update the originator.
        Originator&   orig = const_cast<Originator&> (*org_it);

        orig.set_parameter (*vt_it);
        
        // Evaluate the curve at the given parameter value.
        const Alg_point_2&   p = org_it->curve() (*vt_it);

        p_alg_x = new Algebraic (p.x());
        p_alg_y = new Algebraic (p.y());

        // Update the bounding box.
        Nt_traits                         nt_traits;
        const std::pair<double, double>&  x_bnd = 
                                        nt_traits.double_interval (*p_alg_x);
        const std::pair<double, double>&  y_bnd = 
                                        nt_traits.double_interval (*p_alg_y);

        _bbox.min_x = x_bnd.first;
        _bbox.max_x = x_bnd.second;
        _bbox.min_y = y_bnd.first;
        _bbox.max_y = y_bnd.second;

        return;
      }
    }

    // We should never reach here:
    CGAL_assertion (false);
  }

  // In this case the point is an intersection between two originating
  // curves. Compute the intersections between these curves in the parameter
  // space.
  CGAL_assertion (_origs.size() == 2);

  Orig_iter    org_it1 = _origs.begin();
  Orig_iter    org_it2 = org_it1;

  ++org_it2;

  if (org_it1->curve().id() > org_it2->curve().id())
  {
    // Make sure that org_it1 refers to a curve with smaller ID than org_it2.
    --org_it2;
    ++org_it1;
  }

  Originator&  orig1 = const_cast<Originator&> (*org_it1);
  Originator&  orig2 = const_cast<Originator&> (*org_it2);
  bool         do_ovlp;

  const typename Bezier_cache::Intersection_list&           intr_list =
    cache.get_intersections (orig1.curve().id(),
                             orig1.curve().x_polynomial(),
                             orig1.curve().x_norm(),
                             orig1.curve().y_polynomial(),
                             orig1.curve().y_norm(),
                             orig2.curve().id(),
                             orig2.curve().x_polynomial(),
                             orig2.curve().x_norm(),
                             orig2.curve().y_polynomial(),
                             orig2.curve().y_norm(),
                             do_ovlp);
  typename Bezier_cache::Intersection_iter                  intr_it;
                             
  CGAL_assertion (! do_ovlp);
                    
  // Look for a parameter pair within the ranges of the bounding intervals.
  Nt_traits            nt_traits;
  const Algebraic      s_min = nt_traits.convert (orig1.point_bound().t_min);
  const Algebraic      s_max = nt_traits.convert (orig1.point_bound().t_max);
  const Algebraic      t_min = nt_traits.convert (orig2.point_bound().t_min);
  const Algebraic      t_max = nt_traits.convert (orig2.point_bound().t_max);

  for (intr_it = intr_list.begin(); intr_it != intr_list.end(); ++intr_it)
  {
    if (CGAL::compare (intr_it->s, s_min) != SMALLER &&
        CGAL::compare (intr_it->s, s_max) != LARGER &&
        CGAL::compare (intr_it->t, t_min) != SMALLER &&
        CGAL::compare (intr_it->t, t_max) != LARGER)
    {
      // Update the originators.
      orig1.set_parameter (intr_it->s);
      orig2.set_parameter (intr_it->t);

      // Set the exact point coordinates.
      p_alg_x = new Algebraic (intr_it->x);
      p_alg_y = new Algebraic (intr_it->y);

      // Update the bounding box.
      const std::pair<double, double>&  x_bnd = 
                                          nt_traits.double_interval (*p_alg_x);
      const std::pair<double, double>&  y_bnd = 
                                          nt_traits.double_interval (*p_alg_y);

      _bbox.min_x = x_bnd.first;
      _bbox.max_x = x_bnd.second;
      _bbox.min_y = y_bnd.first;
      _bbox.max_y = y_bnd.second;

      return;
    }
  }

  // We should never reach here:
  CGAL_assertion (false);
}

CGAL_END_NAMESPACE

#endif