1. 计算机字长造成了计算几何的“麻烦”
由于受计算机字长的影响,计算机无法随心所欲地表达实数域内每一个数字,这给数值计算带来诸多麻烦,比如在计算中舍入误差是常常要关注的“麻烦”,它就像噪声对于信号分析与处理问题一样如影随行,不离不弃。对舍入误差非常敏感的数值问题,常常因为计算机的这一局限性带来了对计算结果精度的不确定性。
在计算几何领域,这一问题尤为突出,因为在数值领域这一矛盾虽突出,但不是每个问题都需要关注计算机字长,而在计算几何领域则不然。计算几何的研究对象是几何体,对几何的计算与对数值的计算类似,但无论这些算法是简单还是复杂,从计算结果来看无非涉及两类运算,在CGAL中称为:判断(predictes)和构造(constructions)。所谓判断,主要用于比较和分类,其计算返回值是布尔类型或是枚举类型,而构造则是要利用算法在现有几何对象的基础上产生新的几何对象。通常对于判断来说,其数值精度要求不能有丝毫含糊,要求最高,而对于构造来说其精度要求可与运算速度之间进行平衡选择。
从另一方面来说,计算机字长对舍入误差的影响虽然普遍,但不是对任何类型的数字都有舍入误差。举个来说,对于整型数,如果只做加法、减法或乘法运算,那么计算机是不会产生任何误差的。由于在实际应用中,只需处理整型数的加、减、乘运算的问题占比非常小,所以舍入误差问题看起来才非常普遍。实际上只需要满足两个条件,这个问题就可以认为不需要考虑舍入误差:
1、不会给计算机带来“溢出”的数据类型(如整型,以及分子分母分开存储的有理性的数据)
2、计算在数域内封闭,即计算结果仍然是本数域内的数。
实际上第1条是对参与计算的对象的要求,因为第2条是对计算结果的要求,其本质还是文中最开始所讲的计算机字长的限制。
判断类的问题,其计算通常类似于行列式运算,只涉及加、减、乘,因此根据计算对象的数值类型。在这种条件下,如果是有限精度的整型数可用int,long或double数据类型;如果是任意精度整数类型,可用GMP整数的包装器Gmpz、leda_integer或MP_Float数据类型。注意,除非使用任意精度环类型(支持封闭运算),否则可能会由于溢出而产生不正确的结果。
构造类的问题,通常涉及除法运算,因此其数据类型需要对商运算封闭,具体类型的选择可参考CGAL(2D and 3D Linear Geometry Kernel)。
2. CGAL预定义内核
为了用户方便,CGAL针对笛卡尔坐标计算提供了5种预定义内核:
1、Exact_predicates_inexact_constructions_kernel; 2、Exact_predicates_exact_constructions_kernel; 3、Exact_predicates_exact_constructions_kernel_with_sqrt; 4、Exact_predicates_exact_constructions_kernel_with_kth_root; 5、Exact_predicates_exact_constructions_kernel_with_root_of.
常用的是前三种,特别是第一种和第二种,但是第一种和第二种不会混用,但是第二种和第三种可能会出现混用的情况,因而会有内核转换的问题。
3. CGAL内核转换
内核如何转换,Convert CGAL::Lazy_exact_nt<CGAL::Gmpq> to Double、Cast from Exact_predicates_exact_constructions_kernel to Exact_predicates_exact_constructions_kernel_with_sqrt、CGAL coordinate value transformation not working anymore in 5.2.1等问题都有涉及,本文结合这些帖子,并结合CGAL::Cartesian_converter< K1, K2, NTConverter > Class Template Reference给出较为全面的例子,并扩大了转换的应用范围,将其应用于Surface Mesh和Vector等对象上。
#define CGAL_DO_NOT_USE_BOOST_MP #include <CGAL/Exact_predicates_exact_constructions_kernel_with_sqrt.h> #include <CGAL/Exact_predicates_exact_constructions_kernel.h> #include <CGAL/Cartesian_converter.h> #include <CGAL/Surface_mesh.h> #include <CGAL/Linear_algebraCd.h> typedef CGAL::Exact_predicates_exact_constructions_kernel Epeck; typedef CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt Epecks; typedef Epeck::Point_3 point; typedef Epecks::Point_3 Point; typedef CGAL::Surface_mesh<point> surface_mesh; typedef CGAL::Surface_mesh<Point> Surface_Mesh; typedef boost::graph_traits<surface_mesh>::vertex_descriptor vertex_descriptor; typedef boost::graph_traits<Surface_Mesh>::vertex_descriptor Vertex_Descriptor; typedef CGAL::Linear_algebraCd<Epeck::FT> la_epec; typedef la_epec::Matrix matrix; typedef la_epec::Vector vector; typedef CGAL::Linear_algebraCd<Epecks::FT> LA_epec; typedef LA_epec::Matrix Matrix; typedef LA_epec::Vector Vector; struct Lazy_gmpq_to_Expr_converter : public std::unary_function< CGAL::Lazy_exact_nt<CGAL::Gmpq>, CORE::Expr > { CORE::Expr operator()(const CGAL::Lazy_exact_nt<CGAL::Gmpq> &a) const { return ::CORE::BigRat(exact(a).mpq()); // } }; struct Lazy_Expr_to_gmpq_converter : public std::unary_function<CORE::Expr, CGAL::Lazy_exact_nt<CGAL::Gmpq> > { CGAL::Lazy_exact_nt<CGAL::Gmpq> operator()(const CORE::Expr &a) const { std::stringstream ss; ss << a; CGAL::Lazy_exact_nt<CGAL::Gmpq> b; ss >> b; return b; // } }; typedef CGAL::Cartesian_converter<Epeck, Epecks, Lazy_gmpq_to_Expr_converter> IK_to_EK; typedef CGAL::Cartesian_converter<Epecks, Epeck, Lazy_Expr_to_gmpq_converter> EK_to_IK; void vec_to_Vec(vector &v, Vector &V) { if (v.dimension() != V.dimension()) { std::cerr << "The two vectors have different dimensions!" << std::endl; std::exit(0); return; } IK_to_EK to_sqrt_kernel; int i = 0; for (vector::const_iterator it = v.begin(); it != v.end(); ++it, ++i) { //std::cout << it <<" "<< *it << std::endl; V[i] = to_sqrt_kernel(*it); //std::cout << V[i] << std::endl; } } int main() { //case1:核自身成员来回转核操作 Epeck::Triangle_3 t1( Epeck::Point_3(0., 0., 0.), Epeck::Point_3(1., 0., -1.), Epeck::Point_3(0., 1., 3.) ); Epeck::Line_3 l1( Epeck::Point_3(0.2, 0.25, -7), Epeck::Point_3(0.25, 0.3, 4) ); IK_to_EK to_exact; Epecks::Triangle_3 t2 = to_exact(t1); Epecks::Line_3 l2 = to_exact(l1); const auto inter = CGAL::intersection(t2, l2); // As we are sure that there IS an intersection // and that the intersection IS a point // we do not have to check for this, or put it in a try/catch const Epecks::Point_3& exact_pt = boost::get<Epecks::Point_3>(*inter); EK_to_IK to_inexact; Epeck::Point_3 inexact_pt = to_inexact(exact_pt); std::cout << inexact_pt << std::endl; //case2:非核成员转核操作 point p0(0, 0, 0); point p1(0, 1, 0); point p2(1, 0, 0); surface_mesh sm; vertex_descriptor v0 = sm.add_vertex(p0); vertex_descriptor v1 = sm.add_vertex(p1); vertex_descriptor v2 = sm.add_vertex(p2); sm.add_face(v0, v1, v2); surface_mesh::vertex_iterator it = sm.vertices_begin(); Epecks::Point_3 P0 = to_exact(sm.point(*it++)); Epecks::Point_3 P1 = to_exact(sm.point(*it++)); Epecks::Point_3 P2 = to_exact(sm.point(*it++)); Surface_Mesh SM; Vertex_Descriptor V0 = SM.add_vertex(P0); Vertex_Descriptor V1 = SM.add_vertex(P1); Vertex_Descriptor V2 = SM.add_vertex(P2); SM.add_face(V0, V1, V2); Point P3(0, 0, 1); Point P4(0, 1, 1); Point P5(1, 0, 1); vertex_descriptor V3 = SM.add_vertex(P3); vertex_descriptor V4 = SM.add_vertex(P4); vertex_descriptor V5 = SM.add_vertex(P5); SM.add_face(V3, V4, V5); Surface_Mesh::vertex_iterator jt = SM.vertices_begin()+3; Epeck::Point_3 p3 = to_inexact(SM.point(*jt++)); Epeck::Point_3 p4 = to_inexact(SM.point(*jt++)); Epeck::Point_3 p5 = to_inexact(SM.point(*jt++)); vertex_descriptor v3 = sm.add_vertex(p3); vertex_descriptor v4 = sm.add_vertex(p4); vertex_descriptor v5 = sm.add_vertex(p5); sm.add_face(v3, v4, v5); for (surface_mesh::Vertex_iterator it = sm.vertices_begin(); it != sm.vertices_end(); ++it) { std::cout << sm.point(*it) << std::endl; } std::cout << std::endl; for (Surface_Mesh::Vertex_iterator it = SM.vertices_begin(); it != SM.vertices_end(); ++it) { std::cout << SM.point(*it) << std::endl; } //case3: vector的转核操作 vector vec1(3, 3.0); vector vec2(3, 2.0); Vector vec3(3); Vector vec4(3); vec_to_Vec(vec1, vec3); vec_to_Vec(vec2, vec4); Epecks::FT d = CGAL::sqrt(vec3 * vec4); std::cout << d << std::endl; return 0; }