rational_arc_2.h

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

// Copyright (c) 2005  Tel-Aviv University (Israel).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you may redistribute it under
// the terms of the Q Public License version 1.0.
// See the file LICENSE.QPL distributed with CGAL.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.3-branch/Arrangement_2/include/CGAL/Arr_traits_2/Rational_arc_2.h $
// $Id: Rational_arc_2.h 38368 2007-04-20 11:28:33Z efif $
// 
//
// Author(s)     : Ron Wein <[e-mail]>

#ifndef CGAL_RATIONAL_ARC_2_H
#define CGAL_RATIONAL_ARC_2_H

/*! \file
 * Header file for the _Rational_arc_2 class.
 */

#include <vector>
#include <list>
#include <ostream>
#include <CGAL/Arr_enums.h>

CGAL_BEGIN_NAMESPACE

/*! \class
 * Representation of an segment of a rational function, given as:
 *
 *        D(x)
 *   y = ------              x_l <= x <= x_r
 *        N(x)
 *
 * where D and N are polynomial with integer (or rational) coefficients.
 * The class is templated with two parameters: 
 * Alg_kernel A geometric kernel, where Alg_kernel::FT is the number type
 *            for the coordinates of arrangement vertices, which are algebraic
 *            numbers (defined by Nt_traits::Algebraic).
 * Nt_traits A traits class for performing various operations on the integer,
 *           rational and algebraic types. 
 */

template <class Alg_kernel_, class Nt_traits_>
class _Rational_arc_2
{
public:

  typedef Alg_kernel_                             Alg_kernel;
  typedef Nt_traits_                              Nt_traits;
  typedef _Rational_arc_2<Alg_kernel, Nt_traits>  Self;
  
  typedef typename Alg_kernel::FT                 Algebraic;
  typedef typename Alg_kernel::Point_2            Point_2;

  typedef typename Nt_traits::Integer             Integer;
  typedef typename Nt_traits::Rational            Rational;
  typedef std::vector<Rational>                   Rat_vector;
  typedef typename Nt_traits::Polynomial          Polynomial;

private:

  typedef std::pair<Point_2, unsigned int>        Intersection_point_2;

  enum
  {
    SRC_AT_X_MINUS_INFTY = 1,
    SRC_AT_X_PLUS_INFTY = 2,
    SRC_AT_Y_MINUS_INFTY = 4,
    SRC_AT_Y_PLUS_INFTY = 8,

    SRC_INFO_BITS = SRC_AT_X_MINUS_INFTY + SRC_AT_X_PLUS_INFTY +
                    SRC_AT_Y_MINUS_INFTY + SRC_AT_Y_PLUS_INFTY,

    TRG_AT_X_MINUS_INFTY = 16,
    TRG_AT_X_PLUS_INFTY = 32,
    TRG_AT_Y_MINUS_INFTY = 64,
    TRG_AT_Y_PLUS_INFTY = 128,

    TRG_INFO_BITS = TRG_AT_X_MINUS_INFTY + TRG_AT_X_PLUS_INFTY +
                    TRG_AT_Y_MINUS_INFTY + TRG_AT_Y_PLUS_INFTY,

    IS_DIRECTED_RIGHT = 256,
    IS_CONTINUOUS = 512,
    IS_VALID = 1024
  };
 
  Polynomial        _numer;       // The polynomial in the numerator.
  Polynomial        _denom;       // The polynomial in the denominator.
  Point_2           _ps;          // The source point.
  Point_2           _pt;          // The target point.
  int               _info;        // A set of Boolean flags.  

public:

  /// \name Constrcution methods.
  //@{

  /*!
   * Default constructor.
   */
  _Rational_arc_2 () :
    _info (0)
  {}

