格网DEM生成不规则三角网TIN

#include <iostream>
#include <string>

#include <Vec3.hpp>
#include <threeCGAL.h>
#include <gdal_priv.h>

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Projection_traits_xy_3.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Triangulation_vertex_base_with_info_2.h>
#include <CGAL/Triangulation_face_base_with_info_2.h>
#include <CGAL/boost/graph/graph_traits_Delaunay_triangulation_2.h>
#include <CGAL/boost/graph/copy_face_graph.h>
#include <CGAL/Point_set_3.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/border.h>
#include <CGAL/Polygon_mesh_processing/remesh.h>

using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
using Projection_traits = CGAL::Projection_traits_xy_3<Kernel>;
using Point_2 = Kernel::Point_2;
using Point_3 = Kernel::Point_3;
using Segment_3 = Kernel::Segment_3;
// Triangulated Irregular Network
using TIN = CGAL::Delaunay_triangulation_2<Projection_traits>;

using namespace std;

int main(int argc, char *argv[])
{
    GDALAllRegister();

    string demPath = "D:/Work/DEM2TIN/DEM.tif";
    string tinPath = "D:/Work/DEM2TIN/Tin.ply";

    GDALDataset* img = (GDALDataset *)GDALOpen(demPath.c_str(), GA_ReadOnly);
    if (!img)
    {
        cout << "Can't Open Image!" << endl;
        return 1;
    }

    int imgWidth = img->GetRasterXSize();	//图像宽度
    int imgHeight = img->GetRasterYSize();	//图像高度
    int bandNum = img->GetRasterCount();	//波段数
    //int depth = GDALGetDataTypeSize(img->GetRasterBand(1)->GetRasterDataType()) / 8;	//图像深度
    int depth = sizeof(float);	//图像深度

    double padfTransform[6];
    img->GetGeoTransform(padfTransform);
    double dx = padfTransform[1];
    double startx = padfTransform[0] + 0.5 * dx;
    double dy = -padfTransform[5];
    double starty = padfTransform[3] - imgHeight * dy + 0.5 * dy;

    //申请buf
    int bufWidth = imgWidth;
    int bufHeight = imgHeight;
    size_t imgBufNum = (size_t)bufWidth * bufHeight * bandNum;
    size_t imgBufOffset = (size_t)bufWidth * (bufHeight - 1) * bandNum;
    float *pblock = new float[imgBufNum];

    //读取
    img->RasterIO(GF_Read, 0, 0, bufWidth, bufHeight, pblock + imgBufOffset, bufWidth, bufHeight,
        GDT_Float32, bandNum, nullptr, bandNum*depth, -bufWidth * bandNum*depth, depth);

    CGAL::Point_set_3<Point_3> points;
    double zThreshold = 5;

    //
    for (int yi = 0; yi < imgHeight; yi++)
    {
        for (int xi = 0; xi < imgWidth; xi++)
        {
            //将四个角点的约束加入，保证与DEM范围一致
            if ((xi == 0 && yi == 0) || (xi == imgWidth - 1 && yi == 0) ||
                (xi == imgWidth - 1 && yi == imgHeight - 1) || (xi == 0 && yi == imgHeight - 1))
            {
                double gx1 = startx + dx * xi;
                double gy1 = starty + dy * yi;
                size_t m11 = (size_t)(imgWidth)* yi + xi;
                tinyCG::Vec3d P(gx1, gy1, pblock[m11]);
                points.insert(Point_3(P.x(), P.y(), P.z()));
            }
            else
            {
                double gx0 = startx + dx * (xi - 1);
                double gy0 = starty + dy * (yi - 1);

                double gx1 = startx + dx * xi;
                double gy1 = starty + dy * yi;

                double gx2 = startx + dx * (xi + 1);
                double gy2 = starty + dy * (yi + 1);

                size_t m00 = (size_t)imgWidth * (yi - 1) + xi - 1;
                size_t m01 = (size_t)imgWidth * (yi - 1) + xi;
                size_t m02 = (size_t)imgWidth * (yi - 1) + xi + 1;

                size_t m10 = (size_t)imgWidth* yi + xi - 1;
                size_t m11 = (size_t)imgWidth* yi + xi;
                size_t m12 = (size_t)imgWidth* yi + xi + 1;

                size_t m20 = (size_t)imgWidth * (yi + 1) + xi - 1;
                size_t m21 = (size_t)imgWidth * (yi + 1) + xi;
                size_t m22 = (size_t)imgWidth * (yi + 1) + xi + 1;

                tinyCG::Vec3d P(gx1, gy1, pblock[m11]);

                double zMeanDistance = 0;
                int counter = 0;

                if(m00 < imgBufNum && m22 < imgBufNum)
                {
                    tinyCG::Vec3d A(gx0, gy0, pblock[m00]);
                    tinyCG::Vec3d E(gx2, gy2, pblock[m22]);
                    zMeanDistance = zMeanDistance + tinyCG::threeCGAL::CalDistancePointAndLine(P, A, E);
                    counter++;
                }

                if (m02 < imgBufNum && m20 < imgBufNum)
                {
                    tinyCG::Vec3d C(gx2, gy0, pblock[m02]);
                    tinyCG::Vec3d G(gx0, gy2, pblock[m20]);
                    zMeanDistance = zMeanDistance + tinyCG::threeCGAL::CalDistancePointAndLine(P, C, G);
                    counter++;
                }

                if (m01 < imgBufNum && m21 < imgBufNum)
                {
                    tinyCG::Vec3d B(gx1, gy0, pblock[m01]);
                    tinyCG::Vec3d F(gx1, gy2, pblock[m21]);
                    zMeanDistance = zMeanDistance + tinyCG::threeCGAL::CalDistancePointAndLine(P, B, F);
                    counter++;
                }

                if (m12 < imgBufNum && m10 < imgBufNum)
                {
                    tinyCG::Vec3d D(gx2, gy1, pblock[m12]);
                    tinyCG::Vec3d H(gx0, gy1, pblock[m10]);
                    zMeanDistance = zMeanDistance + tinyCG::threeCGAL::CalDistancePointAndLine(P, D, H);
                    counter++;
                }

                zMeanDistance = zMeanDistance / counter;

                if (zMeanDistance > zThreshold)
                {
                    points.insert(Point_3(P.x(), P.y(), P.z()));
                }
            }
        }
    }


    delete[] pblock;
    pblock = nullptr;

    GDALClose(img);

    // Create DSM
    TIN dsm (points.points().begin(), points.points().end());

    using Mesh = CGAL::Surface_mesh<Point_3>;
    Mesh dsm_mesh;
    CGAL::copy_face_graph (dsm, dsm_mesh);
    std::ofstream dsm_ofile (tinPath, std::ios_base::binary);
    CGAL::set_binary_mode (dsm_ofile);
    CGAL::write_ply (dsm_ofile, dsm_mesh);
    dsm_ofile.close();

    return 0;
}