bezier_cache.h

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

  Point_list                         pts1;

  for (s_it = s_vals.begin(); s_it != s_vals.end(); ++s_it)
  {
    const Algebraic&  x = nt_traits.evaluate_at (polyX_1, *s_it) / denX_1;
    const Algebraic&  y = nt_traits.evaluate_at (polyY_1, *s_it) / denY_1;

    pts1.push_back (My_point_2 (s_it, x, y));
  }

  // Compute all points on (X_2, Y_2) that match a t-value from the list.
  typename Parameter_list::iterator  t_it;
  const Algebraic&                   denX_2 = nt_traits.convert (normX_2);
  const Algebraic&                   denY_2 = nt_traits.convert (normY_2);
  Point_list                         pts2;

  for (t_it = t_vals.begin(); t_it != t_vals.end(); ++t_it)
  {
    const Algebraic&  x = nt_traits.evaluate_at (polyX_2, *t_it) / denX_2;
    const Algebraic&  y = nt_traits.evaluate_at (polyY_2, *t_it) / denY_2;

    pts2.push_back (My_point_2 (t_it, x, y));
  }

  // Go over the points in the pts1 list.
  const bool                x2_simpler = nt_traits.degree(polyX_2) <
                                         nt_traits.degree(polyX_1);
  const bool                y2_simpler = nt_traits.degree(polyY_2) <
                                         nt_traits.degree(polyY_1);
  Point_iter                pit1;
  Point_iter                pit2;
  double                    dx, dy;
  Algebraic                 s, t;
  const Algebraic           one (1);
  unsigned int              k;

  for (pit1 = pts1.begin(); pit1 != pts1.end(); ++pit1)
  {
    // Construct a vector of distances from the current point to all other
    // points in the pts2 list.
    const int                     n_pts2 = pts2.size();
    std::vector<Distance_iter>    dist_vec (n_pts2);

    for (k = 0, pit2 = pts2.begin(); pit2 != pts2.end(); k++, ++pit2)
    {
      // Compute the approximate distance between the teo current points.
      dx = pit1->app_x - pit2->app_x;
      dy = pit1->app_y - pit2->app_y;

      dist_vec[k] = Distance_iter (dx*dx + dy*dy, pit2);
    }
      
    // Sort the vector according to the distances from *pit1.
    std::sort (dist_vec.begin(), dist_vec.end(), 
               Less_distance_iter());

    // Go over the vector entries, starting from the most distant from *pit1
    // to the closest and eliminate pairs of points (we expect that
    // eliminating the distant points is done easily). We stop when we find
    // a pait for *pit1 or when we are left with a single point.
    bool                    found = false;
    const Algebraic&        s = pit1->parameter();
    Algebraic               t;

    for (k = n_pts2 - 1; !found && k > 0; k--)
    {
      pit2 = dist_vec[k].second;

      if (pit1->equals (*pit2))
      {
        // Obtain the parameter value, and try to simplify the representation
        // of the intersection point.
        t = pit2->parameter();

        if (x2_simpler)
          pit1->x = pit2->x;
        if (y2_simpler)
          pit1->y = pit2->y;

        // Remove this point from pts2, as we found a match for it.
        pts2.erase (pit2);
        found = true;
      }
    }

    if (! found)
    {
      // We are left with a single point - pair it with *pit1.
      pit2 = dist_vec[0].second;

      // Obtain the parameter value, and try to simplify the representation
      // of the intersection point.
      t = pit2->parameter();

      if (x2_simpler)
        pit1->x = pit2->x;
      if (y2_simpler)
        pit1->y = pit2->y;

      // Remove this point from pts2, as we found a match for it.
      pts2.erase (pit2);
    }

    // Check if the s- and t-values both lie in the legal range of [0,1].
    // If so, report the the parameter pair as an intersection.
    if (CGAL::sign (s) != NEGATIVE && CGAL::compare (s, one) != LARGER &&
        CGAL::sign (t) != NEGATIVE && CGAL::compare (t, one) != LARGER)
    {
      info.first.push_back (Intersection_point_2 (s, t,
                                                  pit1->x, pit1->y));
    }
  }

  return (info.first);
}

