bezier_bounding_rational_traits.h

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

    }

    std::list<std::pair<Bez_point_bound, Bez_point_bbox> > aux1, aux2;
    vertical_tangency_points(left, t0, t_mid, aux1);
    vertical_tangency_points(right, t_mid, t1, aux2);

    CGAL_assertion(aux1.size() + aux2.size() == 1);
    if (aux1.size() == 1)
    {
      tangency_point = *(aux1.begin());
    }
    else 
    {
      tangency_point = *(aux2.begin());
    }
  }

  ////////////////////////////////////////////////////////////////
  // Slope comparison.
  ///////////////////////////////////////////////////////////////
  Comparison_result compare_slopes_of_bounded_intersection_point(
    const Bez_point_bound& bound1, const Bez_point_bound& bound2)
  {

    // Since we assume the spans do not overlap, we can just compare any vector
    // of the hodograph spans.
    const Control_point_vec& cv1 = bound1.bounding_polyline;
    const Control_point_vec& cv2 = bound2.bounding_polyline;

    // An (expensive) check for the precondition that the spans do not overlap. 
    // Therefor it is commented out.
    /*
    Vector_2 dir1, dir2, leftMost1, leftMost2, rightMost1, rightMost2;
    bool spansOverlap;
    bool cv1SpanOK = _CrvTanAngularSpan(cv1, dir1, leftMost1, rightMost1);
	  bool cv2SpanOK = _CrvTanAngularSpan(cv2, dir2, leftMost2, rightMost2);
    if (cv1SpanOK && cv2SpanOK)
    {
  	  spansOverlap = _AngularSpansOverlap(leftMost1, rightMost1, leftMost2, rightMost2);
    }
    else
    {
      spansOverlap = true;
    }
    CGAL_precondition(!spansOverlap);
    */

    Vector_2 dir1 = cv1.back() - cv1.front();
    Vector_2 dir2 = cv2.back() - cv2.front();

    // RI: use compare_angle_with_x_axis (after doing: if (dir.x()<0) dir=-dir)
    // instead of division.
    NT slope1 = dir1.y()/dir1.x();
    NT slope2 = dir2.y()/dir2.x();
    return CGAL::compare(slope1, slope2);
  }


  // Can be made a template function <Control_point_vec, Kernel>
  void cp_bbox(const Control_point_vec& cp, Bez_point_bbox& bez_bbox) const
  {
   CGAL_assertion (! cp.empty());

    typename Kernel::Compare_x_2       comp_x = kernel.compare_x_2_object();
    typename Kernel::Compare_y_2       comp_y = kernel.compare_y_2_object();

    typename Control_point_vec::const_iterator it = cp.begin();
    typename Control_point_vec::const_iterator min_x = it;
    typename Control_point_vec::const_iterator max_x = it;
    typename Control_point_vec::const_iterator min_y = it;
    typename Control_point_vec::const_iterator max_y = it;

    while (it != cp.end())
    {
      if (comp_x (*it, *min_x) == SMALLER)
        min_x = it;
      else if (comp_x (*it, *max_x) == LARGER)
        max_x = it;

      if (comp_y (*it, *min_y) == SMALLER)
        min_y = it;
      else if (comp_y (*it, *max_y) == LARGER)
        max_y = it;

      ++it;
    }

    bez_bbox.min_x = min_x->x(); 
    bez_bbox.max_x = max_x->x(); 
    bez_bbox.min_y = min_y->y(); 
    bez_bbox.max_y = max_y->y(); 
  }

