hyper_segment_2.h

来自「CGAL is a collaborative effort of severa」· C头文件 代码 · 共 731 行 · 第 1/2 页

   * Compare the slopes of two hyper-segments at their intersection point.
   * \param seg The other segment.
   * \param q A placeholder for the x-coordinate of the intersection point.
   * \return SMALLER if the slope of (*this) is less than seg's at q;
   *         LARGER if it is larger at q;
   *         EQUAL if the two slopes are equal (in case of an overlap).
   * \pre Both (*this) and seg are equal at q.x()
   */
  Comparison_result compare_slopes (const Self& seg,
				    const Point_2& q) const
  {
    // The underlying hyperbola is:
    //   y = sqrt(A*x^2 + B*x + C)
    //
    // If we derive, we obtain:
    //                 2A*x + B            2A*x + B
    //   y' = ------------------------- = ----------
    //         2*sqrt(A*x^2 + B*x + C)        y
    //
    NT                y2 = _A*q.x()*q.x() + _B*q.x() + _C;
    Comparison_result res = _compare (y2, 0);

    CGAL_precondition(_compare (y2,
				seg._A*q.x()*q.x() + seg._B*q.x() + seg._C) 
		      == EQUAL);

    CGAL_assertion(res != SMALLER);
    
    if (res == LARGER)
    {
      // The y-coordinate of the intersection point is greater than 0.
      // It is therefore sufficient to compare the numerators of the 
      // expressions of the two derivatives:
      res = _compare (2*_A*q.x() + _B,
		      2*seg._A*q.x() + seg._B);

      if (res != EQUAL)
	return (res);

      // Compare the second derivatives, where:
      //
      //          2A*y^2 - (2A*x + B)^2
      //   y'' = -----------------------
      //                2*y^3
      //
      return (_compare (2*_A*y2 - (2*_A*q.x() + _B)*(2*_A*q.x() + _B),
			2*seg._A*y2 - (2*seg._A*q.x() + seg._B)*
			(2*seg._A*q.x() + seg._B)));
    }
    else // (res == EQUAL), that is the y-coordinate is 0.
    {
      // Compare the second derivative by x, which in our case equals:
      return (_compare (2*_A*q.x() + _B,
			2*seg._A*q.x() + seg._B));
    }
  }

  /*!
   * Intersect the two hyper-segments.
   * \param seg The second hyper-segment.
   * \param p1 The output intersection point (if one exists).
   * \param p2 The second intersection point (in case of an overlap).
   * \return The number of intersection points computed:
   *         0 - No intersection.
   *         1 - A simple intersection returned as p1.
   *         2 - An overlapping segment, whose end-points are p1 and p2.
   */
  int intersect (const Self& seg,
		 Point_2& p1, Point_2& p2) const
  {
    // First check the case of overlapping hyper-segments.
    if (_has_same_base_hyperbola(seg))
    {
      bool bs = point_is_in_x_range (seg._source); 
      bool bt = point_is_in_x_range (seg._target); 
      
      if (bs && bt)
      {
	// (*this) contains the second segment:
	if (_compare (seg._source.x(), seg._target.x()) == SMALLER)
	{
	  p1 = seg._source;
	  p2 = seg._target;
	}
	else
	{
	  p2 = seg._source;
	  p1 = seg._target;
	}
	return (2);
      }
      else if (bs || bt)
      {
	// Only one of the end-points is contained in (*this):
	const Point_2& p_in  = bs ? seg._source : seg._target;
	const Point_2& p_out = bs ? seg._target : seg._source;
	bool to_right = (_compare(_source.x(), _target.x()) == SMALLER);
	Comparison_result res = _compare (p_out.x(), _source.x());

	if (_compare (p_in.x(), _source.x()) == EQUAL)
	{
	  // The two hyper-segments share just a common end-point.
	  if ((to_right && res == SMALLER) ||
	      (!to_right && res == LARGER))
	  {
	    p1 = p_in;
	    return (1);
	  }
	}
	else if (_compare (p_in.x(), _target.x()) == EQUAL)
	{
	  // The two hyper-segments share just a common end-point.
	  if ((to_right && res == LARGER) ||
	      (!to_right && res == SMALLER))
	  {
	    p1 = p_in;
	    return (1);
	  }
	}

	if (res == SMALLER)
	{
	  p1 = to_right ? _source : _target;
	  p2 = p_in;
	}
	else
	{
	  p1 = p_in;
	  p2 = to_right ? _target : _source;
	}
	return (2);
      }
      else
      {
	// Check if seg contains (*this).
	bs = seg.point_is_in_x_range (_source); 
	bt = seg.point_is_in_x_range (_target); 

	if (bs && bt)
	{
	  if (_compare (_source.x(), _target.x()) == SMALLER)
	  {
	    p1 = _source;
	    p2 = _target;
	  }
	  else
	  {
	    p2 = _source;
	    p1 = _target;
	  }
	  return (2); 
	}
      }

      // The two segments are disjoint:
      return (0);
    }

    // The equations of the two underlying hyperbola are given by:
    //   H1: y = A1*x^2 + B1*x + C1
    //   H2: y = A2*x^2 + B2*x + C2
    //
    // So the following must hold for the x-coordinates of the intersection
    // points:
    //   (A1 - A2)*x^2 + (B1 - B2)*x + (C1 - C2) = 0
    //
    NT     xs[2];
    bool   bs[2];
    int    n_xs;
    int    i;

    n_xs = _solve_quadratic_equation (_A - seg._A,
				      _B - seg._B,
				      _C - seg._C,
				      xs[0], xs[1]);

    // Check the legality of each of the x-coordinates.
    bs[0] = false;
    bs[1] = false;

    for (i = 0; i < n_xs; i++)
    {
      bs[i] = true;

      Comparison_result res1 = _compare(xs[i], _source.x());
      Comparison_result res2 = _compare(xs[i], _target.x());

      if (res1 != EQUAL && res2 != EQUAL && res1 == res2)
	bs[i] = false;
      else
      {
	res1 = _compare(xs[i], seg._source.x());
	res2 = _compare(xs[i], seg._target.x());
      
	if (res1 != EQUAL && res2 != EQUAL && res1 == res2)
	  bs[i] = false;
      }
    }

    // Make sure we return one point at most.
    if (!bs[0] && !bs[1])
      return (0);

    CGAL_assertion (!(bs[0] && bs[1]));

    if (bs[0])
      _get_point_at_x (xs[0], p1);
    else
      _get_point_at_x (xs[1], p1);

    return (1);
  }

  /*!
   * Split a hyper-segment at a given point.
   * \param p The split point.
   * \param seg1 The first resulting hyper-segment.
   * \param seg2 The second resulting hyper-segment.
   * \pre p lies in the interior of the curve (and is not and end-point).
   */
  void split (const Point_2& p,
	      Self& seg1, Self& seg2) const
  {
    CGAL_precondition(point_position(p) == EQUAL);
    CGAL_precondition(_compare(_source.x(), p.x()) != EQUAL);
    CGAL_precondition(_compare(_target.x(), p.x()) != EQUAL);

    seg1 = (*this);
    seg1._target = p;

    seg2 = (*this);
    seg2._source = p;

    return;
  }

  /*
   * Trim the hyper-segment at the given x-coodinates.
   * \param p1 A placeholder for the x-coordinates of the new source.
   * \param p2 A placeholder for the x-coordinates of the new target.
   * \return The trimmed hyper-segment.
   * \pre p1 and p2 are both in the x-range of the hyper-segment.
   */
  Self trim (const Point_2& p1,
	     const Point_2& p2) const
  {
    CGAL_precondition(point_is_in_x_range(p1));
    CGAL_precondition(point_is_in_x_range(p2));

    Self     trimmed (*this);

    _get_point_at_x (p1.x(), trimmed._source);
    _get_point_at_x (p2.x(), trimmed._target);

    return (trimmed);
  }