// ---------------------------------------------------------------------------
// Compute all parameter values of (X_1(s), Y_1(s)) with (X_2(t), Y_2(t)).
//
template<class NtTraits>
bool _Bezier_cache<NtTraits>::_intersection_params
        (const Polynomial& polyX_1, const Integer& normX_1,
         const Polynomial& polyY_1, const Integer& normY_1,
         const Polynomial& polyX_2, const Integer& normX_2,
         const Polynomial& polyY_2, const Integer& normY_2,
         Parameter_list& s_vals) const
{
  // Clear the output parameter list.
  if (! s_vals.empty())
    s_vals.clear();

  // We are looking for the s-values s_0 such that (X_1(s_0), Y_1(s_0)) is
  // an intersection with some point on (X_2(t), Y_2(t)). We therefore have
  // to solve the system of bivariate polynomials in s and t:
  //    I: X_2(t) - X_1(s) = 0
  //   II: Y_2(t) - Y_1(s) = 0
  //
  Integer                  coeff;
  int                      k;

  // Consruct the bivariate polynomial that corresponds to Equation I.
  // Note that we represent a bivariate polynomial as a vector of univariate
  // polynomials, whose i'th entry corresponds to the coefficient of t^i,
  // which is in turn a polynomial it s.
  const int                degX_2 = nt_traits.degree (polyX_2);
  std::vector<Polynomial>  coeffsX_st (degX_2 + 1);

  for (k = degX_2; k >= 0; k--)
  {
    coeff = nt_traits.get_coefficient (polyX_2, k) * normX_1;
    coeffsX_st[k] = nt_traits.construct_polynomial (&coeff, 0);
  }
  coeffsX_st[0] = coeffsX_st[0] - nt_traits.scale (polyX_1, normX_2);

  // Consruct the bivariate polynomial that corresponds to Equation II.
  const int                degY_2 = nt_traits.degree (polyY_2);
  std::vector<Polynomial>  coeffsY_st (degY_2 + 1);
    
  for (k = degY_2; k >= 0; k--)
  {
    coeff = nt_traits.get_coefficient (polyY_2, k) * normY_1;
    coeffsY_st[k] = nt_traits.construct_polynomial (&coeff, 0);
  }
  coeffsY_st[0] = coeffsY_st[0] - nt_traits.scale (polyY_1, normY_2);

  // Compute the resultant of the two bivariate polynomials and obtain
  // a polynomial in s.
  Polynomial            res = _compute_resultant (coeffsX_st, coeffsY_st);

  if (nt_traits.degree (res) < 0)
  {
    // If the resultant is identiaclly zero, then the two curves overlap.
    return (true);
  }

  // Compute the roots of the resultant polynomial and mark that the curves do
  // not overlap.
  nt_traits.compute_polynomial_roots (res, std::back_inserter (s_vals));
  return (false);
}

