conic_arc_2.h

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

    else
    {
      // sin(2*phi) < 0 so PI/2 < phi < PI.
      sin_phi = nt_traits.sqrt((_one - cos_2phi) / _two);
      cos_phi = -nt_traits.sqrt((_one + cos_2phi) / _two);
    }
           
    // Calculate the center (x0, y0) of the conic, given by the formulae:
    //
    //        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 Algebraic  u = nt_traits.convert (or_fact * _u);
    const Algebraic  v = nt_traits.convert (or_fact * _v);
    const Algebraic  det = 4*r*s - t*t;
    Algebraic        x0, y0;

    CGAL_assertion (CGAL::sign (det) == NEGATIVE);
    
    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
    //
    // We store the equation of this line in the extra data structure and also
    // the sign (side of half-plane) our arc occupies with respect to the line.
    _extra_data_P = new Extra_data;

    _extra_data_P->a = cos_phi;
    _extra_data_P->b = sin_phi;
    _extra_data_P->c = - (cos_phi*x0 + sin_phi*y0);

    // Make sure that the two endpoints are located on the same branch
    // of the hyperbola.
    _extra_data_P->side = _sign_of_extra_data (_source.x(), _source.y());

    CGAL_assertion (_extra_data_P->side != ZERO);
    CGAL_assertion (_extra_data_P->side = _sign_of_extra_data(_target.x(),
							      _target.y()));

    return;
  }
  //@}

