conic_arc_2.h

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

      Point_2    tan_ps[2];
      int        n_tan_ps;
      int        i;

      n_tan_ps = vertical_tangency_points (tan_ps);
      for (i = 0; i < n_tan_ps; i++)
      {
        if (CGAL::to_double(tan_ps[i].x()) < x_min)
          x_min = CGAL::to_double(tan_ps[i].x());
        if (CGAL::to_double(tan_ps[i].x()) > x_max)
          x_max = CGAL::to_double(tan_ps[i].x());
      }

      // Go over the horizontal tangency points and try to update the y-points.
      n_tan_ps = horizontal_tangency_points (tan_ps);
      for (i = 0; i < n_tan_ps; i++)
      {
        if (CGAL::to_double(tan_ps[i].y()) < y_min)
          y_min = CGAL::to_double(tan_ps[i].y());
        if (CGAL::to_double(tan_ps[i].y()) > y_max)
          y_max = CGAL::to_double(tan_ps[i].y());
      }
    }

    // Return the resulting bounding box.
    return (Bbox_2 (x_min, y_min, x_max, y_max));
  }
  //@}

  /// \name Modifying functions.
  //@{

  /*!
   * Set the source point of the conic arc.
   * \param ps The new source point.
   * \pre The arc is not a full conic curve.
   *      ps must lie on the supporting conic curve.
   */
  void set_source (const Point_2& ps)
  {
    CGAL_precondition (! is_full_conic());
    CGAL_precondition (_is_on_supporting_conic (ps));
    CGAL_precondition (Alg_kernel().orientation_2_object()
                       (_source, ps, _target) == _orient ||
                       Alg_kernel().orientation_2_object()
                       (ps, _source, _target) == _orient);

    _source = ps;
    return;
  }

  /*!
   * Set the target point of the conic arc.
   * \param pt The new source point.
   * \pre The arc is not a full conic curve.
   *      pt must lie on the supporting conic curve.
   */
  void set_target (const Point_2& pt)
  {
    CGAL_precondition (! is_full_conic());
    CGAL_precondition (_is_on_supporting_conic (pt));
    CGAL_precondition (Alg_kernel().orientation_2_object()
                       (_source, pt, _target) == _orient ||
                       Alg_kernel().orientation_2_object()
                       (_source, _target, pt) == _orient);

    _target = pt;
    return;
  }

  //@}

  /// \name Compute points on the arc.
  //@{

  /*!
   * Calculate the vertical tangency points of the arc.
   * \param vpts The vertical tangency points.
   * \pre The vpts vector should be allocated at the size of 2.
   * \return The number of vertical tangency points.
   */
  int vertical_tangency_points (Point_2* vpts) const
  {
    // No vertical tangency points for line segments:
    if (_orient == COLLINEAR)
      return (0);

    // Calculate the vertical tangency points of the supporting conic.
    Point_2 ps[2];
    int     n;

    n = _conic_vertical_tangency_points (ps);

    // Return only the points that are contained in the arc interior.
    int    m = 0;
    
    for (int i = 0; i < n; i++)
    {
      if (is_full_conic() || _is_strictly_between_endpoints(ps[i]))
      {
        vpts[m] = ps[i];
        m++;
      }
    }

    // Return the number of vertical tangency points found.
    CGAL_assertion (m <= 2);
    return (m);
  }

  /*!
   * Calculate the horizontal tangency points of the arc.
   * \param hpts The horizontal tangency points.
   * \pre The hpts vector should be allocated at the size of 2.
   * \return The number of horizontal tangency points.
   */
  int horizontal_tangency_points (Point_2* hpts) const
  {
    // No horizontal tangency points for line segments:
    if (_orient == COLLINEAR)
      return (0);

    // Calculate the horizontal tangency points of the conic.
    Point_2    ps[2];
    int        n;

    n = _conic_horizontal_tangency_points (ps);

    // Return only the points that are contained in the arc interior.
    int    m = 0;
    
    for (int i = 0; i < n; i++)
    {
      if (is_full_conic() || _is_strictly_between_endpoints(ps[i]))
      {
        hpts[m] = ps[i];
        m++;
      }
    }

    // Return the number of horizontal tangency points found.
    CGAL_assertion (m <= 2);
    return (m);
  }

  /*!
   * Find all points on the arc with a given x-coordinate.
   * \param p A placeholder for the x-coordinate.
   * \param ps The point on the arc at x(p).
   * \pre The vector ps should be allocated at the size of 2.
   * \return The number of points found.
   */
  int get_points_at_x (const Point_2& p,
                       Point_2 *ps) const
  {
    // Get the y coordinates of the points on the conic.
    Algebraic    ys[2];
    int          n;

    n = _conic_get_y_coordinates (p.x(), ys);
    
    // Find all the points that are contained in the arc.
    int   m = 0;
    
    for (int i = 0; i < n; i++)
    {
      ps[m] = Point_2 (p.x(), ys[i]);

      if (is_full_conic() || _is_between_endpoints(ps[m]))
        m++;
    }

    // Return the number of points on the arc.
    CGAL_assertion (m <= 2);
    return (m);
  }

  /*!
   * Find all points on the arc with a given y-coordinate.
   * \param p A placeholder for the y-coordinate.
   * \param ps The point on the arc at x(p).
   * \pre The vector ps should be allocated at the size of 2.
   * \return The number of points found.
   */
  int get_points_at_y (const Point_2& p,
                       Point_2 *ps) const
  {
    // Get the y coordinates of the points on the conic.
    Algebraic    xs[2];
    int          n;

    n = _conic_get_x_coordinates (p.y(), xs);
    
    // Find all the points that are contained in the arc.
    int   m = 0;
    
    for (int i = 0; i < n; i++)
    {
      ps[m] = Point_2 (xs[i], p.y());

      if (is_full_conic() || _is_between_endpoints(ps[m]))
        m++;
    }

    // Return the number of points on the arc.
    CGAL_assertion (m <= 2);
    return (m);
  }
  //@}

