您的位置:首页 > 其它

复杂多边形的三角剖分

2021-06-19 17:56 1301 查看

目录
  • 2. 详论
  • 3. 参考
  • 1. 概述

    1.1. 多边形分类

    需要首先明确的是多边形的分类,第一种是最简单的凸多边形:

    凸多边形的每个内角都是锐角或钝角,这种多边形最普通也最常见。如果至少存在一个角是优角(大于180度小于360度),那么就是凹多边形了:

    以上多边形有一个共同特征就是由单个环线的边界组成。如果存在一个外环和多个内环组成多边形,那么就是带洞多变形了:

    如上图所示的多边形是由一个外环和两个内环组成的,两个内环造成了外环多边形的孔洞,也就是带洞多边形了。

    1.2. 三角剖分

    三角剖分也叫做三角化,或者分格化(tessellation/triangulation),将复杂的多边形剖分成多个三角形。这在图形学上有非常多的好处,便于绘制和计算。这类算法往往与Delaunay三角网算法相关,多边形的边界作为Delaunay三角网的边界约束,从而得到比较好的三角网。

    2. 详论

    我曾经在《通过CGAL将一个多边形剖分成Delaunay三角网》这篇文章中,通过CGAL实现了一个多边形的剖分,不过这个文章介绍的算法内容不支持凹多边形和带洞多边形(这也是很多其他算法实现的问题)。所以我继续翻了CGAL的官方文档,在《2D Triangulation》这一章中确实介绍了带洞多边形的三角剖分的案例。由于带洞多边形最复杂,那么我通过这个案例,来实现一下带洞多边形的三角剖分。

    #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
    #include <CGAL/Constrained_Delaunay_triangulation_2.h>
    #include <CGAL/Triangulation_face_base_with_info_2.h>
    #include <CGAL/Polygon_2.h>
    #include <iostream>
    
    #include <gdal_priv.h>
    #include <ogrsf_frmts.h>
    
    struct FaceInfo2
    {
    FaceInfo2() {}
    int
    56c
    nesting_level;
    bool in_domain() {
    return nesting_level % 2 == 1;
    }
    };
    
    typedef CGAL::Exact_predicates_inexact_constructions_kernel       K;
    typedef CGAL::Triangulation_vertex_base_2<K>                      Vb;
    typedef CGAL::Triangulation_face_base_with_info_2<FaceInfo2, K>    Fbb;
    typedef CGAL::Constrained_triangulation_face_base_2<K, Fbb>        Fb;
    typedef CGAL::Triangulation_data_structure_2<Vb, Fb>               TDS;
    typedef CGAL::Exact_predicates_tag                                Itag;
    typedef CGAL::Constrained_Delaunay_triangulation_2<K, TDS, Itag>  CDT;
    typedef CDT::Point                                                Point;
    typedef CGAL::Polygon_2<K>                                        Polygon_2;
    typedef CDT::Face_handle                                          Face_handle;
    
    void
    mark_domains(CDT& ct,
    Face_handle start,
    int index,
    std::list<CDT::Edge>& border)
    {
    if (start->info().nesting_level != -1) {
    return;
    }
    std::list<Face_handle> queue;
    queue.push_back(start);
    while (!queue.empty()) {
    Face_handle fh = queue.front();
    queue.pop_front();
    if (fh->info().nesting_level == -1) {
    fh->info().nesting_level = index;
    for (int i = 0; i < 3; i++) {
    CDT::Edge e(fh, i);
    Face_handle n = fh->neighbor(i);
    if (n->info().nesting_level == -1) {
    if (ct.is_c
    ad8
    onstrained(e)) border.push_back(e);
    else queue.push_back(n);
    }
    }
    }
    }
    }
    //explore set of facets connected with non constrained edges,
    //and attribute to each such set a nesting level.
    //We start from facets incident to the infinite vertex, with a nesting
    //level of 0. Then we recursively consider the non-explored facets incident
    //to constrained edges bounding the former set and increase the nesting level by 1.
    //Facets in the domain are those with an odd nesting level.
    void
    mark_domains(CDT& cdt)
    {
    for (CDT::Face_handle f : cdt.all_face_handles()) {
    f->info().nesting_level = -1;
    }
    std::list<CDT::Edge> border;
    mark_domains(cdt, cdt.infinite_face(), 0, border);
    while (!border.empty()) {
    CDT::Edge e = border.front();
    border.pop_front();
    Face_handle n = e.first->neighbor(e.second);
    if (n->info().nesting_level == -1) {
    mark_domains(cdt, n, e.first->info().nesting_level + 1, border);
    }
    }
    }
    
    int main()
    {
    //创建三个不相交的嵌套多边形
    Polygon_2 polygon1;
    polygon1.push_back(Point(-0.558868038740926, -0.38960351089588));
    polygon1.push_back(Point(2.77833686440678, 5.37465950363197));
    polygon1.push_back(Point(6.97052814769976, 8.07751967312349));
    polygon1.push_back(Point(13.9207400121065, 5.65046156174335));
    polygon1.push_back(Point(15.5755523607748,-1.98925544794189));
    polygon1.push_back(Point(6.36376361985472, -6.18144673123487));
    
    Polygon_2 polygon2;
    polygon2.push_back(Point(2.17935556413387, 1.4555590039808));
    polygon2.push_back(Point(3.75630057749723, 4.02942327866582));
    polygon2.push_back(Point(5.58700685737883, 4.71820385921534));
    polygon2.push_back(Point(6.54767450919789, 1.76369768475295));
    polygon2.push_back(Point(5.71388749063795, -0.900795613688593));
    polygon2.push_back(Point(3.21252643495814, -0.320769861646896));
    
    Polygon_2 polygon3;
    polygon3.push_back(Point(7.74397762278389, 0.821155837685192));
    polygon3.push_back(Point(9.13966458863422, 4.24693293568146));
    polygon3.push_back(Point(10.1909612642098, 1.83620090375816));
    polygon3.push_back(Point(12.1485481773505, 4.84508449247446));
    polygon3.push_back(Point(11.4416417920497, -2.29648257953892));
    polygon3.push_back(Point(10.1547096547072, 0.712401009177374));
    
    //将多边形插入受约束的三角剖分
    CDT cdt;
    cdt.insert_constraint(polygon1.vertices_begin(), polygon1.vertices_end(), true);
    cdt.insert_constraint(polygon2.vertices_begin(), polygon2.vertices_end(), true);
    cdt.insert_constraint(polygon3.vertices_begin(), polygon3.vertices_end(), true);
    
    //标记由多边形界定的域内的面
    mark_domains(cdt);
    
    //遍历所有的面
    int count = 0;
    for (Face_handle f : cdt.finite_face_handles())
    {
    if (f->info().in_domain()) ++count;
    }
    std::cout << "There are
    ad0
    " << count << " facets in the domain." << std::endl;
    
    //将结果输出成shp文件,方便查看
    {
    GDALAllRegister();
    
    GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("ESRI Shapefile");
    if (!driver)
    {
    printf("Get Driver ESRI Shapefile Error!\n");
    return false;
    }
    
    const char *filePath = "D:/test.shp";
    GDALDataset* dataset = driver->Create(filePath, 0, 0, 0, GDT_Unknown, NULL);
    OGRLayer* poLayer = dataset->CreateLayer("test", NULL, wkbPolygon, NULL);
    
    //创建面要素
    for (Face_handle f : cdt.finite_face_handles())
    {
    if (f->info().in_domain())
    {
    OGRFeature *poFeature = new OGRFeature(poLayer->GetLayerDefn());
    
    OGRLinearRing ogrring;
    for (int i = 0; i < 3; i++)
    {
    ogrring.setPoint(i, f->vertex(i)->point().x(), f->vertex(i)->point().y());
    }
    ogrring.closeRings();
    
    OGRPolygon polygon;
    polygon.addRing(&ogrring);
    poFeature->SetGeometry(&polygon);
    
    if (poLayer->CreateFeature(poFeature) != OGRERR_NONE)
    {
    printf("Failed to create feature in shapefile.\n");
    return false;
    }
    }
    }
    
    //释放
    GDALClose(dataset);
    dataset = nullptr;
    }
    
    return 0;
    }

    在代码的最后,我将生成的三角网输出成shp文件,叠加到原来的多边形中:

    效果似乎不是很明显,那么我将原来的两个内环范围涂黑:

    说明这个算法可以适配于凹多边形以及带洞多边形的三角网剖分,不得不说CGAL这个库真的非常强大。可惜就是这个库太难以使用了,需要一定的计算几何知识和Cpp高级特性的知识,才能运用自如,值得深入学习。

    3. 参考

    CGAL 5.2.2 - 2D Triangulation

    内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
    标签: