conic_x_monotone_arc_2.h

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

  }

  /*!
   * Check whether the given point lies on the arc.
   * \param p The qury point.
   * \param (true) if p lies on the arc; (false) otherwise.
   */
  bool contains_point (const Conic_point_2& p) const
  {
    // First check if p lies on the supporting conic. We first check whether
    // it is one of p's generating conic curves.
    bool       p_on_conic = false;

    if (p.is_generating_conic (_id))
    {
      p_on_conic = true;
    }
    else
    {
      // Check whether p satisfies the supporting conic equation.
      p_on_conic = _is_on_supporting_conic (p.x(), p.y());

      if (p_on_conic)
      {
        // As p lies on the supporting conic of our arc, add its ID to
        // the list of generating conics for p.
        Conic_point_2&  p_non_const = const_cast<Conic_point_2&> (p);
        p_non_const.set_generating_conic (_id);
      }
    }

    if (! p_on_conic)
      return (false);

    // Check if p is between the endpoints of the arc.
    return (_is_between_endpoints (p));
  }
  //@}

  /// \name Constructing points on the arc.
  //@{

  /*!
   * Compute a point on the arc with the same x-coordiante as the given point.
   * \param p The given point.
   * \pre The arc is not vertical and p is in the x-range of the arc.
   * \return A point on the arc with the same x-coordiante as p.
   */
  Point_2 get_point_at_x (const Point_2& p) const
  {
    // Make sure that p is in the x-range of the arc.
    CGAL_precondition ((this->_info & IS_VERTICAL_SEGMENT) == 0);

    CGAL_precondition_code (
      Alg_kernel   ker;
    );

    CGAL_precondition (ker.compare_x_2_object() (p, left()) != SMALLER &&
                       ker.compare_x_2_object() (p, right()) != LARGER);

    if (_is_special_segment())
    {
      // In case of a special segment, the equation of the supported line
      // (a*x + b*y + c) = 0 is stored with the extra data field, and we
      // simply have:
      Algebraic        _y = -(this->_extra_data_P->a*p.x() + 
                              this->_extra_data_P->c) /
                            this->_extra_data_P->b;

      // Return the computed point.
      return (Point_2 (p.x(), _y));
    }

    // Compute the y-coordinate according to the degree of the supporting
    // conic curve.
    Nt_traits        nt_traits;
    Algebraic        y;

    if ((this->_info & DEGREE_MASK) == DEGREE_1)
    {
      // In case of a linear curve, the y-coordinate is a simple linear
      // expression of x(p) (note that v is not 0 as the arc is not vertical):
      //   y = -(u*x(p) + w) / v
      y = -(alg_u*p.x() + alg_w) / alg_v;
    }
    else if (this->_orient == COLLINEAR)
    {
      CGAL_assertion (this->_extra_data_P != NULL);

      // In this case the equation of the supporting line is given by the
      // extra data structure.
      y = -(this->_extra_data_P->a * p.x() +
            this->_extra_data_P->c) / this->_extra_data_P->b;
    }
    else
    {
      CGAL_assertion ((this->_info & DEGREE_MASK) == DEGREE_2);

      // In this case the y-coordinate is one of solutions to the quadratic
      // equation:
      //  s*y^2 + (t*x(p) + v)*y + (r*x(p)^2 + u*x(p) + w) = 0
      Algebraic  A = alg_s;
      Algebraic  B = alg_t*p.x() + alg_v;
      Algebraic  C = (alg_r*p.x() + alg_u)*p.x() + alg_w;

      if (CGAL::sign(this->_s) == ZERO)
      {
        // In this case A is 0 and we have a linear equation.
        CGAL_assertion (CGAL::sign (B) != ZERO);

        y = -C / B;
      }
      else
      {
        // Solve the quadratic equation.
        Algebraic  disc = B*B - 4*A*C;

        CGAL_assertion (CGAL::sign (disc) != NEGATIVE);

        // We take either the root involving -sqrt(disc) or +sqrt(disc)
        // based on the information flags.
        if ((this->_info & PLUS_SQRT_DISC_ROOT) != 0)
        {
          y = (nt_traits.sqrt (disc) - B) / (2*A);
        }
        else

        {
          y = -(B + nt_traits.sqrt (disc)) / (2*A);
        }
      }
    }

    // Return the computed point.
    return (Point_2 (p.x(), y));
  }

  /*!
   * Get a polyline approximating the conic arc.
   * \param n The maximal number of sample points.
   * \param oi An output iterator, whose value-type is pair<double,double>
   *           (representing an approximated point).
   *           In case the arc is a line segment, there are 2 output points,
   *           otherwise the arc is approximated by the polyline defined by
   *           (p_0, p_1, ..., p_n), where p_0 and p_n are the left and right
   *           endpoints of the arc, respectively.
   */
  template <class OutputIterator>
  OutputIterator polyline_approximation (size_t n,
                                         OutputIterator oi) const
  {
    CGAL_precondition (n != 0);

    const double  x_left = CGAL::to_double (left().x());
    const double  y_left = CGAL::to_double (left().y());
    const double  x_right = CGAL::to_double (right().x());
    const double  y_right = CGAL::to_double (right().y());

    if (this->_orient == COLLINEAR)
    {
      // In case of a line segment, return the two endpoints.
      *oi = std::pair<double, double> (x_left, y_left);
      ++oi;
      *oi = std::pair<double, double> (x_right, y_right);
      ++oi;
      return (oi);
    }
    
    // Otherwise, sample (n - 1) equally-spaced points in between.
    const double  app_r = CGAL::to_double (this->_r);
    const double  app_s = CGAL::to_double (this->_s);
    const double  app_t = CGAL::to_double (this->_t);
    const double  app_u = CGAL::to_double (this->_u);
    const double  app_v = CGAL::to_double (this->_v);
    const double  app_w = CGAL::to_double (this->_w);
    const double  x_jump = (x_right - x_left) / n;
    double        x, y;
    const bool    A_is_zero = (CGAL::sign(this->_s) == ZERO);
    double        A = app_s, B, C;
    double        disc;
    size_t        i;

    *oi = std::pair<double, double> (x_left, y_left);   // The left point.
    ++oi;
    for (i = 1; i < n; i++)
    {
      x = x_left + x_jump*i;

      // Solve the quadratic equation: A*x^2 + B*x + C = 0:
      B = app_t*x + app_v;
      C = (app_r*x + app_u)*x + app_w;

      if (A_is_zero)
      {
        y = -C / B;
      }
      else
      {
        disc = B*B - 4*A*C;

        if (disc < 0)
          disc = 0;

        // We take either the root involving -sqrt(disc) or +sqrt(disc)
        // based on the information flags.
        if ((this->_info & PLUS_SQRT_DISC_ROOT) != 0)
        {
          y = (std::sqrt(disc) - B) / (2*A);
        }
        else
        {
          y = -(B + std::sqrt (disc)) / (2*A);
        }
      }

      *oi = std::pair<double, double> (x, y);
      ++oi;
    }
    *oi = std::pair<double, double> (x_right, y_right);   // The right point.
    ++oi;

    return (oi);
  }

  /*!
   * Compare to arcs immediately to the right of their intersection point.
   * \param arc The compared arc.
   * \param p The reference intersection point.
   * \return The relative position of the arcs to the right of p.
   * \pre Both arcs we compare are not vertical segments.
   */
  Comparison_result compare_to_right (const Self& arc,
                                      const Conic_point_2& p) const
  {
    CGAL_precondition ((this->_info & IS_VERTICAL_SEGMENT) == 0 &&
                       (arc._info & IS_VERTICAL_SEGMENT) == 0);

    // In case one arc is facing upwards and another facing downwards, it is
    // clear that the one facing upward is above the one facing downwards.
    if (_has_same_supporting_conic (arc))
    {
      if ((this->_info & FACING_UP) != 0 && (arc._info & FACING_DOWN) != 0)
        return (LARGER);
      else if ((this->_info & FACING_DOWN)!= 0 && (arc._info & FACING_UP) != 0)
        return (SMALLER);

      // In this case the two arcs overlap.
      CGAL_assertion ((this->_info & FACING_MASK) == 
                      (arc._info & FACING_MASK));

      return (EQUAL);
    }

    // Compare the slopes of the two arcs at p, using their first-order
    // partial derivatives.
    Algebraic      slope1_numer, slope1_denom;
    Algebraic      slope2_numer, slope2_denom;

    _derive_by_x_at (p, 1, slope1_numer, slope1_denom);
    arc._derive_by_x_at (p, 1, slope2_numer, slope2_denom);

    // Check if any of the slopes is vertical.
    const bool     is_vertical_slope1 = (CGAL::sign (slope1_denom) == ZERO);
    const bool     is_vertical_slope2 = (CGAL::sign (slope2_denom) == ZERO);

    if (!is_vertical_slope1 && !is_vertical_slope2)
    {
      // The two derivatives at p are well-defined: use them to determine
      // which arc is above the other (the one with a larger slope is below).
      Comparison_result slope_res = CGAL::compare (slope1_numer*slope2_denom,
                                                   slope2_numer*slope1_denom);

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

      // Use the second-order derivative.
      _derive_by_x_at (p, 2, slope1_numer, slope1_denom);
      arc._derive_by_x_at (p, 2, slope2_numer, slope2_denom);

      slope_res = CGAL::compare (slope1_numer*slope2_denom,
                                 slope2_numer*slope1_denom);

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

      // Use the third-order derivative.
      _derive_by_x_at (p, 3, slope1_numer, slope1_denom);
      arc._derive_by_x_at (p, 3, slope2_numer, slope2_denom);
      
      slope_res = CGAL::compare (slope1_numer*slope2_denom,
                                 slope2_numer*slope1_denom);

      // \todo Handle higher-order derivatives:
      CGAL_assertion (slope_res != EQUAL);

      return (slope_res);
    }
    else if (!is_vertical_slope2)
    {
      // The first arc has a vertical slope at p: check whether it is
      // facing upwards or downwards and decide accordingly.
      CGAL_assertion ((this->_info & FACING_MASK) != 0);

      if ((this->_info & FACING_UP) != 0)
        return (LARGER);
      return (SMALLER);
    }
    else if (!is_vertical_slope1)
    {
      // The second arc has a vertical slope at p_int: check whether it is
      // facing upwards or downwards and decide accordingly.
      CGAL_assertion ((arc._info & FACING_MASK) != 0);

      if ((arc._info & FACING_UP) != 0)
        return (SMALLER);
      return (LARGER);
    }

    // The two arcs have vertical slopes at p_int:
    // First check whether one is facing up and one down. In this case the
    // comparison result is trivial.
    if ((this->_info & FACING_UP) != 0 && (arc._info & FACING_DOWN) != 0)
      return (LARGER);
    else if ((this->_info & FACING_DOWN)!= 0 && (arc._info & FACING_UP)!= 0)
      return (SMALLER);

    // Compute the second-order derivative by y and act according to it.
    _derive_by_y_at (p, 2, slope1_numer, slope1_denom);
    arc._derive_by_y_at (p, 2, slope2_numer, slope2_denom);

    Comparison_result slope_res = CGAL::compare (slope1_numer*slope2_denom,
                                                 slope2_numer*slope1_denom);

    // If necessary, use the third-order derivative by y.
    if (slope_res == EQUAL)
    {
      // \todo Check this!
      _derive_by_y_at (p, 3, slope1_numer, slope1_denom);
      arc._derive_by_y_at (p, 3, slope2_numer, slope2_denom);
      
      slope_res = CGAL::compare (slope2_numer*slope1_denom,
                                 slope1_numer*slope2_denom);
    }

    // \todo Handle higher-order derivatives:
    CGAL_assertion(slope_res != EQUAL);

    if ((this->_info & FACING_UP) != 0 && (arc._info & FACING_UP) != 0)
    {
      // Both are facing up.
      return ((slope_res == LARGER) ? SMALLER : LARGER);
    }
    // Both are facing down.
    return (slope_res);
  }

  /*!
   * Compare to arcs immediately to the leftt of their intersection point.
   * \param arc The compared arc.
   * \param p The reference intersection point.
   * \return The relative position of the arcs to the left of p.
   * \pre Both arcs we compare are not vertical segments.
   */
  Comparison_result compare_to_left (const Self& arc,
                                     const Conic_point_2& p) const
  {
    CGAL_precondition ((this->_info & IS_VERTICAL_SEGMENT) == 0 &&
                       (arc._info & IS_VERTICAL_SEGMENT) == 0);

    // In case one arc is facing upwards and another facing downwards, it is
    // clear that the one facing upward is above the one facing downwards.
    if (_has_same_supporting_conic (arc))
    {
      if ((this->_info & FACING_UP) != 0 && (arc._info & FACING_DOWN) != 0)
        return (LARGER);
      else if ((this->_info & FACING_DOWN)!= 0 && (arc._info & FACING_UP)!= 0)
        return (SMALLER);

      // In this case the two arcs overlap.
      CGAL_assertion ((this->_info & FACING_MASK) == 
                      (arc._info & FACING_MASK));

      return (EQUAL);
    }

    // Compare the slopes of the two arcs at p, using their first-order
    // partial derivatives.
    Algebraic      slope1_numer, slope1_denom;
    Algebraic      slope2_numer, slope2_denom;

    _derive_by_x_at (p, 1, slope1_numer, slope1_denom);
    arc._derive_by_x_at (p, 1, slope2_numer, slope2_denom);

    // Check if any of the slopes is vertical.
    const bool     is_vertical_slope1 = (CGAL::sign (slope1_denom) == ZERO);

    const bool     is_vertical_slope2 = (CGAL::sign (slope2_denom) == ZERO);

    if (!is_vertical_slope1 && !is_vertical_slope2)
    {
      // The two derivatives at p are well-defined: use them to determine
      // which arc is above the other (the one with a larger slope is below).
      Comparison_result  slope_res = CGAL::compare(slope2_numer*slope1_denom,
                                                   slope1_numer*slope2_denom);

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

      // Use the second-order derivative.
      _derive_by_x_at (p, 2, slope1_numer, slope1_denom);
      arc._derive_by_x_at (p, 2, slope2_numer, slope2_denom);

      slope_res = CGAL::compare (slope1_numer*slope2_denom,
                                 slope2_numer*slope1_denom);

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

      // Use the third-order derivative.
      _derive_by_x_at (p, 3, slope1_numer, slope1_denom);
      arc._derive_by_x_at (p, 3, slope2_numer, slope2_denom);
      
      slope_res = CGAL::compare (slope2_numer*slope1_denom,
                                 slope1_numer*slope2_denom);

      // \todo Handle higher-order derivatives:
      CGAL_assertion (slope_res != EQUAL);