張興巖 李 琦 梁 棟 蒲 潔
1 河北工業(yè)大學(xué)電子信息工程學(xué)院,天津市西平道 5340 號(hào), 300401 2 河北工業(yè)大學(xué)土木與交通學(xué)院,天津市西平道 5340 號(hào), 300401
三維激光掃描儀在獲取橋梁幾何空間信息方面具有高效率、高精度、非接觸等獨(dú)特優(yōu)勢(shì)[1]。由于采集到的橋梁點(diǎn)云數(shù)據(jù)量較大,不利于開(kāi)展后續(xù)橋梁幾何參數(shù)的測(cè)量分析和模型重建,因此需要將橋梁點(diǎn)云分割為具有語(yǔ)義信息的離散點(diǎn)集。
點(diǎn)云是包含三維坐標(biāo)信息的離散點(diǎn),體素包含若干離散點(diǎn),超體素面片包含若干體素(簡(jiǎn)稱“面片”),融合區(qū)域包含若干面片。點(diǎn)云分割可應(yīng)用于植物葉片分割[2]、地面分割[3]、屋頂分割[4]、建筑物立面分割[5]和室內(nèi)語(yǔ)義分割[6]等領(lǐng)域,主要方法包括隨機(jī)采樣一致性算法[7]、聚類法[8]、區(qū)域生長(zhǎng)法[9]和機(jī)器學(xué)習(xí)法[10]。隨機(jī)采樣一致性算法對(duì)點(diǎn)云的噪聲點(diǎn)和異常點(diǎn)表現(xiàn)出更好的魯棒性,但過(guò)于依賴模型系數(shù);聚類法能根據(jù)點(diǎn)云之間的距離、顏色等屬性分組,空間聚集能力強(qiáng),但對(duì)物體邊界區(qū)分能力較弱;區(qū)域生長(zhǎng)法生長(zhǎng)結(jié)果的好壞取決于種子點(diǎn)的選取和相鄰點(diǎn)區(qū)域的生長(zhǎng)條件,對(duì)于橋梁的細(xì)分割很容易越過(guò)邊界,造成錯(cuò)誤分割;機(jī)器學(xué)習(xí)的分割算法能得到很好的語(yǔ)義分割效果,但訓(xùn)練時(shí)間較長(zhǎng)且前期需要大量的數(shù)據(jù)集。對(duì)于建筑物點(diǎn)云的分割,Xu等[11]提出一種結(jié)合體素結(jié)構(gòu)和區(qū)域增長(zhǎng)的策略,利用平滑度、連續(xù)性和凹凸性作為幾何線索,通過(guò)區(qū)域生長(zhǎng)法對(duì)屋頂體素進(jìn)行融合,但邊界的檢測(cè)仍存在困難;董震等[12]提出一種融合顏色、反射強(qiáng)度和空間距離的多尺度超體素方法,利用圖像分割法分割點(diǎn)云,但需要點(diǎn)云本身帶有多種特征,不適用于只包含空間特征的點(diǎn)云分割。
綜上所述,本文利用橋面在橋梁結(jié)構(gòu)中的位置特點(diǎn)和超體素面片的局部特性,基于橋梁的法線特征、曲率、空間位置和點(diǎn)云數(shù)量占比,提出一種鄰接區(qū)域平面元融合的橋面分割方法。實(shí)驗(yàn)證明,該方法的分割效果較好,可實(shí)現(xiàn)橋梁點(diǎn)云中曲型橋面區(qū)域的分割。
本文方法主要流程見(jiàn)圖1,圖中每個(gè)處理過(guò)程的基本單位都包括點(diǎn)云、體素、超體素面片、融合區(qū)域,其中超體素面片和融合區(qū)域的最小組成單元都是體素。
圖1 橋面分割方法流程
超體素過(guò)分割的最小單元是體素,因此需要先對(duì)點(diǎn)云進(jìn)行體素化處理。體素化是指給所有點(diǎn)云建立一個(gè)軸對(duì)齊包圍盒,以體素分辨率Rvoxel為單位軸向劃分為多個(gè)小包圍盒,用小包圍盒代表其內(nèi)部的點(diǎn)云,該包圍盒即為體素,后續(xù)計(jì)算過(guò)程中用體素的質(zhì)點(diǎn)坐標(biāo)表示體素的空間位置。在3×3×3的空間內(nèi)一個(gè)體素周圍最多有26個(gè)鄰接體素,因此使用26鄰接關(guān)系建立體素間的鄰接圖。
體素可有效減少后續(xù)算法的計(jì)算量,但難以代表一塊區(qū)域的局部特征。超體素過(guò)分割能將具有相似屬性的體素聚類成一塊面片,相較于體素分割,其表示局部特征的能力更強(qiáng)。超體素過(guò)分割的步驟如下:
1)以種子分辨率Rseed為單位建立若干初始種子。
2)以Rsearch為搜索半徑,計(jì)算種子范圍內(nèi)體素的數(shù)量,過(guò)濾掉小于閾值的種子,去除孤立點(diǎn)。
3)以種子為中心逐層向外搜索,使用局部K均值聚類迭代生長(zhǎng)法更新聚類中心,計(jì)算鄰接體素與聚類中心所有體素之間相似性距離的平均值,將平均距離最近的鄰接體素加入聚類。
4)重復(fù)步驟3)直到無(wú)其他鄰接體素或達(dá)到最大搜索范圍,結(jié)束增長(zhǎng),獲得超體素面片。面片中心點(diǎn)由靠近超體素質(zhì)點(diǎn)的體素表示。
生長(zhǎng)過(guò)程中的相似性距離D由空間位置和法線特征的加權(quán)方程決定:
(1)
(2)
Dn=1-cos(ni,nj)
(3)
式中,ws為空間權(quán)值,Ds為體素間的空間距離,wn為法線權(quán)值,Dn為體素間的法線偏差,(xi,yi,zi)為第i個(gè)體素的位置,ni為第i個(gè)體素的法線,(xj,yj,zj)為第j個(gè)體素的位置,nj為第j個(gè)體素的法線,Rseed用于歸一化處理。
面片的屬性特征是后續(xù)融合條件的重要依據(jù),因此本文將面片的法線和曲率作為2個(gè)基本屬性。首先將面片內(nèi)所有體素的法線重定向,對(duì)所有體素法線加權(quán)融合獲得面片法線,然后計(jì)算面片內(nèi)所有體素的曲率平均值獲得面片曲率。同時(shí)為防止融合過(guò)程越過(guò)區(qū)域邊界,還需對(duì)面片進(jìn)行過(guò)濾,過(guò)濾指標(biāo)為面片法線的方向和模長(zhǎng),最終獲得包含橋面的候選面片。
將面片中每個(gè)體素的法向量以y軸正方向?yàn)閰⒖挤较蜻M(jìn)行重定向:
(4)
重定向后利用分層加權(quán)融合的方法計(jì)算面片的法向量。將面片分為中心點(diǎn)和外圍的L個(gè)面片層:
(5)
橋梁不同部位的面片法向量在模長(zhǎng)和方向上具有明顯差異,圖2為面片過(guò)濾的2種指標(biāo)。
圖2 面片過(guò)濾的2種指標(biāo)
第一個(gè)過(guò)濾指標(biāo)是面片法線模長(zhǎng)。由圖2(a)可見(jiàn),雖然橋面部分是曲面,但超體素過(guò)分割后的橋面部分面片曲率幾乎為0(Ⅰ類面片),設(shè)置的法線模長(zhǎng)閾值εm用于過(guò)濾非橋面面片(Ⅱ類面片)。
第二個(gè)過(guò)濾指標(biāo)是面片法線方向。由圖2(b)可見(jiàn),橋面(區(qū)域1)融合遇到邊界時(shí)法線變化不明顯,符合融合條件。越過(guò)區(qū)域1和區(qū)域2的邊界后,2個(gè)區(qū)域合二為一,會(huì)影響分割效果。由于橋面的面片法線和y軸正方向的夾角不超過(guò)橋面的最大坡度,因此設(shè)置法線夾角閾值εθ保留橋面區(qū)域、過(guò)濾邊界,達(dá)到面片過(guò)濾分層的效果。最終獲得的候選面片需滿足公式:
(6)
平面元融合是指將候選面片融合為多個(gè)區(qū)域,是橋面分割的關(guān)鍵步驟。雖然橋面同時(shí)具有凸曲面和凹曲面,但經(jīng)過(guò)過(guò)分割后每個(gè)面片對(duì)于橋面而言相當(dāng)于一塊小平面,稱為“面元”,使用多個(gè)面元表示曲面,對(duì)候選部分進(jìn)行平面元融合。
鄰接區(qū)域平面元融合可以獲取多塊具有幾何意義的橋梁區(qū)域,包括橋面區(qū)域以及橋梁的其他曲面區(qū)域。設(shè)置種子面片Sci和鄰接面片Bcj,輸入為初始面片區(qū)域A、面片法向量n和曲率γ、面片之間的鄰接關(guān)系索引表,輸出為多組融合區(qū)域B。面片融合的具體步驟如下:
1)為提高面片的增長(zhǎng)效率,選取所有曲率最小的面片作為初始種子。
2)從區(qū)域A中取出初始種子,將初始種子加入到一組區(qū)域Rc中,并加入種子隊(duì)列Sc。
3)開(kāi)始區(qū)域增長(zhǎng):
①?gòu)姆N子隊(duì)列Sc中取出隊(duì)頭的種子面片Sci;
②根據(jù)面片鄰接關(guān)系索引表獲得種子面片Sci的所有鄰接面片,存于鄰接區(qū)域Bc中;
③計(jì)算種子面片Sci與其鄰接面片Bcj的相似性度量Sij,若Sij大于相似性度量閾值Sth,則加入該組區(qū)域Rc中;若面片曲率γ大于曲率閾值γth,則將該鄰接面片Bcj加入種子隊(duì)列Sc的隊(duì)尾;
④清空鄰接區(qū)域Bc,返回步驟①,直到種子隊(duì)列Sc為空。
4)將增長(zhǎng)后的區(qū)域Rc添加至區(qū)域B中,清空區(qū)域Rc,返回步驟1),直到區(qū)域A為空。至此,平面元融合結(jié)束,獲得多組融合區(qū)域B。
上述方法中的相似性度量Sij由面片間的法線方向和中心點(diǎn)判決結(jié)果決定。法線方向用于確定面片之間的法線相似性Sn,在滿足法線相似性的情況下,面片夾角可能會(huì)出現(xiàn)銳角或鈍角的情況。由于橋面區(qū)域面片夾角為鈍角,因此中心點(diǎn)判決可用于排除面片夾角為銳角的情況。
根據(jù)余弦定理,面片法線之間的余弦值公式為:
(7)
式中,α為2個(gè)面片法線的夾角,ni為面片Sci的法線向量,nj為面片Bcj的法線向量。由于進(jìn)行了法線重定向,因此法線余弦值取值范圍為[0,1]。
中心點(diǎn)判決過(guò)程如圖3所示。首先計(jì)算種子面片Sci和鄰接面片Bcj之間的交線;然后從種子面片中心點(diǎn)M對(duì)交線作垂線,交點(diǎn)為M0,從鄰接面片中心點(diǎn)P對(duì)交線作垂線,交點(diǎn)為P0;最后根據(jù)向量和向量之間的夾角確定面片中心點(diǎn)判決結(jié)果,融合夾角為鈍角的情況,排除夾角為銳角的情況。判決公式為:
(8)
圖3 相鄰面片的夾角
受到體素分辨率的限制,超體素過(guò)分割后的面片間夾角不會(huì)出現(xiàn)極小銳角的情況。因此,當(dāng)2個(gè)面片法線方向幾乎相同時(shí),無(wú)需進(jìn)行中心點(diǎn)判決。設(shè)置面片法線間的夾角閾值αth=5°,計(jì)算公式為:
(9)
式中,k為鄰接面片中心點(diǎn)判決結(jié)果,Sn為面片間法向量夾角余弦值,α為2個(gè)面片法線的夾角。
為分割出最終的橋面點(diǎn)云,需要將以體素為最小單位的融合區(qū)域映射到以點(diǎn)云為最小單位的融合區(qū)域上,并計(jì)算每組融合區(qū)域點(diǎn)云數(shù)占原始橋梁點(diǎn)云數(shù)的比例,剔除占比不足1%的微小區(qū)域,剩余的融合區(qū)域包括橋面區(qū)域和非橋面區(qū)域。由于點(diǎn)云具有離散性,因此需要使用統(tǒng)計(jì)分布方法從區(qū)域中篩選出橋面點(diǎn)云。分析橋面在橋梁中的結(jié)構(gòu)和位置特點(diǎn)可知,橋面位于橋梁較高處,向地面的投影類似矩形,且投影面積較大。最終使用矩形相似性、相對(duì)區(qū)域高度和投影面積指數(shù)3個(gè)指標(biāo)來(lái)判斷融合區(qū)域是橋面點(diǎn)云的可能性。第n組區(qū)域的3個(gè)指標(biāo)計(jì)算步驟如下:
1) 橋梁點(diǎn)云y軸垂直于地面,y坐標(biāo)最小值為Ymin,統(tǒng)計(jì)融合區(qū)域點(diǎn)云的y坐標(biāo),計(jì)算平均值Ymean_n。
2) 融合區(qū)域的點(diǎn)云向地面xoz投影,計(jì)算x軸方向的最大值Xmax_n和最小值Xmin_n,計(jì)算z軸方向的最大值Zmax_n和最小值Zmin_n。
3) 對(duì)投影區(qū)域進(jìn)行正方形網(wǎng)格劃分,設(shè)置每個(gè)網(wǎng)格邊長(zhǎng)為Gsize,統(tǒng)計(jì)每個(gè)網(wǎng)格中投影點(diǎn)的數(shù)量,生成點(diǎn)云分布矩陣。
4) 將點(diǎn)云分布矩陣轉(zhuǎn)化為灰度圖像,為更直觀地展示投影形狀,使用大津法[13]設(shè)置閾值,進(jìn)行二值化處理,即圖像中大于閾值的像素設(shè)為1,最后生成二值圖。
5) 統(tǒng)計(jì)二值圖中值為1的像素?cái)?shù),記為投影面積Sp_n。使用旋轉(zhuǎn)主軸法[14]計(jì)算二值圖的最小外接矩形,統(tǒng)計(jì)最小外接矩形所占用的像素?cái)?shù),記為外接矩形面積Srect_n。
6) 3個(gè)指標(biāo)的計(jì)算公式為:
(10)
(11)
R′n=1.0-0.02(Rn-1)
(12)
式中,D矩_n為矩形相似性,Hmean_n為相對(duì)區(qū)域高度,Ymean_max為所有融合區(qū)域Ymean_n的最大值,R′n為投影面積指數(shù),Rn為第n組區(qū)域投影面積在所有融合區(qū)域中的名次。由于投影面積不是決定性因素,為減弱該指標(biāo)對(duì)結(jié)果的影響,將其與其他指標(biāo)處于同一度量下,進(jìn)行式(12)的計(jì)算。以上3個(gè)指標(biāo)的取值范圍都是[0,1],指標(biāo)越接近于1,是橋面點(diǎn)云的可能性越大。最終評(píng)價(jià)每一組融合區(qū)域的公式見(jiàn)式(13),通過(guò)比較每一組區(qū)域的Pn,篩選出最大值對(duì)應(yīng)區(qū)域,即橋面點(diǎn)云:
Pn=D矩_n·Hmean_n·R′n
(13)
實(shí)驗(yàn)平臺(tái)為Intel Core i5-9400F CPU@2.90GHz、VMware Workstation 15 Pro、Ubuntu18.04,開(kāi)源點(diǎn)云庫(kù)為PCL 1.11.1。為得到較為真實(shí)的實(shí)驗(yàn)數(shù)據(jù),對(duì)拱形橋、鋼筋橋和斜拉橋的3D模型進(jìn)行采樣,獲得模型點(diǎn)云。由于模型內(nèi)部也會(huì)被采樣,因此需要使用點(diǎn)云軟件Cloud Compare的Hidden Point Removal功能,從橋梁的7個(gè)視點(diǎn)模擬掃描出模型點(diǎn)云數(shù)據(jù)(圖4(a)),最后融合7組結(jié)果并刪除重復(fù)點(diǎn),獲得實(shí)驗(yàn)數(shù)據(jù)。以拱形橋點(diǎn)云作為實(shí)驗(yàn)對(duì)象,結(jié)果如圖4(b)所示。
圖4 拱形橋?qū)嶒?yàn)數(shù)據(jù)
首先對(duì)原始點(diǎn)云數(shù)據(jù)進(jìn)行體素化處理,既能保持橋梁外型輪廓不變,又能降低點(diǎn)云數(shù)據(jù)量、提高后續(xù)處理效率。體素分辨率為0.02 m,得到771 732個(gè)體素。然后進(jìn)行超體素過(guò)分割,種子分辨率為0.2 m,獲得18 627塊面片。接著以ey=(0,1,0)為基準(zhǔn)向量對(duì)面片的法向量作重定向處理,重定向后法向量與y軸的夾角范圍為[0,π/2]??紤]到曲型橋面的坡度,為面片法線偏離y軸預(yù)留π/6的閾值。同時(shí),為使橋面在候選面片中有更高的占比,設(shè)置法線模長(zhǎng)閾值為0.8,過(guò)濾法線模長(zhǎng)小于閾值的面片。圖5為面片過(guò)濾前后的效果對(duì)比,過(guò)濾后的候選區(qū)域剩余14 667塊面片,法線方向垂直于y軸的面片被過(guò)濾,部分橋墩和護(hù)欄也被過(guò)濾,分割后獲得較好的分層效果。
圖5 面片過(guò)濾
最后按照§1.3的步驟進(jìn)行鄰接面片的平面元融合,設(shè)置生長(zhǎng)條件的相似性度量閾值為0.8,加入種子隊(duì)列條件的曲率閾值為0.08,融合后共計(jì)4 546組區(qū)域,部分融合區(qū)域具有一定的語(yǔ)義信息,便于后續(xù)橋面的獲取。
圖6為融合區(qū)域的篩選過(guò)程,借助PCL點(diǎn)云庫(kù)中的getLabeledCloud函數(shù)將以體素為單位的融合區(qū)域映射到以點(diǎn)云為單位的融合區(qū)域上,計(jì)算每組區(qū)域點(diǎn)云數(shù)占橋梁點(diǎn)云數(shù)的比例,剔除不足1%的區(qū)域。圖6(b)中大部分微小區(qū)域被剔除,剩下5組融合區(qū)域,分別計(jì)算對(duì)應(yīng)的評(píng)價(jià)指標(biāo),得到Pn。詳細(xì)數(shù)據(jù)如表1所示,雖然1號(hào)區(qū)域在投影面積方面不占優(yōu)勢(shì),但另外2項(xiàng)指標(biāo)得分較高,因此1號(hào)區(qū)域的Pn最高,被認(rèn)為是橋面點(diǎn)云。結(jié)合圖6(c)可知,1號(hào)區(qū)域確實(shí)是橋面點(diǎn)云。
圖7 鋼筋橋和斜拉橋的分割實(shí)驗(yàn)
圖6 篩選過(guò)程
表1 5組區(qū)域的評(píng)價(jià)指標(biāo)
為進(jìn)一步驗(yàn)證本文方法的可行性和穩(wěn)定性,對(duì)鋼筋橋和斜拉橋進(jìn)行橋面分割實(shí)驗(yàn)(圖7),記錄下3種橋梁(1號(hào)為拱形橋,2號(hào)為鋼筋橋,3號(hào)為斜拉橋)的調(diào)試參數(shù)和結(jié)果(表2、表3)。由圖7可見(jiàn),在面片過(guò)濾階段,大部分鋼筋、繩索和橋墩等非橋面區(qū)域被過(guò)濾掉,使得后續(xù)融合區(qū)域更加準(zhǔn)確,最后統(tǒng)計(jì)各區(qū)域的點(diǎn)云分布,篩選出橋面點(diǎn)云。
表2 調(diào)試參數(shù)
表3 橋面點(diǎn)云的數(shù)量
由表2可見(jiàn),體素分辨率Rvoxel和種子分辨率Rseed的設(shè)置主要與橋梁點(diǎn)云的密度和結(jié)構(gòu)尺寸有關(guān),二者之間的關(guān)系可設(shè)置為Rseed≈(6~10)Rvoxel,Rvoxel和Rseed是影響分割時(shí)間的主要參數(shù),其幾何意義具有參考性,其他參數(shù)則相差不大。表3為橋面點(diǎn)云的數(shù)量,為直觀地展示本文方法的分割效果,使用Cloud Compare軟件的Closest Point Set功能和Remove Duplicate Points功能計(jì)算算法分割橋面的結(jié)果以及人工分割橋面結(jié)果的交并比。由表可見(jiàn),交并比均在95%以上。因此,本文方法對(duì)以上3種類型的橋面分割均具有較高的準(zhǔn)確性。
由于在融合橋面區(qū)域的過(guò)程中會(huì)出現(xiàn)法線變化不明顯的情況,從而導(dǎo)致融合的區(qū)域越界。為驗(yàn)證本文方法在特殊邊界融合區(qū)域上的優(yōu)勢(shì),采用區(qū)域生長(zhǎng)法、超體素區(qū)域生長(zhǎng)法[9]和本文方法進(jìn)行比較,并截取S3DIS(stanford large-scale 3D indoor spaces)數(shù)據(jù)集的墻面制作點(diǎn)云數(shù)據(jù),包括3塊區(qū)域(R1、R2、R3)和2處邊界(B1、B2),如圖8(a)所示。實(shí)驗(yàn)對(duì)象的數(shù)據(jù)無(wú)顏色信息,在MATLAB中使用漸變色展示其結(jié)構(gòu)特征。T為處理數(shù)據(jù)所用的時(shí)間。
圖8 對(duì)比實(shí)驗(yàn)
分析邊界B2的越界情況可知,在點(diǎn)或面片過(guò)濾前,3種方法都會(huì)將區(qū)域R1、邊界B2、區(qū)域R3融合到一起,這是由于邊界B2點(diǎn)與點(diǎn)、面與面之間的法線方向變化不明顯所致。在點(diǎn)或面片過(guò)濾后,點(diǎn)法線或面片法線過(guò)濾掉部分邊界B2,邊界B2被截?cái)?,區(qū)域R1和區(qū)域R3融合為單獨(dú)區(qū)域。
分析邊界B1的越界情況可知,在點(diǎn)或面片過(guò)濾后,仍然不能截?cái)噙吔鏐1,這是因?yàn)樵撨吔绲狞c(diǎn)法線或面片法線與y軸的夾角在閾值εθ范圍內(nèi),無(wú)法被過(guò)濾,最終導(dǎo)致區(qū)域R1和區(qū)域R2融合在一起。區(qū)域生長(zhǎng)法依靠點(diǎn)法線變化判斷生長(zhǎng)條件,在融合區(qū)域時(shí)無(wú)法分辨邊界B1的情況(圖8(b))。從圖8(a)的鄰接關(guān)系圖可以看出,區(qū)域R1和邊界B1交界處的面片具有鄰接關(guān)系,超體素區(qū)域生長(zhǎng)法的生長(zhǎng)條件只依靠相鄰面片中心點(diǎn)法線的變化,因此面片融合時(shí)同樣無(wú)法分辨邊界B1的情況(圖8(c))。本文方法在法線變化的基礎(chǔ)上加入中心點(diǎn)判決準(zhǔn)則,防止相鄰面片間夾角為銳角的情況出現(xiàn)(圖8(d)),區(qū)域R1和區(qū)域R2融合為單獨(dú)區(qū)域。
綜上所述,相較于2種傳統(tǒng)方法,本文方法使3塊區(qū)域單獨(dú)融合的同時(shí)不會(huì)越過(guò)邊界,在跨邊界融合方面具有較好的抑制作用。在速度方面,以面片為單位的處理速度要優(yōu)于以點(diǎn)為單位的處理速度。
本文提出一種鄰接區(qū)域平面元融合的橋面分割方法,主要用于分割建筑物點(diǎn)云中的大型曲面,可極大減少人工干預(yù)。相較于傳統(tǒng)分割方法,本文方法在融合過(guò)程中加入面片過(guò)濾和面片判決準(zhǔn)則,可進(jìn)一步避免融合區(qū)域越過(guò)邊界的情況出現(xiàn)。
雖然本文方法的閾值參數(shù)具有較好的通用性,但超體素過(guò)分割種子點(diǎn)的選擇仍然會(huì)影響過(guò)分割的效果,進(jìn)而導(dǎo)致后續(xù)融合區(qū)域的邊界較為粗糙。如何得到融合區(qū)域的清晰邊界,將是后續(xù)的研究重點(diǎn)。