  /*!
   * Constructor of a whole polynomial curve.
   * \param pcoeffs The rational coefficients of the polynomial p(x).
   */
  _Rational_arc_2 (const Rat_vector& pcoeffs) :
    _info (0)
  {
    // Mark that the endpoints of the polynomial are unbounded (the source is
    // at x = -oo and the target is at x = +oo).
    _info = (_info | SRC_AT_X_MINUS_INFTY);
    _info = (_info | TRG_AT_X_PLUS_INFTY);
    _info = (_info | IS_DIRECTED_RIGHT);

    // Set the numerator polynomial.
    Nt_traits    nt_traits;
    Integer      p_factor;

    nt_traits.construct_polynomial (&(pcoeffs[0]),
                                    pcoeffs.size() - 1,
                                    _numer,
                                    p_factor);

    // Define the denominator to be a constant polynomial.
    Integer      denom_coeffs[1];

    denom_coeffs [0] = p_factor;
    _denom = nt_traits.construct_polynomial (denom_coeffs, 0);
    
    // Check whether the end points lie at y = -oo or at y = +oo.
    const int    deg_num = nt_traits.degree (_numer);
    CGAL::Sign   lead_sign;

    if (deg_num > 0)
    {
      // Check if the degree is even or odd and check the sign of the leading
      // coefficient of the polynomial.
      lead_sign = CGAL::sign (nt_traits.get_coefficient (_numer, deg_num));
      CGAL_assertion (lead_sign != CGAL::ZERO);

      if (deg_num % 2 == 0)
      {
        // Polynomial of an even degree.
        if (lead_sign == CGAL::NEGATIVE)
          _info = (_info | SRC_AT_Y_MINUS_INFTY | TRG_AT_Y_MINUS_INFTY);
        else
          _info = (_info | SRC_AT_Y_PLUS_INFTY | TRG_AT_Y_PLUS_INFTY);
      }
      else
      {
        // Polynomial of an odd degree.
        if (lead_sign == CGAL::NEGATIVE)
          _info = (_info | SRC_AT_Y_PLUS_INFTY | TRG_AT_Y_MINUS_INFTY);
        else
          _info = (_info | SRC_AT_Y_MINUS_INFTY | TRG_AT_Y_PLUS_INFTY);
      }
    }
    else
    {
      // In the case of a constant polynomial it is possible to set a finite
      // y-coordinate for the source and target points.
      Integer lead_coeff = (deg_num < 0) ? Integer(0) : 
        nt_traits.get_coefficient (_numer, deg_num);

      _ps = _pt = Point_2 (Algebraic(), nt_traits.convert (lead_coeff));
    }

    // Mark that the arc is continuous and valid.
    _info = (_info | IS_CONTINUOUS);
    _info = (_info | IS_VALID);
  }

