俗话说好记性不如烂笔头，平时项目中遇到的问题极有可能在接下来的项目中再次偶遇，最近在搞一个基于钻孔数据构建地层模型的项目，所以把项目过程中遇到的问题做一下延伸并记录下来。
在Cesium中进行地形抬升以及填挖方量计算时大多需要再所选的研究区域构建不规则三角网，js构建不规则三角网有现成的库，如delaunator、d3js；推荐使用delaunator。有个问题，因为Cesium地形获取的坐标是小数点位数非常多的浮点数，导致构建出来的三角网不符合实际需求；因为浮点数计算的问题可能导致一条直线上的三个点会构成一个三角形。
针对以上问题，我的解决方式可分为以下几个步骤:

1、研究区域细分

在Cesium中，可以通过拾取的顶点获取研究区域的矩形包围盒；对矩形包围盒的东西、南北按照一定的距离进行插值细分；会得到一大堆坐标点，注意这里要转换为度单位的经纬度坐标，然后判断插值的点是否在研究区域内，如果在研究区域内，就把点保存下来

// pickedpositions为鼠标点选的坐标信息
// excateOptions.lerpMeter 为插值距离（单位为米）
// meshpositions研究区域内坐标点集合（笛卡尔坐标）
 // meshpoints 研究区域内坐标点集合（经纬度坐标，单位为角度）用于构建不规则三角网
const rectangle = Rectangle.fromCartesianArray(pickedpositions)
// 西北角
const southwest = Cartographic.toCartesian(Rectangle.southwest(rectangle))
// 东南角
const southeast = Cartographic.toCartesian(Rectangle.southeast(rectangle))
// 西北角
const northwest = Cartographic.toCartesian(Rectangle.northwest(rectangle))
const distanceH = Cartesian3.distance(southwest, southeast)
const distanceV = Cartesian3.distance(southwest, northwest)
const spanH = Math.ceil(distanceH / excateOptions.lerpMeter)
const spanV = Math.ceil(distanceV / excateOptions.lerpMeter)


// 构建turf多边形
const tmps = pickedpositions.reduce((prev, curr) => {
    const {longitude, latitude} = LonLatHeight.fromCartesian3(curr)
    return (prev.push([longitude, latitude]), prev)
}, [] as any [])
const  polygonTurf = Turf.polygon([[...tmps, tmps[0]]])

// 构建网格
const meshpositions: Cartesian3[] = []
const meshpoints: Array<[number, number]> = []
for (let i = 0; i <= spanH; i++) {
    const offsetH = Cartesian3.lerp(southwest, southeast, i / spanH, new Cartesian3())
    const { longitude } = Cartographic.fromCartesian(offsetH)
    for (let j = 0; j <= spanV; j++) {
        const offsetV = Cartesian3.lerp(southwest, northwest, j / spanV, new Cartesian3())
        const {latitude} = Cartographic.fromCartesian(offsetV)
        const position = Cartesian3.fromRadians(longitude, latitude)
        const lon = CMath.toDegrees(longitude)
        const lat = CMath.toDegrees(latitude)
        const pointTurf = Turf.point([lon, lat])
        if (!Turf.booleanPointInPolygon(pointTurf, polygonTurf))  continue
        meshpoints.push([lon, lat])
        meshpositions.push(position)
    }
} 
                    

2、研究边细分

对点选的坐标集合，两点构成直线，并以上述距离进行线段插值细分。

const boundarypositions: Cartesian3[] =  []
const boundarypoints:Array<[number, number]> = []
{
    
    for (let i = 0; i < pickedpositions.length; i++) {
        const start = pickedpositions[i]
        const end = pickedpositions[(i+1)%pickedpositions.length]
        const distance = Cartesian3.distance(start, end)
        const span = Math.ceil(distance / excateOptions.lerpMeter)

        for (let j = 0; j < span; j++) {
            const position = Cartesian3.lerp(start, end, j / span, new Cartesian3())
            const { longitude, latitude} = LonLatHeight.fromCartesian3(position)
            boundarypoints.push([longitude, latitude])
            boundarypositions.push(position)
        }

    }   
}

最后将研究区域细分结果以及研究边细分结果进行合并，并构建三角网。

const combinedpoints = meshpoints.concat(boundarypoints)
const combinedpositions = meshpositions.concat(boundarypositions)
const coords = combinedpositions.reduce((prev, curr) => {
    const {longitude, latitude} = LonLatHeight.fromCartesian3(curr)
    return prev.push([longitude, latitude]), prev
}, [] as any [])

const delaunay = Delaunay.from(combinedpoints)
const triangles = delaunay.triangles

实现效果如下：

地面三角网

可以看到图中的三角网是不符合需要的，红线范围内的三角形明显是不规范的，这是因为浮点数计算导致的误差。接下来还需要去除这种不规则的三角形。

3、去除狭长三角形

对于上述不规则三角形，他们的三个点基本上在一条直线上；实际情况这三个点就应该在一条直线上，且不应该构成三角形。但是由于经纬度小数点位太多，导致计算结果出现误差。我的处理方法是遍历三角网中所有三角形，取出三角形中的两个顶点构成直线计算他们的直线方程的斜率，并对斜率取1位有效数字，判断这三条直线的斜率是否相等，如果相等则舍弃这个三角形，反之把三角形保存下来。

// 考虑计算量可能比较大，这里使用webworker子线程计算斜率
self.onmessage = function(e) {
    const indices = []
    const {triangles, coords} = e.data
    for (let i = 0; i < triangles.length; i += 3) {
        const [a, b, c] = triangles.slice(i, i +3)
        const [x1, y1] = coords[a]
        const [x2, y2] = coords[b]
        const [x3, y3] = coords[c]
        let {k: k1} = computeLinearQquation([x1, y1], [x2, y2])
        let {k: k2} = computeLinearQquation([x1, y1], [x3, y3])
        let {k: k3} = computeLinearQquation([x2, y2], [x3, y3])
        k1 = Number(k1.toFixed(1))
        k2 = Number(k2.toFixed(1))
        k3 = Number(k3.toFixed(1))
        if (k1==k2&&k2==k3) continue
        indices.push(a, b, c)
    }
    self.postMessage(indices)
}   

/**
 * 计算线性方程的斜率和截距
 *
 * @param p1 数组类型，包含两个元素，分别为点1的横坐标和纵坐标
 * @param p2 数组类型，包含两个元素，分别为点2的横坐标和纵坐标
 * @returns 返回包含斜率和截距的对象，其中k为斜率，b为截距
 */
function computeLinearQquation(p1, p2){
    const x1 = p1[0], y1 = p1[1]
    const x2 = p2[0], y2 = p2[1]
    // 斜率
    let k = (y2 - y1) / (x2 - x1)
    if (k === Infinity || k === -Infinity) {
        return {k: x1, b: 0}
    }
    // 截距
    let b = y1 - k * x1
    return {k, b}
}

最后符合需求的不规则三角网索引，使用primitive进行绘制如下：

索引重构之后的三角网
自此，地面不规则三角网构建完毕，后续地形抬升以及填挖方量可以基于地表面三角网构建广义三菱柱来进行计算。
完整源代码已经上传至Gitee，有需要的可以拉下来试下（https://gitee.com/AnZunYi/cesium-features），代码写的有点那啥，请轻点喷。