bezier_cache.h

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

// Copyright (c) 2006  Tel-Aviv University (Israel).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you may redistribute it under
// the terms of the Q Public License version 1.0.
// See the file LICENSE.QPL distributed with CGAL.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.3-branch/Arrangement_2/include/CGAL/Arr_traits_2/Bezier_cache.h $
// $Id: Bezier_cache.h 37177 2007-03-17 08:48:10Z afabri $
// 
//
// Author(s)     : Ron Wein     <[e-mail]>

#ifndef CGAL_BEZIER_CACHE_H
#define CGAL_BEZIER_CACHE_H

/*! \file
 * Header file for the _Bezier_cache class.
 */

#include <set>
#include <map>
#include <ostream>

CGAL_BEGIN_NAMESPACE

/*! \class _Bezier_cache
 * Stores all cached intersection points and vertical tangency points.
 */
template <class Nt_traits_> class _Bezier_cache
{
public:

  typedef Nt_traits_                              Nt_traits;

  typedef typename Nt_traits::Integer             Integer;
  typedef typename Nt_traits::Polynomial          Polynomial;
  typedef typename Nt_traits::Algebraic           Algebraic;

  typedef _Bezier_cache<Nt_traits>                Self;

  /// \name Type definitions for the vertical tangency-point mapping.
  //@{
  typedef size_t                                     Curve_id;
  typedef std::list<Algebraic>                       Vertical_tangency_list;
  typedef
    typename Vertical_tangency_list::const_iterator  Vertical_tangency_iter;
  //@}

  /// \name Type definitions for the intersection-point mapping.
  //@{

  /*! \struct Intersection_point_2
   * Representation of an intersection point (in both parameter and physical
   * spaces).
   */
  struct Intersection_point_2
  {
    Algebraic          s;      // The parameter for the first curve.
    Algebraic          t;      // The parameter for the second curve.
    Algebraic          x;      // The x-coordinate.
    Algebraic          y;      // The y-coordinate.

    /*! Constructor. */
    Intersection_point_2 (const Algebraic& _s, const Algebraic& _t,
                          const Algebraic& _x, const Algebraic& _y) :
      s(_s), t(_t),
      x(_x), y(_y)
    {}
  };

  typedef std::pair<Curve_id, Curve_id>              Curve_pair;
  typedef std::pair<Algebraic, Algebraic>            Parameter_pair;
  typedef std::list<Intersection_point_2>            Intersection_list;
  typedef
    typename Intersection_list::const_iterator       Intersection_iter;

private:

  typedef std::map<Curve_id,
                   Vertical_tangency_list>           Vert_tang_map;
  typedef typename Vert_tang_map::value_type         Vert_tang_map_entry;
  typedef typename Vert_tang_map::iterator           Vert_tang_map_iterator;

  /*! \struct Less_curve_pair
   * An auxiliary functor for comparing pair of curve IDs.
   */
  struct Less_curve_pair
  {
    bool operator() (const Curve_pair& cp1, const Curve_pair& cp2) const
    {
      // Compare the pairs of IDs lexicographically.
      return (cp1.first < cp2.first ||
              (cp1.first == cp2.first && cp1.second < cp2.second));
    }
  };

  typedef std::list<Algebraic>                       Parameter_list;
  typedef std::pair<Intersection_list, bool>         Intersection_info;
  typedef std::map<Curve_pair,
                   Intersection_info,
                   Less_curve_pair>                  Intersect_map;
  typedef typename Intersect_map::value_type         Intersect_map_entry;
  typedef typename Intersect_map::iterator           Intersect_map_iterator;

  /*! \struct My_point_2
   * Representation of exact and approximate point coordinates.
   */
  struct My_point_2
  {
    typename Parameter_list::const_iterator  prm_it;
    Algebraic                                x;
    Algebraic                                y;
    double                                   app_x;
    double                                   app_y;

    /*! Default constructor. */
    My_point_2 () :
      app_x (0),
      app_y (0)
    {}

    /*! Constructor. */
    My_point_2 (typename Parameter_list::const_iterator it,
                const Algebraic& _x, const Algebraic& _y) :
      prm_it (it),
      x (_x),
      y (_y),
      app_x (CGAL::to_double(_x)),
      app_y (CGAL::to_double(_y))
    {}

    /*! Get the parameter value. */
    const Algebraic& parameter () const
    {
      return (*prm_it);
    }

    /*! Check if the given two points are equal. */
    bool equals (const My_point_2& other) const
    {
      return (CGAL::compare (x, other.x) == EQUAL &&
              CGAL::compare (y, other.y) == EQUAL);
    }
  };

  /*! \struct Less_distance_iter
   * An auxiliary functor for comparing approximate distances.
   */
  typedef std::list<My_point_2>                      Point_list;
  typedef typename Point_list::iterator              Point_iter;
  typedef std::pair<double, Point_iter>              Distance_iter;

  struct Less_distance_iter
  {
    bool operator() (const Distance_iter& dit1,
                     const Distance_iter& dit2) const
    {
      return (dit1.first < dit2.first);
    }
  };

  // Data members:
  Nt_traits         nt_traits;        /*! Number-type traits. */
  Vert_tang_map     vert_tang_map;    /*! Maps curves to their vertical
                                          tangency parameters. */
  Intersect_map     intersect_map;    /*! Maps curve pairs to their
                                          intersection parameter pairs. */

  // Copy constructor and assignment operator - not supported.
  _Bezier_cache (const Self& );
  Self& operator= (const Self& );

public:

  /*! Constructor. */
  _Bezier_cache ()
  {}

  /*!
   * Get the vertical tangency parameters of a given curve (X(t), Y(t)).
   * \param id The curve ID.
   * \param polyX The coefficients of X(t).
   * \param normX The normalizing factor of X(t).
   * \return A list of parameters 0 < t_1 < ...< t_k < 1 such that X'(t_i) = 0.
   */
  const Vertical_tangency_list&
  get_vertical_tangencies (const Curve_id& id,
                           const Polynomial& polyX, const Integer& normX);

  /*!
   * Get the intersection parameters of two given curve (X_1(s), Y_1(s))
   * and (X_2(t), Y_2(t)).
   * \param id1 The ID of the first curve.
   * \param polyX_1 The coefficients of X_1.
   * \param normX_1 The normalizing factor of X_1.
   * \param polyY_1 The coefficients of Y_1.
   * \param normY_1 The normalizing factor of Y_1.
   * \param id2 The ID of the second curve.
   * \param polyX_2 The coefficients of X_2.
   * \param normX_2 The normalizing factor of X_2.
   * \param polyY_2 The coefficients of Y_2.
   * \param normY_2 The normalizing factor of Y_2.
   * \param do_ovlp Output: Do the two curves overlap.
   * \pre id1 < id2 (swap their order if this is not the case).
   * \return A list of parameter pairs (s_1, t_1), ..., (s_k, t_k) such that
   *         X_1(s_i) = X_2(t_i) and Y_1(s_i) = Y_2(t_i).
   *         Each pair is also associated with the physical point coordinates.
   */
  const Intersection_list&
  get_intersections (const Curve_id& id1,
                     const Polynomial& polyX_1, const Integer& normX_1,
                     const Polynomial& polyY_1, const Integer& normY_1,
                     const Curve_id& id2,
                     const Polynomial& polyX_2, const Integer& normX_2,
                     const Polynomial& polyY_2, const Integer& normY_2,
                     bool& do_ovlp);

private:

  /*!
   * Compute all parameter values of (X_1(s), Y_1(s)) with (X_2(t), Y_2(t)).
   * \param polyX_1 The coefficients of X_1.
   * \param normX_1 The normalizing factor of X_1.
   * \param polyY_1 The coefficients of Y_1.
   * \param normY_1 The normalizing factor of Y_1.
   * \param polyX_2 The coefficients of X_2.
   * \param normX_2 The normalizing factor of X_2.
   * \param polyY_2 The coefficients of Y_2.
   * \param normY_2 The normalizing factor of Y_2.
   * \param s_vals Output: A list of s-values such that (X_1(s), Y_1(s))
   *                       is an intersection point with (X_2(t), Y_2(t)).
   *                       The list is given in ascending order.
   *                       Note that the values are not bounded to [0,1].
   * \return Do the two curves overlap (if they do, s_vals is empty).
   */
  bool _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;

  /*!
   * Compute the resultant of two bivariate polynomials in x and y with
   * respect to y. The bivariate polynomials are given as vectors of
   * polynomials, where bp1[i] is a coefficient of y^i, which is in turn a
   * polynomial in x.
   * \param bp1 The first bivariate polynomial.
   * \param bp2 The second bivariate polynomial.
   * \return The resultant polynomial (a polynomial in x).
   */
  Polynomial _compute_resultant (const std::vector<Polynomial>& bp1,
                                 const std::vector<Polynomial>& bp2) const;
};

