bezier_bounding_rational_traits.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_bounding_rational_traits.h $
// $Id: Bezier_bounding_rational_traits.h 39309 2007-07-05 10:23:02Z golubevs $
// 
//
// Author(s)     : Iddo Hanniel <[e-mail]>

#ifndef CGAL_BEZIER_BOUNDING_RATIONAL_TRAITS_H
#define CGAL_BEZIER_BOUNDING_RATIONAL_TRAITS_H

/*! \file
 * Definition of the Bezier_bounding_rational_traits<Kernel> class.
 */

#include <CGAL/functional.h>
#include <CGAL/Cartesian.h>
#include <CGAL/Polygon_2_algorithms.h> 
#include <CGAL/Arr_traits_2/de_Casteljau_2.h> 

#include <deque>
#include <list>
#include <math.h>

CGAL_BEGIN_NAMESPACE

template <class _NT> 
struct _Bez_point_bound ;

template <class _NT> 
struct _Bez_point_bbox ;

/*!\ class Bezier_bounding_rational_traits
 * A traits class for performing bounding operations on Bezier curves,
 * assuming the NT is rational.
 */
template <typename _Kernel>
class Bezier_bounding_rational_traits
{
public:
  typedef _Kernel Kernel; // Used in polygon and point/vector operations
  typedef typename Kernel::FT NT; // The number type used to represent control points and t-parameters.
  typedef typename Kernel::Point_2 Point_2; // The control point type.
  typedef typename Kernel::Vector_2 Vector_2;
  typedef typename Kernel::Line_2 Line_2;
  typedef std::deque<Point_2>                   Control_point_vec;

  typedef  _Bez_point_bound<NT>                 Bez_point_bound;
  typedef  _Bez_point_bbox<NT>                  Bez_point_bbox;

  // Auxilary structure for output of intersection points.
  struct Bound_pair
  {
    Bez_point_bound bound1;
    Bez_point_bound bound2;
    Bez_point_bbox bbox;
  
    Bound_pair(const Bez_point_bound& b1, const Bez_point_bound& b2,
      const Bez_point_bbox& bb) : bound1(b1), bound2(b2), bbox(bb) {}
    Bound_pair() : bound1(), bound2(), bbox() {}
  };

  typedef std::list<Bound_pair> Bound_pair_lst;

  //TODO - add here typedefs for vertical tangency output (pair<..>).

  Kernel  kernel;
  NT      _can_refine_bound;  

  // a ctr that enables to change can_refine_bound.
  // iddo: important! we need to make the arr traits receive this traits as a reference,
  // since it now contains state.
  // The line below currently causes a runtime bug, so for now.
  //Bezier_bounding_rational_traits(double can_refine_bound = std::pow(2,-53))
                                  //0.00000000000000000000001)
  Bezier_bounding_rational_traits (double can_refine_bound =
                                   0.00000000000000011)
    : _can_refine_bound(can_refine_bound)
  {}
  
  // iddo: can_refine() function can be based on Euclidean or parametric space.
  // For rationals it is a simple parametric-space check.
  bool can_refine(const Control_point_vec& , const NT& l, const NT& r)
  {
    // iddo: might be made more efficient if needed (based on denom...)
    if (r-l < _can_refine_bound)
      return false;
    return true;  
  }

  // Can be made a template function <Control_point_vec, NT>
  /*
  void DeCasteljau(const Control_point_vec& cp, const NT& t0,
    Control_point_vec& lcp, Control_point_vec& rcp) const
  {
    const NT t1 = NT(1) - t0; 

    lcp.push_back(cp.front());
    rcp.push_front(cp.back());

    Control_point_vec aux(cp.begin(), cp.end());
    unsigned int i=0;
    for (; i<cp.size()-1; ++i)
    {
      Control_point_vec aux1(cp.size()-i-1); //in every step of the subdivision we reduce by 1 the size of the auxilary vector
      unsigned int j=0;
      for (; j<cp.size()-i-1; ++j)
      {
        aux1[j] = Point_2(t1*(aux[j].x()) + t0*(aux[j+1].x()), t1*(aux[j].y()) + t0*(aux[j+1].y())); 
      }
      lcp.push_back(aux1.front());
      rcp.push_front(aux1.back());

      aux = aux1;
    }
  }
  */

  ////////////////////////////////////////////////////////////////
  //Intersection points
  ///////////////////////////////////////////////////////////////
  void CrvCrvInter(const Control_point_vec& cp1, const Control_point_vec& cp2, 
                  Bound_pair_lst& intersection_pairs)
  {
    //iddo: in order to support interval, we will need to convert Rational to
    //      interval do we need to do it outside the call to the function or
    //      here (needs conversion from Rat_Point...).
    
    std::set<NT> endpoint_intersections; // Auxilary variable needed in the recursion
    // to identify already existing points due to endpoint intersections.

    // Call to recursive function.
    _CrvCrvInterAux (cp1, 0, 1,
                     cp2, 0, 1,
                     true,
                     endpoint_intersections,
                     intersection_pairs);
  }

