bezier_curve_2.h

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

   * \param oi Output: An output iterator for the samples. The value-type of
   *                   this iterator must be std::pair<double, double>.
   * \return A past-the-end iterator for the samples.
   */
  template <class OutputIterator>
  OutputIterator sample (const double& t_start, const double& t_end,
                         unsigned int n_samples,
                         OutputIterator oi) const 
  {
    // Convert the coordinates of the control points to doubles.
    typedef Simple_cartesian<double>                        App_kernel;
    typedef App_kernel::Point_2                             App_point_2;

    const unsigned int             n_pts = number_of_control_points();
    std::vector<App_point_2>       app_ctrl_pts (n_pts);
    unsigned int                   k;

    for (k = 0; k < n_pts; k++)
    {
      const Rat_point_2&   pt = control_point(k);

      app_ctrl_pts[k] = App_point_2 (CGAL::to_double (pt.x()),
                                     CGAL::to_double (pt.y()));
    }

    // Sample the approximated curve.
    const unsigned int   n = (n_samples >= 2) ? n_samples : 2; 
    const double         delta_t = (t_end - t_start) / (n - 1);
    App_point_2          p;

    for (k = 0; k < n; k++)
    {
      p = point_on_Bezier_curve_2 (app_ctrl_pts.begin(), app_ctrl_pts.end(),
                                   t_start + k * delta_t);

      *oi = std::make_pair (p.x(), p.y());
      ++oi;
    }

    return (oi);
  }

  /*!
   * Compute all parameter values t such that the x-coordinate of B(t) is x0.
   * Note that the function does not return only values between 0 and 1, so
   * the output t-values may belong to the imaginary continuation of the curve.
   * \param x0 The given x-coordinate.
   * \param oi Output: An output iterator for the t-values.
   * \return A past-the-end iterator.
   */
  template <class OutputIterator>
  OutputIterator get_t_at_x (const Rational& x0,
                             OutputIterator oi) const
  {
    return (_solve_t_values (_rep().x_polynomial(), _rep().x_norm(), x0,
                             oi));
  }

  /*!
   * Compute all parameter values t such that the y-coordinate of B(t) is y0.
   * Note that the function does not return only values between 0 and 1, so
   * the output t-values may belong to the imaginary continuation of the curve.
   * \param y0 The given y-coordinate.
   * \param oi Output: An output iterator for the t-values.
   * \return A past-the-end iterator.
   */
  template <class OutputIterator>
  OutputIterator get_t_at_y (const Rational& y0,
                             OutputIterator oi) const
  {
    return (_solve_t_values (_rep().y_polynomial(), _rep.y_norm(), y0,
                             oi));
  }

  /*!
   * Check if the two curves have the same support.
   */
  bool has_same_support (const Self& bc) const;

  /*!
   * Get the bounding box of the curve.
   */
  const Bbox_2& bbox () const
  {
    return (_rep()._bbox);
  }

private:

  // Get the representation.
  inline const Bcv_rep& _rep () const
  {
    return (*(this->ptr()));
  }

  inline Bcv_rep& _rep ()
  {
    return (*(this->ptr()));
  }

  /*!
   * Compute all parameter values t, such that P(t) = val.
   * \param poly The polynomial.
   * \param norm Its normalizing factor.
   * \param val The required value.
   * \param oi Output: An output iterator for the t-values.
   * \return A past-the-end iterator.
   */
  template <class OutputIterator>
  OutputIterator _solve_t_values (const Polynomial& poly,
                                  const Integer& norm,
                                  const Rational& val,
                                  OutputIterator oi) const
  {
    // Construct the polynomial P(t) - val = 0:
    Nt_traits             nt_traits;
    const Integer         numer = nt_traits.numerator (val);
    const Integer         denom = nt_traits.denominator (val);
    const int             deg = nt_traits.degree (poly);
    Integer              *coeffs = new Integer [deg + 1];
    int                   k;

    for (k = 1; k <= deg; k++)
      coeffs[k] = nt_traits.get_coefficient (poly, k) * denom;
    coeffs[0] = nt_traits.get_coefficient (poly, 0) * denom -
                numer * norm;

    // Solve the polynomial and obtain the t-values.
    OutputIterator  end = nt_traits.compute_polynomial_roots
                            (nt_traits.construct_polynomial (coeffs, deg),
                             oi);

    delete[] coeffs;
    return (end);
  }
};