protected:

  /// \name Auxiliary functions.
  //@{

  /*!
   * Evaluate the sign of (a*x + b*y + c) stored with the extra data field
   * at a given point.
   * \param px The x-coordinate of query point.
   * \param py The y-coordinate of query point.
   * \return The sign of (a*x + b*y + c).
   */
  Sign _sign_of_extra_data (const Algebraic& px,
			    const Algebraic& py) const
  {
    CGAL_assertion (_extra_data_P != NULL);

    if (_extra_data_P == NULL)
      return (ZERO);

    Algebraic         val = (_extra_data_P->a*px + _extra_data_P->b*py + 
                             _extra_data_P->c);

    return (CGAL::sign (val));
  }

  /*!
   * Check whether the given point lies on the supporting conic of the arc.
   * \param p The query point.
   * \return (true) if p lies on the supporting conic; (false) otherwise.
   */
  bool _is_on_supporting_conic (const Point_2& p) const
  {
    // Check whether p satisfies the conic equation.
    // The point must satisfy: r*x^2 + s*y^2 + t*xy + u*x + v*y + w = 0.
    Nt_traits        nt_traits;
    const Algebraic  val = (nt_traits.convert(_r)*p.x() + 
			    nt_traits.convert(_t)*p.y() +
			    nt_traits.convert(_u)) * p.x() +
                           (nt_traits.convert(_s)*p.y() + 
			    nt_traits.convert(_v)) * p.y() +
                           nt_traits.convert(_w);

    return (CGAL::sign (val) == ZERO);
  }
 
  /*!
   * 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
  {
    CGAL_precondition (! is_full_conic());

    // Check if p is one of the endpoints.
    Alg_kernel                         ker;

    if (ker.equal_2_object() (p, _source) || 
        ker.equal_2_object() (p, _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);

    // Check if we have extra data available.
    if (_extra_data_P != NULL)
    {
      if (_extra_data_P->side != ZERO)
      {
        // In case of a hyperbolic arc, make sure the point is located on the
        // same branch as the arc.
        if (_sign_of_extra_data(p.x(), p.y()) != _extra_data_P->side)
          return (false);
      }
      else
      {
        // In case we have a segment of a line pair, make sure that p really
        // satisfies the equation of the line.
        if (_sign_of_extra_data(p.x(), p.y()) != ZERO)
          return (false);
      }
    }

    // Act according to the conic degree.
    Alg_kernel                         ker;

    if (_orient == COLLINEAR)
    {
      Comparison_result  res1;
      Comparison_result  res2;

      if (ker.compare_x_2_object() (_source, _target) == EQUAL)
      {
        // 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.
        res1 = ker.compare_y_2_object() (p, _source);
        res2 = ker.compare_y_2_object() (p, _target);
      }
      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.
        res1 = ker.compare_x_2_object() (p, _source);
        res2 = ker.compare_x_2_object() (p, _target);
      }

      // If p is not in the (open) x-range (or y-range) of the segment, it
      // cannot be contained in the segment.
      if (res1 == EQUAL || res2 == EQUAL || res1 == res2)
	return (false);

      // Perform an orientation test: This is crucial for segment of line
      // pairs, as we want to make sure that p lies on the same line as the
      // source and the target.
      return (ker.orientation_2_object()(_source, p, _target) == COLLINEAR);
    }
    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.
      if (_orient == COUNTERCLOCKWISE)
        return (ker.orientation_2_object()(_source, p, _target) == LEFT_TURN);
      else
        return (ker.orientation_2_object()(_source, p, _target) == RIGHT_TURN);
    }
  }
  
  /*!
   * 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
  {
    // In case the base conic is of degree 1 (and not 2), the arc has no
    // vertical tangency points.
    if (CGAL::sign (_s) == ZERO)
      return (0);

    // 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 Integer _two = 2;
    const Integer _four = 4;
    Algebraic     xs[2];
    Algebraic    *xs_end;
    int           n_xs;
    Nt_traits     nt_traits;
    
    xs_end = nt_traits.solve_quadratic_equation (_t*_t - _four*_r*_s,
						 _two*_t*_v - _four*_s*_u,
						 _v*_v - _four*_s*_w,
						 xs);
    n_xs = xs_end - xs;

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

    if (CGAL::sign (_t) == ZERO)
    {
      // The two vertical tangency points have the same y coordinate:
      ys[0] = nt_traits.convert (-_v) /nt_traits.convert (_two*_s);
      n_ys = 1;
    }
    else
    {
      ys_end =
        nt_traits.solve_quadratic_equation (_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);
      n_ys = ys_end - 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]);
        n++;
      }
      else
      {
        for (j = 0; j < n_ys; j++)
        {
          if (CGAL::compare (ys[j], 
                             -(nt_traits.convert(_t) * xs[i] + 
                               nt_traits.convert(_v)) /
                             nt_traits.convert(_two*_s)) == EQUAL)
          {
            ps[n] = Point_2 (xs[i], ys[j]);
            n++;
            break;
          }
        }
      }
    }

    CGAL_assertion (n <= 2);
    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 Integer _zero = 0;

    // In case the base conic is of degree 1 (and not 2), the arc has no
    // vertical tangency points.
    if (CGAL::sign (_r) == ZERO)
      return (0);

    // 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 Integer _two = 2;
    const Integer _four = 4;
    int           n;
    Algebraic     ys[2];
    Algebraic    *ys_end;
    Nt_traits     nt_traits;

    ys_end = nt_traits.solve_quadratic_equation (_t*_t - _four*_r*_s,
						 _two*_t*_u - _four*_r*_v,
						 _u*_u - _four*_r*_w,
						 ys);
    n = ys_end - ys;

    // Compute the x coordinates and construct the horizontal tangency points.
    Algebraic     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 = -(nt_traits.convert(_t)*ys[i] + nt_traits.convert(_u)) / 
        nt_traits.convert(_two*_r);

      ps[i] = Point_2 (x, ys[i]);
    }
      
    CGAL_assertion (n <= 2);
    return (n);
  }

  /*!
   * Find the y coordinates of the underlying conic at a given x coordinate.
   * \param x The x coordinate.
   * \param ys The output y coordinates. 
   * \pre The vector ys 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 Algebraic& x,
                                Algebraic *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
    Nt_traits     nt_traits;
    Algebraic     A = nt_traits.convert(_s);
    Algebraic     B = nt_traits.convert(_t)*x + nt_traits.convert(_v);
    Algebraic     C = (nt_traits.convert(_r)*x + 
                       nt_traits.convert(_u))*x + nt_traits.convert(_w);

    return (_solve_quadratic_equation (A, B, C, ys[0], ys[1]));
  }

  /*!
   * Find the x coordinates of the underlying conic at a given y coordinate.
   * \param y The y coordinate.
   * \param xs The output x coordinates. 
   * \pre The vector xs 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 Algebraic& y,
                                Algebraic *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
    Nt_traits     nt_traits;
    Algebraic     A = nt_traits.convert(_r);
    Algebraic     B = nt_traits.convert(_t)*y + nt_traits.convert(_u);
    Algebraic     C = (nt_traits.convert(_s)*y + 
                       nt_traits.convert(_v))*y + nt_traits.convert(_w);

    return (_solve_quadratic_equation (A, B, C, xs[0], xs[1]));
  }

  /*!
   * Solve the given quadratic equation: Ax^2 + B*x + C = 0.
   * \param x_minus The root obtained from taking -sqrt(discriminant).
   * \param x_plus The root obtained from taking -sqrt(discriminant).
   * \return The number of disticnt solutions to the equation.
   */
  int _solve_quadratic_equation (const Algebraic& A,
                                 const Algebraic& B,
                                 const Algebraic& C,
                                 Algebraic& x_minus, Algebraic& x_plus) const
  {
    // Check if we actually have a linear equation.
    if (CGAL::sign(A) == ZERO)
    {
      if (CGAL::sign(B) == ZERO)
	return (0);

      x_minus = x_plus = -C / B;
      return (1);
    }

    // Compute the discriminant and act according to its sign.
    const Algebraic  disc = B*B - 4*A*C;
    Sign             sign_disc = CGAL::sign (disc);

    if (sign_disc == NEGATIVE)
    {
      // No real-valued solutions:
      return (0);
    }
    else if (sign_disc == ZERO)
    {
      // One distinct solution:
      x_minus = x_plus = -B / (2*A);
      return (1);
    }

    // Compute the two distinct solutions:
    Algebraic     _2A = 2*A;
    Nt_traits     nt_traits;
    Algebraic     sqrt_disc = nt_traits.sqrt (disc);

    x_minus = -(B + sqrt_disc) / _2A;
    x_plus = (sqrt_disc - B) / _2A;
    return (2);
  }
  //@}

};

