CGAL预定义内核转换的问题

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;
}