/*!
 * Exporter for Bezier curves.
 */
template <class Rat_kernel, class Alg_kernel, class Nt_traits, 
          class Bounding_traits>
std::ostream& 
operator<< (std::ostream& os, 
            const _Bezier_curve_2<Rat_kernel, Alg_kernel, Nt_traits,
                                  Bounding_traits> & bc)
{
  const unsigned int  n = bc.number_of_control_points();
  unsigned int        k;

  os << n;
  for (k = 0; k < n; k++)
    os << "  " << bc.control_point(k);

  return (os);
}

/*!
 * Importer for Bezier curves.
 */
template <class Rat_kernel, class Alg_kernel, class Nt_traits, 
          class Bounding_traits>
std::istream& 
operator>> (std::istream& is, 
            _Bezier_curve_2<Rat_kernel, Alg_kernel, Nt_traits,
                            Bounding_traits> & bc)
{
  // Read the number of control points.
  unsigned int  n;

  is >> n;

  // Read the control points.
  std::vector<typename Rat_kernel::Point_2>   ctrl_pts (n);
  unsigned int                                k;

  for (k = 0; k < n; k++)
    is >> ctrl_pts[k];

  // Set the Bezier curve.
  bc = _Bezier_curve_2<Rat_kernel, Alg_kernel, Nt_traits,
                       Bounding_traits> (ctrl_pts.begin(),
                                         ctrl_pts.end());

  return (is);
}

// ---------------------------------------------------------------------------
// Construct the representation of X(t) and Y(t).
//
template <class RatKer, class AlgKer, class NtTrt, class BndTrt>
void _Bezier_curve_2_rep<RatKer, AlgKer, NtTrt,
                         BndTrt>::_construct_polynomials () const
{
  const int        n = _ctrl_pts.size() - 1;
  Rational        *coeffsX = new Rational [n + 1];
  Rational        *coeffsY = new Rational [n + 1];
  const Rational   rat_zero = Rational (0);
  int              j, k;
  
  CGAL_precondition_msg (n > 0,
                         "There must be at least 2 control points.");
  
  for (j = 0; j <= n; j++)
    coeffsX[j] = coeffsY[j] = rat_zero;
  
  // Compute the rational coefficients, given by the formula:
  //
  //                     n
  //                   *****
  //                   *   *      / n \   k        n-k
  //   (X(t), Y(t)) =   *    p_k (     ) t  (1 - t)
  //                   *   *      \ k /
  //                   *****
  //                    k=0
  //
  Rational               px, py;
  Integer                n_over_k_j;
  bool                   even_exp;
  
  typename Control_point_vec::const_iterator pts_begin = _ctrl_pts.begin();  
  typename Control_point_vec::const_iterator pts_end = _ctrl_pts.end();  

  for (k = 0; pts_begin != pts_end; ++pts_begin, k++)
  {
    px = pts_begin->x();
    py = pts_begin->y();
    
    // By simplifying (1 - t)^(n-k) we obtain that the k'th expression of
    // the above sum is given by:
    //
    //     n-k
    //    *****
    //    *   *     j            n!         j+k
    //     *    (-1)  p_k ---------------- t
    //    *   *            j! k! (n-k-j)!
    //    *****
    //     j=0
    //
    even_exp = true; 
    for (j = 0; j <= n - k; j++)
    {
      n_over_k_j = _choose (n, k, j);
      
      if (even_exp)
      {
        // We should add the current values to the coefficients of the
        // monomial t^(n_j).
        coeffsX[j + k] += px * n_over_k_j;
        coeffsY[j + k] += py * n_over_k_j;
      }
      else
      {
        // We should subtract the current values from the coefficients of the
        // monomial t^(n_j).
        coeffsX[j + k] -= px * n_over_k_j;
        coeffsY[j + k] -= py * n_over_k_j;
      }
      
      // As we increment j, negate the "even" flag for the exponent (n-j).
      even_exp = !even_exp;
    } // loop on j.
  } // loop on k.
  
  // Convert the rational polynomials to polynomials with rational
  // coefficients (plus normalizing factors).
  Nt_traits        nt_traits;
  p_polyX = new Polynomial;
  p_normX = new Integer;
  p_polyY = new Polynomial;
  p_normY = new Integer;
  
  nt_traits.construct_polynomial (coeffsX, n,
                                  *p_polyX, *p_normX);
  delete[] coeffsX;
  
  nt_traits.construct_polynomial (coeffsY, n,
                                  *p_polyY, *p_normY);
  delete[] coeffsY;

  CGAL_assertion (nt_traits.degree (*p_polyX) > 0);
  CGAL_assertion (nt_traits.degree (*p_polyY) > 0);

  return;
}

