一、問題簡述
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比對不僅面積值
相差小，且基本重合，可知思路正确。但可供測試的超大頂點模型構件數量有限，
且測試構件複雜度不足，仍需擷取更多複雜構件進行進一步測試并優化方案。