超大顶点模型构件投影优化

一、问题简述

BIM审核项目中 在建筑模型数据入库后 经常有模型构件进行投影的场景 模型投影调用超图Supermap iObjects Java提供的GeoModel3D.converToRegion()接口。而当对某些复杂的模型构件 节点数目过多 进行投影操作时 运算效率较低 有时超过30分钟。经过项目组的优化 有一定的优化效果 但是并没有彻底改善这个问题。

如下图所示的来源于某桥梁模型的构件 以三角形网格表达 该模型共有173926个节点和346459个面。

成功运行构件投影测试脚本 在windows虚拟机上运行时间大约为16分钟 合960s。

二、问题分析2.1现象观察

从这个图像中 我们有下述两点观察

①真正对投影面积有较大影响的底面结构比较简单 可以清晰地看到网格。

②所占数据较大的 其实是图中的圆柱体结构。如果利用超图接口进行计算 输入数据是进行过坐标转换到WGS84的三维坐标数据。耗时的方法是GeoModel3D.convertToRegion()方法。该方法执行时间超过16分钟 通过该方法计算后并转换回投影坐标系中计算出来的面积为5339370351平方毫米。
针对这一问题 首先可以尝试寻找超图接口的替代方法 除此之外 也可尝试探索对原有网格进行简化后再进行算法的方法。

2.2方案选定

经过多方资料查询 尝试从网格简化和shapely库地理计算的两个角度出发解决该问题。

2.2.1基于Shapely的自行实现

经过讨论 技术团队认为本问题可以抽象为两个步骤

①第一步将三维模型中的所有节点都投影至XY平面 并完整的保留三角形网格 此时三角形网格会在同一个XY坐标发生重叠。

②第二步针对超多数量的三角形网格进行一元聚合(unary_union) 即在三角形重合的地方进行合并。合并之后即可以调用聚合后的二维图形的面积接口 计算出投影面积。

具体实现上 第一步比较简单 直接忽略z轴坐标 并按照三角形网格的相同接口组成polygon即可以 第二步利用到的则是python shapely库的unary_union接口 该接口是GEOS库中UnaryUnionOp类的简单封装。

面积的计算结果为5340510605平方毫米 与超图的计算结果仅有1.14平方米的误差。鉴于超图的计算过程中从原始数据进行了坐标转换后又转换回来 怀疑是这个过程中有精度损失。

在同一台机器的相同的计算环境下 计算时间则大幅度减少到158s 仅仅是原来的16.5%。

2.2.2网格简化算法

经过调研 学术界在网格简化问题上影响力比较大的是卡内基梅隆大学在20世纪末提出的基于Quadric Error Metrics的算法。16年发布的开源软件MeshLab基于VCG Library对该算法有一个实现(功能名称为Quadric Edge Collapse Decimation) 因此可以借助该软件对简化后再进行计算的思路进行验证。而且MeshLib提供了cmdline的接口 从而使得该算法能力能够简单地被嵌入系统 或者可以直接包装VCG Library的相关接口 接入该算法能力。

利用该算法将原有34万 面的三维模型分别简化为16万、8万、4万、2万、1万个面 从图像上的直观上看几乎没有区别。技术团队也针对简化后的模型计算了投影面积 并测量了计算时间 总结如下表

顶点数量 ???面数量 ???面积/平方毫米 ?计算时长/秒 ?

原始数据?????173926 ????346459 ? ?5340510605 ? ? 158 ? ?????已验证

1/2简化??????80630 ????159992 ? ?5340510605 ? ? 119 ?? ?????已验证

1/4简化?? ? ??40584 ?????80000 ? ?5340510605 ? ? ?64 ???????已验证

1/8简化?? ? ??20518 ?????39999 ? ?5340510604 ? ? ?24 ???????已验证

1/16简化?? ??10391 ?????20000 ? ?5340510885 ? ? ?10 ????????已验证

1/32简化?? ? ??5289 ??????9999 ? ?5340510877 ? ???4 ?? ?? ????已验证

从结果中可以看出 即使简化为1万个面 计算出来的面积与原始数据的偏差只有200平方毫米。

2.2.3初步结论

