rational_arc_2.h

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

   * Compare the x-coordinate of a vertical asymptotes of the two arcs.
   */
  Comparison_result compare_ends (Curve_end ind1,
                                  const Self& arc, Curve_end ind2) const
  {
    // Get the x-coordinates of the first vertical asymptote.
    Algebraic                x1;

    if (ind1 == MIN_END)
    {
      CGAL_assertion (left_boundary_in_x() == NO_BOUNDARY &&
                      left_boundary_in_y() != NO_BOUNDARY);
      if ((_info & IS_DIRECTED_RIGHT) != 0)
        x1 = _ps.x();
      else
        x1 = _pt.x();
    }
    else
    {
      CGAL_assertion (right_boundary_in_x() == NO_BOUNDARY &&
                      right_boundary_in_y() != NO_BOUNDARY);
      if ((_info & IS_DIRECTED_RIGHT) != 0)
        x1 = _pt.x();
      else
        x1 = _ps.x();
    }

    // Get the x-coordinates of the second vertical asymptote.
    Algebraic                x2;

    if (ind2 == MIN_END)
    {
      CGAL_assertion (arc.left_boundary_in_x() == NO_BOUNDARY &&
                      arc.left_boundary_in_y() != NO_BOUNDARY);
      if ((arc._info & IS_DIRECTED_RIGHT) != 0)
        x2 = arc._ps.x();
      else
        x2 = arc._pt.x();
    }
    else
    {
      CGAL_assertion (arc.right_boundary_in_x() == NO_BOUNDARY &&
                      arc.right_boundary_in_y() != NO_BOUNDARY);
      if ((arc._info & IS_DIRECTED_RIGHT) != 0)
        x2 = arc._pt.x();
      else
        x2 = arc._ps.x();
    }

    // Compare the x-coordinates. In case they are not equal we are done.
    const Comparison_result  res = CGAL::compare (x1, x2);

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

    // If the x-coordinates of the asymptote are equal, but one arc is
    // defined to the left of the vertical asymptote and the other to its
    // right, we can easily determine the comparison result.
    if (ind1 == MAX_END && ind2 == MIN_END)
      return (SMALLER);
    else if (ind1 == MIN_END && ind2 == MAX_END)
      return (LARGER);

    // Determine the multiplicities of the pole for the two arcs.
    Nt_traits         nt_traits;
    int               mult1 = 1;
    Polynomial        p_der = nt_traits.derive (_denom);

    while (CGAL::sign (nt_traits.evaluate_at (p_der, x1)) == CGAL::ZERO)
    {
      mult1++;
      p_der = nt_traits.derive (p_der);
    }

    int               mult2 = 1;

    p_der = nt_traits.derive (arc._denom);
    while (CGAL::sign (nt_traits.evaluate_at (p_der, x1)) == CGAL::ZERO)
    {
      mult2++;
      p_der = nt_traits.derive (p_der);
    }

    if (mult1 != mult2)
    {
      if (ind1 == MIN_END)
      {
        // In case we compare to the right of the pole, the curve with larger
        // multiplicity is to the left.
        return ((mult1 > mult2) ? SMALLER : LARGER);
      }
      else
      {
        // In case we compare to the left of the pole, the curve with larger
        // multiplicity is to the right.
        return ((mult1 < mult2) ? SMALLER : LARGER);
      }
    }

    // The pole multiplicities are the same. If we denote the two rational
    // functions we have as P1(x)/Q1(x) and P2(X)/Q2(x), then we can divide
    // Q1(x) and Q2(x). As Q1(x') = Q2(x') = 0 (where x' is our current pole),
    // the zeros cancel themselves and we consider the value of our function
    // at x'.
    int            deg_q1 = nt_traits.degree (_denom);
    Rat_vector     coeffs_q1 (deg_q1 + 1);
    int            deg_q2 = nt_traits.degree (arc._denom);
    Rat_vector     coeffs_q2 (deg_q2 + 1);
    int            k;

    for (k = 0; k <= deg_q1; k++)
      coeffs_q1[k] = Rational (nt_traits.get_coefficient (_denom, k), 1);

    for (k = 0; k <= deg_q2; k++)
      coeffs_q2[k] = Rational (nt_traits.get_coefficient (arc._denom, k), 1);

    Polynomial         norm_q1, norm_q2;
    nt_traits.construct_polynomials (&(coeffs_q1[0]), deg_q1,
                                     &(coeffs_q2[0]), deg_q2,
                                     norm_q1, norm_q2);

    const Algebraic    val_q1 = CGAL::abs (nt_traits.evaluate_at (norm_q1, x1));
    const Algebraic    val_q2 = CGAL::abs (nt_traits.evaluate_at (norm_q2, x1));
    Comparison_result  val_res = CGAL::compare (val_q1, val_q2);
    
    if (val_res != EQUAL)
    {
      if (ind1 == MIN_END)
        // Swap the comparison result.
        val_res = (val_res == LARGER ? SMALLER : LARGER);

      return (val_res);
    }

    // Both arcs are defined to the same side (left or right) of the vertical
    // asymptote. If one is defined at y = -oo and the other at y = +oo, we
    // preform a "lexicographic" comparison.
    const Boundary_type  inf_y1 = 
      (ind1 == MIN_END ? left_boundary_in_y() : right_boundary_in_y());
    const Boundary_type  inf_y2 = 
      (ind2 == MIN_END ? arc.left_boundary_in_y() : arc.right_boundary_in_y());

    if (inf_y1 == MINUS_INFINITY && inf_y2 == PLUS_INFINITY)
      return (SMALLER);
    else if (inf_y1 == PLUS_INFINITY && inf_y2 == MINUS_INFINITY)
      return (LARGER);

    // If we reached here, the denominators Q1 and Q2 are equal. We therefore
    // evaluate P1(x') and P2(x') and find the relative x-positions of the
    // curve ends by considering their comparison result.
    const Algebraic    val_p1 = nt_traits.evaluate_at (_numer, x1);
    const Algebraic    val_p2 = nt_traits.evaluate_at (arc._numer, x1);
    
    val_res = CGAL::compare (val_p1, val_p2);

    if (val_res == EQUAL)
    {
      CGAL_assertion (_has_same_base (arc));
      return (EQUAL);
    }

    if (ind1 == MAX_END)
      // Swap the comparison result.
      val_res = (val_res == LARGER ? SMALLER : LARGER);

    return (val_res);
  }

  /*!
   * Compare the slopes of the arc with another given arc at their given 
   * intersection point.
   * \param cv The given arc.
   * \param p The intersection point.
   * \param mult Output: The mutiplicity of the intersection point.
   * \return SMALLER if (*this) slope is less than cv's;
   *         EQUAL if the two slopes are equal;
   *         LARGER if (*this) slope is greater than cv's.
   */
  Comparison_result compare_slopes (const Self& arc,
				    const Point_2& p,
				    unsigned int& mult) const
  {
    // Make sure that p is in the x-range of both arcs.
    CGAL_precondition (is_valid());
    CGAL_precondition (arc.is_valid());
    CGAL_precondition (_is_in_true_x_range (p.x()) &&
		       arc._is_in_true_x_range (p.x()));

    // Check the case of overlapping arcs.
    if (_has_same_base (arc))
      return (EQUAL);

    // The intersection point may have a maximal multiplicity value of:
    //   max (deg(p1) + deg(q2), deg(q1) + deg(p2)).
    const Algebraic&  _x = p.x();
    Nt_traits         nt_traits;
    Polynomial        pnum1 = this->_numer;
    Polynomial        pden1 = this->_denom;
    const bool        simple_poly1 = (nt_traits.degree (pden1) <= 0);
    Polynomial        pnum2 = arc._numer;
    Polynomial        pden2 = arc._denom;
    const bool        simple_poly2 = (nt_traits.degree (pden2) <= 0);
    int               max_mult;
    Algebraic         d1, d2;
    Comparison_result res;

    max_mult = nt_traits.degree (pnum1) + nt_traits.degree (pden2);

    if (nt_traits.degree (pden1) + nt_traits.degree (pnum2) > max_mult)
      max_mult = nt_traits.degree (pden1) + nt_traits.degree (pnum2);

    for (mult = 1; mult <= static_cast<unsigned int>(max_mult); mult++)
    {
      // Compute the current derivative. Use the equation:
      //
      // (p(x) / q(x))' = (p'(x)*q(x) - p(x)*q'(x)) / q^2(x)
      if (simple_poly1)
      {
	pnum1 = nt_traits.derive(pnum1);
      }
      else
      {
	pnum1 = nt_traits.derive(pnum1)*pden1 - pnum1*nt_traits.derive(pden1);
	pden1 *= pden1;
      }

      if (simple_poly2)
      {
	pnum2 = nt_traits.derive(pnum2);
      }
      else
      {
	pnum2 = nt_traits.derive(pnum2)*pden2 - pnum2*nt_traits.derive(pden2);
	pden2 *= pden2;
      }

      // Compute the two derivative values and compare them. 
      d1 = nt_traits.evaluate_at (pnum1, _x) / 
	   nt_traits.evaluate_at (pden1, _x);
      d2 = nt_traits.evaluate_at (pnum2, _x) / 
	   nt_traits.evaluate_at (pden2, _x);

      res = CGAL::compare (d1, d2);

      // Stop here in case the derivatives are not equal.
      if (res != EQUAL)
	return (res);
    }

    // If we reached here, there is an overlap.
    return (EQUAL);
  }

  /*!
   * Compare the two arcs at x = -oo.
   * \param arc The given arc.
   * \pre Both arcs have a left end which is unbounded in x.
   * \return SMALLER if (*this) lies below the other arc;
   *         EQUAL if the two supporting functions are equal;
   *         LARGER if (*this) lies above the other arc.
   */
  Comparison_result compare_at_minus_infinity (const Self& arc) const
  {
    CGAL_precondition (left_boundary_in_x() == MINUS_INFINITY &&
                       arc.left_boundary_in_x() == MINUS_INFINITY);

    // Check for easy cases, when the infinity at y of both ends is different.
    const Boundary_type  inf_y1 = left_boundary_in_y();
    const Boundary_type  inf_y2 = arc.left_boundary_in_y();

    if (inf_y1 != inf_y2)
    {
      if (inf_y1 == MINUS_INFINITY || inf_y2 == PLUS_INFINITY)
        return (SMALLER);
      else if (inf_y1 == PLUS_INFINITY || inf_y2 == MINUS_INFINITY)
        return (LARGER);
      else
        CGAL_assertion (false);
    }

    // First compare the signs of the two denominator polynomials at x = -oo.
    const CGAL::Sign  sign1 = _sign_at_minus_infinity (_denom);
    const CGAL::Sign  sign2 = _sign_at_minus_infinity (arc._denom);
   
    CGAL_assertion (sign1 != CGAL::ZERO && sign2 != CGAL::ZERO);

    const bool        flip_res = (sign1 != sign2);

    // We wish to compare the two following functions:
    //
    //   y = p1(x)/q1(x)    and     y = p2(x)/q2(x)
    //
    // It is clear that we should look at the sign of the polynomial
    // p1(x)*q2(x) - p2(x)*q1(x) at x = -oo.
    const CGAL::Sign  sign_ip = _sign_at_minus_infinity (_numer*arc._denom - 
                                                         arc._numer*_denom);

    if (sign_ip == CGAL::ZERO)
    {
      CGAL_assertion (_has_same_base (arc));
      return (EQUAL);
    }

    if (sign_ip == CGAL::NEGATIVE)
      return (flip_res ? LARGER : SMALLER);
    else
      return (flip_res ? SMALLER : LARGER);
  }

  /*!
   * Compare the two arcs at x = +oo.
   * \param arc The given arc.
   * \pre Both arcs are have a right end which is unbounded in x.
   * \return SMALLER if (*this) lies below the other arc;
   *         EQUAL if the two supporting functions are equal;
   *         LARGER if (*this) lies above the other arc.
   */
  Comparison_result compare_at_plus_infinity (const Self& arc) const
  {
    CGAL_precondition (right_boundary_in_x() == PLUS_INFINITY &&
                       arc.right_boundary_in_x() == PLUS_INFINITY);

    // Check for easy cases, when the infinity at y of both ends is different.
    const Boundary_type  inf_y1 = right_boundary_in_y();
    const Boundary_type  inf_y2 = arc.right_boundary_in_y();

    if (inf_y1 != inf_y2)
    {
      if (inf_y1 == MINUS_INFINITY || inf_y2 == PLUS_INFINITY)
        return (SMALLER);
      else if (inf_y1 == PLUS_INFINITY || inf_y2 == MINUS_INFINITY)
        return (LARGER);
      else
        CGAL_assertion (false);
    }

    // First compare the signs of the two denominator polynomials at x = +oo.
    const CGAL::Sign  sign1 = _sign_at_plus_infinity (_denom);
    const CGAL::Sign  sign2 = _sign_at_plus_infinity (arc._denom);
   
    CGAL_assertion (sign1 != CGAL::ZERO && sign2 != CGAL::ZERO);

    const bool        flip_res = (sign1 != sign2);

    // We wish to compare the two following functions:
    //
    //   y = p1(x)/q1(x)    and     y = p2(x)/q2(x)
    //
    // It is clear that we should look at the sign of the polynomial
    // p1(x)*q2(x) - p2(x)*q1(x) at x = +oo.
    const CGAL::Sign  sign_ip = _sign_at_plus_infinity (_numer*arc._denom - 
                                                        arc._numer*_denom);

    if (sign_ip == CGAL::ZERO)
    {
      CGAL_assertion (_has_same_base (arc));
      return (EQUAL);
    }

    if (sign_ip == CGAL::NEGATIVE)
      return (flip_res ? LARGER : SMALLER);
    else
      return (flip_res ? SMALLER : LARGER);
  }

  /*!
   * Check whether the two arcs are equal (have the same graph).
   * \param arc The compared arc.
   * \return (true) if the two arcs have the same graph; (false) otherwise.
   */
  bool equals (const Self& arc) const
  {
    // The two arc must have the same base rational function.
    CGAL_precondition (is_valid());
    CGAL_precondition (arc.is_valid());
    if (! _has_same_base (arc))
      return (false);

    // Check that the arc endpoints are the same.
    Boundary_type inf1 = left_boundary_in_x ();
    Boundary_type inf2 = arc.left_boundary_in_x ();

    if (inf1 != inf2)
      return (false);

    if (inf1 == NO_BOUNDARY && CGAL::compare(left().x(), arc.left().x()) != EQUAL)
      return (false);

    inf1 = right_boundary_in_x ();
    inf2 = arc.right_boundary_in_x ();

    if (inf1 != inf2)
      return (false);

    if (inf1 == NO_BOUNDARY && CGAL::compare(right().x(), arc.right().x()) != EQUAL)
      return (false);

    // If we reached here, the two arc are equal:
    return (true);
  }

  /*!
   * Check whether it is possible to merge the arc with the given arc.
   * \param arc The query arc.
   * \return (true) if it is possible to merge the two arcs;
   *         (false) otherwise.
   */
  bool can_merge_with (const Self& arc) const
  {
    // In order to merge the two arcs, they should have the same base rational
    // function.
    CGAL_precondition (is_valid());
    CGAL_precondition (arc.is_valid());
    if (! _has_same_base (arc))
      return (false);

    // Check if the left endpoint of one curve is the right endpoint of the
    // other.
    Alg_kernel   ker;

    return ((right_boundary_in_x() == NO_BOUNDARY &&
             right_boundary_in_y() == NO_BOUNDARY &&