  void refine_intersection_point (const Bound_pair& intersection_point,
                                  Bound_pair& refined_point)
  {
    // Assumes we do not need to check span, intersection_point has all the
    //necessary data.
    CGAL_precondition(intersection_point.bound1.can_refine);
    CGAL_precondition(intersection_point.bound2.can_refine);
    // The following precondition makes sure the point is not rational already
    // (in which case it would not need refinement).
    //CGAL_precondition(intersection_point.bbox.min_x != intersection_point.bbox.max_x);

    const Control_point_vec& cv1 = intersection_point.bound1.bounding_polyline;
    const NT& l1 = intersection_point.bound1.t_min;
    const NT& r1 = intersection_point.bound1.t_max;
    const Control_point_vec& cv2 = intersection_point.bound2.bounding_polyline;
    const NT& l2 = intersection_point.bound2.t_min;
    const NT& r2 = intersection_point.bound2.t_max;

    /* Subdivide the two curves and recurse. */
    // This is probably not needed as it will be handled
    // in the recursion.
    bool can_refine1 = can_refine(cv1, l1, r1);
    bool can_refine2 = can_refine(cv2, l2, r2);
    if (!can_refine1 || !can_refine2) 
    {
      // Failed in the bounded approach - stop the subdivision and indicate
      // that the subdivision has failed.
      refined_point = intersection_point;
      refined_point.bound1.can_refine = false;
      refined_point.bound2.can_refine = false;

      return;
    }

    Control_point_vec cv1a, cv1b;
    
    bisect_control_polygon_2 (cv1.begin(), cv1.end(),
                              std::back_inserter(cv1a),
                              std::front_inserter(cv1b));
    //DeCasteljau(cv1, 0.5, cv1a, cv1b);
    NT t_mid1 = NT(0.5) * (l1 + r1);

    Control_point_vec cv2a, cv2b; 

    bisect_control_polygon_2 (cv2.begin(), cv2.end(),
                              std::back_inserter(cv2a),
                              std::front_inserter(cv2b));
    //DeCasteljau(cv2, 0.5, cv2a, cv2b);
    NT t_mid2 = NT(0.5) * (l2 + r2);

    // Iddo: Catch situations that we cannot refine, in roder to prevent
    //       having multiple points as the refined point.
    bool can_refine_pt_cv1a = can_refine (cv1a, l1, t_mid1);
    bool can_refine_pt_cv1b = can_refine (cv1b, t_mid1, r1);
    bool can_refine_pt_cv2a = can_refine (cv2a, l2, t_mid2);
    bool can_refine_pt_cv2b = can_refine (cv2b, t_mid2, r2);
    
    if (!can_refine_pt_cv1a || !can_refine_pt_cv1b ||
        !can_refine_pt_cv2a || !can_refine_pt_cv2b)
    {
      // Construct a Bez_point_bound that is inconclusive:
      // all parameters found so far + can_refine = false, with the bounding
      // box of cv1, and add to output.
      Bez_point_bbox bbox1;

      cp_bbox (cv1, bbox1);
      refined_point = 
        Bound_pair (Bez_point_bound (cv1, l1, r1,
                                     Bez_point_bound::INTERSECTION_PT, false),
                    Bez_point_bound (cv2, l2, r2,
                                     Bez_point_bound::INTERSECTION_PT, false),
                    bbox1);

      return;
    }

    //Recursion:
    Bound_pair_lst intersection_pairs;
    std::set<NT> endpoint_intersections;
    bool should_check_span = false; //this precondition is hard to check efficiently.

    _CrvCrvInterAux (cv1a, l1, t_mid1, cv2a, l2, t_mid2,
                     should_check_span, endpoint_intersections,
                     intersection_pairs);
    _CrvCrvInterAux (cv1a, l1, t_mid1, cv2b, t_mid2, r2,
                     should_check_span, endpoint_intersections,
                     intersection_pairs);
    _CrvCrvInterAux (cv1b, t_mid1, r1, cv2a, l2, t_mid2,
                     should_check_span, endpoint_intersections,
                     intersection_pairs);
    _CrvCrvInterAux (cv1b, t_mid1, r1, cv2b, t_mid2, r2,
                     should_check_span, endpoint_intersections,
                     intersection_pairs);

    refined_point = intersection_pairs.front(); 

    CGAL_assertion (intersection_pairs.size() == 1);

    return;
  }

  ////////////////////////////////////////////////////////////////
  // Vertical tangency points
  ///////////////////////////////////////////////////////////////

