function_objects.h

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

    typedef Arity_tag< 2 >   Arity;

    Construct_base_vector_3() {}
    Construct_base_vector_3(const Construct_orthogonal_vector_3& co_)
      : co(co_)
    {}

    Vector_3
    operator()( const Plane_3& h, int index ) const
    {
      if (index == 1) {
	// point():
	// a() != RT0 : Point_3( -d(), RT0, RT0, a() );
	// b() != RT0 : Point_3( RT0, -d(), RT0, b() );
	//            : Point_3( RT0, RT0, -d(), c() );
	// point1():
	// a() != RT0 : Point_3( -b()-d(), a(), RT0, a() );
	// b() != RT0 : Point_3( RT0, -c()-d(), b(), b() );
	//            : Point_3( c(), RT0, -a()-d(), c() );
	
	const RT RT0(0);
	if ( h.a() != RT0 )
	  {
	    return Vector_3( -h.b(), h.a(), RT0, h.a() );
	  }
	if ( h.b() != RT0 )
	  {
	    return Vector_3( RT0, -h.c(), h.b(), h.b() );
	  }
	CGAL_kernel_assertion ( h.c() != RT(0) );
	return Vector_3( h.c(), RT0, -h.a(), h.c() );
      } else {
	Vector_3 a = co(h);
	Vector_3 b = this->operator()(h, 1);
	return Vector_3(a.hy()*b.hz() - a.hz()*b.hy(),
			a.hz()*b.hx() - a.hx()*b.hz(),
			a.hx()*b.hy() - a.hy()*b.hx(),
			a.hw()*b.hw() );
      }
    }
  };

  template <typename K>
  class Construct_bbox_2
  {
    typedef typename K::Point_2          Point_2;
    typedef typename K::Segment_2        Segment_2;
    typedef typename K::Iso_rectangle_2  Iso_rectangle_2;
    typedef typename K::Triangle_2       Triangle_2;
    typedef typename K::Circle_2         Circle_2;
  public:
    typedef Bbox_2           result_type;
    typedef Arity_tag< 1 >   Arity;

    Bbox_2
    operator()( const Point_2& p) const
    {
      Interval_nt<> ihx = CGAL_NTS to_interval(p.hx());
      Interval_nt<> ihy = CGAL_NTS to_interval(p.hy());
      Interval_nt<> ihw = CGAL_NTS to_interval(p.hw());

      Interval_nt<> ix = ihx/ihw;
      Interval_nt<> iy = ihy/ihw;

      return Bbox_2(ix.inf(), iy.inf(), ix.sup(), iy.sup());
    }

    Bbox_2
    operator()( const Segment_2& s) const
    { return s.source().bbox() + s.target().bbox(); }

    Bbox_2
    operator()( const Triangle_2& t) const
    {
      typename K::Construct_bbox_2 construct_bbox_2;
      return construct_bbox_2(t.vertex(0))
	+ construct_bbox_2(t.vertex(1))
	+ construct_bbox_2(t.vertex(2));
    }

    Bbox_2
    operator()( const Iso_rectangle_2& r) const
    {
      typename K::Construct_bbox_2 construct_bbox_2;
      return construct_bbox_2((r.min)()) + construct_bbox_2((r.max)());
    }

    Bbox_2
    operator()( const Circle_2& c) const
    {
      typename K::Construct_bbox_2 construct_bbox_2;
      Bbox_2 b = construct_bbox_2(c.center());

      Interval_nt<> x (b.xmin(), b.xmax());
      Interval_nt<> y (b.ymin(), b.ymax());

      Interval_nt<> sqr = CGAL_NTS to_interval(c.squared_radius());
      Interval_nt<> r = CGAL::sqrt(sqr);
      Interval_nt<> minx = x-r;
      Interval_nt<> maxx = x+r;
      Interval_nt<> miny = y-r;
      Interval_nt<> maxy = y+r;

      return Bbox_2(minx.inf(), miny.inf(), maxx.sup(), maxy.sup()); }
  };


  template <typename K>
  class Construct_bbox_3
  {
    typedef typename K::Point_3          Point_3;
    typedef typename K::Segment_3        Segment_3;
    typedef typename K::Triangle_3       Triangle_3;
    typedef typename K::Tetrahedron_3    Tetrahedron_3;
    typedef typename K::Iso_cuboid_3     Iso_cuboid_3;
    typedef typename K::Sphere_3         Sphere_3;
  public:
    typedef Bbox_3           result_type;
    typedef Arity_tag< 1 >   Arity;

    Bbox_3
    operator()(const Point_3& p) const
    {
       Interval_nt<> ihx = CGAL_NTS to_interval(p.hx());
       Interval_nt<> ihy = CGAL_NTS to_interval(p.hy());
       Interval_nt<> ihz = CGAL_NTS to_interval(p.hz());
       Interval_nt<> ihw = CGAL_NTS to_interval(p.hw());

       Interval_nt<> ix = ihx/ihw;
       Interval_nt<> iy = ihy/ihw;
       Interval_nt<> iz = ihz/ihw;

       return Bbox_3(ix.inf(), iy.inf(), iz.inf(),
                     ix.sup(), iy.sup(), iz.sup());
    }

    Bbox_3
    operator()(const Segment_3& s) const
    { return s.source().bbox() + s.target().bbox(); }

    Bbox_3
    operator()(const Triangle_3& t) const
    {
      typename K::Construct_bbox_3 construct_bbox;
      return construct_bbox(t.vertex(0))
	   + construct_bbox(t.vertex(1))
	   + construct_bbox(t.vertex(2));
    }

    Bbox_3
    operator()(const Iso_cuboid_3& r) const
    {
      typename K::Construct_bbox_3 construct_bbox;
      return construct_bbox((r.min)()) + construct_bbox((r.max)());
    }

    Bbox_3
    operator()(const Tetrahedron_3& t) const
    {
      typename K::Construct_bbox_3 construct_bbox_3;
      return construct_bbox_3(t.vertex(0)) + construct_bbox_3(t.vertex(1))
           + construct_bbox_3(t.vertex(2)) + construct_bbox_3(t.vertex(3));
    }

    Bbox_3
    operator()(const Sphere_3& s) const
    {
      Bbox_3 b = s.center().bbox();

      Interval_nt<> x (b.xmin(), b.xmax());
      Interval_nt<> y (b.ymin(), b.ymax());
      Interval_nt<> z (b.zmin(), b.zmax());

      Interval_nt<> sqr = CGAL_NTS to_interval(s.squared_radius());
      Interval_nt<> r = CGAL::sqrt(sqr);
      Interval_nt<> minx = x-r;
      Interval_nt<> maxx = x+r;
      Interval_nt<> miny = y-r;
      Interval_nt<> maxy = y+r;
      Interval_nt<> minz = z-r;
      Interval_nt<> maxz = z+r;

      return Bbox_3(minx.inf(), miny.inf(), minz.inf(),
		    maxx.sup(), maxy.sup(), maxz.sup());
    }
  };


  template <typename K>
  class Construct_bisector_2
  {
    typedef typename K::RT      RT;
    typedef typename K::FT      FT;
    typedef typename K::Point_2 Point_2;
    typedef typename K::Line_2  Line_2;
  public:
    typedef Line_2           result_type;
    typedef Arity_tag< 2 >   Arity;

    Line_2
    operator()(const Point_2& p, const Point_2& q) const
    {
      // Bisector equation is based on equation
      // ( X - p.x())^2 + (Y - p.y())^2 == ( X - q.x())^2 + (Y - q.y())
      // and x() = hx()/hw() ...

      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 a = RT(2) * ( phx*phw*qhw*qhw - qhx*qhw*phw*phw );
      RT b = RT(2) * ( phy*phw*qhw*qhw - qhy*qhw*phw*phw );
      RT c = qhx*qhx*phw*phw + qhy*qhy*phw*phw
	   - phx*phx*qhw*qhw - phy*phy*qhw*qhw;

      return Line_2( a, b, c );
    }

    Line_2
    operator()(const Line_2& p, const Line_2& q) const
    {
      RT a, b, c;
      bisector_of_linesC2(p.a(), p.b(), p.c(),
                          q.a(), q.b(), q.c(),
                          a, b, c);
      return Line_2(a, b, c);
    }
  };

  template <typename K>
  class Construct_bisector_3
  {
    typedef typename K::RT      RT;
    typedef typename K::FT      FT;
    typedef typename K::Point_3 Point_3;
    typedef typename K::Plane_3 Plane_3;
  public:
    typedef Plane_3          result_type;
    typedef Arity_tag< 2 >   Arity;

    Plane_3
    operator()(const Point_3& p, const Point_3& q) const
    {
      // Bisector equation is based on equation
      // ( X - p.x())^2 + (Y - p.y())^2 == ( X - q.x())^2 + (Y - q.y())
      // and x() = hx()/hw() ...

      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();

      RT a = RT(2) * ( phx*phw*qhw*qhw - qhx*qhw*phw*phw );
      RT b = RT(2) * ( phy*phw*qhw*qhw - qhy*qhw*phw*phw );
      RT c = RT(2) * ( phz*phw*qhw*qhw - qhz*qhw*phw*phw );
      RT d = qhx*qhx*phw*phw + qhy*qhy*phw*phw + qhz*qhz*phw*phw
	   - phx*phx*qhw*qhw - phy*phy*qhw*qhw - phz*phz*qhw*qhw;

      return Plane_3( a, b, c, d );
    }

    Plane_3
    operator()(const Plane_3& p, const Plane_3& q) const
    {
      RT a, b, c, d;
      bisector_of_planesC3(p.a(), p.b(), p.c(), p.d(),
	                   q.a(), q.b(), q.c(), q.d(),
			   a, b, c, d);
      return Plane_3(a, b, c, d);
    }
  };


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

    Point_2
    operator()(const Point_2& p, const Point_2& q, const Point_2& r) const
    {
      typedef typename K::RT  RT;
      const RT phw(p.hw());
      const RT qhw(q.hw());
      const RT rhw(r.hw());
      RT hx(p.hx()*qhw*rhw + q.hx()*phw*rhw + r.hx()*phw*qhw);
      RT hy(p.hy()*qhw*rhw + q.hy()*phw*rhw + r.hy()*phw*qhw);
      RT hw( phw*qhw*rhw * 3);
      return Point_2(hx, hy, hw);
    }

    Point_2
    operator()(const Triangle_2& t) const
    {
      return this->operator()(t.vertex(0), t.vertex(1), t.vertex(2));
    }

    Point_2
    operator()(const Point_2& p, const Point_2& q,
               const Point_2& r, const Point_2& s) const
    {
      typedef typename K::RT  RT;
      const RT phw(p.hw());
      const RT qhw(q.hw());
      const RT rhw(r.hw());
      const RT shw(s.hw());
      RT hx(p.hx()*qhw*rhw*shw + q.hx()*phw*rhw*shw + r.hx()*phw*qhw*shw
	    + s.hx()*phw*qhw*rhw);
      RT hy(p.hy()*qhw*rhw*shw + q.hy()*phw*rhw*shw + r.hy()*phw*qhw*shw
	    + s.hy()*phw*qhw*rhw);
      RT hw( phw*qhw*rhw*shw * 4);
      return Point_2(hx, hy, hw);
    }
  };

  template <typename K>
  class Construct_centroid_3
  {
    typedef typename K::RT             RT;
    typedef typename K::Point_3        Point_3;
    typedef typename K::Triangle_3     Triangle_3;
    typedef typename K::Tetrahedron_3  Tetrahedron_3;
  public:
    typedef Point_3          result_type;
    typedef Arity_tag< 3 >   Arity;

    Point_3
    operator()(const Point_3& p, const Point_3& q, const Point_3& r) const
    {
      const RT& phw = p.hw();
      const RT& qhw = q.hw();
      const RT& rhw = r.hw();
      RT hx(p.hx()*qhw*rhw + q.hx()*phw*rhw + r.hx()*phw*qhw);
      RT hy(p.hy()*qhw*rhw + q.hy()*phw*rhw + r.hy()*phw*qhw);
      RT hz(p.hz()*qhw*rhw + q.hz()*phw*rhw + r.hz()*phw*qhw);
      RT hw( phw*qhw*rhw * RT(3));
      return Point_3(hx, hy, hz, hw);
    }

    Point_3
    operator()(const Point_3& p, const Point_3& q,
               const Point_3& r, const Point_3& s) const
    {
      const RT& phw = p.hw();
      const RT& qhw = q.hw();
      const RT& rhw = r.hw();
      const RT& shw = s.hw();
      RT hx(p.hx()*qhw*rhw*shw + q.hx()*phw*rhw*shw + r.hx()*phw*qhw*shw
	    + s.hx()*phw*qhw*rhw);
      RT hy(p.hy()*qhw*rhw*shw + q.hy()*phw*rhw*shw + r.hy()*phw*qhw*shw
	    + s.hy()*phw*qhw*rhw);
      RT hz(p.hz()*qhw*rhw*shw + q.hz()*phw*rhw*shw + r.hz()*phw*qhw*shw
	    + s.hz()*phw*qhw*rhw);
      RT hw( phw*qhw*rhw*shw * RT(4));
      return Point_3(hx, hy, hz, hw);
    }

    Point_3
    operator()(const Triangle_3& t) const
    {
      return this->operator()(t.vertex(0), t.vertex(1), t.vertex(2));
    }

    Point_3
    operator()(const Tetrahedron_3& t) const
    {
      return this->operator()(t.vertex(0), t.vertex(1),
                              t.vertex(2), t.vertex(3));
    }
  };

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

    Point_2
    operator()(const Point_2& p, const Point_2& q) const
    {
      typename K::Construct_midpoint_2 construct_midpoint_2;
      return construct_midpoint_2(p, q);
    }

    Point_2
    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();