// ---------------------------------------------------------------------------
// Compute the value of n! / (j! k! (n-k-j)!).
//
template <class RatKer, class AlgKer, class NtTrt, class BndTrt>
typename _Bezier_curve_2_rep<RatKer, AlgKer, NtTrt, BndTrt>::Integer
_Bezier_curve_2_rep<RatKer, AlgKer, NtTrt, BndTrt>::_choose (int n,
                                                             int j,
                                                             int k) const
{
  Integer   reduced_fact = 1;
  Integer   j_fact = 1, k_fact = 1;
  int       i;
  
  for (i = n - k - j + 1; i <= n; i++)
    reduced_fact *= Integer (i);
  
  for (i = 2; i <= j; i++)
    j_fact *= Integer (i);
  
  for (i = 2; i <= k; i++)
    k_fact *= Integer (i);
  
  return (reduced_fact / (j_fact * k_fact));
}

// ---------------------------------------------------------------------------
// Compute a point on the Bezier curve B(t) given a rational t-value.
//
template <class RatKer, class AlgKer, class NtTrt, class BndTrt>
typename _Bezier_curve_2<RatKer, AlgKer, NtTrt, BndTrt>::Rat_point_2
_Bezier_curve_2<RatKer, AlgKer, NtTrt, BndTrt>::operator()
        (const Rational& t) const 
{
  // Check for extermal t values (either 0 or 1).
  const CGAL::Sign   sign_t = CGAL::sign (t);

  CGAL_precondition (sign_t != NEGATIVE);

  if (sign_t == ZERO)
  {
    // If t is 0, simply return the first control point.
    return (_rep()._ctrl_pts[0]);
  }

  Comparison_result  res = CGAL::compare (t, Rational(1));

  CGAL_precondition (res != LARGER);

  if (res == EQUAL)
  {
    // If t is 1, simply return the first control point.
    return (_rep()._ctrl_pts[_rep()._ctrl_pts.size() - 1]);
  }

  // Evaluate the point for 0 < t < 1.
  if (! _rep().has_polynomials())
  {
    // Use de Casteljau's algorithm to evaluate the polynomial at t.
    return (point_on_Bezier_curve_2 (_rep()._ctrl_pts.begin(),
                                     _rep()._ctrl_pts.end(),
                                     t));
  }

  // Compute the x and y coordinates using the X(t), Y(t) polynomials.
  Rational           x, y;
  Nt_traits          nt_traits;
    
  x = nt_traits.evaluate_at (_rep().x_polynomial(), t) /
    Rational (_rep().x_norm(), 1);
  y = nt_traits.evaluate_at (_rep().y_polynomial(), t) /
    Rational (_rep().y_norm(), 1);
  
  // Return the point.
  return (Rat_point_2 (x, y));
}

