如何防止CGAL生成3D Alpha形状的简并四面体？

Question

我正在尝试在3D中检索alpha形状内的四面体以进行一些FEM计算。我使用的代码如下所示：

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Triangulation_vertex_base_with_info_3.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Alpha_shape_3.h>
#include <CGAL/Alpha_shape_vertex_base_3.h>
#include <CGAL/Alpha_shape_cell_base_3.h>


typedef CGAL::Exact_predicates_inexact_constructions_kernel                 Kernel;
typedef Kernel::FT                                                          FT;
typedef CGAL::Triangulation_vertex_base_with_info_3<std::size_t, Kernel>    Vb3;
typedef CGAL::Alpha_shape_vertex_base_3<Kernel, Vb3>                        asVb3;
typedef CGAL::Alpha_shape_cell_base_3<Kernel>                               asCb3;
typedef CGAL::Triangulation_data_structure_3<asVb3,asCb3>                   asTds3;
typedef CGAL::Delaunay_triangulation_3<Kernel,asTds3, CGAL::Fast_location>  asTriangulation_3;
typedef CGAL::Alpha_shape_3<asTriangulation_3>                              Alpha_shape_3;
typedef Kernel::Point_3                                                     Point_3;
typedef Kernel::Tetrahedron_3                                               Tetrahedron_3;

void triangulateAlphaShape3D()
{
    ...

    const Alpha_shape_3 as(pointsList.begin(), pointsList.end(),
                           FT(alpha),
                           Alpha_shape_3::GENERAL);

    for(Alpha_shape_3::Finite_cells_iterator fit = as.finite_cells_begin() ;
        fit != as.finite_cells_end() ; ++fit)
    {
        if(as.classify(fit) == Alpha_shape_3::INTERIOR)
        {
            const Alpha_shape_3::Cell_handle cell{fit};

            const std::vector<std::size_t> element{cell->vertex(0)->info(),
                                                   cell->vertex(1)->info(),
                                                   cell->vertex(2)->info(),
                                                   cell->vertex(3)->info()};
        }
    }

    ...
}

我的点列表由在盒子中形成立方体的点组成。 CGAL的输出是（在gmsh中已Vizualized）：

Output of CGAL vizualized in gmsh

然而，图中似乎存在某种交叉/重叠的元素，当在调试时通过gdb在alpha形状构造函数处使用断点执行该程序时，程序会抛出以下内容：文件Uncertain.h：

T make_certain() const
  {
    if (is_certain())
      return _i;
    CGAL_PROFILER(std::string("Uncertain_conversion_exception thrown for CGAL::Uncertain< ")
          + typeid(T).name() + " >");
    throw Uncertain_conversion_exception(
                  "Undecidable conversion of CGAL::Uncertain<T>");
  }

调用堆栈为：

#0 0xa31c08 libstdc++-6!.cxa_throw() (libstdc++-6.dll:??)
#1 0x62060a4a   CGAL::Uncertain<bool>::make_certain(this=0x22d0ac) (CGAL/Uncertain.h:125)
#2 0x62060ae5   CGAL::Uncertain<bool>::operator bool(this=0x22d0ac) (CGAL/Uncertain.h:133)
#3 0x620240bd   CGAL::coplanar_orientationC3<CGAL::Interval_nt<false> >(px=..., py=..., pz=..., qx=..., qy=..., qz=..., rx=..., ry=..., rz=...) (CGAL/predicates/kernel_ftC3.h:209)
#4 0x62056e68   CGAL::CartesianKernelFunctors::Coplanar_orientation_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >::operator() (this=0x22d5af, p=..., q=..., r=...) (CGAL/Cartesian/function_objects.h:3649)
#5 0x620518ea   CGAL::Filtered_predicate<CGAL::CartesianKernelFunctors::Coplanar_orientation_3<CGAL::Simple_cartesian<CGAL::Mpzf> >, CGAL::CartesianKernelFunctors::Coplanar_orientation_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >, CGAL::Cartesian_converter<CGAL::Type_equality_wrapper<CGAL::Cartesian_base_no_ref_count<double, CGAL::Epick>, CGAL::Epick>, CGAL::Simple_cartesian<CGAL::Mpzf>, CGAL::NT_converter<double, CGAL::Mpzf> >, CGAL::Cartesian_converter<CGAL::Type_equality_wrapper<CGAL::Cartesian_base_no_ref_cou (CGAL/Filtered_predicate.h:167)
#6 0x6205081f   CGAL::Triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_bas (CGAL/Triangulation_3.h:589)
#7 0x62059ec2   CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_v (CGAL/Delaunay_triangulation_3.h:1509)
#8 0x6205950e   CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_v (CGAL/Delaunay_triangulation_3.h:1545)
#9 0x6205988c   CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_v (CGAL/Delaunay_triangulation_3.h:775)
#10 0x620599c9  CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_v (CGAL/Delaunay_triangulation_3.h:857)
#11 0x6204e578  CGAL::Triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_bas (CGAL/Triangulation_3.h:1323)
#12 0x6201d924  CGAL::Triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_bas (CGAL/Triangulation_3.h:3696)
#13 0x620254cc  CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_v (CGAL/Delaunay_triangulation_3.h:1165)
#14 0x620257bd  CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_v (CGAL/Delaunay_triangulation_3.h:1150)
#15 0x62025636  CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_v (CGAL/Delaunay_triangulation_3.h:537)
#16 0x62027f9c  CGAL::Triangulation_hierarchy_3<CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CG (CGAL/Triangulation_hierarchy_3.h:279)
#17 0x6202816f  CGAL::Triangulation_hierarchy_3<CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Triangulation_hierarchy_vertex_base_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CG (CGAL/Triangulation_hierarchy_3.h:300)
#18 0x62017b27  CGAL::Alpha_shape_3<CGAL::Delaunay_triangulation_3<CGAL::Epick, CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Triangulation_vertex_base_with_info_3<unsigned long long, CGAL::Epick, CGAL::Triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<void> > >, CGAL::Boolean_tag<false>, CGAL::Boolean_tag<false> >, CGAL::Alpha_shape_cell_base_3<CGAL::Epick, CGAL::Default, CGAL::Boolean_tag<false>, CGAL::Boolean_tag<false> >, CGAL::Sequential_tag>, CGAL (CGAL/Alpha_shape_3.h:269)