// ---------------------------------------------------------------------------
// Get the vertical tangency parameters of a given curve (X(t), Y(t)).
//
template<class NtTraits>
const typename _Bezier_cache<NtTraits>::Vertical_tangency_list&
_Bezier_cache<NtTraits>::get_vertical_tangencies
        (const Curve_id& id,
         const Polynomial& polyX, const Integer& /* normX */)
{
  // Try to find the curve ID in the map.
  Vert_tang_map_iterator    map_iter = vert_tang_map.find (id);

  if (map_iter != vert_tang_map.end())
  {
    // Found in the map: return the cached parameters' list.
    return (map_iter->second);
  }

  // We need to compute the vertical tangency parameters: find all t-values
  // such that X'(t) = 0, and store them in the cache.
  Vertical_tangency_list&  vert_tang_list = vert_tang_map[id];
  const Polynomial&        polyX_der = nt_traits.derive (polyX);
  
  nt_traits.compute_polynomial_roots (polyX_der, 0, 1,
                                      std::back_inserter(vert_tang_list));
  return (vert_tang_list);
}

// ---------------------------------------------------------------------------
// Get the intersection parameters of two given curve (X_1(s), Y_1(s))
// and (X_2(t), Y_2(t)).
//
template<class NtTraits>
const typename _Bezier_cache<NtTraits>::Intersection_list&
_Bezier_cache<NtTraits>::get_intersections
        (const Curve_id& id1,
         const Polynomial& polyX_1, const Integer& normX_1,
         const Polynomial& polyY_1, const Integer& normY_1,
         const Curve_id& id2,
         const Polynomial& polyX_2, const Integer& normX_2,
         const Polynomial& polyY_2, const Integer& normY_2,
         bool& do_ovlp)
{
  CGAL_precondition (id1 < id2);

  // Construct the pair of curve IDs, and try to find it in the map.
  Curve_pair                curve_pair (id1, id2);
  Intersect_map_iterator    map_iter = intersect_map.find (curve_pair);

  if (map_iter != intersect_map.end())
  {
    // Found in the map: return the cached information.
    do_ovlp = map_iter->second.second;
    return (map_iter->second.first);
  }

  // We need to compute the intersection-parameter pairs.
  Intersection_info&      info = intersect_map[curve_pair];

  // Compute s-values and t-values such that (X_1(s), Y_1(s)) and
  // (X_2(t), Y_2(t)) are the intersection points.
  Parameter_list          s_vals;

  do_ovlp = _intersection_params (polyX_1, normX_1, polyY_1, normY_1,
                                  polyX_2, normX_2, polyY_2, normY_2,
                                  s_vals);

  if (do_ovlp)
  {
    // Update the cache and return an empty list of intersection parameters.
    info.second = true;
    return (info.first);
  }

  Parameter_list          t_vals;

  do_ovlp = _intersection_params (polyX_2, normX_2, polyY_2, normY_2,
                                  polyX_1, normX_1, polyY_1, normY_1,
                                  t_vals);

  CGAL_assertion (! do_ovlp);

  // We should now have the same number of s-values and t-values, and we have
  // to pair them together. 
  CGAL_assertion (s_vals.size() == s_vals.size());

  // Compute all points on (X_1, Y_1) that match an s-value from the list.
  typename Parameter_list::iterator  s_it;
  const Algebraic&                   denX_1 = nt_traits.convert (normX_1);
  const Algebraic&                   denY_1 = nt_traits.convert (normY_1);