①针对该问题 通过采用比较成熟的开源计算方案 计算效率可以得到大幅度提高 从16分钟缩短到3分钟之内。

②网格简化算法能进一步加速计算 最快4s可以得到结果 误差微乎其微。

三、接口设计思路

希望以新接口替代超图的GeoModel3D.convertToRegion() 所以设计入参 GeoModel3D return GeoRegion

实现路径 GeoModel3D-- 骨架(顶点串及序列)-- meshlab 接口接受的数

据结构(off)-- 简化网格-- 节点降维-- 重合处三角合并-- 转换GeoJSON ??-- GeoRegion

三、代码实现3.1 获取三维模型构件点集

获取超图 Model 对象的顶点集合与序列并输出 off 格式文件的方法参考以下代码

/**

* 将 Model 对象转换为 off 格式并在指定位置输出

*

* param model Model 对象

* param originalOffPath ?off 输出路径

*/

private void outputOffFile(Model model, String originalOffPath) {

List Double verticesList new ArrayList ();

List Integer vertexIndexesList new ArrayList ();

//获取 GeoModel3D 顶点点串与序列

for (int index 0; index model.getSkeletonCount(-1); index ) {

SkeletonID id new SkeletonID(-1, index);

Skeleton skeleton model.getSkeleton(id);

double[] vertices skeleton.getVertices();

List Double tempVerticeList Doubles.asList(vertices);

verticesList.addAll(tempVerticeList);

int[] vertexIndexes skeleton.getVertexIndexes();

List Integer tempIndexlist Ints.asList(vertexIndexes);

vertexIndexesList.addAll(tempIndexlist);

}

//输出 originalOff.off

writeOff(verticesList, vertexIndexesList, originalOffPath);

}

注意 此处的代码 getSkeletonCount()和 SkeletonId()方法中的 LOD 层级取值-1 代表当前 LOD 层级 此种取值是约定俗成的。

3.2 Windows环境调用meshlabserver方法

我们应用meshlab的Quadric Edge Collapse Decimation方法对模型骨架转化得到的off文件进行模型简化。为将这部分工作做代码实现 我们通过Meshlab的.mlx脚本 保存对原数据的操作 然后通过调用meshlabserver进行批量处理。

3.2.1 meshlab方法参数研究

在 做 三 维 模 型 简 化 时 主 要 使 用MeshLab的Quadric Edge Collapse Decimation方法 其参数较多 经过试验 建议以下参数配比

其中 需要注意的参数有以下几个

①Target number of faces 目标面数。简化目标面的数量 该参数决定了简化的程度 软件中针对特定模型进行手动操作尚可 代码实现需要频繁改变脚本参数 不适用。

②Percentage reduction (0..1) 简化百分比。默认为0 当设置不为0时 简化面数是以该参数为主 替代Target Number of Faces 如此可以跳过动态编辑mlx脚本 减少运行时间。且根据经验 最好设置为0.5的整数次幂。

③Preserve Topology 保留拓扑关系。因为模型可能存在复杂的拓扑关系 简化后的模型进行投影后 希望保留其应有的岛洞关系 因此此处勾选。

其余参数 默认值即可。

--

3.2.2 meshlabserver命令行执行

如图为meshlabserver在Windows上的cmd调用实现

cmd 语法为

meshlabserver 路径 -i 原始 off 路径 -o 简化输出 off 路径 -m vc vn -s mlx

脚本路径

如此得到简化off后 python读取其中的顶点位置与序列信息 进行降维与投影合并。

3.3 python库导入

为实现三维骨架点串的节点降维、重合处三角合并、聚合后二维图形面积获取等阶段性目标 需要导入若干python库 首先安装python 3.8版本即可。

导入shapely库 使用unary_union()做投影面合并 导入geojson库 实现

shapely.Polygon转化为GeoJSON文件。使用pip install命令安装以上库即可。

①?shapely

②?geojson

3.4 off的降维合并与输出

简化后的off python读取其点串信息 降维构件二维面 并调用shapely

库的unary_union()方法实现面的合并 代码如下

?

def get_projection_area(vertices, triangles):

polygons []

for i1, i2, i3 in triangles:

p1 vertices[i1]

