一、问题简述
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比对不仅面积值
相差小，且基本重合，可知思路正确。但可供测试的超大顶点模型构件数量有限，
且测试构件复杂度不足，仍需获取更多复杂构件进行进一步测试并优化方案。
