function_objects.h

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

      // FT num_x = det3x3_by_formula(psy,psz,ps2,
      //                              qsy,qsz,qs2,
      //                              rsy,rsz,0);
      // FT num_y = det3x3_by_formula(psx,psz,ps2,
      //                              qsx,qsz,qs2,
      //                              rsx,rsz,0);
      // FT num_z = det3x3_by_formula(psx,psy,ps2,
      //                              qsx,qsy,qs2,
      //                              rsx,rsy,0);

      FT num_x = ps2 * det2x2_by_formula(qsy,qsz,rsy,rsz)
	       - qs2 * det2x2_by_formula(psy,psz,rsy,rsz);
      FT num_y = ps2 * det2x2_by_formula(qsx,qsz,rsx,rsz)
	       - qs2 * det2x2_by_formula(psx,psz,rsx,rsz);
      FT num_z = ps2 * det2x2_by_formula(qsx,qsy,rsx,rsy)
	       - qs2 * det2x2_by_formula(psx,psy,rsx,rsy);

      FT den   = det3x3_by_formula(psx,psy,psz,
                                   qsx,qsy,qsz,
                                   rsx,rsy,rsz);

      CGAL_kernel_assertion( den != 0 );
      FT inv = 1 / (2 * den);

      FT x = s.x() + num_x*inv;
      FT y = s.y() - num_y*inv;
      FT z = s.z() + num_z*inv;
      return construct_point_3(x, y, z);
    }

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

    Point_3
    operator()(const Point_3& p, const Point_3& q,
	       const Point_3& r, const Point_3& s) const
    {
      typename K::Construct_point_3 construct_point_3;
      // Translate p to origin to simplify the expression.
      FT qpx = q.x()-p.x();
      FT qpy = q.y()-p.y();
      FT qpz = q.z()-p.z();
      FT qp2 = CGAL_NTS square(qpx) + CGAL_NTS square(qpy) + CGAL_NTS square(qpz);
      FT rpx = r.x()-p.x();
      FT rpy = r.y()-p.y();
      FT rpz = r.z()-p.z();
      FT rp2 = CGAL_NTS square(rpx) + CGAL_NTS square(rpy) + CGAL_NTS square(rpz);
      FT spx = s.x()-p.x();
      FT spy = s.y()-p.y();
      FT spz = s.z()-p.z();
      FT sp2 = CGAL_NTS square(spx) + CGAL_NTS square(spy) + CGAL_NTS square(spz);

      FT num_x = det3x3_by_formula(qpy,qpz,qp2,
                                   rpy,rpz,rp2,
                                   spy,spz,sp2);
      FT num_y = det3x3_by_formula(qpx,qpz,qp2,
                                   rpx,rpz,rp2,
                                   spx,spz,sp2);
      FT num_z = det3x3_by_formula(qpx,qpy,qp2,
                                   rpx,rpy,rp2,
                                   spx,spy,sp2);
      FT den   = det3x3_by_formula(qpx,qpy,qpz,
                                   rpx,rpy,rpz,
                                   spx,spy,spz);
      CGAL_kernel_assertion( ! CGAL_NTS is_zero(den) );
      FT inv = 1 / (2 * den);

      FT x = p.x() + num_x*inv;
      FT y = p.y() - num_y*inv;
      FT z = p.z() + num_z*inv;
      return construct_point_3(x, y, z);
    }

    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_cross_product_vector_3
  {
    typedef typename K::Vector_3  Vector_3;
  public:
    typedef Vector_3              result_type;
    typedef Arity_tag< 2 >        Arity;

    Vector_3
    operator()(const Vector_3& v, const Vector_3& w) const
    {
      return Vector_3(v.y() * w.z() - v.z() * w.y(),
		      v.z() * w.x() - v.x() * w.z(),
		      v.x() * w.y() - v.y() * w.x());
    }
  };

  template <typename K>
  class Construct_lifted_point_3
  {
    typedef typename K::Point_2                    Point_2;
    typedef typename K::Point_3                    Point_3;
    typedef typename K::Plane_3                    Plane_3;
    typedef typename K::Construct_base_vector_3    Construct_base_vector_3;
    typedef typename K::Construct_point_on_3       Construct_point_on_3;
    typedef typename K::Construct_scaled_vector_3  Construct_scaled_vector_3;
    typedef typename K::Construct_translated_point_3
    Construct_translated_point_3;
    Construct_base_vector_3 cb;
    Construct_point_on_3 cp;
    Construct_scaled_vector_3 cs;
    Construct_translated_point_3 ct;
  public:
    typedef Point_3          result_type;
    typedef Arity_tag< 2 >   Arity;

    Construct_lifted_point_3() {}
    Construct_lifted_point_3(const Construct_base_vector_3& cb_,
			     const Construct_point_on_3& cp_,
			     const Construct_scaled_vector_3& cs_,
			     const Construct_translated_point_3& ct_)
      : cb(cb_), cp(cp_), cs(cs_), ct(ct_)
    {}

    Point_3
    operator()(const Plane_3& h, const Point_2& p) const
    {
      return ct(ct(cp(h), cs(cb(h,1), p.x())), cs(cb(h,2), p.y()));
    }
  };

  template <typename K>
  class Construct_direction_2
  {
    typedef typename K::Direction_2     Direction_2;
    typedef typename Direction_2::Rep   Rep;
    typedef typename K::Point_2         Point_2;
    typedef typename K::Vector_2        Vector_2;
    typedef typename K::Line_2          Line_2;
    typedef typename K::Ray_2           Ray_2;
    typedef typename K::Segment_2       Segment_2;
    typedef typename K::RT              RT;

  public:
    typedef Direction_2                 result_type;
    typedef Arity_tag< 1 >              Arity;

    Rep // Direction_2
    operator()(Return_base_tag, const RT& x, const RT& y) const
    { return Rep(x, y); }

    Rep // Direction_2
    operator()(Return_base_tag, const Vector_2& v) const
    {
      return Rep(v.x(),v.y()); }

    Rep // Direction_2
    operator()(Return_base_tag, const Line_2& l) const
    { return Rep(l.b(), -l.a()); }

    Rep // Direction_2
    operator()(Return_base_tag, const Point_2& p, const Point_2& q) const
    {
      return Rep(q.x() - p.x(), q.y() - p.y());
    }

    Rep // Direction_2
    operator()(Return_base_tag, const Ray_2& r) const
    {
      return this->operator()(Return_base_tag(), r.source(), r.second_point());
    }

    Rep // Direction_2
    operator()(Return_base_tag, const Segment_2& s) const
    {
      return this->operator()(Return_base_tag(), s.source(), s.target());
    }


    Direction_2
    operator()(const RT& x, const RT& y) const
    { return this->operator()(Return_base_tag(), x, y); }

    Direction_2
    operator()(const Vector_2& v) const
    {
      return this->operator()(Return_base_tag(), v); }

    Direction_2
    operator()(const Line_2& l) const
    { return this->operator()(Return_base_tag(), l); }

    Direction_2
    operator()(const Point_2& p, const Point_2& q) const
    {
      return this->operator()(Return_base_tag(), p, q);
    }

    Direction_2
    operator()(const Ray_2& r) const
    {
      return this->operator()(Return_base_tag(), r);
    }

    Direction_2
    operator()(const Segment_2& s) const
    {
      return this->operator()(Return_base_tag(), s);
    }
  };

  template <typename K>
  class Construct_direction_3
  {
    typedef typename K::Direction_3     Direction_3;
    typedef typename K::Vector_3        Vector_3;
    typedef typename K::Line_3          Line_3;
    typedef typename K::Ray_3           Ray_3;
    typedef typename K::Segment_3       Segment_3;
    typedef typename K::RT              RT;
    typedef typename Direction_3::Rep   Rep;
  public:
    typedef Direction_3       result_type;
    typedef Arity_tag< 1 >    Arity;

    Rep // Direction_3
    operator()(Return_base_tag, const RT& x, const RT& y, const RT& z) const
    { return Rep(x, y, z); }

    Rep // Direction_3
    operator()(Return_base_tag, const Vector_3& v) const
    { return Rep(v); }

    Rep // Direction_3
    operator()(Return_base_tag, const Line_3& l) const
    { return Rep(l); }

    Rep // Direction_3
    operator()(Return_base_tag, const Ray_3& r) const
    { return Rep(r); }

    Rep // Direction_3
    operator()(Return_base_tag, const Segment_3& s) const
    { return Rep(s); }


    Direction_3
    operator()(const RT& x, const RT& y, const RT& z) const
    { return this->operator()(Return_base_tag(), x, y, z); }

    Direction_3
    operator()(const Vector_3& v) const
    { return this->operator()(Return_base_tag(), v); }

    Direction_3
    operator()(const Line_3& l) const
    { return this->operator()(Return_base_tag(), l); }

    Direction_3
    operator()(const Ray_3& r) const
    { return this->operator()(Return_base_tag(), r); }

    Direction_3
    operator()(const Segment_3& s) const
    { return this->operator()(Return_base_tag(), s); }
  };

  template <typename K>
  class Construct_equidistant_line_3
  {
    typedef typename K::FT          FT;
    typedef typename K::Point_3     Point_3;
    typedef typename K::Vector_3    Vector_3;
    typedef typename K::Line_3      Line_3;
    typedef typename Line_3::Rep    Rep;
  public:
    typedef Line_3           result_type;
    typedef Arity_tag< 3 >   Arity;

    Line_3
    operator()( const Point_3& p, const Point_3& q, const Point_3& s) const
    {
      CGAL_kernel_precondition(! collinear(p, q, s));

      // Translate s to origin to simplify the expression.
      FT psx = p.x()-s.x();
      FT psy = p.y()-s.y();
      FT psz = p.z()-s.z();
      FT ps2 = CGAL_NTS square(psx) + CGAL_NTS square(psy) + CGAL_NTS square(psz);
      FT qsx = q.x()-s.x();
      FT qsy = q.y()-s.y();
      FT qsz = q.z()-s.z();
      FT qs2 = CGAL_NTS square(qsx) + CGAL_NTS square(qsy) + CGAL_NTS square(qsz);
      FT rsx = psy*qsz-psz*qsy;
      FT rsy = psz*qsx-psx*qsz;
      FT rsz = psx*qsy-psy*qsx;

      // The following determinants can be developped and simplified.
      //
      // FT num_x = det3x3_by_formula(psy,psz,ps2,
      //                              qsy,qsz,qs2,
      //                              rsy,rsz,0);
      // FT num_y = det3x3_by_formula(psx,psz,ps2,
      //                              qsx,qsz,qs2,
      //                              rsx,rsz,0);
      // FT num_z = det3x3_by_formula(psx,psy,ps2,
      //                              qsx,qsy,qs2,
      //                              rsx,rsy,0);

      FT num_x = ps2 * det2x2_by_formula(qsy,qsz,rsy,rsz)
	       - qs2 * det2x2_by_formula(psy,psz,rsy,rsz);
      FT num_y = ps2 * det2x2_by_formula(qsx,qsz,rsx,rsz)
	       - qs2 * det2x2_by_formula(psx,psz,rsx,rsz);
      FT num_z = ps2 * det2x2_by_formula(qsx,qsy,rsx,rsy)
	       - qs2 * det2x2_by_formula(psx,psy,rsx,rsy);

      FT den   = det3x3_by_formula(psx,psy,psz,
                                   qsx,qsy,qsz,
                                   rsx,rsy,rsz);

      CGAL_kernel_assertion( den != 0 );
      FT inv = 1 / (2 * den);

      FT x = s.x() + num_x*inv;
      FT y = s.y() - num_y*inv;
      FT z = s.z() + num_z*inv;
      return Rep(Point_3(x, y, z), Vector_3(rsx, rsy, rsz));
    }

  };

  template <typename K>
  class Construct_iso_rectangle_2
  {
    typedef typename K::RT               RT;
    typedef typename K::FT               FT;
    typedef typename K::Point_2          Point_2;
    typedef typename K::Iso_rectangle_2  Iso_rectangle_2;
    typedef typename Iso_rectangle_2::Rep     Rep;

  public:
    typedef Iso_rectangle_2              result_type;
    typedef Arity_tag< 2 >               Arity;

    Rep // Iso_rectangle_2
    operator()(Return_base_tag, const Point_2& p, const Point_2& q, int) const
    {
      // I have to remove the assertions, because of Cartesian_converter.
      // CGAL_kernel_assertion(p.x()<=q.x());
      // CGAL_kernel_assertion(p.y()<=q.y());
      return Rep(p, q, 0);
    }

    Rep // Iso_rectangle_2
    operator()(Return_base_tag, const Point_2& p, const Point_2& q) const
    {
      FT minx, maxx, miny, maxy;
      if (p.x() < q.x()) { minx = p.x(); maxx = q.x(); }
      else               { minx = q.x(); maxx = p.x(); }
      if (p.y() < q.y()) { miny = p.y(); maxy = q.y(); }
      else               { miny = q.y(); maxy = p.y(); }

      return Rep(Point_2(minx, miny),
	         Point_2(maxx, maxy), 0);
    }

    Rep // Iso_rectangle_2
    operator()(Return_base_tag, const Point_2 &left,   const Point_2 &right,
               const Point_2 &bottom, const Point_2 &top) const
    {
      CGAL_kernel_assertion_code(typename K::Less_x_2 less_x;)
      CGAL_kernel_assertion_code(typename K::Less_y_2 less_y;)
      CGAL_kernel_assertion(!less_x(right, left));
      CGAL_kernel_assertion(!less_y(top, bottom));
      return Rep(Point_2(left.x(), bottom.y()),
		 Point_2(right.x(), top.y()), 0);
    }

    Rep // Iso_rectangle_2
    operator()(Return_base_tag, const RT& min_hx, const RT& min_hy,
	       const RT& max_hx, const RT& max_hy) const
    {
      CGAL_kernel_precondition(min_hx <= max_hx);
      CGAL_kernel_precondition(min_hy <= max_hy);
      return Rep(Point_2(min_hx, min_hy),
		 Point_2(max_hx, max_hy), 0);
    }

    Rep // Iso_rectangle_2
    operator()(Return_base_tag, const RT& min_hx, const RT& min_hy,
	       const RT& max_hx, const RT& max_hy, const RT& hw) const
    {
      if (hw == 1)
	return Rep(Point_2(min_hx, min_hy),
		   Point_2(max_hx, max_hy), 0);
      return Rep(Point_2(min_hx/hw, min_hy/hw),
		 Point_2(max_hx/hw, max_hy/hw), 0);
    }


    Iso_rectangle_2
    operator()(const Point_2& p, const Point_2& q, int i) const
    {
      return this->operator()(Return_base_tag(), p, q, i);
    }

    Iso_rectangle_2
    operator()(const Point_2& p, const Point_2& q) const
    {
      return this->operator()(Return_base_tag(), p, q);
    }

    Iso_rectangle_2
    operator()(const Point_2 &left,   const Point_2 &right,
               const Point_2 &bottom, const Point_2 &top) const
    {
      return this->operator()(Return_base_tag(), left, right, bottom, top);
    }

    Iso_rectangle_2
    operator()(const RT& min_hx, const RT& min_hy,
	       const RT& max_hx, const RT& max_hy) const
    {
      return this->operator()(Return_base_tag(), min_hx, min_hy, max_hx, max_hy);
    }

    Iso_rectangle_2
    operator()(const RT& min_hx, const RT& min_hy,
	       const RT& max_hx, const RT& max_hy, const RT& hw) const
    {
      return this->operator()(Return_base_tag(), min_hx, min_hy, max_hx, max_hy, hw);
    }
  };

  template <typename K>
  class Construct_line_2
  {
    typedef typename K::RT                        RT;
    typedef typename K::FT                        FT;
    typedef typename K::Point_2                   Point_2;
    typedef typename K::Direction_2               Direction_2;
    typedef typename K::Vector_2                  Vector_2;
    typedef typename K::Segment_2                 Segment_2;
    typedef typename K::Ray_2                     Ray_2;
    typedef typename K::Line_2                    Line_2;
    typedef typename Line_2::Rep                  Rep;
    typedef typename K::Construct_point_on_2      Construct_point_on_2;
    Construct_point_on_2 c;
  public:
    typedef Line_2            result_type;
    typedef Arity_tag< 2 >    Arity;

    Construct_line_2() {}
    Construct_line_2(const Construct_point_on_2& c_) : c(c_) {}

    Rep // Line_2
    operator()(Return_base_tag, const RT& a, const RT