private:

  /*
  void _construct_bboxes (const Control_point_vec& cv1,
                          const Control_point_vec& cv2,
                          Bez_point_bbox& bbox1, Bez_point_bbox& bbox2) const
  {
    const double   diff_bound = 0.1875; // radians, a bit more than 10 degrees.
    const double   _pi = 3.14159265;

    // Compute the angles theta1 and theta2 the lines that connect the first
    // and last control points in each curve form with the x-axis.
    const double   dx1 = CGAL::to_double (cv1.back().x()) -
                         CGAL::to_double (cv1.front().x());
    const double   dy1 = CGAL::to_double (cv1.back().y()) -
                         CGAL::to_double (cv1.front().y());
    const double   theta1 = atan2 (dy1, dx1);

    const double   dx2 = CGAL::to_double (cv2.back().x()) -
                         CGAL::to_double (cv2.front().x());
    const double   dy2 = CGAL::to_double (cv2.back().y()) -
                         CGAL::to_double (cv2.front().y());
    const double   theta2 = atan2 (dy2, dx2);

    // Check if the two angles are close, and compute half their average.
    const double   diff_theta = fabs (theta2 - theta1);

    double         half_theta;

    if (diff_theta < diff_bound)
    {
      // The two angles are close:
      half_theta = (theta1 + theta2) / 4;
    }
    else if (fabs (diff_theta - _pi) < diff_bound)
    {
      // We have theta2 ~= theta1 +/- PI, thus we take:
      if (theta2 > 0)
        half_theta = (theta1 + theta2 - _pi) / 4;
      else
        half_theta = (theta1 + theta2 + _pi) / 4;
    }
    else
    {
      // In this case the angles are not close, so we compute their
      // axis-alligned bounding boxes.
      cp_bbox(cv1, bbox1);
      cp_bbox(cv2, bbox2);
      return;
    }

    // Check if theta/2 is close to 0 or to +/- PI/2.
    if (fabs(half_theta) < 0.0625 || _pi/2 - fabs(half_theta) < 0.0625)
    {
      // In this case, using the axis-alligned bounding box is fine.
      cp_bbox(cv1, bbox1);
      cp_bbox(cv2, bbox2);
      return;
    }

    // Compute a rational tau ~= tan(theta/2) and use it to find a rotation
    // angle close to theta with rational sine and cosine.
    const double   tan_half_theta = tan (half_theta);
    const int      denom = 1000;
    const int      numer = static_cast<int> (denom * tan_half_theta + 0.5);
    const NT       tau (numer, denom);
    const NT       sqr_tau = tau*tau;
    const NT       sin_theta = 2*tau / (1 + sqr_tau);
    const NT       cos_theta = (1 - sqr_tau) / (1 + sqr_tau);

    // Rotate both control polygons by the minus the approximated angle. 
    Control_point_vec                          rot_cv1;
    Control_point_vec                          rot_cv2;
    typename Control_point_vec::const_iterator it;

    for (it = cv1.begin(); it != cv1.end(); ++it)
    {
      rot_cv1.push_back (Point_2 (cos_theta * it->x() + sin_theta * it->y(),
                                  -sin_theta * it->x() + cos_theta * it->y()));
    }

    for (it = cv2.begin(); it != cv2.end(); ++it)
    {
      rot_cv2.push_back (Point_2 (cos_theta * it->x() + sin_theta * it->y(),
                                  -sin_theta * it->x() + cos_theta * it->y()));
    }

    // Compute the axis-alligned bounding boxes of the rotated curves, which
    // are equivalent of the rotated skewed bounding boxes of the curves.
    cp_bbox (rot_cv1, bbox1);
    cp_bbox (rot_cv2, bbox2);

    return;
  }
  */

  // Can be made a template function <Control_point_vec, Kernel>
  bool _is_monotone(const Control_point_vec& cp, bool check_x) const
  {
    CGAL_precondition(cp.size() > 1);
    Comparison_result                          curr_res;
    Comparison_result                          res = EQUAL;

    typename Control_point_vec::const_iterator cp_iter = cp.begin();
    typename Control_point_vec::const_iterator cp_iter_end = cp.end();
    typename Control_point_vec::const_iterator cp_iter_next = cp_iter; 
    ++cp_iter_next;

    while (cp_iter_next != cp_iter_end)
    {
      curr_res = check_x ? compare_x(*cp_iter, *cp_iter_next) :
                           compare_y(*cp_iter, *cp_iter_next);
      if (curr_res != EQUAL)
      {
        res = curr_res;
        ++cp_iter; 
        ++cp_iter_next;
        break;
      }
      ++cp_iter;
      ++cp_iter_next;
    }
    if (cp_iter_next == cp_iter_end)
      // iso-parallel segment.
      // In our context, we consider a horizontal segment as non-y-monotone
      // but a vertical segment is considered x-monotone (is this right?).
    {
      if (check_x)
        return true; // Vertical segment when checking x-monoticity
      return false;  // Horizontal segment when checking y-monoticity
    }

    while (cp_iter_next != cp_iter_end)
    {
      curr_res = check_x ? compare_x(*cp_iter, *cp_iter_next) :
                           compare_y(*cp_iter, *cp_iter_next);
      if (curr_res == EQUAL)
      {
        ++cp_iter;
        ++cp_iter_next;
        continue;
      }
      if (res != curr_res)
        break;

      ++cp_iter;
      ++cp_iter_next;
    }
    if (cp_iter_next != cp_iter_end)
      return false;
    return true;
  }

  inline bool _is_x_monotone(const Control_point_vec& cp) const
  {
    return _is_monotone(cp, true);
  }

  inline bool _is_y_monotone(const Control_point_vec& cp) const
  {
    return _is_monotone(cp, false);
  }


  // Can be made a template function <Control_point_vec, Kernel>
  // Returns true if the span is less than 90 degrees (we don't want to work with larger spans).
  bool _CrvTanAngularSpan(const Control_point_vec& cp,
				Vector_2& dir,
				Vector_2& leftMost,
				Vector_2& rightMost)
  {
    Vector_2 v;
    const Point_2& O(ORIGIN);

    dir = cp.back() - cp.front();

	  // Initialize LeftMost and RightMost
	  leftMost = dir;
	  rightMost = dir;

    typename Control_point_vec::const_iterator iterCP0 = cp.begin();
    typename Control_point_vec::const_iterator iterCP1 = iterCP0;
    ++iterCP1;
    for (; iterCP1!=cp.end(); ++iterCP1, ++iterCP0) 
    {
      v = (*iterCP1) - (*iterCP0);

      // Test for 90 degree angle, we can loosen it to test for 180 degree between
	    // left and right but we don't really want to work with larger spans.
      // TODO - try get rid of these ORIGINs
      if (angle(ORIGIN+v, O, ORIGIN+dir) != ACUTE)
	      return false;

	    if (left_turn(ORIGIN+v, O, ORIGIN+dir)) 
      {
		    if (left_turn(ORIGIN+v, O, ORIGIN+leftMost)) 
        {
			    leftMost = v;
		    }
	    }
	    else 
      {
		    if (right_turn(ORIGIN+v, O, ORIGIN+rightMost)) 
        {
			    rightMost = v;		
		    }
	    }
    } //for

    return true;
  }

  // TODO - there is already a kernel function that does this ccw?
  // Can be made a template function <Kernel>
  inline bool _VecIsBetween(const Vector_2& vec,
							      const Vector_2& leftVec,
							      const Vector_2& rightVec)
  {
    const Point_2& O(ORIGIN);
	  return 	(right_turn(ORIGIN+rightVec, O, ORIGIN+vec) && left_turn(ORIGIN+leftVec, O, ORIGIN+vec));
  }

  // Can be made a template function <Kernel>
  bool _AngularSpansOverlap(
				  const Vector_2& leftMost1,
				  const Vector_2& rightMost1,
				  const Vector_2& leftMost2,
				  const Vector_2& rightMost2
				  )
  {
	  Vector_2 negVec;

	  // Check if any of span2 vectors is between span1.
	  if (_VecIsBetween(leftMost2, leftMost1, rightMost1))
		  return true;

	  negVec = -leftMost2;
	  if (_VecIsBetween(negVec, leftMost1, rightMost1))
		  return true;

	  if (_VecIsBetween(rightMost2, leftMost1, rightMost1))
		  return true;

	  negVec = -rightMost2;
	  if (_VecIsBetween(negVec, leftMost1, rightMost1))
		  return true;


	  // Note: If we got here, the only possibility for intersection
	  // is that span1 is totally between span2

	  // Check if span1 is between span2
	  if (_VecIsBetween(leftMost1, leftMost2, rightMost2))
		  return true;

	  negVec = -leftMost1;
	  if (_VecIsBetween(negVec, leftMost2, rightMost2))
		  return true;


	  // The rest of the symmetric conditions below are redundant (see note above)
	  return false;

    /*
	  if (_VecIsBetween(rightMost1, leftMost2, rightMost2))
		  return true;

	  negVec = -rightMost1;
	  if (_VecIsBetween(negVec, leftMost2, rightMost2))
		  return true;

    return false;
    */			
  }

  // Constructing the curve skewed bbox - actually returning
  // two parallel lines (parallel to the first-last line) that bound
  // the curve from above and below.
  // Using less_signed_distance_to_line_2 enables not to construct
  // new points (only take maximum and minimum signed dist points
  // and add the first-last vector).
  // Can be made a template function <Control_point_vec, Kernel>
  void _cvSkewBbox(const Control_point_vec& cp,
				Line_2& la,
				Line_2& lb)
  {
    // Functor for comparing two points signed distance from a given line (used for min/max_elem).
    // Binder taken from <CGAL/functional.h>, if this doesn't compile, construct your own functor.
    typedef CGALi::Binder<typename Kernel::Less_signed_distance_to_line_2, Arity_tag< 3 >, Line_2, 1 > Less_signed_dist_to_line;
    // Explanation: Take the 3-ary kernel functor Less_signed_distance_to_line_2 and bind its
    // first parameter, the line, to be our given line l.

    Line_2 l(cp.front(), cp.back());
    Vector_2 v(cp.front(), cp.back());
    typename Kernel::Less_signed_distance_to_line_2 lsdtl2 = kernel.less_signed_distance_to_line_2_object();

    typename Control_point_vec::const_iterator aux;
    aux = std::min_element(cp.begin(), cp.end(), Less_signed_dist_to_line(lsdtl2, l));
    la = Line_2(*aux, v);

    aux = std::max_element(cp.begin(), cp.end(), Less_signed_dist_to_line(lsdtl2, l));
    lb = Line_2(*aux, v);