conic_arc_2_core.h

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

    //        t*v - 2*s*u                t*u - 2*r*v
    //  x0 = -------------   ,     y0 = -------------
    //        4*r*s - t^2                4*r*s - t^2
    //
    // The denominator (4*r*s - t^2) must be negative for hyperbolas.
    const CoNT  u = or_fact * CoNT(_u);
    const CoNT  v = or_fact * CoNT(_v);
    const CoNT  det = 4*r*s - t*t;
    CoNT        x0, y0;

    CGAL_assertion (det < _zero);
    
    x0 = (t*v - _two*s*u) / det;
    y0 = (t*u - _two*r*v) / det;
    
    // The axis separating the two branches of the hyperbola is now given by:
    // 
    //  cos(phi)*x + sin(phi)*y - (cos(phi)*x0 + sin(phi)*y0) = 0
    //
    _data.hyper_P = new Hyperbolic_arc_data;

    _data.hyper_P->a = cos_phi;
    _data.hyper_P->b = sin_phi;
    _data.hyper_P->c = - (cos_phi*x0 + sin_phi*y0);

    // Make sure that the two endpoints are located on the same branch
    // of the hyperbola.
    _data.hyper_P->side = _hyperbolic_arc_side(_source);

    CGAL_assertion (_data.hyper_P->side = _hyperbolic_arc_side(_target));

    return;
  }

  /*!
   * Find on which branch of the hyperbola is the given point located.
   * The point is assumed to be on the hyperbola.
   * \param p The query point.
   * \return The branch ID (either -1 or 1).
   */
  int _hyperbolic_arc_side (const Point_2& p) const
  {
    if (_data.hyper_P == NULL)
      return (0);

    CoNT       val;

    val = _data.hyper_P->a*p.x() + _data.hyper_P->b*p.y() + _data.hyper_P->c;
    return ((val > 0) ? 1 : -1);
  }
 
  /*!
   * Check whether the given point is between the source and the target.
   * The point is assumed to be on the conic's boundary.
   * \param p The query point.
   * \return (true) if the point is between the two endpoints, 
   *         (false) if it is not.
   */
  bool _is_between_endpoints (const Point_2& p) const
  {
    if (p.equals(_source) || p.equals(_target))
      return (true);
    else
      return (_is_strictly_between_endpoints(p));
  }

  /*!
   * Check whether the given point is strictly between the source and the
   * target (but not any of them).
   * The point is assumed to be on the conic's boundary.
   * \param p The query point.
   * \return (true) if the point is strictly between the two endpoints, 
   *         (false) if it is not.
   */
  bool _is_strictly_between_endpoints (const Point_2& p) const
  {
    // In case this is a full conic, any point on its boundary is between
    // its end points.
    if (is_full_conic())
      return (true);

    // In case of a hyperbolic arc, make sure the point is located on the
    // same branch as the arc.
    if ((_info & IS_HYPERBOLA) != 0)
    {
      if (_hyperbolic_arc_side(p) != _data.hyper_P->side)
	return (false);
    }

    // Act according to the conic degree.
    if ((_info & DEGREE_MASK) == DEGREE_1)
    {
      if (is_vertical_segment())
      {
	// In case of a vertical segment - just check whether the y coordinate
	// of p is between those of the source's and of the target's.
	Comparison_result r1 = compare_y (p, _source);
	Comparison_result r2 = compare_y (p, _target);

	return ((r1 == SMALLER && r2 == LARGER) ||
		(r1 == LARGER && r2 == SMALLER));
      }
      else
      {
	// Otherwise, since the segment is x-monotone, just check whether the
	// x coordinate of p is between those of the source's and of the 
	// target's.
	Comparison_result r1 = compare_x (p, _source);
	Comparison_result r2 = compare_x (p, _target);

	return ((r1 == SMALLER && r2 == LARGER) ||
		(r1 == LARGER && r2 == SMALLER));
      }
    }
    else
    {
      // In case of a conic of degree 2, make a decision based on the conic's
      // orientation and whether (source,p,target) is a right or a left turn.
      static Kernel                  ker;
      typename Kernel::Orientation_2 orient_f = ker.orientation_2_object();

      if (_orient == CGAL::COUNTERCLOCKWISE)
	return (orient_f(_source, p, _target) == LEFT_TURN);
      else
	return (orient_f(_source, p, _target) == RIGHT_TURN);
    }
  }

  /*!
   * Check whether the underlying conic contains a point on its boundary.
   * \param q The query point.
   * \return (true) if the underlying conic contains the point on its boundary.
   */
  bool _conic_has_on_boundary (const Point_2& q) const
  {
    return (_conic_has_on_boundary (q.x(), q.y()));
  }

  /*!
   * Check whether the underlying conic contains (x,y) on its boundary.
   * \param x The x coordinate of the query point.
   * \param y The y coordinate of the query point.
   * \return (true) if the underlying conic contains the point on its boundary.
   */
  bool _conic_has_on_boundary (const CoNT& x, const CoNT& y) const
  {
    const CoNT _zero = 0;
    CoNT       val;

    // The point must satisfy: r*x^2 + s*y^2 + t*xy + u*x + v*y + w = 0.
    val = CoNT(_r)*x*x + CoNT(_s)*y*y + CoNT(_t)*x*y +
          CoNT(_u)*x + CoNT(_v)*y + CoNT(_w);

    return (val == _zero);
  }

  /*!
   * Find the y coordinates of the underlying conic at a given x coordinate.
   * \param x The x coordinate.
   * \param ys The output y coordinates. 
   *           This area must be allocated at the size of 2.
   * \return The number of y coordinates computed (either 0, 1 or 2).
   */
  int _conic_get_y_coordinates (const CoNT& x,
                                CoNT *ys) const
  {
    // Solve the quadratic equation for a given x and find the y values:
    //  s*y^2 + (t*x + v)*y + (r*x^2 + u*x + w) = 0
    return (solve_quadratic_eq<CoNT,CoNT> 
	    (CoNT(_s),
	     CoNT(_t)*x + CoNT(_v),
	     (CoNT(_r)*x + CoNT(_u))*x + CoNT(_w),
	     ys));
  }

  /*!
   * Find the x coordinates of the underlying conic at a given y coordinate.
   * \param y The y coordinate.
   * \param xs The output x coordinates. 
   *           This area must be allocated at the size of 2.
   * \return The number of x coordinates computed (either 0, 1 or 2).
   */
  int _conic_get_x_coordinates (const CoNT& y,
                                CoNT *xs) const
  {
    // Solve the quadratic equation for a given y and find the x values:
    //  r*x^2 + (t*y + u)*x + (s*y^2 + v*y + w) = 0
    return (solve_quadratic_eq<CoNT,CoNT> 
	    (CoNT(_r),
	     CoNT(_t)*y + CoNT(_u),
	     (CoNT(_s)*y + CoNT(_v))*y + CoNT(_w),
	     xs));
  }
  
  /*!
   * Find the vertical tangency points of the undelying conic.
   * \param ps The output points of vertical tangency.
   *           This area must be allocated at the size of 2.
   * \return The number of vertical tangency points.
   */
  int _conic_vertical_tangency_points (Point_2* ps) const
  {
    const CfNT _zero = 0;

    // In case the base conic is of degree 1 (and not 2), the arc has no
    // vertical tangency points.
    if ((_info & DEGREE_MASK) == DEGREE_1 || _s == _zero)
      return (0);

    // Special treatment for circles, where the vertical tangency points
    // are simply (x0-r,y0) and (x0+r,y0).
    if ((_info & IS_CIRCLE) != 0)
    {
      ps[0] = Point_2 (CoNT(_data.circ_P->x0 - _data.circ_P->r), 
		       CoNT(_data.circ_P->y0),
		       _conic_id);

      ps[1] = Point_2 (CoNT(_data.circ_P->x0 + _data.circ_P->r), 
                       CoNT(_data.circ_P->y0),
		       _conic_id);

      return (2);
    }

    // We are interested in the x coordinates where the quadratic equation:
    //  s*y^2 + (t*x + v)*y + (r*x^2 + u*x + w) = 0
    // has a single solution (obviously if s = 0, there are no such points).
    // We therefore demand that the discriminant of this equation is zero:
    //  (t*x + v)^2 - 4*s*(r*x^2 + u*x + w) = 0
    const CfNT _two = 2;
    const CfNT _four = 4;
    CoNT       xs[2];
    int        n_xs;

    n_xs = solve_quadratic_eq<CfNT,CoNT> (_t*_t - _four*_r*_s,
					  _two*_t*_v - _four*_s*_u,
					  _v*_v - _four*_s*_w,
					  xs);

    // Find the y-coordinates of the vertical tangency points.
    CoNT     ys[2];
    int      n_ys;

    if (_t == _zero)
    {
      // The two vertical tangency points have the same y coordinate:
      ys[0] = CoNT(-_v) / CoNT(_two*_s);
      n_ys = 1;
    }
    else
    {
      n_ys = solve_quadratic_eq<CfNT,CoNT> (_four*_r*_s*_s - _s*_t*_t,
					    _four*_r*_s*_v - _two*_s*_t*_u,
					    _r*_v*_v - _t*_u*_v + _t*_t*_w,
					    ys);
    }

    // Pair the x and y coordinates and obtain the vertical tangency points.
    int   n = 0;
    int   i, j;

    for (i = 0; i < n_xs; i++)
    {
      if (n_ys == 1)
      {
        ps[n] = Point_2 (xs[i], ys[0],
			 _conic_id);
	n++;
      }
      else
      {
	for (j = 0; j < n_ys; j++)
	{
	  if (ys[j] == -(CoNT(_t)*xs[i] + CoNT(_v)) / CoNT(_two*_s))
	  {
	    ps[n] = Point_2 (xs[i], ys[j],
			     _conic_id);
	    n++;
	    break;
	  }
	}
      }
    }

    return (n);
  }

  /*!
   * Find the horizontal tangency points of the undelying conic.
   * \param ps The output points of horizontal tangency.
   *           This area must be allocated at the size of 2.
   * \return The number of horizontal tangency points.
   */
  int _conic_horizontal_tangency_points (Point_2* ps) const
  {
    const CfNT _zero = 0;

    // In case the base conic is of degree 1 (and not 2), the arc has no
    // vertical tangency points.
    if ((_info & DEGREE_MASK) == DEGREE_1 || _r == _zero)
      return (0);

    // Special treatment for circles, where the horizontal tangency points
    // are simply (x0,y0-r) and (x0,y0+r).
    if ((_info & IS_CIRCLE) != 0)
    {
      ps[0] = Point_2 (CoNT(_data.circ_P->x0), 
                       CoNT(_data.circ_P->y0 - _data.circ_P->r),
		       _conic_id);

      ps[1] = Point_2 (CoNT(_data.circ_P->x0), 
                       CoNT(_data.circ_P->y0 + _data.circ_P->r),
		       _conic_id);

      return (2);
    }

    // We are interested in the y coordinates were the quadratic equation:
    //  r*x^2 + (t*y + u)*x + (s*y^2 + v*y + w) = 0
    // has a single solution (obviously if r = 0, there are no such points).
    // We therefore demand that the discriminant of this equation is zero:
    //  (t*y + u)^2 - 4*r*(s*y^2 + v*y + w) = 0
    const CfNT _two = 2;
    const CfNT _four = 4;
    int        n;
    CoNT       ys[2];

    n = solve_quadratic_eq<CfNT,CoNT> (_t*_t - _four*_r*_s,
				       _two*_t*_u - _four*_r*_v,
				       _u*_u - _four*_r*_w,
				       ys);

    // Compute the x coordinates and construct the horizontal tangency points.
    CoNT       x;
    int        i;

    for (i = 0; i < n; i++)
    {
      // Having computed y, x is the simgle solution to the quadratic equation
      // above, and since its discriminant is 0, x is simply given by:
      x = -(CoNT(_t)*ys[i] + CoNT(_u)) / CoNT(_two*_r);

      ps[i] = Point_2 (x, ys[i],
		       _conic_id);
    }
      
    return (n);
  }
  
  /*!
   * Set the facing information for the (x-monotone) arc: It is facing up if
   * it lies above the line segments that connect its two endpoints, and facing
   * down if it lies below it.
   */
  void _set_facing ()
  {
    // Check whether the arc (which is x-monotone of degree 2) lies above or 
    // below the segement that contects its two end-points (x1,y1) and (x2,y2).
    // To do that, we find the y coordinate of a point on the arc whose x
    // coordinate is (x1+x2)/2 and compare it to (y1+y2)/2.
    const CoNT   _two = 2;
    const CoNT   x_mid = (_source.x() + _target.x()) / _two;
    const CoNT   y_mid = (_source.y() + _target.y()) / _two;
    Point_2      p_mid (x_mid, y_mid);
    Point_2      ps[2];
    int          n_ps;

    n_ps = get_points_at_x (p_mid, ps);

    CGAL_assertion (n_ps == 1);

    Comparison_result res = ps[0].compare_y (p_mid);

    if (res == LARGER)
    {
      // The arc is above the connecting segment, so it is facing upwards.
      _info = _info | FACING_UP;
    }
    else if (res == SMALLER)
    {
      // The arc is below the connecting segment, so it is facing downwards.
      _info = _info | FACING_DOWN;
    }
    
    CGAL_assertion(res != EQUAL);
    return;
  }

  /*!
   * Calculate all x coordinates of intersection points between the two
   * base conics of (*this) and the given arc.
   * \param arc The arc whose underlying conic we intersect.
   * \param xs The output x coordinates.
   *           This area must be allocated to the size of 4.
   * \return The number of unique x coordinates.
   */
  int _x_coordinates_of_intersections_with (const Self& arc,
					    CoNT* xs) const
  {
    const CfNT _zero = 0;
    int        n_roots;         // The number of distinct x values.

    // Check whether both arcs are line segments.
    if ((_info & DEGREE_MASK) == DEGREE_1 && 
	(arc._info & DEGREE_MASK) == DEGREE_1)
    {
      // The two conics are: u*x + v*y + w = 0
      //                and: u'*x + v'*y + w' = 0
      // There's a single solution for x, which is:
      const CfNT denom = _v*arc._w - _w*arc._v;
      const CfNT numer = _u*arc._v - _v*arc._u;

      if (numer == _zero)
      {
	n_roots = 0;
      }
      else
      {
	xs[0] = CoNT(denom) / CoNT(numer);
	n_roots = 1;
      }
    }
    // Check whether the second arc is really a line segment.
    else if ((arc._info & DEGREE_MASK) == DEGREE_1)
    {
      // The two conics are: r*x^2 + s*y^2 + t*xy + u*x + v*y + w = 0
      //                and: a*x + b*y + c = 0
      // There are therefore 2 possible x-values, the solutions for:
      const CfNT  _two = 2;
      const CfNT& a = arc._u;
      const CfNT& b = arc._v;
      const CfNT& c = arc._w;

      n_roots = solve_quadratic_eq<CfNT, CoNT> 
	(a*a*_s + b*b*_r - a*b*_t,
	 _two*a*c*_s + b*b*_u - a*b*_v - b*c*_t,
	 c*c*_s + b*b*_w - b*c*_v,
	 xs);
    }
    // Check if the two arcs are circular arcs.
    else if ((_info & IS_CIRCLE) != 0 && (arc._info & IS_CIRCLE) != 0)
    {
      // Special treatment for two circles.
      // The two curves are: r*x^2 + r*y^2 + u*x + v*y + w = 0
      //                and: r'*x^2 + r'*y^2 + u'*x + v'*y + w' = 0
      //
      // Thus, r'*C1-r*C2 is a line whose equation is: a*x + b*y + = 0, where:
      const CfNT  _two = 2;
      const CfNT a = arc._r*_u - _r*arc._u;
      const CfNT b = arc._r*_v - _r*arc._v;
      const CfNT c = arc._r*_w - _r*arc._w;

      if (b == _zero)
      {
	// The line a*x + c = 0 connects both intersection points of the two
	// circles, so the both have an x-coordinate of -c/a.
	if