p2 vertices[i2]

p3 vertices[i3]

poly Polygon(((p1[0], p1[1]), (p2[0], p2[1]), (p3[0], p3[1])))

if not poly.is_valid:

continue

polygons.append(poly)

print( polygon contructed )

unioned unary_union(polygons)

print( unioned )

print( Area: %f % unioned.area)

return unioned

求得合并投影面后 转换为GeoJSON格式 以供构建超图GeoRegion使用。

代码如下

def write_geojson(polygon, output_path):

#p wkt.loads((str)(polygon))

# using geojson module to convert from WKT back into GeoJSON format

geojson_out geojson.Feature(geometry polygon, properties {})

with open(output_path, w ) as outfile:

json.dump(geojson_out, outfile, indent 3)

outfile.close()

3.5 超图GeoRegion构建与坐标转换

由于3.1中得到Model及骨架中的点集不含坐标系信息 所以在此基础上进行一系列计算得到的结果GeoRegion也是原点在o处的模型的投影面。因此合并投影面GeoJSON在转换回超图的GeoRegion对象后 需要根据其原始的空间矩阵做矩阵转换 恢复位置、比例、姿态 赋予空间信息。矩阵转换代码如下

/**
?* 对Model骨架点串简化投影合并得到的GeoRegion进行转化 使其恢复位置、比例、姿态
?*
?* param geoRegion ???????网格简化与unary_union()后获得的结果GeoRegion对象
?* param transformParams 矩阵转换所需的参数
?* return ??????????????????转换后的GeoRegion
?*/
private GeoRegion transformGeoRegion(GeoRegion geoRegion, TransformParams transformParams) {
????GeoRegion resGeoRegion new GeoRegion();
????if (null geoRegion || geoRegion.isEmpty()) {
????????return resGeoRegion;
????}
????try {
????????for (int i 0; i geoRegion.getPartCount(); i ) {
????????????Point2Ds point2Ds geoRegion.getPart(i);
????????????Point2Ds point2DsNew new Point2Ds();
????????????for (int j 0; j point2Ds.getCount(); j ) {
????????????????Point2D point2D point2Ds.getItem(j);
????????????????Double x point2D.getX();
????????????????Double y point2D.getY();
????????????????// 1 缩放
????????????????Double xRes1 x * transformParams.getScaleX();
????????????????Double yRes1 y * transformParams.getScaleY();
????????????????// 2 旋转
????????????????//绕x轴
????????????????Double yRes2 yRes1 * Math.cos(transformParams.getRotateX());
????????????????//绕y轴
????????????????Double xRes2 xRes1 * Math.cos(transformParams.getRotateY());
????????????????//绕z轴
????????????????Double xRes3 xRes2 * Math.cos(transformParams.getRotateZ()) - yRes2 * Math.sin(transformParams.getRotateZ());
????????????????Double yRes3 yRes2 * Math.cos(transformParams.getRotateZ()) - xRes2 * Math.sin(transformParams.getRotateZ());

????????????????// 3 平移
????????????????Double xRes4 xRes3 transformParams.getOffsetX();
????????????????Double yRes4 yRes3 transformParams.getOffsetY();

????????????????Point2D point2dNew new Point2D(xRes4, yRes4);
????????????????point2DsNew.add(point2dNew);
????????????}
????????????resGeoRegion.insertPart(i, point2DsNew);
????????}
????} catch (Exception e) {
????????e.getMessage();
????}
????return resGeoRegion;
}

四、测试

将目前的方法与超图convertToRegion()方法计算结果对比 在表中取形状复杂的构件进行测试 分别比较各自投影面的面积相似度与空间覆盖率 结果较为理想。

测试结果

最终结果与超图convertToRegion()方法获取的GeoRegion比对不仅面积值

相差小 且基本重合 可知思路正确。但可供测试的超大顶点模型构件数量有限

且测试构件复杂度不足 仍需获取更多复杂构件进行进一步测试并优化方案。

本文转自网络，原文链接：https://developer.aliyun.com/article/784555
本站部分内容转载于网络，版权归原作者所有，转载之目的在于传播更多优秀技术内容，如有侵权请联系QQ/微信：153890879删除，谢谢！