// ---------------------------------------------------------------------------
// Compute the resultant of two bivariate polynomials.
//
template<class NtTraits>
typename _Bezier_cache<NtTraits>::Polynomial
_Bezier_cache<NtTraits>::_compute_resultant
        (const std::vector<Polynomial>& bp1,
         const std::vector<Polynomial>& bp2) const
{
  // Create the Sylvester matrix of polynomial coefficients. Also prepare
  // the exp_fact vector, that represents the normalization factor (see
  // below).
  const int        m = bp1.size() - 1;
  const int        n = bp2.size() - 1;
  const int        dim = m + n;
  const Integer    zero = 0;
  const Polynomial zero_poly = nt_traits.construct_polynomial (&zero, 0);
  int              i, j, k;

  std::vector<std::vector<Polynomial> >  mat (dim);
  std::vector <int>                      exp_fact (dim);
  
  for (i = 0; i < dim; i++)
  {
    mat[i].resize (dim);
    exp_fact[i] = 0;
    
    for (j = 0; j < dim; j++)
      mat[i][j] = zero_poly;
  }
  
  // Initialize it with copies of the two bivariate polynomials.
  for (i = 0; i < n; i++)
    for (j = m; j >= 0; j--)
      mat[i][i + j] = bp1[j];
  
  for (i = 0; i < m; i++)
    for (j = n; j >= 0; j--)
      mat[n + i][i + j] = bp2[j];
  
  // Perform Gaussian elimination on the Sylvester matrix. The goal is to
  // reach an upper-triangular matrix, whose diagonal elements are mat[0][0]
  // to mat[dim-1][dim-1], such that the determinant of the original matrix
  // is given by:
  //
  //              dim-1
  //             *******
  //              *   *  mat[i][i]
  //              *   *
  //               i=0
  //      ---------------------------------
  //         dim-1
  //        *******            exp_fact[i]
  //         *   *  (mat[i][i])
  //         *   *
  //          i=0
  //
  bool             found_row;
  Polynomial       value;
  int              sign_fact = 1;

  for (i = 0; i < dim; i++)
  {
    // Check if the current diagonal value is a zero polynomial.
    if (nt_traits.degree (mat[i][i]) < 0)
    {
      // If the current diagonal value is a zero polynomial, try to replace
      // the current i'th row with a row with a higher index k, such that
      // mat[k][i] is not a zero polynomial.
      
      found_row = false;
      for (k = i + 1; k < dim; k++)
      {
        if (nt_traits.degree (mat[k][i]) <= 0)
        {
          found_row = true;
          break;
        }
      }
      
      if (found_row)
      {
        // Swap the i'th and the k'th rows (note that we start from the i'th
        // column, because the first i entries in every row with index i or
        // higher should be zero by now).
        for (j = i; j < dim; j++)
        {
          value = mat[i][j];
          mat[i][j] = mat[k][j];
          mat[k][j] = value;
        }
        
        // Swapping two rows should change the sign of the determinant.
        // We therefore swap the sign of the normalization factor.
        sign_fact = -sign_fact;
      }
      else
      {
        // In case we could not find a non-zero value, the matrix is
        // singular and its determinant is a zero polynomial.
        return (mat[i][i]);
      }
    }
    
    // Zero the whole i'th column of the following rows.
    for (k = i + 1; k < dim; k++)
    {
      if (nt_traits.degree (mat[k][i]) >= 0)
      {
        value = mat[k][i];
        mat[k][i] = zero_poly;
        
        for (j = i + 1; j < dim; j++)
        {
          mat[k][j] = mat[k][j] * mat[i][i] - mat[i][j] * value; 
        }
        
        // We multiplied the current row by the i'th diagonal entry, thus
        // multipling the determinant value by it. We therefore increment
        // the exponent of mat[i][i] in the normalization factor.
        exp_fact[i] = exp_fact[i] + 1;
      }
    }
  }

  // Now the determinant is simply the product of all diagonal items,
  // divided by the normalizing factor.
  const Integer    sgn (sign_fact);
  Polynomial       det_factor = nt_traits.construct_polynomial (&sgn, 0);
  Polynomial       diag_prod = mat[dim - 1][dim - 1];
  
  CGAL_assertion (exp_fact [dim - 1] == 0);
  for (i = dim - 2; i >= 0; i--)
  {
    // Try to avoid unnecessary multiplications by ignoring the current
    // diagonal item if its exponent in the normalization factor is greater
    // than 0.
    if (exp_fact[i] > 0)
    {
      exp_fact[i] = exp_fact[i] - 1;
    }
    else
    {
      diag_prod *= mat[i][i];
    }
    
    for (j = 0; j < exp_fact[i]; j++)
      det_factor *= mat[i][i];
  }
  
  // In case of a trivial normalization factor, just return the product
  // of diagonal elements.
  if (nt_traits.degree(det_factor) == 0)
    return (diag_prod);
  
  // Divide the product of diagonal elements by the normalization factor
  // and obtain the determinant (note that we should have no remainder).
  Polynomial       det, rem;
  
  det = nt_traits.divide (diag_prod, det_factor, rem);

  CGAL_assertion (nt_traits.degree(rem) < 0);
  return (det);
}

CGAL_END_NAMESPACE

#endif