function_objects.h

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

		       const Compare_xy_2& cxy_)
      : cp(cp_), cxy(cxy_)
    {}

    result_type
    operator()( const Ray_2& r, const Point_2& p) const
    {
      const Point_2 & source = cp(r,0);
      return p == source || Direction_2(p - source) == r.direction();
    } // FIXME

    result_type
    operator()( const Segment_2& s, const Point_2& p) const
    {
      return co(cp(s,0), p, cp(s,1));
    }
  };

  template <typename K>
  class Collinear_2
  {
    typedef typename K::Point_2        Point_2;
    typedef typename K::Orientation_2  Orientation_2;
    Orientation_2 o;
  public:
    typedef typename K::Bool_type      result_type;
    typedef Arity_tag< 3 >             Arity;

    Collinear_2() {}
    Collinear_2(const Orientation_2 o_) : o(o_) {}

    result_type
    operator()(const Point_2& p, const Point_2& q, const Point_2& r) const
    {
      typedef typename K::RT RT;

      const RT& phx = p.hx();
      const RT& phy = p.hy();
      const RT& phw = p.hw();
      const RT& qhx = q.hx();
      const RT& qhy = q.hy();
      const RT& qhw = q.hw();
      const RT& rhx = r.hx();
      const RT& rhy = r.hy();
      const RT& rhw = r.hw();

      // | A B |
      // | C D |

      RT  A = phx*rhw - phw*rhx;
      RT  B = phy*rhw - phw*rhy;
      RT  C = qhx*rhw - qhw*rhx;
      RT  D = qhy*rhw - qhw*rhy;

      RT  det =  A*D - B*C;

      /*
	RT det_old =   p.hx() * (q.hy()*r.hw() - q.hw()*r.hy() )
	+ p.hy() * (q.hw()*r.hx() - q.hx()*r.hw() )
	+ p.hw() * (q.hx()*r.hy() - q.hy()*r.hx() );
	
	if ( !(CGAL_NTS sign(det) == CGAL_NTS sign(det_old)) )
	{
	std::cerr << "det: " << det << " det_old: " << det_old << flush;
	}
      */

      return CGAL_NTS is_zero(det);
    }
  };

  template <typename K>
  class Compare_angle_with_x_axis_2
  {
    typedef typename K::Point_2            Point_2;
    typedef typename K::Vector_2           Vector_2;
    typedef typename K::Direction_2        Direction_2;
  public:
    typedef typename K::Comparison_result  result_type;
    typedef Arity_tag< 2 >                 Arity;

    result_type
    operator()(const Direction_2& d1, const Direction_2& d2) const
    {
      typedef typename K::RT  RT;
      CGAL_kernel_precondition(
          static_cast<int>(COUNTERCLOCKWISE) == static_cast<int>(LARGER)
       && static_cast<int>(COLLINEAR)        == static_cast<int>(EQUAL)
       && static_cast<int>(CLOCKWISE)        == static_cast<int>(SMALLER) );

      const RT RT0(0);

      Vector_2 dirvec1(d1.x(), d1.y());      // Added
      Point_2   p1 = CGAL::ORIGIN + dirvec1; // Added
      Vector_2 dirvec2(d2.x(), d2.y());      // Added
      Point_2   p2 = ORIGIN + dirvec2;       // Added
      //  Point_2   p1 = ORIGIN + d1.vector(); // Commented out
      //  Point_2   p2 = ORIGIN + d2.vector(); // Commented out

      CGAL_kernel_precondition( RT0 < p1.hw() );
      CGAL_kernel_precondition( RT0 < p2.hw() );

      int       x_sign1 = static_cast<int>(CGAL_NTS sign( p1.hx() ));
      int       x_sign2 = static_cast<int>(CGAL_NTS sign( p2.hx() ));
      int       y_sign1 = static_cast<int>(CGAL_NTS sign( p1.hy() ));
      int       y_sign2 = static_cast<int>(CGAL_NTS sign( p2.hy() ));

      if ( y_sign1 * y_sign2 < 0)
	{
	  return (0 < y_sign1 ) ? SMALLER : LARGER;
	}

      Point_2   origin( RT0  , RT0   );

      if ( 0 < y_sign1 * y_sign2 )
	{
	  return enum_cast<Comparison_result>(orientation(origin, p2, p1));

	  // Precondition on the enums:
	  // COUNTERCLOCKWISE == LARGER   ( ==  1 )
	  // COLLINEAR        == EQUAL    ( ==  0 )
	  // CLOCKWISE        == SMALLER  ( == -1 )
	}

      // ( y_sign1 * y_sign2 == 0 )

      bool b1 = (y_sign1 == 0) && (x_sign1 >= 0);
      bool b2 = (y_sign2 == 0) && (x_sign2 >= 0);

      if ( b1 ) { return  b2 ? EQUAL : SMALLER; }
      if ( b2 ) { return  b1 ? EQUAL : LARGER; }
      if ( y_sign1 == y_sign2 )  // == 0
	  return EQUAL;
      else
	  return (orientation(origin, p1, p2) == COUNTERCLOCKWISE) ?
	    SMALLER : LARGER;
    }
  };

  template <typename K>
  class Compare_distance_2
  {
    typedef typename K::Point_2            Point_2;
  public:
    typedef typename K::Comparison_result  result_type;
    typedef Arity_tag< 3 >                 Arity;

    result_type
    operator()(const Point_2& p, const Point_2& q, const Point_2& r) const
    {
      typedef typename K::RT RT;

      const RT & phx = p.hx();
      const RT & phy = p.hy();
      const RT & phw = p.hw();
      const RT & qhx = q.hx();
      const RT & qhy = q.hy();
      const RT & qhw = q.hw();
      const RT & rhx = r.hx();
      const RT & rhy = r.hy();
      const RT & rhw = r.hw();

      RT dosd =   // difference of squared distances

	//            phx * phx   *   qhw * qhw * rhw * rhw
	//   -RT(2) * phx * qhx   *   phw * qhw * rhw * rhw
	//   +        qhx * qhx   *   phw * phw * rhw * rhw
	//
	//   +        phy * phy   *   qhw * qhw * rhw * rhw
	//   -RT(2) * phy * qhy   *   phw * qhw * rhw * rhw
	//   +        qhy * qhy   *   phw * phw * rhw * rhw
	//
	// - (        phx * phx   *   qhw * qhw * rhw * rhw
	//   -RT(2) * phx * rhx   *   phw * qhw * qhw * rhw
	//   +        rhx * rhx   *   phw * phw * qhw * qhw
	//
	//   +        phy * phy   *   qhw * qhw * rhw * rhw
	//   -RT(2) * phy * rhy   *   phw * qhw * qhw * rhw
	//   +        rhy * rhy   *   phw * phw * qhw * qhw
	
	rhw*rhw * (         phw * ( qhx*qhx + qhy*qhy )
			    - 2 * qhw * ( phx*qhx + phy*qhy )
			    )
	- qhw*qhw * (         phw * ( rhx*rhx + rhy*rhy )
			      - 2 * rhw * ( phx*rhx + phy*rhy )
			      );

      return enum_cast<Comparison_result>(CGAL_NTS sign(dosd));
    }
  };

  template <typename K>
  class Compare_squared_distance_2
  {
    typedef typename K::Point_2            Point_2;
    typedef typename K::FT                 FT;
  public:
    typedef typename K::Comparison_result  result_type;
    typedef Arity_tag< 3 >                 Arity;

    result_type
    operator()(const Point_2& p, const Point_2& q, const FT& d2) const
    {
      return CGAL_NTS compare(squared_distance(p, q), d2);
    }
  };

  template <typename K>
  class Compare_distance_3
  {
    typedef typename K::Point_3            Point_3;
  public:
    typedef typename K::Comparison_result  result_type;
    typedef Arity_tag< 3 >                 Arity;

    result_type
    operator()(const Point_3& p, const Point_3& q, const Point_3& r) const
    {
      typedef typename K::RT RT;

      const RT & phx = p.hx();
      const RT & phy = p.hy();
      const RT & phz = p.hz();
      const RT & phw = p.hw();
      const RT & qhx = q.hx();
      const RT & qhy = q.hy();
      const RT & qhz = q.hz();
      const RT & qhw = q.hw();
      const RT & rhx = r.hx();
      const RT & rhy = r.hy();
      const RT & rhz = r.hz();
      const RT & rhw = r.hw();

      RT dosd =   // difference of squared distances

	rhw*rhw * (         phw * ( qhx*qhx + qhy*qhy + qhz*qhz )
			    - 2 * qhw * ( phx*qhx + phy*qhy + phz*qhz )
			    )
	- qhw*qhw * (         phw * ( rhx*rhx + rhy*rhy + rhz*rhz )
			      - 2 * rhw * ( phx*rhx + phy*rhy + phz*rhz )
			      );

      return enum_cast<Comparison_result>(CGAL_NTS sign(dosd));
    }
  };

  template <typename K>
  class Compare_squared_distance_3
  {
    typedef typename K::Point_3            Point_3;
    typedef typename K::FT                 FT;
  public:
    typedef typename K::Comparison_result  result_type;
    typedef Arity_tag< 3 >                 Arity;

    result_type
    operator()(const Point_3& p, const Point_3& q, const FT& d2) const
    {
      return CGAL_NTS compare(squared_distance(p, q), d2);
    }
  };

  template <typename K>
  class Compare_slope_2
  {
    typedef typename K::Line_2             Line_2;
    typedef typename K::Segment_2          Segment_2;
  public:
    typedef typename K::Comparison_result  result_type;
    typedef Arity_tag< 2 >                 Arity;

    result_type
    operator()(const Line_2& l1, const Line_2& l2) const
    {
      typedef typename K::RT RT;
      if (l1.is_horizontal())
	return l2.is_vertical() ?
	  SMALLER : enum_cast<Comparison_result>(CGAL_NTS sign(l2.a() * l2.b()));
      if (l2.is_horizontal())
	return l1.is_vertical() ?
	  LARGER : enum_cast<Comparison_result>(- CGAL_NTS sign(l1.a() * l1.b()));
      if (l1.is_vertical()) return l2.is_vertical() ? EQUAL : LARGER;
      if (l2.is_vertical()) return SMALLER;
      int l1_sign = CGAL_NTS sign(-l1.a() * l1.b());
      int l2_sign = CGAL_NTS sign(-l2.a() * l2.b());

      if (l1_sign < l2_sign) return SMALLER;
      if (l1_sign > l2_sign) return LARGER;

      if (l1_sign > 0)
	return CGAL_NTS compare( CGAL_NTS abs(l1.a() * l2.b()),
				 CGAL_NTS abs(l2.a() * l1.b()) );

      return CGAL_NTS compare( CGAL_NTS abs(l2.a() * l1.b()),
			       CGAL_NTS abs(l1.a() * l2.b()) );
    } // FIXME

    result_type
    operator()(const Segment_2& s1, const Segment_2& s2) const
    {
      typedef typename K::FT        FT;

      typename K::Comparison_result cmp_y1 = compare_y(s1.source(), s1.target());
      if (cmp_y1 == EQUAL) // horizontal
	{
	  typename K::Comparison_result cmp_x2 = compare_x(s2.source(), s2.target());

	  if (cmp_x2 == EQUAL) return SMALLER;
	  FT s_hw = s2.source().hw();
	  FT t_hw = s2.target().hw();
	  return enum_cast<Comparison_result>(
	    - CGAL_NTS sign((s2.source().hy()*t_hw - s2.target().hy()*s_hw) *
			    (s2.source().hx()*t_hw - s2.target().hx()*s_hw)) );
	}

      typename K::Comparison_result cmp_y2 = compare_y(s2.source(), s2.target());
      if (cmp_y2 == EQUAL)
	{
	  typename K::Comparison_result cmp_x1 = compare_x(s1.source(), s1.target());

	  if (cmp_x1 == EQUAL) return LARGER;
	  FT s_hw = s1.source().hw();
	  FT t_hw = s1.target().hw();
	  return enum_cast<Comparison_result>(
	    CGAL_NTS sign((s1.source().hy()*t_hw - s1.target().hy()*s_hw) *
			  (s1.source().hx()*t_hw - s1.target().hx()*s_hw)) );
	}

      typename K::Comparison_result cmp_x1 = compare_x(s1.source(), s1.target());
      typename K::Comparison_result cmp_x2 = compare_x(s2.source(), s2.target());
      if (cmp_x1 == EQUAL)
	return cmp_x2 == EQUAL ? EQUAL : LARGER;

      if (cmp_x2 == EQUAL) return SMALLER;

      FT s1_s_hw = s1.source().hw();
      FT s1_t_hw = s1.target().hw();
      FT s2_s_hw = s2.source().hw();
      FT s2_t_hw = s2.target().hw();
      FT s1_xdiff = s1.source().hx()*s1_t_hw - s1.target().hx()*s1_s_hw;
      FT s1_ydiff = s1.source().hy()*s1_t_hw - s1.target().hy()*s1_s_hw;
      FT s2_xdiff = s2.source().hx()*s2_t_hw - s2.target().hx()*s2_s_hw;
      FT s2_ydiff = s2.source().hy()*s2_t_hw - s2.target().hy()*s2_s_hw;
      typename K::Sign s1_sign = CGAL_NTS sign(s1_ydiff * s1_xdiff);
      typename K::Sign s2_sign = CGAL_NTS sign(s2_ydiff * s2_xdiff);

      if (s1_sign < s2_sign) return SMALLER;
      if (s1_sign > s2_sign) return LARGER;

      if (s1_sign > 0)
	return enum_cast<Comparison_result>(
	 CGAL_NTS sign ( CGAL_NTS abs(s1_ydiff * s2_xdiff) -
			 CGAL_NTS abs(s2_ydiff * s1_xdiff)) );

      return enum_cast<Comparison_result>(
       CGAL_NTS sign ( CGAL_NTS abs(s2_ydiff * s1_xdiff) -
		       CGAL_NTS abs(s1_ydiff * s2_xdiff)) );
    }
  };

  template <typename K>
  class Compare_x_at_y_2
  {
    typedef typename K::Point_2            Point_2;
    typedef typename K::Line_2             Line_2;
  public:
    typedef typename K::Comparison_result  result_type;
    typedef Arity_tag< 3 >                 Arity;

    result_type
    operator()( const Point_2& p, const Line_2& h) const
    {
      typedef typename K::RT RT;
      CGAL_kernel_precondition( ! h.is_horizontal() );
      typename K::Oriented_side ors = h.oriented_side( p );
      if ( h.a() < RT(0) )
	  ors = opposite( ors );
      if ( ors == ON_POSITIVE_SIDE )
	  return LARGER;
      return ( ors == ON_NEGATIVE_SIDE ) ? SMALLER : EQUAL;
    } // FIXME

    result_type
    operator()( const Point_2& p, const Line_2& h1, const Line_2& h2) const
    { return CGAL_NTS compare(h1.x_at_y( p.y() ), h2.x_at_y( p.y() )); }
    // FIXME

    result_type
    operator()( const Line_2& l1, const Line_2& l2, const Line_2& h) const
    { return compare_x_at_y( gp_linear_intersection( l1, l2 ), h); }
    // FIXME

    result_type
    operator()( const Line_2& l1, const Line_2& l2,
	        const Line_2& h1, const Line_2& h2) const
    { return compare_x_at_y( gp_linear_intersection( l1, l2 ), h1, h2 ); }
    // FIXME
  };

  template <typename K>
  class Compare_xyz_3
  {
    typedef typename K::Point_3            Point_3;
  public:
    typedef typename K::Comparison_result  result_type;
    typedef Arity_tag< 2 >                 Arity;

    result_type
    operator()( const Point_3& p, const Point_3& q) const
    {
      typedef typename K::RT RT;
      RT pV = p.hx()*q.hw();
      RT qV = q.hx()*p.hw();
      if ( pV < qV )
	{
	  return SMALLER;
	}
      if ( qV < pV )    //   ( pV > qV )
	{
	  return LARGER;
	}
      // same x
      pV = p.hy()*q.hw();
      qV = q.hy()*p.hw();
      if ( pV < qV )
	{
	  return SMALLER;
	}
      if ( qV < pV )    //   ( pV > qV )
	{
	  return LARGER;
	}
      // same x and y
      pV = p.hz()*q.hw();
      qV = q.hz()*p.hw();
      return CGAL_NTS compare(pV, qV);
    }
  };

  template <typename K>
  class Compare_xy_2
  {
    typedef typename K::Point_2            Point_2;
  public:
    typedef typename K::Comparison_result  result_type;
    typedef Arity_tag< 2 >                 Arity;

    result_type
    operator()( const Point_2& p, const Point_2& q) const
    {
      typedef typename K::RT RT;

      const RT& phx = p.hx();
      const RT& phy = p.hy();
      const RT& phw = p.hw();
      const RT& qhx = q.hx();
      const RT& qhy = q.hy();
      const RT& qhw = q.hw();

      RT pV = phx*qhw;
      RT qV = qhx*phw;
      if ( pV == qV )