  // Vertical tangency points (cp is a Bez curve between t0-t1).
  // Output are pairs of <bound,bbox> representing the vertical tangency points.
  void vertical_tangency_points
      (const Control_point_vec& cp,
       const NT& t0, const NT& t1,
       std::list<std::pair<Bez_point_bound, Bez_point_bbox> >& tangency_points)
  {
    //TODO - handle special case of degree two curves.

    bool can_refine_pt = can_refine(cp, t0, t1);
    if (!can_refine_pt)
    {
      // Construct a Bez_point_bound that is inconclusive:
      // all parameters found so far + can_refine = false,
      // and add to output.
      Bez_point_bbox bbox;
      cp_bbox(cp, bbox);
      tangency_points.push_back
        (std::make_pair(Bez_point_bound(cp, t0, t1,
                                        Bez_point_bound::VERTICAL_TANGENCY_PT,
                                        false),
                        bbox));
      return;
    }

    // Check for x-monoticity (if exists return)
    bool is_x_monotone = _is_x_monotone(cp);
    if (is_x_monotone)
      return;

    // Check for y-monoticity (else recurse)
    bool is_y_monotone = _is_y_monotone(cp);
    if (!is_y_monotone)
    {
      // Recurse
      Control_point_vec left, right;

      bisect_control_polygon_2 (cp.begin(), cp.end(),
                                std::back_inserter(left),
                                std::front_inserter(right));
      //DeCasteljau(cp, 0.5, left, right);
      NT t_mid = NT(0.5) * (t0 + t1);
      // TODO - handle the case where t_mid is a vertical (rational) tangency
      // point, by checking whether the first vector of right is vertical.
      // Ron: This happens in the Bezier.dat input file!
      CGAL_assertion(CGAL::compare(right[0].x(), right[1].x()) != EQUAL);

      vertical_tangency_points(left, t0, t_mid, tangency_points);
      vertical_tangency_points(right, t_mid, t1, tangency_points);
      return;
    }

    // Check for convexity (since it is y-monotone and not x-monotone)
    // if convex - add to list, else recurse.
    if (is_convex_2 (cp.begin(), cp.end(), kernel))
    {
      Bez_point_bbox bbox;
      typename Control_point_vec::const_iterator res;
      res = bottom_vertex_2 (cp.begin(), cp.end(), kernel);
      bbox.min_y = res -> y();
      res = top_vertex_2 (cp.begin(), cp.end(), kernel);
      bbox.max_y = res -> y();
      res = left_vertex_2 (cp.begin(), cp.end(), kernel);
      bbox.min_x = res -> x();
      res = right_vertex_2 (cp.begin(), cp.end(), kernel);
      bbox.max_x = res -> x();

      tangency_points.push_back(std::make_pair(
        Bez_point_bound(cp, t0, t1, Bez_point_bound::VERTICAL_TANGENCY_PT, true),
        bbox));
    }
    else
    {
      // Recurse
      Control_point_vec left, right;

      bisect_control_polygon_2 (cp.begin(), cp.end(),
                                std::back_inserter(left),
                                std::front_inserter(right));
      //DeCasteljau(cp, 0.5, left, right);
      NT t_mid = NT(0.5) * (t0 + t1);

      // TODO - handle the case where t_mid is a vertical (rational) tangency
      // point, by checking whether the first vector of right is vertical.
      CGAL_assertion(CGAL::compare(right[0].x(), right[1].x()) != EQUAL);

      vertical_tangency_points(left, t0, t_mid, tangency_points);
      vertical_tangency_points(right, t_mid, t1, tangency_points);
      return;
    }

    // Debug print
    //std::cout << "number of tangency points is " << tangency_points.size() << std::endl;
  }

  // Refinement of vertical tangency point.
  // Assumes the tangency point is isolated (only a single point in the interval).
  void refine_tangency_point(const Control_point_vec& cp, const NT& t0, const NT& t1,
    std::pair<Bez_point_bound, Bez_point_bbox>& tangency_point)
  {
    // TODO - a better implmentation that just checks which of the two subdivided
    // curves is x-monotone (since they should both be y-monotone and convex).

    Control_point_vec left, right;
    /*
    // iddo: the can_refine check is probably not needed here (the check
    // in the recursion will handle it).
    bool can_refine_pt = can_refine(cp, left, right);
    if (!can_refine)
    {
      CGAL_assertion(false); 
    }
    */

    bisect_control_polygon_2 (cp.begin(), cp.end(),
                              std::back_inserter(left),
                              std::front_inserter(right));
    //DeCasteljau(cp, 0.5, left, right);
    NT t_mid = NT(0.5) * (t0 + t1);

    // TODO - handle the case where t_mid is a vertical (rational) tangency point,
    // by checking whether the first vector of right is vertical.
    CGAL_assertion(CGAL::compare(right[0].x(), right[1].x()) != EQUAL);

    bool can_refine_pt_left = can_refine(left, t0, t_mid);
    bool can_refine_pt_right = can_refine(right, t_mid, t1);

    if (!can_refine_pt_left || !can_refine_pt_right)
    {
      // Construct a Bez_point_bound that is inconclusive:
      // all parameters found so far + can_refine = false,
      // and add to output.
      Bez_point_bbox bbox;
      cp_bbox(cp, bbox);
      tangency_point = 
        std::make_pair(Bez_point_bound(cp, t0, t1,
                                       Bez_point_bound::VERTICAL_TANGENCY_PT,
                                       false),
                       bbox);
      return;