private:

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

  /*!
   * Set the properties of a conic arc (for the usage of the constructors).
   * \param rat_coeffs A vector of size 6, storing the rational coefficients
   *                   of x^2, y^2, xy, x, y and the free coefficient resp.
   */
  void _set (const Rational* rat_coeffs)
  {
    _info = IS_VALID;

    // Convert the coefficients vector to an equivalent vector of integer
    // coefficients.
    Nt_traits         nt_traits;
    Integer           int_coeffs[6];

    nt_traits.convert_coefficients (rat_coeffs, rat_coeffs + 6,
                                    int_coeffs);

    // Check the orientation of conic curve, and negate the conic coefficients
    // if its given orientation.
    typename Rat_kernel::Conic_2   temp_conic (rat_coeffs[0], rat_coeffs[1],
                                               rat_coeffs[2], rat_coeffs[3],
                                               rat_coeffs[4], rat_coeffs[5]);

    if (_orient == temp_conic.orientation())
    {
      _r = int_coeffs[0];
      _s = int_coeffs[1];
      _t = int_coeffs[2];
      _u = int_coeffs[3];
      _v = int_coeffs[4];
      _w = int_coeffs[5];
    }
    else
    {
      _r = -int_coeffs[0];
      _s = -int_coeffs[1];
      _t = -int_coeffs[2];
      _u = -int_coeffs[3];
      _v = -int_coeffs[4];
      _w = -int_coeffs[5];
    }

    // Make sure both endpoint lie on the supporting conic.
    if (! _is_on_supporting_conic (_source) ||
	! _is_on_supporting_conic (_target))
    {
      _info = 0;          // Invalid arc.
      return;
    }

    _extra_data_P = NULL;

    // Check whether we have a degree 2 curve.
    if ((CGAL::sign (_r) != ZERO ||
	 CGAL::sign (_s) != ZERO ||
	 CGAL::sign (_t) != ZERO))
    {
      if (_orient == COLLINEAR)
      {
        // We have a segment of a line pair with rational coefficients.
        // Compose the equation of the underlying line
        // (with algebraic coefficients).
        const Algebraic        x1 = _source.x(), y1 = _source.y();
        const Algebraic        x2 = _target.x(), y2 = _target.y();
        
        // The supporting line is A*x + B*y + C = 0, where:
        //
        //  A = y2 - y1,    B = x1 - x2,    C = x2*y1 - x1*y2 
        //
        // We use the extra dat field to store the equation of this line.
        _extra_data_P = new Extra_data;
        _extra_data_P->a = y2 - y1;
        _extra_data_P->b = x1 - x2;
        _extra_data_P->c = x2*y1 - x1*y2;
        _extra_data_P->side = ZERO;

        // Make sure the midpoint is on the line pair (thus making sure that
        // the two points are not taken from different lines).
        Alg_kernel       ker;
        Point_2          p_mid = ker.construct_midpoint_2_object() (_source,
                                                                    _target);
        
        if (CGAL::sign ((nt_traits.convert(_r)*p_mid.x() + 
                         nt_traits.convert(_t)*p_mid.y() + 
                         nt_traits.convert(_u)) * p_mid.x() +
                        (nt_traits.convert(_s)*p_mid.y() +
                         nt_traits.convert(_v)) * p_mid.y() +
                        nt_traits.convert(_w)) != ZERO)
        {
          _info = 0;          // Invalid arc.
          return;
        }
      }
      else 
      {
        // The sign of (4rs - t^2) detetmines the conic type:
        // - if it is possitive, the conic is an ellipse,
        // - if it is negative, the conic is a hyperbola,
        // - if it is zero, the conic is a parabola.
        CGAL::Sign   sign_conic = CGAL::sign (4*_r*_s - _t*_t);

        if (sign_conic == NEGATIVE)
          // Build the extra hyperbolic data
          _build_hyperbolic_arc_data ();

        if (sign_conic != POSITIVE)
        {
          // In case of a non-degenerate parabola or a hyperbola, make sure 
          // the arc is not infinite.
          Alg_kernel       ker;
          Point_2          p_mid = ker.construct_midpoint_2_object() (_source,
                                                                      _target);
          Point_2          ps[2];

          bool  finite_at_x = (get_points_at_x(p_mid, ps) > 0);
          bool  finite_at_y = (get_points_at_y(p_mid, ps) > 0);
      
          if (! finite_at_x && ! finite_at_y)
          {
            _info = 0;          // Invalid arc.
            return;
          }
        }
      }
    }

    // Mark that this arc valid and is not a full conic curve.
    _info = IS_VALID;

    return;
  }

  /*!
   * Set the properties of a conic arc that is really a full curve
   * (that is, an ellipse).
   * \param rat_coeffs A vector of size 6, storing the rational coefficients
   *                   of x^2, y^2, xy, x, y and the free coefficient resp.
   * \param comp_orient Should we compute the orientation of the given curve.
   */
  void _set_full (const Rational* rat_coeffs,
                  const bool& comp_orient)
  {
    // Convert the coefficients vector to an equivalent vector of integer
    // coefficients.
    Nt_traits         nt_traits;
    Integer           int_coeffs[6];

    nt_traits.convert_coefficients (rat_coeffs, rat_coeffs + 6,
                                    int_coeffs);

    // Check the orientation of conic curve, and negate the conic coefficients
    // if its given orientation.
    typename Rat_kernel::Conic_2   temp_conic (rat_coeffs[0], rat_coeffs[1],
                                               rat_coeffs[2], rat_coeffs[3],
                                               rat_coeffs[4], rat_coeffs[5]);
    const Orientation              temp_orient = temp_conic.orientation();

    if (comp_orient)
      _orient = temp_orient;

    if (_orient == temp_orient)
    {
      _r = int_coeffs[0];
      _s = int_coeffs[1];
      _t = int_coeffs[2];
      _u = int_coeffs[3];
      _v = int_coeffs[4];
      _w = int_coeffs[5];
    }
    else
    {
      _r = -int_coeffs[0];
      _s = -int_coeffs[1];
      _t = -int_coeffs[2];
      _u = -int_coeffs[3];
      _v = -int_coeffs[4];
      _w = -int_coeffs[5];
    }

    // Make sure the conic is a non-degenerate ellipse:
    // The coefficients should satisfy (4rs - t^2) > 0.
    const bool  is_ellipse = (CGAL::sign (4*_r*_s - _t*_t) == POSITIVE);
    CGAL_assertion (is_ellipse);

    // We do not have to store any extra data with the arc.
    _extra_data_P = NULL;

    // Mark that this arc is a full conic curve.
    if (is_ellipse)
      _info = IS_VALID | IS_FULL_CONIC;
    else
      _info = 0;

    return;
  }

  /*!
   * Build the data for hyperbolic arc, contaning the characterization of the
   * hyperbolic branch the arc is placed on.
   */
  void _build_hyperbolic_arc_data ()
  {
    // Let phi be the rotation angle of the conic from its canonic form.
    // We can write:
    // 
    //                          t
    //  sin(2*phi) = -----------------------
    //                sqrt((r - s)^2 + t^2)
    // 
    //                        r - s
    //  cos(2*phi) = -----------------------
    //                sqrt((r - s)^2 + t^2)
    //
    Nt_traits        nt_traits;
    const int        or_fact = (_orient == CLOCKWISE) ? -1 : 1;
    const Algebraic  r = nt_traits.convert (or_fact * _r);
    const Algebraic  s = nt_traits.convert (or_fact * _s);
    const Algebraic  t = nt_traits.convert (or_fact * _t);
    const Algebraic  cos_2phi = (r - s) / nt_traits.sqrt((r-s)*(r-s) + t*t);
    const Algebraic  _zero = 0;
    const Algebraic  _one = 1;
    const Algebraic  _two = 2;
    Algebraic        sin_phi;
    Algebraic        cos_phi;

    // Calculate sin(phi) and cos(phi) according to the half-angle formulae:
    // 
    //  sin(phi)^2 = 0.5 * (1 - cos(2*phi))
    //  cos(phi)^2 = 0.5 * (1 + cos(2*phi))
    Sign             sign_t = CGAL::sign (t);

    if (sign_t == ZERO)
    {
      // sin(2*phi) == 0, so phi = 0 or phi = PI/2
      if (CGAL::sign (cos_2phi) == POSITIVE)
      {
        // phi = 0.
        sin_phi = _zero;
        cos_phi = _one;
      }
      else
      {
        // phi = PI/2.
        sin_phi = _one;
        cos_phi = _zero;
      }
    }
    else if (sign_t == POSITIVE)
    {
      // sin(2*phi) > 0 so 0 < phi < PI/2.
      sin_phi = nt_traits.sqrt((_one - cos_2phi) / _two);
      cos_phi = nt_traits.sqrt((_one + cos_2phi) / _two);
    }