#ifdef CGAL_EXPANDED_CIRCUMCENTER_COMPUTATION
      RT vvx =
	( qhy*qhw*phw*phw - phy*phw*qhw*qhw )
	*( phx*phx*rhw*rhw + phy*phy*rhw*rhw -
	   rhx*rhx*phw*phw - rhy*rhy*phw*phw )
	-  ( rhy*rhw*phw*phw - phy*phw*rhw*rhw )
	*( phx*phx*qhw*qhw + phy*phy*qhw*qhw -
	   qhx*qhx*phw*phw - qhy*qhy*phw*phw );

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

      RT vvw = RT(2) *
	(  ( qhx*qhw*phw*phw - phx*phw*qhw*qhw )
	   *( rhy*rhw*phw*phw - phy*phw*rhw*rhw )
	   -  ( rhx*rhw*phw*phw - phx*phw*rhw*rhw )
	   *( qhy*qhw*phw*phw - phy*phw*qhw*qhw ) );
#endif // CGAL_EXPANDED_CIRCUMCENTER_COMPUTATION

      RT qy_py = ( qhy*qhw*phw*phw - phy*phw*qhw*qhw );
      RT qx_px = ( qhx*qhw*phw*phw - phx*phw*qhw*qhw );
      RT rx_px = ( rhx*rhw*phw*phw - phx*phw*rhw*rhw );
      RT ry_py = ( rhy*rhw*phw*phw - phy*phw*rhw*rhw );

      RT px2_py2_rx2_ry_2 =
	phx*phx*rhw*rhw + phy*phy*rhw*rhw -
	rhx*rhx*phw*phw - rhy*rhy*phw*phw ;
      RT px2_py2_qx2_qy_2 =
	phx*phx*qhw*qhw + phy*phy*qhw*qhw -
	qhx*qhx*phw*phw - qhy*qhy*phw*phw ;

      RT vvx = qy_py * px2_py2_rx2_ry_2 - ry_py * px2_py2_qx2_qy_2;
      RT vvy = rx_px * px2_py2_qx2_qy_2 - qx_px * px2_py2_rx2_ry_2;
      RT vvw = RT(2) * ( qx_px * ry_py - rx_px * qy_py );

      return Point_2( vvx, vvy, vvw );
    }

    Point_2
    operator()(const Triangle_2&