cesium计算河道填挖方体积(基于点云数据与BIM模型)
·
背景：
时间为23年7月份某一天，接到领导通知，说要做一个三维的小项目。听到是个小项目，以为是常见的三维需求，沉浸在可以丰富技术的喜悦中，没成想是苦逼的开始....
需求：
计算某个河道的填挖方体积，甲方可以提供的只有河道的BIM模型、若干点云数据与一条河道中心线的数据 。以下第一张图为河道设计BIM模型，就是河道要改造成的样子，称之为设计数据。第二张图为河道开始改造之前提取的点云模型，称之为原始数据。第三张图是河道改造之后的点云数据模型，称之为验收数据。由原始数据和验收数据要计算出河道实际的填方与挖方的体积(根据谁的点云数据在上面判断是填方还是挖方)，以及原始数据与设计数据对比、验收数据与设计数据对比。
解决思路：
1.梳理目前已有的数据，有河道中心线数据、BIM模型、若干点云数据
2.根据传统水利行业计算填挖方的计算方式，是在河道上设置两个断面，以这两个断面的填方或者挖方面积作为底面，来计算断面间的填挖方体积。例如填方体积，断面1填方面积为s1，断面2填方面积为s2，长度为len。则计算填方体积公式为 。断面面积计算如下图：
3.现在思路差不多了，接下来就是如何根据现有的数据去实现。首先有河道的BIM模型与河道的中心线数据，可以计算河道的长度，然后每隔若干米在中心线取一个点，再根据这个点位前后的中心点确定断面的方向向量。同时为了数据计算的准确度，需要找出河道的拐点，由于中心点数据提供的比较密集，所以暂时不需要处理特殊拐点。
下面是根据中心线找出所有100m处点位的点，并且根据向量与点坐标生成河道的断面代码
/**
* 计算所有点之间距离，找出所有100米的位置，返回每个位置的法向量、坐标点、距离、位于哪两个点
* @params points 点列表
*/
const GeneratePolygons = () => {
let distancetotal = 0 //累计距离
let directionpolyline = null //每个平面法向量线
//遍历点数组，计算当前点与下一点的距离与方向
for (let i = 0; i < centerLine.length - 1; i++) {
let currentPoint = centerLine[i] //当前点
let nextPoint = centerLine[i + 1] //下一点
//添加点位
// viewer.entities.add({
// position: Cesium.Cartesian3.fromDegrees(...centerLine[i], 20),
// point: {
// color: Cesium.Color.GREEN, //点位颜色
// pixelSize: 5 //像素点大小
// },
// });
let startCartesian = Cesium.Cartesian3.fromDegrees(currentPoint[0], currentPoint[1]);
let endCartesian = Cesium.Cartesian3.fromDegrees(nextPoint[0], nextPoint[1]);
directionpolyline = [currentPoint[0], currentPoint[1], nextPoint[0], nextPoint[1]]
// 计算直线的方向向量
let direction = Cesium.Cartesian3.subtract(endCartesian, startCartesian, new Cesium
.Cartesian3());
//归一向量用于计算两点距离
let normalizedDirection = Cesium.Cartesian3.normalize(direction, new Cesium.Cartesian3());
//计算两点之间是否大于100米
let distance = Cesium.Cartesian3.distance(startCartesian, endCartesian).toFixed(2)
//距离总和
distancetotal += Number(distance)
// 判断存在几个100点
let num = Math.floor(distancetotal / 100)
//余数为多少
let remainder = distancetotal % 100
if (i == 0) {
GeneratePlane(startCartesian, normalizedDirection)
GeneratePolyline(directionpolyline)
}
if (num) {
let flag = true;
//100米位置点距离上一个点的距离
let distancePois = 100 - distancetotal + Number(distance)
for (let i = 0; i < num; i++) {
if (i > 0) {
distancePois = distancePois + 100
}
//计算100米位置偏移量
let offset = Cesium.Cartesian3.multiplyByScalar(normalizedDirection, distancePois,
new Cesium
.Cartesian3())
//中心坐标
let pointOnLine = Cesium.Cartesian3.add(startCartesian, offset, new Cesium
.Cartesian3());
GeneratePlane(pointOnLine, normalizedDirection)
if (flag) {
//这里代码只执行一次
GeneratePolyline(directionpolyline)
flag = false;
}
}
//距离总和改为整除100的余数
distancetotal = remainder
}
}
}
/**
* @param {directionpolyline} 法向量线
*/
const GeneratePolyline = (directionpolyline) => {
let arrowMaterial = new Cesium.PolylineArrowMaterialProperty(Cesium.Color.BLUE);
viewer.entities.add({
polyline: {
// fromDegrees返回给定的经度和纬度值数组（以度为单位），该数组由Cartesian3位置组成。
// Cesium.Cartesian3.fromDegreesArray([经度1, 纬度1, 经度2, 纬度2,])
// Cesium.Cartesian3.fromDegreesArrayHeights([经度1, 纬度1, 高度1, 经度2, 纬度2, 高度2])
positions: Cesium.Cartesian3.fromDegreesArray(directionpolyline),
// 宽度
width: 5,
// 线的颜色
material: arrowMaterial,
// 线的顺序,仅当`clampToGround`为true并且支持地形上的折线时才有效。
zIndex: 10,
// 显示在距相机的距离处的属性，多少区间内是可以显示的
distanceDisplayCondition: new Cesium.DistanceDisplayCondition(0, 1500),
// 是否显示
show: true,
clampToGround: true,
}
});
}
/**
* 创建平面
* @param {point} 平面的中心点
* @param {normal} 两点间法向量
*/
const GeneratePlane = (point, normal) => {
let plane2 = new Cesium.Plane.fromPointNormal(point, normal);
let plane1 = new Cesium.Plane(normal, 0);
let plane = new Cesium.Plane(Cesium.Cartesian3.UNIT_X, 0.0)
let entity = viewer.entities.add({ //Cesium.Math.toRadians(1)
position: point,
orientation: Cesium.Transforms.headingPitchRollQuaternion(
point,
new Cesium.HeadingPitchRoll(Cesium.Math.toRadians(-65), 0, 0)//平面手动纠偏
),
plane: {
plane: new Cesium.CallbackProperty(() => plane, false),
// palne: planeOrientation,
dimensions: new Cesium.Cartesian2(45.0, 25.0), // Size of the plane
outline: true,
material: Cesium.Color.LIGHTGRAY.withAlpha(0.7), // 墙体材质与透明度
outlineColor: Cesium.Color.YELLOW.withAlpha(0.3), // 轮廓颜色
rotation: Cesium.Math.toRadians(90)
},
});
}
创建完平面之后发现，平面与向量并不是垂直的，使用cesium提供的好几种方案都无效，最后只能手动纠偏。下图为实现的效果：
4.能够生成断面了，接下来就是最关键的一步了，如何在断面上计算填方挖方的面积呢？这里也是最让我头疼的地方，在ceisum以及它的生态中并没有直接提供计算垂直断面面积的api，只能一步步来实现了。如上面图2所示，计算面积需要获取设计断面线，也就是模型与断面的交线。现在已知断面双端的点坐标，最终决定运用cesium提供的viewer.scene.sampleHeight方法对两个点之间的模型进行高度提取，这个方法支持输入一个位置，会给你返回当前位置的高度。
使用时需要注意：1.不要加地形，会影响你的判断，不加地形会返回null，加地形之后就无法分辨哪些点是模型上的点了 2.其他的entity、BIM、primitives都会影响你的返回值，例如已经加载了这个断面，回返回断面的高度。 3.这个方法支持屏蔽某些entity，接收一个entity数组参数，在这个数组内的entity会被屏蔽掉，不会被取到高度。
具体实现的方案就是，对断面两端点之间得连线进行遍历，每隔1m取一个点，并且测量点的高度，只保存有返回值的点，这样能保证全部点都在BIM上面，最后返回模型点的集合。
这里比较简单，代码就不贴了。效果如下图所示
5.模型点已经提取完毕了，下面还需要提取原始数据的点，那么如何在点云数据提取断面相交的点呢?研究了很久，发现只能遍历点云的每个点数据，然后计算这个点位到平面的距离，点云有几百万个点，这样显卡肯定是吃不消的。这时只能有请后端大哥出场了，使用python去计算这些点，然后返回每个断面提取到的点的数组。这时我们就能拿到全部断面上原始数据的点坐标了。拿到坐标之后，发现只需要在河道正上方的点，所以需要除掉河道两边多余的原始点。这里也比较简单就不贴了，如果有需要的可以给我私信。
注意注意~， 加载原始数据与验收数据点的时候，一定要使用primitive形式去加载，不要一味entity，不理解的查一下它俩的原理，点位太多会导致内存溢出！！
6.现在我们拿到了同一平面下的原始数据的线与设计断面的线，接下来就是计算这两条线的焦点，并且需要判断出是填方还是挖方面积(设计线在上则为填方，反之为挖方)。话不多说，上代码，具体逻辑如下，看不懂就问问gtp~~
// 主函数：计算两个线段集合的交点，并计算填方和挖方面积
//line1设计数据点集合，line2原始数据点集合 viewer视图，其余参数为记录断面entity以及展示各项数据所用
//判断填挖方是判断交线最低点，
export default async function intersection(line1, line2, viewer, type, entitylist, station) {
line1.unshift(line2[0]);
line1.push(line2[line2.length - 1]);
let line1index, line2index;//记录交点位置
let fillArea = 0,
cutArea = 0; //累加填挖方面积
let beforeinter; //记录交点坐标
const line1Clone = lodash.cloneDeep(line1);
const line2Clone = lodash.cloneDeep(line2);
// 迭代所有线段，找到交点并计算面积
for (let i = 0; i < line1.length - 1; i++) {
for (let j = 0; j < line2.length - 1; j++) {
const intersection = getLineIntersection3D(line1[i], line1[i + 1], line2[j], line2[j + 1]);
if (intersection) {
if (i > 0 && beforeinter) {
//考虑特殊情况，在BIM模型相邻两点间可能有多个交点，此时需要相邻交点作为BIM线高度
let linefirst = [beforeinter, ...line1Clone.slice(line1index, i + 1), intersection]
let lineSecond = [beforeinter, ...line2Clone.slice(line2index, j + 1).reverse(),
intersection
]
// 判断两条线的最低点
let minHeightFirst = Number(Cesium.Cartographic.fromCartesian(linefirst[1]).height)
let minHeightSecond = Number(Cesium.Cartographic.fromCartesian(lineSecond[1]).height)
if (linefirst.length > 2) {
for (let i = 1; i < linefirst.length - 1; i++) {
let height = Number(Cesium.Cartographic.fromCartesian(linefirst[i]).height)
let difference = minHeightFirst - height
if(difference < 0){
minHeightFirst = height
}
}
} else {
let difference = minHeightFirst - Number(Cesium.Cartographic.fromCartesian(linefirst[0]).height)
if(difference < 0){
minHeightFirst = Number(Cesium.Cartographic.fromCartesian(linefirst[0]).height)
}
}
if (lineSecond.length > 2) {
for (let i = 1; i < lineSecond.length - 1; i++) {
let height = Number(Cesium.Cartographic.fromCartesian(lineSecond[i]).height)
let difference = minHeightSecond - height
if(difference < 0){
minHeightSecond = height
}
}
} else {
let difference = minHeightSecond - Number(Cesium.Cartographic.fromCartesian(lineSecond[0]).height)
if(difference < 0){
minHeightSecond = Number(Cesium.Cartographic.fromCartesian(lineSecond[0]).height)
}
}
const arr = [beforeinter, ...line1Clone.slice(line1index, i + 1), intersection, ...
line2Clone.slice(
line2index, j + 1).reverse()
];
let center = Cesium.BoundingSphere.fromPoints(arr).center;
let result = await computeArea(arr)
let differenceAB = Number(minHeightFirst) - Number(minHeightSecond)
if (differenceAB < 0) {
// 计算填方面积
fillArea += result.area;
let label = result.area;
} else if (differenceAB > 0) {
// 计算挖方面积
cutArea += result.area;
let label = result.area;
}
if (i == line1.length - 2 && j == line2.length - 2) {
return ({
fillarea: fillArea,
cutarea: cutArea
});
}
}
beforeinter = intersection;
line1index = i + 1;
line2index = j + 1;
}
}
}
}
// 计算3D空间中两条线段的交点
function getLineIntersection3D(P1, P2, Q1, Q2) {
const u = Cesium.Cartesian3.subtract(P2, P1, new Cesium.Cartesian3());
const v = Cesium.Cartesian3.subtract(Q2, Q1, new Cesium.Cartesian3());
const w = Cesium.Cartesian3.subtract(P1, Q1, new Cesium.Cartesian3());
const a = Cesium.Cartesian3.dot(u, u);
const b = Cesium.Cartesian3.dot(u, v);
const c = Cesium.Cartesian3.dot(v, v);
const d = Cesium.Cartesian3.dot(u, w);
const e = Cesium.Cartesian3.dot(v, w);
const denominator = a * c - b * b;
// 检查是否平行
if (Math.abs(denominator) < 1e-12) {
console.log('存在平行线段')
return null; // 线段平行，无交点
}
// 计算交点参数
const sc = (b * e - c * d) / denominator;
const tc = (a * e - b * d) / denominator;
// 检查交点是否在线段范围内
if (sc < 0 || sc > 1 || tc < 0 || tc > 1) {
if (Cesium.Cartesian3.equals(P1, Q1)) {
return Cesium.Cartesian3.clone(P1); // 或者返回Q1/Q2/P2，取决于哪个点重合
} else if (Cesium.Cartesian3.equals(P2, Q2)) {
return Cesium.Cartesian3.clone(P2);
}
return null; // 交点不在线段范围内
}
// 计算交点位置
return Cesium.Cartesian3.add(P1, Cesium.Cartesian3.multiplyByScalar(u, sc, new Cesium.Cartesian3()), new Cesium
.Cartesian3());
}
7.接下来就是计算两条线围城面积了，这里是最头疼的一步，不仅要判断是填方还是挖方，还要找出交点，最重要是如何计算面积。期间我实验了好多方法，要么没法实现，要么有缺陷，要么计算不准确。
最终通过两个方案实现了面积的计算，第一就是通过earcut.js对围成的面进行三角剖分，它的返回值是一个以三个顶点为一组的数组，格式[顶点1，顶点2，顶点3，顶点4，顶点5，顶点6]，每三个相邻点为一个三角形的顶点，原理是把一个平面划分为若干个三角形，并返回这些三角形的顶点坐标。如下图所示：
通过这个方法可以拿到线段围成断面的若干三角形顶点，只要计算出这些三角形的面积并加起来，就是断面面积的和。
通过顶点计算面积的代码如下（具体含义自己去查，用法如下）:
// 计算多边形的面积，采用三角剖分法
function calculatePolygonArea(positions) {
let area = 0;
const arr = positions.flatMap(cartographicToWgs84);
const indices = earcut(arr, null, 3); //耳切法
//earcut方法有时会获取不到，则需要使用其他方法计算
if (indices.length === 0) {
const lengthline = [];
const degrees = [];
const positionsCopy = lodash.cloneDeep(positions);
// 计算每条边的长度和夹角
for (let i = 0; i < positions.length - 1; i++) {
lengthline.push(computeLengthFromPoint(positions[i], positions[i + 1]));
degrees.push(computeDegreesFromLine(positionsCopy[i], positionsCopy[i + 1], positionsCopy[i + 2]));
}
// 计算多边形面积
return calculateVertices(lengthline, degrees);
}
// 循环通过三角剖分索引创建三角形，并计算每个三角形的面积
for (let i = 0; i < indices.length; i += 3) {
area += calculateTriangleArea(positions[indices[i]], positions[indices[i + 1]], positions[indices[i + 2]]);
}
return Number(area.toFixed(2));
}
// 根据边长和角度计算多边形的顶点坐标
function calculateVertices(sides, angles) {
const vertices = [{
x: 0,
y: 0
}]; // 初始点位于原点
let currentAngle = 0;
for (let i = 0; i < sides.length; i++) {
const lastVertex = vertices[vertices.length - 1];
vertices.push({
x: lastVertex.x + sides[i] * Math.cos(currentAngle),
y: lastVertex.y + sides[i] * Math.sin(currentAngle),
});
currentAngle -= angles[i];
}
return calculateArea(vertices);
}
// 计算多边形面积（使用向量叉积）
function calculateArea(vertices) {
return Math.abs(vertices.reduce((acc, vertex, i, arr) => {
const nextVertex = arr[(i + 1) % arr.length];
return acc + vertex.x * nextVertex.y - nextVertex.x * vertex.y;
}, 0) / 2);
}
// 计算两点之间的距离
function computeLengthFromPoint(startCartesian, endCartesian) {
return Cesium.Cartesian3.distance(startCartesian, endCartesian).toFixed(2);
}
// 计算两个向量之间的夹角
function computeDegreesFromLine(Cartesian1, Cartesian2, Cartesian4) {
const direction1 = Cesium.Cartesian3.subtract(Cartesian1, Cartesian2, new Cesium.Cartesian3());
const direction2 = Cesium.Cartesian3.subtract(Cartesian2, Cartesian4, new Cesium.Cartesian3());
return Cesium.Math.toDegrees(Cesium.Cartesian3.angleBetween(direction1, direction2));
}
// 计算三角形的面积
function calculateTriangleArea(p1, p2, p3) {
const v1 = Cesium.Cartesian3.subtract(p2, p1, new Cesium.Cartesian3());
const v2 = Cesium.Cartesian3.subtract(p3, p1, new Cesium.Cartesian3());
const crossProduct = Cesium.Cartesian3.cross(v1, v2, new Cesium.Cartesian3());
return Cesium.Cartesian3.magnitude(crossProduct) / 2;
}
// 将 Cartesian3 转换为 WGS84 坐标
function cartographicToWgs84(cartographic) {
const cartographics = Cesium.Cartographic.fromCartesian(cartographic);
return [
Cesium.Math.toDegrees(cartographics.longitude),
Cesium.Math.toDegrees(cartographics.latitude),
cartographics.height
];
}
最后也没用到这个方法，因为一个河道有大量断面，计算时间往往需要页面卡住几分钟，但是有了思路就好办了，让后端用python去实现了。为了后期优化，对计算数据做了保存，如果能拿到数据，就不需要在进行计算了，就是说，只有第一次导入数据时需要计算一遍。这样可以提前操作，用起来就不会卡顿了。
结语
最终历经千辛万苦终于是实现了，后续又做了许多算法的优化，以及联动效果。回头看，轻舟已过万重山~
最后看一下实现的效果吧
创作不易，请勿搬运！！！如果对您有帮助，请三连支持一下(*^_^*)~~
DAMO开发者矩阵，由阿里巴巴达摩院和中国互联网协会联合发起，致力于探讨最前沿的技术趋势与应用成果，搭建高质量的交流与分享平台，推动技术创新与产业应用链接，围绕“人工智能与新型计算”构建开放共享的开发者生态。
更多推荐
19
0- 0
已为社区贡献1条内容
所有评论(0)