它们是某种简并的四面体吗？如何防止CGAL生成它们？

我目前使用CGAL 4.11。

更新：我切换到CGAL 5.0.2，并且抛出的问题消失了。但是我仍然有这个奇怪的问题。测试代码：

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Triangulation_vertex_base_with_info_3.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Alpha_shape_3.h>
#include <CGAL/Alpha_shape_vertex_base_3.h>
#include <CGAL/Alpha_shape_cell_base_3.h>

#include <gmsh.h>

#include <algorithm>
#include <fstream>
#include <random>

typedef CGAL::Exact_predicates_inexact_constructions_kernel                 Kernel;
typedef Kernel::FT                                                          FT;
typedef CGAL::Triangulation_vertex_base_with_info_3<std::size_t, Kernel>    Vb3;
typedef CGAL::Alpha_shape_vertex_base_3<Kernel, Vb3>                        asVb3;
typedef CGAL::Alpha_shape_cell_base_3<Kernel>                               asCb3;
typedef CGAL::Triangulation_data_structure_3<asVb3,asCb3>                   asTds3;
typedef CGAL::Delaunay_triangulation_3<Kernel,asTds3, CGAL::Fast_location>  asTriangulation_3;
typedef CGAL::Alpha_shape_3<asTriangulation_3>                              Alpha_shape_3;
typedef Kernel::Point_3                                                     Point_3;
typedef Kernel::Tetrahedron_3                                               Tetrahedron_3;

