rational_arc_2.h

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

  /*!
   * Print the rational arc.
   */
  std::ostream& print (std::ostream& os) const
  {    
    // Print y as a rational function of x.
    os << "y = (";
    _print_polynomial (os, _numer, 'x');
    os << ") / (";
    _print_polynomial (os, _denom, 'x');
    os << ") on ";

    // Print the definition range.
    Boundary_type      inf_x = source_boundary_in_x();
    if (inf_x == MINUS_INFINITY) 
      os << "(-oo";
    else if (inf_x == PLUS_INFINITY)
      os << "(+oo";
    else if (source_boundary_in_y() != NO_BOUNDARY)
      os << '(' << source_x();
    else
      os << '[' << source().x();
 
    os << ", ";
    inf_x = target_boundary_in_x();

    if (inf_x == MINUS_INFINITY) 
      os << "-oo)";
    else if (inf_x == PLUS_INFINITY)
      os << "+oo)";
    else if (target_boundary_in_y() != NO_BOUNDARY)
      os << target_x() << ')';
    else
      os << target().x() << ']';

    return (os);
  }

  //@}

private:

  /// \name Auxiliary (private) functions.
  //@{

  /*!
   * Check if the given x-value is in the x-range of the arc.
   * \param x The x-value.
   * \param eq_src Output: Is this value equal to the x-coordinate of the
   *                       source point.
   * \param eq_trg Output: Is this value equal to the x-coordinate of the
   *                       target point.
   */
  bool _is_in_x_range (const Algebraic& x,
                       bool& eq_src, bool& eq_trg) const
  {
    Comparison_result  res1;

    eq_src = eq_trg = false;
    if ((_info & IS_DIRECTED_RIGHT) != 0)
    {
      // Compare to the left endpoint (the source in this case).
      if ((_info & SRC_AT_X_MINUS_INFTY) != 0)
      {
        res1 = LARGER;
      }
      else
      {
        res1 = CGAL::compare (x, _ps.x());

        if (res1 == SMALLER)
          return (false);
        
        if (res1 == EQUAL)
        {
          eq_src = true;
          return (true);
        }
      }
      
      // Compare to the right endpoint (the target in this case).
      if ((_info & TRG_AT_X_PLUS_INFTY) != 0)
        return (true);

      const Comparison_result  res2 = CGAL::compare (x, _pt.x());
      
      if (res2 == LARGER)
          return (false);
       
      if (res2 == EQUAL)
        eq_trg = true;

      return (true);      
    }
    
    // Compare to the left endpoint (the target in this case).
    if ((_info & TRG_AT_X_MINUS_INFTY) != 0)
    {
      res1 = LARGER;
    }
    else
    {
      res1 = CGAL::compare (x, _pt.x());

      if (res1 == SMALLER)
        return (false);
        
      if (res1 == EQUAL)
      {
        eq_trg = true;
        return (true);
      }
    }
      
    // Compare to the right endpoint (the source in this case).
    if ((_info & SRC_AT_X_PLUS_INFTY) != 0)
      return (true);

    const Comparison_result  res2 = CGAL::compare (x, _ps.x());
      
    if (res2 == LARGER)
      return (false);
       
    if (res2 == EQUAL)
      eq_src = true;
    
    return (true);
  }

  /*!
   * Check if the given x-value is in the x-range of the arc, excluding its
   * open ends.
   */
  bool _is_in_true_x_range (const Algebraic& x) const
  {
    bool          eq_src, eq_trg;
    const bool    is_in_x_range_closure = _is_in_x_range (x, eq_src, eq_trg);

    if (! is_in_x_range_closure)
      return (false);

    // Check if we have a vertical asymptote at the source point.
    if (eq_src &&
        (_info & (SRC_AT_Y_MINUS_INFTY | SRC_AT_Y_PLUS_INFTY)) != 0)
      return (false);

    // Check if we have a vertical asymptote at the target point.
    if (eq_trg &&
        (_info & (TRG_AT_Y_MINUS_INFTY | TRG_AT_Y_PLUS_INFTY)) != 0)
      return (false);

    // If we reached here, the value is in the true x-range of the arc.
    return (true);
  }

  /*!
   * Check if the underlying rational fucntion is the same in the given arc.
   * \param arc The given arc.
   * \return (true) if arc's underlying rational fucntion is the same
   *         as of *this; (false) otherwise.
   */
  bool _has_same_base (const Self& arc) const
  {
    // p1(x)/q1(x) == p2(x)/q2(x) if and only if p1*q2 = p2*q1:
    return (_numer * arc._denom == _denom * arc._numer);
  }

  /*!
   * Compute the sign of the given polynomial at x = -oo.
   */
  CGAL::Sign _sign_at_minus_infinity (const Polynomial& poly) const
  {
    // Get the degree.
    Nt_traits    nt_traits;
    const int    degree = nt_traits.degree (poly);

    if (degree < 0)
      return (CGAL::ZERO);

    // Get the leading coefficient. Its sign is the sign of the polynomial
    // at x = -oo if the degree is even, and the opposite sign if it is odd.
    const CGAL::Sign  lead_sign = 
      CGAL::sign (nt_traits.get_coefficient (poly, degree));

    CGAL_assertion (lead_sign != CGAL::ZERO);

    if (degree % 2 == 0)
      return (lead_sign);
    else
      return ((lead_sign == CGAL::POSITIVE) ? CGAL::NEGATIVE : CGAL::POSITIVE);
  }

  /*!
   * Compute the sign of the given polynomial at x = +oo.
   */
  CGAL::Sign _sign_at_plus_infinity (const Polynomial& poly) const
  {
    // Get the degree.
    Nt_traits    nt_traits;
    const int    degree = nt_traits.degree (poly);

    if (degree < 0)
      return (CGAL::ZERO);

    // Get the leading coefficient. Its sign is the sign of the polynomial
    // at x = +oo.
    return (CGAL::sign (nt_traits.get_coefficient (poly, degree)));
  }

  /*!
   * Compute infinity type of the rational function P(x)/Q(x) at x = -oo.
   * \param y Output: The value of the horizontal asymptote (if exists).
   * \return The infinity type for the y-coordinate at x = -oo.
   */
  Boundary_type _analyze_at_minus_infinity (const Polynomial& P,
                                            const Polynomial& Q,
                                            Algebraic& y) const
  {
    // Get the degree of the polynomials.
    Nt_traits    nt_traits;
    const int    deg_p = nt_traits.degree (P);
    const int    deg_q = nt_traits.degree (Q);

    if (deg_p < 0 || deg_p < deg_q)
    {
      // We have a zero polynomial or a zero asymptote.
      y = 0;
      return (NO_BOUNDARY);
    }

    // Get the leading coefficients.
    Integer      p_lead = nt_traits.get_coefficient (P, deg_p);
    Integer      q_lead = nt_traits.get_coefficient (Q, deg_q);
    
    if (deg_p == deg_q)
    {
      // We have a horizontal asymptote.
      y = p_lead / q_lead;
      return (NO_BOUNDARY);
    }

    // We have a tendency to infinity.
    const int    def_diff = deg_p - deg_q;

    if (CGAL::sign (p_lead) == CGAL::sign (q_lead))
      return ((def_diff % 2 == 0) ? PLUS_INFINITY : MINUS_INFINITY);
    else
      return ((def_diff % 2 == 0) ? MINUS_INFINITY : PLUS_INFINITY);
  }

  /*!
   * Compute infinity type of the rational function P(x)/Q(x) at x = +oo.
   * \param y Output: The value of the horizontal asymptote (if exists).
   * \return The infinity type for the y-coordinate at x = +oo.
   */
  Boundary_type _analyze_at_plus_infinity (const Polynomial& P,
                                           const Polynomial& Q,
                                           Algebraic& y) const
  {
    // Get the degree of the polynomials.
    Nt_traits    nt_traits;
    const int    deg_p = nt_traits.degree (P);
    const int    deg_q = nt_traits.degree (Q);

    if (deg_p < 0 || deg_p < deg_q)
    {
      // We have a zero polynomial or a zero asymptote.
      y = 0;
      return (NO_BOUNDARY);
    }

    // Get the leading coefficients.
    Integer      p_lead = nt_traits.get_coefficient (P, deg_p);
    Integer      q_lead = nt_traits.get_coefficient (Q, deg_q);
    
    if (deg_p == deg_q)
    {
      // We have a horizontal asymptote.
      y = p_lead / q_lead;
      return (NO_BOUNDARY);
    }

    // We have a tendency to infinity.
    if (CGAL::sign (p_lead) == CGAL::sign (q_lead))
      return (PLUS_INFINITY);
    else
      return (MINUS_INFINITY);
  }

  /*!
   * Compute all zeros of the denominator polynomial that lie within the
   * x-range of the arc.
   */
  template <class OutputIterator>
  OutputIterator _denominator_roots (OutputIterator oi,
                                     bool& root_at_ps, bool& root_at_pt) const
  {
    Nt_traits         nt_traits;

    root_at_ps = root_at_pt = false;

    if (nt_traits.degree (_denom) <= 0)
      return (oi);

    // Compute the roots of the denominator polynomial.
    std::list<Algebraic>                           q_roots;
    bool                                           eq_src, eq_trg;
    typename std::list<Algebraic>::const_iterator  x_iter;

    nt_traits.compute_polynomial_roots (_denom,
                                        std::back_inserter (q_roots));

    // Go over the roots and check whether they lie in the x-range of the arc.
    for (x_iter = q_roots.begin(); x_iter != q_roots.end(); ++x_iter)
    {
      if (_is_in_x_range (*x_iter, eq_src, eq_trg))
      {
        if (eq_src)
        {
          root_at_ps = true;
        }
        else if (eq_trg)
        {
          root_at_pt = true;
        }
        else
        {
          // The root lies in the interior of the arc.
          *oi = *x_iter;
          ++oi;
        }
      }
    }

    return (oi);
  }

  /*!
   * Determine the signs of the rational functions infinitisimally to the left
   * and to the right of the given pole.
   * \param x0 The x-coordinate of the pole.
   * \pre x0 lies in the interior of the arc.
   * \return The signs to the left and to the right of x0.
   */
  std::pair<CGAL::Sign, CGAL::Sign>
  _analyze_near_pole (const Algebraic& x0) const
  {
    // Note that as the rational function is always normalized, the numerator
    // value is non-zero at the pole x0.
    Nt_traits         nt_traits;
    const Algebraic   numer_at_x0 = nt_traits.evaluate_at (_numer, x0);
    const CGAL::Sign  numer_sign = CGAL::sign (numer_at_x0);

    CGAL_assertion (numer_sign != CGAL::ZERO);

    // Determine the multiplicity of the pole and the sign of the first
    // non-zero derivative of the denominator polynomial at x0.
    int               mult = 1;
    Polynomial        p_der = nt_traits.derive (_denom);
    CGAL::Sign        der_sign;

    while ((der_sign = CGAL::sign (nt_traits.evaluate_at (p_der,
                                                          x0))) == CGAL::ZERO)
    {
      mult++;
      p_der = nt_traits.derive (p_der);
    }

    // Determine the tendency of the rational function to the left and to the
    // right of the pole (to y = -oo or to y = +oo).
    CGAL::Sign        sign_left, sign_right;

    if (mult % 2 == 1)
    {
      // Odd pole multiplicity: different signs from both sides of the pole.
      if (der_sign == numer_sign)
      {
        sign_left = CGAL::NEGATIVE;
        sign_right = CGAL::POSITIVE;
      }
      else
      {
        sign_left = CGAL::POSITIVE;
        sign_right = CGAL::NEGATIVE;
      }
    }
    else
    {
      // Even pole multiplicity: equal signs from both sides of the pole.
      if (der_sign == numer_sign)
      {
        sign_left = CGAL::POSITIVE;
        sign_right = CGAL::POSITIVE;
      }
      else
      {
        sign_left = CGAL::NEGATIVE;
        sign_right = CGAL::NEGATIVE;
      }
    }

    return (std::make_pair (sign_left, sign_right));
  }

  /*!
   * Split the arc into two at a given pole. The function returns the sub-arc
   * to the left of the pole and sets (*this) to be the right sub-arc.
   * \param x0 The x-coordinate of the pole.
   * \pre x0 lies in the interior of the arc.
   * \return The sub-arc to the left of the pole.
   */
  Self _split_at_pole (const Algebraic& x0)
  {
    // Analyze the behaviour of the function near the given pole.
    const std::pair<CGAL::Sign, CGAL::Sign>  signs = _analyze_near_pole (x0);
    const CGAL::Sign    sign_left = signs.first;
    const CGAL::Sign    sign_right = signs.second;

    // Create a fictitious point that represents the x-coordinate of the pole.
    Point_2    p0 (x0, 0);

    // Make a copy of the current ar