/*!
 * Exporter for conic arcs.
 */
template <class Rat_kernel, class Alg_kernel, class Nt_traits>
std::ostream& 
operator<< (std::ostream& os, 
            const _Conic_arc_2<Rat_kernel, Alg_kernel, Nt_traits> & arc)
{
  os << "{" << CGAL::to_double(arc.r()) << "*x^2 + "
     << CGAL::to_double(arc.s()) << "*y^2 + "
     << CGAL::to_double(arc.t()) << "*xy + " 
     << CGAL::to_double(arc.u()) << "*x + "
     << CGAL::to_double(arc.v()) << "*y + "
     << CGAL::to_double(arc.w()) << "}";

  if (arc.is_full_conic())
  {
    os << " - Full curve";
  }
  else
  {
    os << " : (" << CGAL::to_double(arc.source().x()) << "," 
       << CGAL::to_double(arc.source().y()) << ") ";

    if (arc.orientation() == CLOCKWISE)
      os << "--cw-->";
    else if (arc.orientation() == COUNTERCLOCKWISE)
      os << "--ccw-->";
    else
      os << "--l-->";

    os << " (" << CGAL::to_double(arc.target().x()) << "," 
       << CGAL::to_double(arc.target().y()) << ")";
  }

  return (os);
}

CGAL_END_NAMESPACE

#endif