// ---------------------------------------------------------------------------
// Compute a point on the Bezier curve B(t) given an algebraic t-value.
//
template <class RatKer, class AlgKer, class NtTrt, class BndTrt>
typename _Bezier_curve_2<RatKer, AlgKer, NtTrt, BndTrt>::Alg_point_2
_Bezier_curve_2<RatKer, AlgKer, NtTrt, BndTrt>::operator()
          (const Algebraic& t) const
{
  // Check for extermal t values (either 0 or 1).
  Nt_traits          nt_traits;
  const CGAL::Sign   sign_t = CGAL::sign (t);

  CGAL_precondition (sign_t != NEGATIVE);

  if (sign_t == ZERO)
  {
    // If t is 0, simply return the first control point.
    const Rat_point_2&  p_0 = _rep()._ctrl_pts[0];
    
    return (Alg_point_2 (nt_traits.convert (p_0.x()),
                         nt_traits.convert (p_0.y())));
  }

  Comparison_result  res = CGAL::compare (t, Algebraic(1));

  CGAL_precondition (res != LARGER);

  if (res == EQUAL)
  {
    // If t is 1, simply return the first control point.
    const Rat_point_2&  p_n = _rep()._ctrl_pts[_rep()._ctrl_pts.size() - 1];
    
    return (Alg_point_2 (nt_traits.convert (p_n.x()),
                         nt_traits.convert (p_n.y())));
  }

  // IDDO: Check if it is worth evaluating this using de Casteljau 
  //       (since the parameter value is algebraic).
  // IDDO: Check if there is a way to know that an Algebraic is a rational
  //       number...

  // The t-value is between 0 and 1: Compute the x and y coordinates.
  const Algebraic    x = nt_traits.evaluate_at (_rep().x_polynomial(), t) /
                         nt_traits.convert (_rep().x_norm());
  const Algebraic    y = nt_traits.evaluate_at (_rep().y_polynomial(), t) /
                         nt_traits.convert (_rep().y_norm());

  return (Alg_point_2 (x, y));
}

// ---------------------------------------------------------------------------
// Check if the two curves have the same support.
//
template <class RatKer, class AlgKer, class NtTrt, class BndTrt>
bool _Bezier_curve_2<RatKer, AlgKer, NtTrt, BndTrt>::has_same_support
        (const Self& bc) const
{
  // If one curve is of degree d1 and the other of degree d2, there can be
  // at most d1*d2 intersection points between them.
  const int      deg1 = number_of_control_points() - 1;
  const int      deg2 = bc.number_of_control_points() - 1;
  const int      n_samples = deg1*deg2;
  Rat_point_2    p1;
  int            k;
  
  for (k = 0; k <= n_samples; k++)
  {
    // Compute p1 = B1(k/n_samples), where B1 is (*this) curve.
    if (k == 0)
      p1 = (_rep()._ctrl_pts[0]);
    else if (k == 1)
      p1 = (_rep()._ctrl_pts[_rep()._ctrl_pts.size() - 1]);
    else
      p1 = this->operator() (Rational (k, n_samples));
    
    // Get all t-values such that the x-coordinate of B2(t) equals x1,
    // and check if there exists a t-value such that the y-coordinate of
    // b2(t) equals the y-coordinate of p1.
    std::list<Algebraic>                           t_vals;
    typename std::list<Algebraic>::const_iterator  t_iter;
    Nt_traits                             nt_traits;
    const Algebraic&                      y1 = nt_traits.convert (p1.y());
    bool                                  eq_y = false;
    
    bc.get_t_at_x (p1.x(), std::back_inserter(t_vals));
    
    for (t_iter = t_vals.begin(); t_iter != t_vals.end(); ++t_iter)
    {
      const Alg_point_2&  p2 = bc (*t_iter, false);
      
      if (CGAL::compare (y1, p2.y()) == CGAL::EQUAL)
      {
        eq_y = true;
        break;
      }
    }
    
    // If we found a point on B1 which is not of B2, the two curves do not
    // have the same support.
    if (! eq_y)
      return (false);
  }
  
  // If we reached here, we found (d1*d2 + 1) common points of B1 and B2.
  // This means they have the same support.
  return (true);
}

CGAL_END_NAMESPACE

#endif