protected:

  /*!
   * Compare two values.
   */
  inline Comparison_result _compare (const NT& val1, const NT& val2) const
  {
    return (CGAL_NTS compare(val1, val2));
  }

  /*!
   * Check whether the underlying hyperbola contains a point with the given
   * x-coordinate and compute it if it does.
   * \param x The x-coordinate.
   * \param p The output point.
   * \return (true) if the underlying hyperbola contains a point with the 
   *         given x-coordinate.
   */
  bool _get_point_at_x (const NT& x,
                        Point_2& p) const
  {
    // Try to compute the y-coordinate.
    NT              val = _A*x*x + _B*x + _C;
    static const NT _zero = 0;

    if (_compare(val, _zero) == SMALLER)
      return (false);

    p = Point_2(x, CGAL::sqrt(val));
    return (true);
  }

  /*!
   * Check whether the two hyper-segments lie on the same hyperbola.
   * \param seg The compared hyper-segment.
   * \return (true) if the two hyper-segments lie on the same hyperbola.
   */
  bool _has_same_base_hyperbola (const Self& seg) const
  {
    return ((_compare (_A, seg._A) == EQUAL) &&
	    (_compare (_B, seg._B) == EQUAL) &&
	    (_compare (_C, seg._C) == EQUAL));
  }

  /*!
   * Solve the equation a*x^2 + b*x + c = 0.
   * \return The number of solutions found.
   */
  int _solve_quadratic_equation (const NT& a, const NT& b, const NT& c,
				 NT& x1, NT& x2) const
  {
    static const NT _zero = 0;
    static const NT _two = 2;
    static const NT _four = 4;

    if (_compare (a, _zero) == EQUAL)
    {
      // If this is a linear equation:
      if (_compare (b, _zero) == EQUAL)
	return (0);

      x1 = -c / b;
      return (1);
    }

    // Act according to the sign of (b^2 - 4ac):
    const NT          disc = b*b - _four*a*c;
    Comparison_result res = _compare (disc, _zero);

    if (res == SMALLER)
    {
      return (0);
    }
    else if (res == EQUAL)
    {
      x1 = (-b) / (_two*a);
      return (1);
    }
    
    const NT          sqrt_disc = CGAL::sqrt (disc);
    
    x1 = (sqrt_disc - b) / (_two*a);
    x2 = -(sqrt_disc + b) / (_two*a);
    return (2);
  }
};

#ifndef NO_OSTREAM_INSERT_CONIC_ARC_2
template <class Kernel>
std::ostream& operator<< (std::ostream& os, 
			  const Hyper_segment_2<Kernel>& seg)
{
  os << "{ y^2 = " << CGAL::to_double(seg.A()) << "*x^2 + "
     << CGAL::to_double(seg.B()) << "*x + " 
     << CGAL::to_double(seg.C()) << " } : "
     << "(" << CGAL::to_double(seg.source().x()) << "," 
     << CGAL::to_double(seg.source().y()) << ") -> "
     << "(" << CGAL::to_double(seg.target().x()) << "," 
     << CGAL::to_double(seg.target().y()) << ")";

  return (os);
}
#endif // NO_OSTREAM_INSERT_CONIC_ARC_2

CGAL_END_NAMESPACE

#endif