  /*!
   * Constructor of a polynomial ray, defined by y = p(x), for x_s <= x if the
   * ray is directed to the right, or for x_s >= x if it is directed to the
   * left.
   * \param pcoeffs The rational coefficients of the polynomial p(x).
   * \param x_s The x-coordinate of the source point.
   * \param dir_right Is the ray directed to the right (to +oo)
   *                  or to the left (to -oo).
   */
  _Rational_arc_2 (const Rat_vector& pcoeffs,
                   const Algebraic& x_s, bool dir_right) :
    _info (0)
  {
    // Mark that the target points of the polynomial is unbounded.
    if (dir_right)
    {
      _info = (_info | TRG_AT_X_PLUS_INFTY);
      _info = (_info | IS_DIRECTED_RIGHT);
    }
    else
    {
      _info = (_info | TRG_AT_X_MINUS_INFTY);
    }

    // Set the numerator polynomial.
    Nt_traits    nt_traits;
    Integer      p_factor;

    nt_traits.construct_polynomial (&(pcoeffs[0]),
                                    pcoeffs.size() - 1,
                                    _numer,
                                    p_factor);

    // Define the denominator to be a constant polynomial.
    Integer      denom_coeffs[1];

    denom_coeffs [0] = p_factor;
    _denom = nt_traits.construct_polynomial (denom_coeffs, 0);

    // Set the source point.
    _ps = Point_2 (x_s, nt_traits.evaluate_at (_numer, x_s));
    
    // Check whether the target point lies at y = -oo or at y = +oo.
    const int    deg_num = nt_traits.degree (_numer);
    CGAL::Sign   lead_sign;

    if (deg_num > 0)
    {
      // Check if the degree is even or odd and check the sign of the leading
      // coefficient of the polynomial.
      lead_sign = CGAL::sign (nt_traits.get_coefficient (_numer, deg_num));
      CGAL_assertion (lead_sign != CGAL::ZERO);

      if (dir_right)
      {
        // The target is at x= +oo, thus:
        if (lead_sign == CGAL::POSITIVE)
          _info = (_info | TRG_AT_Y_PLUS_INFTY);
        else
          _info = (_info | TRG_AT_Y_MINUS_INFTY);
      }
      else
      {
        // The target is at x= -oo, thus:
        if ((deg_num % 2 == 0 && lead_sign == CGAL::POSITIVE) ||
            (deg_num % 2 == 1 && lead_sign == CGAL::NEGATIVE))
          _info = (_info | TRG_AT_Y_PLUS_INFTY);
        else
          _info = (_info | TRG_AT_Y_MINUS_INFTY);
      }
    }
    else
    {
      // In the case of a constant polynomial it is possible to set a finite
      // y-coordinate for the target point.
      Integer lead_coeff = (deg_num < 0) ? Integer(0) : 
        nt_traits.get_coefficient (_numer, deg_num);

      _pt = Point_2 (Algebraic(), nt_traits.convert (lead_coeff));
    }

    // Mark that the arc is continuous and valid.
    _info = (_info | IS_CONTINUOUS);
    _info = (_info | IS_VALID);
  }

  /*!
   * Constructor of a polynomial arc, defined by y = p(x), x_min <= x <= x_max.
   * \param pcoeffs The rational coefficients of the polynomial p(x).
   * \param x_s The x-coordinate of the source point.
   * \param x_t The x-coordinate of the target point.
   * \pre The two x-coordinate must not be equal.
   */
  _Rational_arc_2 (const Rat_vector& pcoeffs,
		   const Algebraic& x_s, const Algebraic& x_t) :
    _info (0)
  {
    // Compare the x-coordinates and determine the direction.
    Comparison_result   x_res = CGAL::compare (x_s, x_t);

    CGAL_precondition (x_res != EQUAL);
    if (x_res == EQUAL)
      return;

    if (x_res == SMALLER)
      _info = (_info | IS_DIRECTED_RIGHT);

    // Set the numerator polynomial.
    Nt_traits    nt_traits;
    Integer      p_factor;

    nt_traits.construct_polynomial (&(pcoeffs[0]),
                                    pcoeffs.size() - 1,
                                    _numer,
                                    p_factor);

    // Define the denominator to be a constant polynomial.
    Integer      denom_coeffs[1];

    denom_coeffs [0] = p_factor;
    _denom = nt_traits.construct_polynomial (denom_coeffs, 0);
    
    // Set the endpoints.
    _ps = Point_2 (x_s, nt_traits.evaluate_at (_numer, x_s));
    _pt = Point_2 (x_t, nt_traits.evaluate_at (_numer, x_t));

    // Mark that the arc is continuous and valid.
    _info = (_info | IS_CONTINUOUS);
    _info = (_info | IS_VALID);
  }    

