conic_x_monotone_arc_2.h

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

    // Clear the _info bits.
    this->_info = Conic_arc_2::IS_VALID;

    // Check if the arc is directed right (the target is lexicographically
    // greater than the source point), or to the left.
    Alg_kernel         ker;
    Comparison_result  dir_res = ker.compare_xy_2_object() (this->_source, 
							    this->_target);

    CGAL_assertion (dir_res != EQUAL);

    if (dir_res == SMALLER)
      this->_info = (this->_info | IS_DIRECTED_RIGHT);

    // Compute the degree of the underlying conic.
    if (CGAL::sign (this->_r) != ZERO ||
        CGAL::sign (this->_s) != ZERO ||
        CGAL::sign (this->_t) != ZERO)
    {
      this->_info = (this->_info | DEGREE_2);
      
      if (this->_orient == COLLINEAR)
      {
        this->_info = (this->_info | IS_SPECIAL_SEGMENT);
        
        if (ker.compare_x_2_object() (this->_source, this->_target) == EQUAL)
        {
          // The arc is a vertical segment:
          this->_info = (this->_info | IS_VERTICAL_SEGMENT);
        }
        
        return;
      }
    }
    else
    {
      CGAL_assertion (CGAL::sign (this->_u) != ZERO ||
                      CGAL::sign (this->_v) != ZERO);

      if (CGAL::sign (this->_v) == ZERO)
      {

        // The supporting curve is of the form: _u*x + _w = 0
        this->_info = (this->_info | IS_VERTICAL_SEGMENT);
      }

      this->_info = (this->_info | DEGREE_1);

      return;
    }

    if (this->_orient == COLLINEAR)
      return;

    // Compute a midpoint between the source and the target and get the y-value
    // of the arc at its x-coordiante.
    Point_2          p_mid = ker.construct_midpoint_2_object() (this->_source,
                                                                this->_target);
    Algebraic        ys[2];
    int              n_ys;

    n_ys = _conic_get_y_coordinates (p_mid.x(), ys);

    CGAL_assertion (n_ys != 0);

    // Check which solution lies on the x-monotone arc.
    Point_2          p_arc_mid (p_mid.x(), ys[0]);

    if (_is_strictly_between_endpoints (p_arc_mid))
    {
      // Mark that we should use the -sqrt(disc) root for points on this
      // x-monotone arc.
      this->_info = (this->_info & ~PLUS_SQRT_DISC_ROOT);
    }
    else
    {
      CGAL_assertion (n_ys == 2);
      p_arc_mid = Point_2 (p_mid.x(), ys[1]);

      CGAL_assertion (_is_strictly_between_endpoints (p_arc_mid));

      // Mark that we should use the +sqrt(disc) root for points on this
      // x-monotone arc.
      this->_info = (this->_info | PLUS_SQRT_DISC_ROOT);
    }

    // Check whether the conic is facing up or facing down:
    // Check whether the arc (which is x-monotone of degree 2) lies above or
    // below the segement that contects its two end-points (x1,y1) and (x2,y2).
    // To do that, we find the y coordinate of a point on the arc whose x
    // coordinate is (x1+x2)/2 and compare it to (y1+y2)/2.
    Comparison_result res = ker.compare_y_2_object() (p_arc_mid, p_mid);

    if (res == LARGER)
    {
      // The arc is above the connecting segment, so it is facing upwards.
      this->_info = (this->_info | FACING_UP);
    }
    else if (res == SMALLER)
    {
      // The arc is below the connecting segment, so it is facing downwards.
      this->_info = (this->_info | FACING_DOWN);
    }

    return;
  }

  /*!
   * Check if the arc is a special segment connecting two algebraic endpoints
   * (and has no undelying integer conic coefficients).
   */
  bool _is_special_segment () const
  {
    return ((this->_info & IS_SPECIAL_SEGMENT) != 0);
  }

  /*!
   * Check whether the given point lies on the supporting conic of the arc.
   * \param px The x-coordinate of query point.
   * \param py The y-coordinate of query point.
   * \return (true) if p lies on the supporting conic; (false) otherwise.
   */
  bool _is_on_supporting_conic (const Algebraic& px,
                                const Algebraic& py) const
  {
    CGAL::Sign       _sign;

    if (! _is_special_segment())
    {
      // Check whether p satisfies the conic equation.
      // The point must satisfy: r*x^2 + s*y^2 + t*xy + u*x + v*y + w = 0.
      _sign = CGAL::sign ((alg_r*px + alg_t*py + alg_u) * px +
                          (alg_s*py + alg_v) * py +
                          alg_w);
    }
    else
    {
      // Check whether p satisfies the equation of the line stored with the
      // extra data.
      _sign = _sign_of_extra_data (px, py);
    }

    return (_sign == ZERO);
  }

  /*!
   * Check whether the two arcs have the same supporting conic.
   * \param arc The compared arc.
   * \return (true) if the two supporting conics are the same.
   */
  bool _has_same_supporting_conic (const Self& arc) const
  {
    // Check if the two arcs originate from the same conic:
    if (_id == arc._id && _id.is_valid() && arc._id.is_valid())
      return (true);

    // In case both arcs are collinear, check if they have the same
    // supporting lines.
    if (this->_orient == COLLINEAR && arc._orient == COLLINEAR)
    {
      // Construct the two supporting lines and compare them.
      Alg_kernel                             ker;
      typename Alg_kernel::Construct_line_2  construct_line =
                                                 ker.construct_line_2_object();
      typename Alg_kernel::Line_2          l1 = construct_line (this->_source,
                                                                this->_target);
      typename Alg_kernel::Line_2          l2 = construct_line (arc._source,
                                                                arc._target);
      typename Alg_kernel::Equal_2         equal = ker.equal_2_object();

      if (equal (l1, l2))
        return (true);
      
      // Try to compare l1 with the opposite of l2.
      l2 = construct_line (arc._target, arc._source);

      return (equal (l1, l2));
    }
    else if (this->_orient == COLLINEAR || arc._orient == COLLINEAR)
    {
      // Only one arc is collinear, so the supporting curves cannot be the
      // same:
      return (false);
    }

    // Check whether the coefficients of the two supporting conics are equal
    // up to a constant factor.
    Integer        factor1 = 1;
    Integer        factor2 = 1;

    if (CGAL::sign (this->_r) != ZERO)
      factor1 = this->_r;
    else if (CGAL::sign (this->_s) != ZERO)
      factor1 = this->_s;
    else if (CGAL::sign (this->_t) != ZERO)
      factor1 = this->_t;
    else if (CGAL::sign (this->_u) != ZERO)
      factor1 = this->_u;
    else if (CGAL::sign (this->_v) != ZERO)
      factor1 = this->_v;
    else if (CGAL::sign (this->_w) != ZERO)
      factor1 = this->_w;

    if (CGAL::sign (arc._r) != ZERO)
      factor2 = arc._r;
    else if (CGAL::sign (arc._s) != ZERO)
      factor2 = arc._s;
    else if (CGAL::sign (arc._t) != ZERO)
      factor2 = arc._t;
    else if (CGAL::sign (arc._u) != ZERO)

      factor2 = arc._u;
    else if (CGAL::sign (arc._v) != ZERO)
      factor2 = arc._v;
    else if (CGAL::sign (arc._w) != ZERO)
      factor2 = arc._w;

    return (CGAL::compare  (this->_r * factor2, arc._r * factor1) == EQUAL &&
            CGAL::compare  (this->_s * factor2, arc._s * factor1) == EQUAL &&
            CGAL::compare  (this->_t * factor2, arc._t * factor1) == EQUAL &&
            CGAL::compare  (this->_u * factor2, arc._u * factor1) == EQUAL &&
            CGAL::compare  (this->_v * factor2, arc._v * factor1) == EQUAL &&
            CGAL::compare  (this->_w * factor2, arc._w * factor1) == EQUAL);
  }

  /*!
   * Get the i'th order derivative by x of the conic at the point p=(x,y).
   * \param p The point where we derive.
   * \param i The order of the derivatives (either 1, 2 or 3).
   * \param slope_numer The numerator of the slope.
   * \param slope_denom The denominator of the slope.
   * \todo Allow higher order derivatives.
   */
  void _derive_by_x_at (const Point_2& p, const unsigned int& i,
                        Algebraic& slope_numer, Algebraic& slope_denom) const
  {
    if (_is_special_segment())
    {
      // Special treatment for special segments, given by (a*x + b*y + c = 0),
      // so their first-order derivative by x is simply -a/b. The higher-order
      // derivatives are all 0.
      if (i == 1)
      {
        if (CGAL::sign (this->_extra_data_P->b) != NEGATIVE)
        {          
          slope_numer = - this->_extra_data_P->a;
          slope_denom = this->_extra_data_P->b;
        }
        else
        {
          slope_numer = this->_extra_data_P->a;
          slope_denom = - this->_extra_data_P->b;
        }
      }
      else
      {
        slope_numer = 0;
        slope_denom = 1;
      }

      return;
    }

    // The derivative by x of the conic
    //   C: {r*x^2 + s*y^2 + t*xy + u*x + v*y + w = 0}
    // at the point p=(x,y) is given by:
    //
    //           2r*x + t*y + u       alpha
    //   y' = - ---------------- = - -------
    //           2s*y + t*x + v       beta
    //
    const Algebraic  _two = 2;
    const Algebraic  sl_numer = _two*alg_r*p.x() + alg_t*p.y() + alg_u;
    const Algebraic  sl_denom = _two*alg_s*p.y() + alg_t*p.x() + alg_v;
    
    if (i == 1)
    {
      // Make sure that the denominator is always positive.
      if (CGAL::sign (sl_denom) != NEGATIVE)
      {
        slope_numer = -sl_numer;
        slope_denom = sl_denom;
      }
      else
      {
        slope_numer = sl_numer;
        slope_denom = -sl_denom;
      }

      return;
    }

    // The second-order derivative is given by:
    //
    //             s*alpha^2 - t*alpha*beta + r*beta^2     gamma
    //   y'' = -2 ------------------------------------- = -------
    //                           beta^3                    delta
    //
    const Algebraic  sl2_numer = alg_s * sl_numer*sl_numer -
                                 alg_t * sl_numer*sl_denom +
                                 alg_r * sl_denom*sl_denom;
    const Algebraic  sl2_denom = sl_denom*sl_denom*sl_denom;

    if (i == 2)
    {
      // Make sure that the denominator is always positive.
      if (CGAL::sign (sl_denom) != NEGATIVE)
      {
        slope_numer = -_two *sl2_numer;
        slope_denom = sl2_denom;
      }
      else
      {
        slope_numer = _two *sl2_numer;
        slope_denom = -sl2_denom;
      }

      return;
    }

    // The third-order derivative is given by:
    //
    //              (2s*alpha - t*beta) * gamma
    //   y''' = -6 ------------------------------
    //                    beta^2 * delta
    //
    const Algebraic  sl3_numer = (_two * alg_s * sl_numer -
                                  alg_t * sl_denom) * sl2_numer;
    const Algebraic  sl3_denom = sl_denom*sl_denom * sl2_denom;

    if (i == 3)
    {
      // Make sure that the denominator is always positive.
      if (CGAL::sign (sl_denom) != NEGATIVE)
      {
        slope_numer = -6 * sl3_numer;
        slope_denom = sl3_denom;
      }
      else
      {
        slope_numer = 6 * sl3_numer;
        slope_denom = -sl2_denom;
      }

      return;
    }

    // \todo Handle higher-order derivatives as well.
    CGAL_assertion (false);
    return;
  }

  /*!
   * Get the i'th order derivative by y of the conic at the point p=(x,y).
   * \param p The point where we derive.
   * \param i The order of the derivatives (either 1, 2 or 3).
   * \param slope_numer The numerator of the slope.
   * \param slope_denom The denominator of the slope.
   * \todo Allow higher order derivatives.
   */
  void _derive_by_y_at (const Point_2& p, const int& i,
                        Algebraic& slope_numer, Algebraic& slope_denom) const
  {
    if (_is_special_segment())
    {
      // Special treatment for special segments, given by (a*x + b*y + c = 0),
      // so their first-order derivative by x is simply -b/a. The higher-order
      // derivatives are all 0.
      if (i == 1)
      {
        if (CGAL::sign (this->_extra_data_P->a) != NEGATIVE)
        {          
          slope_numer = - this->_extra_data_P->b;
          slope_denom = this->_extra_data_P->a;
        }
        else
        {
          slope_numer = this->_extra_data_P->b;
          slope_denom = - this->_extra_data_P->a;
        }
      }
      else
      {
        slope_numer = 0;
        slope_denom = 1;
      }

      return;
    }

    // The derivative by y of the conic
    //   C: {r*x^2 + s*y^2 + t*xy + u*x + v*y + w = 0}
    // at the point p=(x,y) is given by:
    //
    //           2s*y + t*x + v     alpha
    //   x' = - ---------------- = -------
    //           2r*x + t*y + u      beta
    //
    const Algebraic  _two = 2;
    const Algebraic  sl_numer = _two*alg_s*p.y() + alg_t*p.x() + alg_v;
    const Algebraic  sl_denom = _two*alg_r*p.x() + alg_t*p.y() + alg_u;

    if (i == 1)
    {
      // Make sure that the denominator is always positive.
      if (CGAL::sign (sl_denom) != NEGATIVE)
      {
        slope_numer = -sl_numer;
        slope_denom = sl_denom;
      }
      else
      {
        slope_numer = sl_numer;
        slope_denom = -sl_denom;
      }


      return;
    }

    // The second-order derivative is given by:
    //
    //             r*alpha^2 - t*alpha*beta + s*beta^2
    //   x'' = -2 -------------------------------------
    //                           beta^3
    //
    const Algebraic  sl2_numer = alg_r * sl_numer*sl_numer -