int main(int argc, char **argv)
{
    std::random_device rd;
    std::mt19937 mt(rd());
    std::uniform_real_distribution<double> dist(-0.000000001, 0.000000001);
    //std::uniform_real_distribution<double> dist(0,0);

    std::vector<std::vector<std::size_t>> elementsList;

    const double alpha = 0.9;

    // We have to construct an intermediate representation for CGAL. We also reset
    // nodes properties.
    std::vector<std::pair<Point_3, std::size_t>> pointsList;

    std::size_t counter = 0;

    /********************************************************************************
                                     Box Nodes
     *******************************************************************************/
    for(double y = 1; y < 10 ; y += 1)
    {
        for(double z = 1; z <= 10 ; z += 1)
        {
            std::pair<Point_3, std::size_t> point = std::make_pair(Point_3(10, y + dist(mt), z + dist(mt)), counter);
            pointsList.push_back(point);
            counter++;
        }
    }

    for(double x = 0; x <= 10 ; x += 1)
    {
        for(double z = 1; z <= 10 ; z += 1)
        {
            std::pair<Point_3, std::size_t> point = std::make_pair(Point_3(x+ dist(mt), 0, z+ dist(mt)), counter);
            pointsList.push_back(point);
            counter++;
        }
    }

    for(double y = 1; y <= 10 ; y += 1)
    {
        for(double z = 1; z <= 10 ; z += 1)
        {
            std::pair<Point_3, std::size_t> point = std::make_pair(Point_3(0, y+ dist(mt), z+ dist(mt)), counter);
            pointsList.push_back(point);
            counter++;
        }
    }

    for(double x = 1; x <= 10 ; x += 1)
    {
        for(double z = 1; z <= 10 ; z += 1)
        {
            std::pair<Point_3, std::size_t> point = std::make_pair(Point_3(x+ dist(mt), 10, z+ dist(mt)), counter);
            pointsList.push_back(point);
            counter++;
        }
    }

    for(double x = 0; x <= 10 ; x += 1)
    {
        for(double y = 0; y <= 10 ; y += 1)
        {
            std::pair<Point_3, std::size_t> point = std::make_pair(Point_3(x+ dist(mt), y+ dist(mt), 0), counter);
            pointsList.push_back(point);
            counter++;
        }
    }


    /********************************************************************************
                                     Cube nodes
     *******************************************************************************/
    for(double x = 1; x <= 5 ; x += 1)
    {
        for(double y = 1; y <= 5 ; y += 1)
        {
            for(double z = 1; z <= 5 ; z += 1)
            {
                std::pair<Point_3, std::size_t> point = std::make_pair(Point_3(x+ dist(mt), y+ dist(mt), z+ dist(mt)), counter);
                pointsList.push_back(point);
                counter++;
            }
        }
    }

    std::random_shuffle(pointsList.begin(), pointsList.end());

    for(std::size_t n = 0; n < pointsList.size() ; ++n)
    {
        pointsList[n].second = n;
    }

    for(std::size_t n = 0 ; n < pointsList.size() ; ++n)
    {
        for(std::size_t n2 = 0 ; n2 < pointsList.size() ; ++n2)
        {
            if(n != n2)
            {
                if(pointsList[n].first == pointsList[n2].first)
                {
                    std::cerr << "Duplicate node found: " << "(" << pointsList[n].first.x() << ", "
                                                                 << pointsList[n].first.y() << ", "
                                                                 << pointsList[n].first.z() << ")";
                    std::cerr << std::endl;
                    std::cerr << "Counter: " << pointsList[n].second << " vs " << pointsList[n2].second << std::endl;
                    throw std::runtime_error("Duplicates nodes!");
                }
            }
        }
    }

    std::cout << "Number of point: " << pointsList.size() << std::endl;

    const Alpha_shape_3 as(pointsList.begin(), pointsList.end(),
                           FT(alpha),
                           Alpha_shape_3::GENERAL);

    // We check for each triangle which one will be kept (alpha shape), then we
    // perfom operations on the remaining elements
    for(Alpha_shape_3::Finite_cells_iterator fit = as.finite_cells_begin() ;
        fit != as.finite_cells_end() ; ++fit)
    {
        // If true, the elements are fluid elements
        if(as.classify(fit) == Alpha_shape_3::INTERIOR)
        {
            const Alpha_shape_3::Cell_handle cell{fit};

            const std::vector<std::size_t> element{cell->vertex(0)->info(),
                                                   cell->vertex(1)->info(),
                                                   cell->vertex(2)->info(),
                                                   cell->vertex(3)->info()};

            elementsList.push_back(element);
        }
    }

    std::cout << "Number of Element: " << elementsList.size() << std::endl;

    gmsh::initialize();
    gmsh::model::add("theModel");
    gmsh::model::setCurrent("theModel");
    gmsh::model::addDiscreteEntity(0, 1);
    gmsh::model::addDiscreteEntity(3, 2);
    gmsh::view::add("Nothing", 1);

    std::vector<std::size_t> nodesTags(pointsList.size());
    std::vector<double> nodesCoord(3*pointsList.size());

    for(std::size_t n = 0 ; n < pointsList.size() ; ++n)
    {
        nodesTags[n] = n + 1;
        nodesCoord[3*n] = pointsList[n].first.x();
        nodesCoord[3*n + 1] = pointsList[n].first.y();
        nodesCoord[3*n + 2] = pointsList[n].first.z();
    }

    gmsh::model::mesh::addNodes(0, 1, nodesTags, nodesCoord);

    std::vector<std::size_t> elementTags(elementsList.size());
    std::vector<std::size_t> nodesTagsPerElement(4*elementsList.size());
    for(std::size_t elm = 0 ; elm < elementsList.size() ; ++elm)
    {
        elementTags[elm] = elm + 1;
        for(std::size_t n = 0 ; n < 4 ; ++n)
            nodesTagsPerElement[4*elm + n] = pointsList[elementsList[elm][n]].second + 1;
    }

    gmsh::model::mesh::addElementsByType(2, 4, elementTags, nodesTagsPerElement); //Tetrahedron 4
    gmsh::model::mesh::addElementsByType(1, 15, nodesTags, nodesTags); //Point

    std::vector<std::vector<double>> dataNothing(pointsList.size());

    for(std::size_t n = 0 ; n < pointsList.size() ; ++n)
    {
        const std::vector<double> nothingVec{0};
        dataNothing[n] = nothingVec;
    }

    gmsh::view::addModelData(1, 0, "theModel", "NodeData", nodesTags, dataNothing, 0, 1);

    gmsh::view::write(1, "mesh.msh", false);

    gmsh::model::remove();

    gmsh::finalize();

    return 0;
}

[如果我稍微干扰点位置（在我的程序中，我使用gmsh生成初始状态，并且节点位置并非完全“完整”），则会出现那些奇怪的元素，而如果节点位置是完全整数，则元素数量为正确（在这种情况下为1017）。

Answer 1

您的代码基本上是将点放在规则的网格上，然后稍稍扰动。结果，您最终得到了通常所说的条子（在您最喜欢的搜索引擎中查找“ sets tet mesh”）。我的建议是使用常规网格生成初始网格并将扰动应用于输入（由于顶点上有ID，因此可以轻松地将扰动存储在循环中）。