  /*!
   * Constructor of a polynomial function, defined by y = p(x)/q(x) for any x.
   * \param pcoeffs The rational coefficients of the polynomial p(x).
   * \param qcoeffs The rational coefficients of the polynomial q(x).
   */
  _Rational_arc_2 (const Rat_vector& pcoeffs, const Rat_vector& qcoeffs) :
    _info (0)
  {
    // Mark that the endpoints of the rational functions are unbounded (the
    // source is at x = -oo and the target is at x = +oo).
    _info = (_info | SRC_AT_X_MINUS_INFTY);
    _info = (_info | TRG_AT_X_PLUS_INFTY);
    _info = (_info | IS_DIRECTED_RIGHT);

    // Set the numerator and denominator polynomials.
    Nt_traits    nt_traits;
    const bool   valid = 
      nt_traits.construct_polynomials (&(pcoeffs[0]),
					pcoeffs.size() - 1,
                                       &(qcoeffs[0]),
					qcoeffs.size() - 1,
                                       _numer, _denom);

    CGAL_precondition_msg (valid,
                           "The denominator polynomial must not be 0.");
    if (! valid)
      return;

    // Analyze the bahaviour of the rational function at x = -oo (the source).
    Algebraic           y0;
    const Boundary_type inf_s = _analyze_at_minus_infinity (_numer, _denom,
                                                            y0);

    if (inf_s == MINUS_INFINITY)
      _info = (_info | SRC_AT_Y_MINUS_INFTY);
    else if (inf_s == PLUS_INFINITY)
      _info = (_info | SRC_AT_Y_PLUS_INFTY);
    else // if (inf_s == NO_BOUNDARY)
      _ps = Point_2 (0, y0);

    // Analyze the bahaviour of the rational function at x = +oo (the target).
    const Boundary_type inf_t = _analyze_at_plus_infinity (_numer, _denom,
                                                           y0);

    if (inf_t == MINUS_INFINITY)
      _info = (_info | TRG_AT_Y_MINUS_INFTY);
    else if (inf_t == PLUS_INFINITY)
      _info = (_info | TRG_AT_Y_PLUS_INFTY);
    else // if (inf_t == NO_BOUNDARY)
      _pt = Point_2 (0, y0);

    // Mark that the arc is valid. As it may have poles, we do not mark it
    // as continuous.
    _info = (_info | IS_VALID);
  }

  /*!
   * Constructor of a ray of a rational function, defined by y = p(x)/q(x),
   * for x_s <= x if the ray is directed to the right, or for x_s >= x if it
   * is directed to the left.
   * \param pcoeffs The rational coefficients of the polynomial p(x).
   * \param qcoeffs The rational coefficients of the polynomial q(x).
   * \param x_s The x-coordinate of the source point.
   * \param dir_right Is the ray directed to the right (to +oo)
   *                  or to the left (to -oo).
   */
  _Rational_arc_2 (const Rat_vector& pcoeffs, const Rat_vector& qcoeffs,
                   const Algebraic& x_s, bool dir_right) :
    _info (0)
  {
    // Mark that the target points of the polynomial is unbounded.
    if (dir_right)
    {
      _info = (_info | TRG_AT_X_PLUS_INFTY);
      _info = (_info | IS_DIRECTED_RIGHT);
    }
    else
    {
      _info = (_info | TRG_AT_X_MINUS_INFTY);
    }

    // Set the numerator and denominator polynomials.
    Nt_traits    nt_traits;
    const bool   valid = 
      nt_traits.construct_polynomials (&(pcoeffs[0]),
					pcoeffs.size() - 1,
                                       &(qcoeffs[0]),
					qcoeffs.size() - 1,
                                       _numer, _denom);

    CGAL_precondition_msg (valid,
                           "The denominator polynomial must not be 0.");
    if (! valid)
      return;

    // The source point has a bounded x-coordinate. Set the source point
    // and check if it lies next to a pole.
    if (CGAL::sign (nt_traits.evaluate_at (_denom, x_s)) != CGAL::ZERO)
    {
      // We have a nomral endpoint.
      _ps = Point_2 (x_s, nt_traits.evaluate_at (_numer, x_s) /
                     nt_traits.evaluate_at (_denom, x_s));
    }
    else
    {
      // The y-coodinate is unbounded, but we can set its sign.
      _ps = Point_2 (x_s, 0);
      
      std::pair<CGAL::Sign, CGAL::Sign>  signs = _analyze_near_pole (x_s);