rational_arc_2.h

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

      const CGAL::Sign                   sign_s = 
        (dir_right ? signs.second : signs.first);

      if (sign_s == CGAL::NEGATIVE)
        _info = (_info | SRC_AT_Y_MINUS_INFTY);
      else
        _info = (_info | SRC_AT_Y_PLUS_INFTY);
    }

    // Set the properties of the target.
    Algebraic           y0;

    if (dir_right)
    {
      // The target point is at x = +oo.
      const Boundary_type inf_t = _analyze_at_plus_infinity (_numer, _denom,
                                                             y0);

      if (inf_t == MINUS_INFINITY)
        _info = (_info | TRG_AT_Y_MINUS_INFTY);
      else if (inf_t == PLUS_INFINITY)
        _info = (_info | TRG_AT_Y_PLUS_INFTY);
      else // if (inf_t == NO_BOUNDARY)
        _pt = Point_2 (0, y0);
    }
    else
    {
      // The target point is at x = -oo.
      const Boundary_type inf_t = _analyze_at_minus_infinity (_numer, _denom,
                                                              y0);

      if (inf_t == MINUS_INFINITY)
        _info = (_info | TRG_AT_Y_MINUS_INFTY);
      else if (inf_t == PLUS_INFINITY)
        _info = (_info | TRG_AT_Y_PLUS_INFTY);
      else // if (inf_t == NO_BOUNDARY)
        _pt = Point_2 (0, y0);
    }

    // Mark that the arc is valid. As it may have poles, we do not mark it
    // as continuous.
    _info = (_info | IS_VALID);
  }
  
  /*!
   * Constructor of a bounded rational arc, defined by y = p(x)/q(x), 
   * where: x_min <= x <= x_max.
   * \param pcoeffs The rational coefficients of the polynomial p(x).
   * \param qcoeffs The rational coefficients of the polynomial q(x).
   * \param x_s The x-coordinate of the source point.
   * \param x_t The x-coordinate of the target point.
   * \pre The two x-coordinate must not be equal,
   *      and q(x) != 0 for all x_min <= x <= x_max.
   */
  _Rational_arc_2 (const Rat_vector& pcoeffs, const Rat_vector& qcoeffs,
		   const Algebraic& x_s, const Algebraic& x_t) :
    _info (0)
  {
    // Compare the x-coordinates and determine the direction.
    Comparison_result   x_res = CGAL::compare (x_s, x_t);

    CGAL_precondition (x_res != EQUAL);
    if (x_res == EQUAL)
      return;

    if (x_res == SMALLER)
      _info = (_info | IS_DIRECTED_RIGHT);

    // Set the numerator and denominator polynomials.
    Nt_traits    nt_traits;
    const bool   valid = 
      nt_traits.construct_polynomials (&(pcoeffs[0]),
					pcoeffs.size() - 1,
                                       &(qcoeffs[0]),
					qcoeffs.size() - 1,
                                       _numer, _denom);

    CGAL_precondition_msg (valid,
                           "The denominator polynomial must not be 0.");
    if (! valid)
      return;

    // Set the source point and check if it lies next to a pole.
    if (CGAL::sign (nt_traits.evaluate_at (_denom, x_s)) != CGAL::ZERO)
    {
      // We have a nomral endpoint.
      _ps = Point_2 (x_s, nt_traits.evaluate_at (_numer, x_s) /
                          nt_traits.evaluate_at (_denom, x_s));
    }
    else
    {
      // The y-coodinate is unbounded, but we can set its sign.
      _ps = Point_2 (x_s, 0);

      std::pair<CGAL::Sign, CGAL::Sign>  signs = _analyze_near_pole (x_s);
      const CGAL::Sign                   sign_s =
        ((_info & IS_DIRECTED_RIGHT) != 0) ? signs.second : signs.first;

      if (sign_s == CGAL::NEGATIVE)
        _info = (_info | SRC_AT_Y_MINUS_INFTY);
      else
        _info = (_info | SRC_AT_Y_PLUS_INFTY);
    }

    // Set the target point and check if it lies next to a pole.
    if (CGAL::sign (nt_traits.evaluate_at (_denom, x_t)) != CGAL::ZERO)
    {
      // We have a nomral endpoint.
      _pt = Point_2 (x_t, nt_traits.evaluate_at (_numer, x_t) /
                     nt_traits.evaluate_at (_denom, x_t));
    }
    else
    {
      // The y-coodinate is unbounded, but we can set its sign.
      _pt = Point_2 (x_t, 0);
      
      std::pair<CGAL::Sign, CGAL::Sign>  signs = _analyze_near_pole (x_t);
      const CGAL::Sign                   sign_t =
        ((_info & IS_DIRECTED_RIGHT) != 0) ? signs.first : signs.second;

      if (sign_t == CGAL::NEGATIVE)
        _info = (_info | TRG_AT_Y_MINUS_INFTY);
      else
        _info = (_info | TRG_AT_Y_PLUS_INFTY);
    }

    // Mark that the arc is valid. As it may have poles, we do not mark it
    // as continuous.
    _info = (_info | IS_VALID);
  }

  /*!
   * Break an arc of a rational function into continuous sub-arcs, splitting
   * it at its poles.
   */
  template <class OutputIterator>
  OutputIterator make_continuous (OutputIterator oi) const
  {
    // Compute the roots of the denominator polynomial.
    std::list<Algebraic>                 q_roots;
    bool                                 root_at_ps, root_at_pt;

    if ((_info & IS_CONTINUOUS) == 0)
      _denominator_roots (std::back_inserter (q_roots),
                          root_at_ps, root_at_pt);

    // Check the case of a continuous arc:
    if (q_roots.empty())
    {
      Self    arc = *this;

      arc._info = (arc._info | IS_CONTINUOUS);
      *oi = arc;
      ++oi;
      return (oi);
    }

    // Split the arc accordingly.
    typename std::list<Algebraic>::const_iterator iter;
    Self                                          arc = *this;

    for (iter = q_roots.begin(); iter != q_roots.end(); ++iter)
    {
      *oi = arc._split_at_pole (*iter);
      ++oi;
    }

    // Add the final x-monotone sub-arc.
    arc._info = (arc._info | IS_CONTINUOUS);
    *oi = arc;
    ++oi;

    return (oi);
  }
  //@}

  /// \name Accessing the arc properties.
  //@{

  /*! Get the numerator polynomial of the underlying rational function. */
  const Polynomial& numerator () const
  {
    return (_numer);
  }

  /*! Get the denominator polynomial of the underlying rational function. */
  const Polynomial& denominator () const
  {
    return (_denom);
  }

  /*! Check if the x-coordinate of the source point is infinite. */
  Boundary_type source_boundary_in_x () const
  {
    if ((_info & SRC_AT_X_MINUS_INFTY) != 0)
      return (MINUS_INFINITY);
    else if ((_info & SRC_AT_X_PLUS_INFTY) != 0)
      return (PLUS_INFINITY);
    else
      return (NO_BOUNDARY);
  }

  /*! Check if the y-coordinate of the source point is infinite. */
  Boundary_type source_boundary_in_y () const
  {
    if ((_info & SRC_AT_Y_MINUS_INFTY) != 0)
      return (MINUS_INFINITY);
    else if ((_info & SRC_AT_Y_PLUS_INFTY) != 0)
      return (PLUS_INFINITY);
    else
      return (NO_BOUNDARY);
  }

  /*! Check if the x-coordinate of the target point is infinite. */
  Boundary_type target_boundary_in_x () const
  {
    if ((_info & TRG_AT_X_MINUS_INFTY) != 0)
      return (MINUS_INFINITY);
    else if ((_info & TRG_AT_X_PLUS_INFTY) != 0)
      return (PLUS_INFINITY);
    else
      return (NO_BOUNDARY);
  }

  /*! Check if the y-coordinate of the target point is infinite. */
  Boundary_type target_boundary_in_y () const
  {
    if ((_info & TRG_AT_Y_MINUS_INFTY) != 0)
      return (MINUS_INFINITY);
    else if ((_info & TRG_AT_Y_PLUS_INFTY) != 0)
      return (PLUS_INFINITY);
    else
      return (NO_BOUNDARY);
  }

  /*! Get the source point. */
  const Point_2& source () const
  {
    CGAL_precondition ((_info & IS_VALID) != 0 &&
                       source_boundary_in_x() == NO_BOUNDARY &&
                       source_boundary_in_y() == NO_BOUNDARY);
    return (_ps);
  }

  /*! Get the x-coordinate of the source point. */
  Algebraic source_x () const
  {
    CGAL_precondition ((_info & IS_VALID) != 0 &&
                       source_boundary_in_x() == NO_BOUNDARY);
    return (_ps.x());
  }

  /*! Get the y-coordinate of the source point. */
  Algebraic source_y () const
  {
    CGAL_precondition ((_info & IS_VALID) != 0 &&
                       source_boundary_in_y() == NO_BOUNDARY);
    return (_ps.y());
  }

  /*! Get the target point. */
  const Point_2& target () const
  {
    CGAL_precondition ((_info & IS_VALID) != 0 &&
                       target_boundary_in_x() == NO_BOUNDARY &&
                       target_boundary_in_y() == NO_BOUNDARY);
    return (_pt);
  }

  /*! Get the x-coordinate of the target point. */
  Algebraic target_x () const
  {
    CGAL_precondition ((_info & IS_VALID) != 0 &&
                       target_boundary_in_x() == NO_BOUNDARY);
    return (_pt.x());
  }

  /*! Get the y-coordinate of the target point. */
  Algebraic target_y () const
  {
    CGAL_precondition ((_info & IS_VALID) != 0 &&
                       target_boundary_in_y() == NO_BOUNDARY);
    return (_pt.y());
  }

  /*! Check if the x-coordinate of the left point is infinite. */
  Boundary_type left_boundary_in_x () const
  {
    if ((_info & IS_DIRECTED_RIGHT) != 0)
      return (source_boundary_in_x());
    else
      return (target_boundary_in_x());
  }

  /*! Check if the y-coordinate of the left point is infinite. */
  Boundary_type left_boundary_in_y () const
  {
    if ((_info & IS_DIRECTED_RIGHT) != 0)
      return (source_boundary_in_y());
    else
      return (target_boundary_in_y());
  }

  /*! Check if the x-coordinate of the right point is infinite. */
  Boundary_type right_boundary_in_x () const
  {
    if ((_info & IS_DIRECTED_RIGHT) != 0)
      return (target_boundary_in_x());
    else
      return (source_boundary_in_x());
  }

  /*! Check if the y-coordinate of the right point is infinite. */
  Boundary_type right_boundary_in_y () const
  {
    if ((_info & IS_DIRECTED_RIGHT) != 0)
      return (target_boundary_in_y());
    else
      return (source_boundary_in_y());
  }

  /*! Get the left endpoint. */
  const Point_2& left () const
  {
    CGAL_precondition (left_boundary_in_x() == NO_BOUNDARY &&
                       left_boundary_in_y() == NO_BOUNDARY);
    return ((_info & IS_DIRECTED_RIGHT) ? _ps : _pt);
  }

  /*! Get the right endpoint. */
  const Point_2& right () const
  {
    CGAL_precondition (right_boundary_in_x() == NO_BOUNDARY &&
                       right_boundary_in_y() == NO_BOUNDARY);
    return ((_info & IS_DIRECTED_RIGHT) ? _pt : _ps);
  }

  /*! Check if the arc is valid. */
  bool is_valid () const
  {
    return ((_info & IS_VALID) != 0);
  }

  /*! Check if the arc is continuous. */
  bool is_continuous () const
  {
    return ((_info & IS_CONTINUOUS) != 0);
  }

  /*! Check if the arc is directed right. */
  bool is_directed_right () const
  {
    return ((_info & IS_DIRECTED_RIGHT) != 0);
  }
  //@}

  /// \name Predicates.
  //@{

  /*!
   * Get the relative position of the point with respect to the rational arc.
   * \param p The query point.
   * \pre p is in the x-range of the arc.
   * \return SMALLER if the point is below the arc;
   *         LARGER if the point is above the arc;
   *         EQUAL if p lies on the arc.
   */
  Comparison_result point_position (const Point_2& p) const
  {
    // Make sure that p is in the x-range of the arc and check whether it
    // has the same x-coordinate as one of the endpoints.
    CGAL_precondition (is_continuous());
    CGAL_precondition (_is_in_true_x_range (p.x()));

    // Evaluate the rational function at x(p), which lies at the interior
    // of the x-range.
    Nt_traits   nt_traits;
    Algebraic   y = nt_traits.evaluate_at (_numer, p.x()) /
                    nt_traits.evaluate_at (_denom, p.x());
    
    // Compare the resulting y-coordinate with y(p):
    return (CGAL::compare (p.y(), y));
  }

  /*!
   * Compare the x-coordinate of a vertical asymptote of the arc (one of its
   * ends) and the given point.
   */
  Comparison_result compare_end (const Point_2& p,
                                 Curve_end ind) const
  {
    Algebraic                x0;

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

    // Compare the x-coordinates.
    const Comparison_result  res = CGAL::compare (p.x(), x0);

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

    return ((ind == MIN_END) ? SMALLER : LARGER);
  }

  /*!