工業(yè)CT重建模型與原始CAD模型粗配準的改進

段黎明，王 侃，陳 中

（1．重慶大學ICT研究中心，重慶 400030；2．重慶大學 機械工程學院，重慶 400030）

0 引言

在使用計算機輔助設計（Computer Aided Design，CAD）軟件設計產(chǎn)品，建立其原始CAD模型并加工成型后，為判斷加工成型后的產(chǎn)品在形狀和尺寸方面是否達到設計要求，需要分析加工成型后的產(chǎn)品和原始CAD模型之間的差別。所采用的方法是：首先對產(chǎn)品進行數(shù)字化掃描，再根據(jù)掃描獲取的數(shù)據(jù)重建其三維數(shù)字化模型，然后與其原始三維CAD模型進行比對分析。

目前對產(chǎn)品進行數(shù)字化掃描常用的技術(shù)有三坐標機測量、超聲波測量和工業(yè)計算機斷層掃描（Computed Tomography，CT）等。其中工業(yè)CT掃描技術(shù)能以圖像形式準確反映工件的內(nèi)外部情況，且不受工件復雜程度的影響，目前已成為檢測復雜產(chǎn)品（特別是含有內(nèi)腔的產(chǎn)品）的重要手段［1］。由于工業(yè)CT對產(chǎn)品進行掃描時有其獨立的坐標系，導致重建的三維數(shù)字化模型坐標系不同于產(chǎn)品原始三維CAD模型的坐標系，在對這兩個模型進行比對分析時，必須首先進行配準。

目前，配準領域多以迭代最近點算法（Iterative Closest Point，ICP）及其改進算法為主［2］。該算法對兩模型之間的初始位置要求較高，且很容易陷入局部極值，過早的收斂會導致配準失敗［3－4］。而近年提出的利用遺傳算法進行配準［5－6］的方法，雖然有較好的全局搜索能力，但耗時較長。無論采取哪種算法，都有必要進行粗配準，以減少兩模型之間的平移和角度偏差，在較短時間內(nèi)盡量使兩模型對齊，為精配準提供良好的初始條件。

常用的粗配準方法主要有三點對齊法、最小包圍盒法、遺傳算法和主成分分析（Principal Component Analysis，PCA）法等。三點對齊法［7］不僅引入了人工操作的誤差，不能保證配準結(jié)果，還無法實現(xiàn)自動配準；最小包圍盒方法［8］僅利用模型的外部輪廓數(shù)據(jù)，對內(nèi)部結(jié)構(gòu)不敏感，在配準由工業(yè)CT圖像重建的模型時，容易忽略工件內(nèi)部關(guān)鍵結(jié)構(gòu)，導致粗配準成功率較低；遺傳算法［9］耗時較長，適用于精配準領域的多次迭代運算而粗配準要求在較短時間內(nèi)完成的情形。

PCA方法［10］是一種正交化線性變換方法，它將數(shù)據(jù)變換到一個新的坐標系統(tǒng)中，構(gòu)成較少的不相關(guān)的新變量以代替原始變量，而每個新變量都包含盡量多的原始變量信息，新變量叫作原始變量的主成分或特征向量［11］。對于一個三維剛體，存在特定的表征其質(zhì)量分布的三個主軸，即三個主方向，而主軸就是慣量矩陣的特征向量。該方法利用兩模型點云具有相近的點分布特征進行配準，其經(jīng)典算法的思路為：計算數(shù)據(jù)模型的質(zhì)心及其慣量矩陣的特征向量，并利用質(zhì)心對齊和主軸對齊，得到平移向量和旋轉(zhuǎn)矩陣，從而完成配準。

文獻［12］將PCA方法具體化為力矩主軸法進行粗配準，該方法的優(yōu)點是運算量小，缺點是對數(shù)據(jù)缺失比較敏感。工業(yè)CT掃描圖像重建的三維模型數(shù)據(jù)完整，適宜采用該方法進行粗配準。文獻［13］使用PCA方法對點云數(shù)據(jù)進行粗配準，提到其反向問題，即由于PCA方法得到的主軸存在兩個方向，粗配準的結(jié)果可能是其中任意一個，不符合配準預期，易導致配準失敗。該文獻采用人工旋轉(zhuǎn)模型將其配準到正確方向。文獻［14］使用最小包圍盒來糾正PCA的反向問題，但同樣需要多次翻轉(zhuǎn)計算，時間空間效率較差。本文針對PCA方法的反向問題，以工業(yè)CT切片圖像三維重建模型的特點為基礎，提出基于模型形心和質(zhì)心矢量關(guān)系的PCA修正方法，避免主軸反向的發(fā)生。

1 重建模型與CAD模型配準問題的數(shù)學描述

利用工業(yè)CT對產(chǎn)品進行掃描時，需要根據(jù)產(chǎn)品的內(nèi)外結(jié)構(gòu)情況來確定裝夾方式，掃描坐標系不一定與原始CAD模型的坐標系一致，這樣根據(jù)掃描數(shù)據(jù)重建的三維數(shù)字化模型與其原始的三維CAD模型常常不會處于同一坐標系中。因此，需要對三維重建模型進行空間變換，使其與CAD模型在空間上對齊。

工業(yè)CT三維重建模型一般為立體光刻成型（Stereo Lithography，STL）格式的表面模型，它通過大量的三角面片擬合實體表面來表現(xiàn)模型的空間特征，為方便配準計算，將CAD模型離散為表面點云模型。本文假設兩模型均已轉(zhuǎn)化為不包含冗余數(shù)據(jù)的點云模型。

設CAD點云模型為點集矩陣P，工業(yè)CT掃描重建點云模型為點集矩陣Q，存在一種剛體變換，使得Q經(jīng)過該變換得到新點集矩陣Q′后，Q′與P在空間上能夠重合，即Q′上任意一點到P的距離最小。配準問題可依此描述為：求得一個空間剛體變換S，且滿足distance（S（Q），P）最?。?5］。根據(jù)計算機圖形學理論［16］，剛體變換S可由一個旋轉(zhuǎn)矩陣R和一個平移向量T來描述，

則新點集Q′＝Q·R＋T。

2 主成分分析配準方法的改進

經(jīng)典的PCA配準方法具體步驟如下：

（1）求兩模型的質(zhì)心GP和GQ。

（2）根據(jù)力學公式分別求出兩模型的慣量矩陣IP和IQ。

（3）計算主方向，主方向即為兩模型慣量矩陣IP和IQ的特征向量。由于IP和IQ均為實對稱矩陣，故其有三個正交的特征向量，即三個正交的主方向。

（4）計算旋轉(zhuǎn)矩陣R。將兩組特征向量按特征值的大小分別排列為三階方陣并單位化，得到EP和EQ，則可通過R·EQ＝EP求得旋轉(zhuǎn)矩陣R。

（5）計算兩模型間的平移向量T，T即為兩模型質(zhì)心間的矢量。

（6）將旋轉(zhuǎn)矩陣R和平移矩陣T應用到Q上，得到Q′，完成配準。

該經(jīng)典方法存在一個致命的缺陷，即計算出的兩模型慣量矩陣的特征向量沒有確定的正負號，導致模型的第二和第三主方向符號不確定，這將對后續(xù)旋轉(zhuǎn)矩陣R的求解產(chǎn)生致命影響。在某些情況下主方向會出現(xiàn)反向，求出R并配準后，Q′與P關(guān)于某軸線對稱，直接導致配準結(jié)果與預期嚴重不符。因此，有必要在PCA配準后判斷是否反向，以及出現(xiàn)反向的主方向向量，從而修正配準結(jié)果。

2．1 判斷是否反向

形心和質(zhì)心作為點集的固有特征，在正確配準時，兩點之間的相對位置關(guān)系不會隨空間剛體的變換而變化。當配準反向時，點云的形心與質(zhì)心的相對位置關(guān)系與正常情況反向，因此該位置關(guān)系便可作為一種是否反向的判斷標準。圖1所示為點集模型中各點和矢量在aGb平面的投影示意圖，其中兩段弧分別表示點集P和Q′，a和b是P的兩個主方向，G是P和Q′配準后共同的質(zhì)心，PMP和PMQ′分別是P和Q′的形心MP和MQ′在aGb平面的投影點，點集的形心和質(zhì)心的位置關(guān)系在aGb平面內(nèi)用矢量GPMP和GPMQ′表示。利用PCA配準正常時（如圖1a），由于配準誤差，GPMP與GPMQ′之間有一定夾角，但兩者在a和b方向上不存在方向相反的分量。只有在PCA配準產(chǎn)生反向時（如圖1b），兩矢量在a軸上的分量會近似地大小相等、方向相反，在b軸上的分量也一樣。因此，可以通過對比兩點集形心與質(zhì)心的位置關(guān)系，判斷出是否有反向情況發(fā)生，方法為首先計算由點集形心和質(zhì)心構(gòu)成的兩個矢量GMP和GMQ′，將它們分別投影在三個主方向a，b，c上，然后判斷兩矢量在三個方向上的分量，有如下三種情況：

（1）若兩矢量在其中一個主方向上的兩個投影分量同向，在另兩個主方向上的分量分別反向，則判定模型發(fā)生了反向，且同向的投影分量所在主方向為旋轉(zhuǎn)軸線。

（2）若兩矢量在兩個主方向上的投影分量同向，在另一個主方向上的兩投影分量反向，則判定配準模型不存在反向問題。

（3）若兩矢量在三個主方向上的投影分量分別同向，則判定配準模型不存在反向問題。

2．2 修正反向的點集

因為PCA配準反向相當于點集Q′繞一個主方向（用c表示）旋轉(zhuǎn)了180°，找出該主方向并逆向旋轉(zhuǎn)180°即可完成修正。如圖2所示，在由點集形心和質(zhì)心構(gòu)成的兩矢量GMP和GMQ′確定的平面內(nèi)，這兩個矢量關(guān)于c對稱，因此可將兩者在平面內(nèi)的對稱軸s作為逆向旋轉(zhuǎn)的軸。

將Q′繞s軸旋轉(zhuǎn)180°，即可完成對反向點集的修正。

2．3 算法步驟

假設已由PCA的前述步驟得到Q′及兩點集質(zhì)心。

步驟1 計算P和Q′的形心MP和MQ′。

步驟2 計算P的形心與質(zhì)心的距離，若小于一定閾值，則認為兩點重合，無法完成后續(xù)配準運算，退出；否則，繼續(xù)步驟3。

步驟3 計算兩點集質(zhì)心與形心的位置關(guān)系，用向量表示為u＝GMP，v＝GMQ′。

步驟4 求得兩向量在三個主方向上的分量，ua，ub，uc，va，vb，vc，分別判斷ua與va、ub與vb、uc與vc的方向，若有兩個主方向上的分量同向，則認為不存在反向情況，結(jié)束；若有一個主方向上分量同向，另兩方向上分量反向，則認為存在反向，轉(zhuǎn)步驟5。

步驟5 計算u和v之間的對稱軸矢量s，將Q′繞s旋轉(zhuǎn)180°，轉(zhuǎn)步驟3。

由此得到的模型點集Q′即為修正后的配準結(jié)果。

2．4 方法適用范圍

本文通過計算點云模型的形心與質(zhì)心位置關(guān)系來判斷是否有反向發(fā)生，該方法能夠正確配準的前提為點云模型的形心與質(zhì)心不重合。由于工業(yè)CT自身的貫穿性特點，其掃描的產(chǎn)品大多包含內(nèi)腔，且不規(guī)則，質(zhì)量在空間上分布不均勻，重建出的點云模型形心與質(zhì)心一般不重合，適合采用本文方法。形心與質(zhì)心重合的產(chǎn)品使用本方法計算時有可能產(chǎn)生錯誤，不屬于本文方法適用的范圍。

3 實驗驗證

在Windows XP系統(tǒng)中使用VC＋＋6．0開發(fā)了實驗軟件，計算機為主頻2．0GHz CPU，1 000 MB內(nèi)存，圖形開發(fā)包使用Open Inventor，版本為SIM公司的Coin3D。以下配準實驗均在該實驗軟件下完成。

3．1 實驗一

為驗證本文方法在PCA方法反向問題上的有效性，選取GE公司公布的Engine Block模型進行配準。對Engine Block切片進行三維重建得到點云模型，包括985 993個點，為方便觀察，導入軟件后對點云模型進行表面著色渲染，如圖3的灰色模型所示；對模型進行一定的空間變換后生成一新點云模型，如圖3的黑色模型所示。該發(fā)動機模型為腔體類工件，外形輪廓近似長方體，但內(nèi)部結(jié)構(gòu)復雜，點云的空間分布也不均勻，符合本文方法的適用范圍。首先對其使用經(jīng)典的PCA方法進行配準，配準結(jié)果如圖4所示，在灰色模型的上側(cè)和燃燒室部分存在明顯不匹配的情況，黑色模型由于反向問題而產(chǎn)生上下部分和前后部分的顛倒。使用本文方法配準后的結(jié)果如圖5所示，圖中模型被扭轉(zhuǎn)過來，上部和燃燒室部分都已基本對齊，說明已經(jīng)不存在方向性的配準錯誤，解決了PCA方法配準的反向問題。

3．2 實驗二

為驗證本文方法能滿足產(chǎn)品工業(yè)CT掃描重建模型與原始CAD模型粗配準的要求，選取某型號化油器進行工業(yè)CT掃描，將切片圖像進行三維重建得到其點云模型，包括241 964個點，如圖6中的黑色模型所示；將該化油器的原始CAD模型離散化為點云模型，包括258 044個點，如圖6中的灰色模型所示?；推髂Ｐ徒Y(jié)構(gòu)復雜，油路等表面特征無規(guī)則，點云的空間分布不均勻，滿足本文方法適用范圍。使用經(jīng)典PCA方法配準的結(jié)果如圖7所示，重建模型背面的法蘭盤被變換到了CAD模型正面的進氣孔部位。使用本文方法配準，結(jié)果如圖8所示，兩模型的進氣孔及上部的凸緣等其他部分已基本對齊，不存在方向性的配準錯誤，從而得到了正確的粗配準結(jié)果。

利用Metro［17］對兩點云模型進行點采樣，并求出兩模型的最大誤差為9．584 792mm，平均誤差為1．982 866mm，標準誤差為2．471 730mm。實驗表明，本文方法在完成模型方位正確配準的前提下，能夠?qū)⑴錅收`差保持在較小范圍，為后續(xù)精配準提供了較好的初始位置。

3．3 實驗三

為了對比本文方法與已有修正PCA反向的方法，通過在本文實驗軟件平臺上編寫程序，實現(xiàn)了文獻［14］使用的修正方法，并對如圖9所示的工件進行粗配準實驗。該工件包括51 770個點云。文獻［14］通過計算兩模型包圍盒的重合比例來判斷是否反向，若反向則通過翻轉(zhuǎn)模型繼續(xù)判斷，直到該比例超過閾值k。但該方法的比例閾值k的選取直接影響了配準的正確性，圖10所示為k分別取60%，80%和95%時該方法的配準結(jié)果?？梢钥闯?，當閾值k較低，而兩模型空間姿態(tài)相差較大，且兩包圍盒重合度較高時，該方法可能會產(chǎn)生誤判，沒有修正反向問題；當閾值k較高時，又可能出現(xiàn)無法達到閾值的情況，導致配準結(jié)果不確定。使用本文方法對兩模型進行配準，結(jié)果修正了反向問題，如圖11所示。表1所示為本文方法與文獻［14］方法在配準耗時和方向正確性兩方面的對比，本文方法一次配準0．266s得到正確結(jié)果，相對文獻［14］的方法有較少的耗時和較高的可靠性。

表1 本文方法與文獻［14］方法粗配準實驗對比

4 結(jié)束語

本文根據(jù)工業(yè)CT掃描對象的特點，通過計算點云模型形心與質(zhì)心的位置關(guān)系，判斷并修正PCA方法配準的反向模型。通過對照實驗，驗證了本文方法能夠解決PCA方法配準時的反向問題，在完成模型方位正確配準的前提下，能夠?qū)⑴錅收`差保持在較小范圍，為后續(xù)精配準提供較好的初始位置。但本文只研究了不規(guī)則工件模型的粗配準問題，規(guī)則工件模型的粗配準問題還有待于進一步研究。

［1］DUAN Liming，LIU Yuanbao，WU Zhifang，et al．Method of reconstructing 3－D CAD model based on industrial computed tomography［J］．Computer Integrated Manufacturing Systems，2009，15（3）：479－486（in Chinese）．［段黎明，劉元寶，吳志芳，等．基于工業(yè)計算機斷層成像技術(shù)的三維CAD模型重構(gòu)方法［J］．計算機集成制造系統(tǒng)，2009，15（3）：479－486．］

［2］RUSINKIEWIEZ S，LEVOY M．Efficient variants of the ICP algorithm［C］／／Proceedings of the 3rd International Conference on 3DDigital Imaging and Modeling．Washington，D．C．，USA：IEEE Computer Society，2001：145－152．

［3］EZRA E，SHARIR M，EFRAT A．On the performance of the ICP algorithm［J］．Computational Geometry：Theory and Applications，2008，41（1／2）：77－93．

［4］POTTMANN H，HUANG Qixing，YANG Yongliang，et al．Geometry and convergence analysis of algorithms for registration of 3Dshapes［J］．International Journal of Computer Vision，2006，67（3）：277－296．

［5］LUCIANO S，OLGA B，KIM B．Robust range image registration using genetic algorithms［M］．Singapore：World Scientific Publishing Co．，2005．

［6］JING Shikai，CHENG Yunyong，ZHANG Dinghua，et al．Tolerance zone constrained alignment method for turbine blade model［J］．Computer Integrated Manufacturing Systems，2010，16（4）：883－886（in Chinese）．［敬石開，程云勇，張定華，等．一種區(qū)域公差約束的葉片模型配準方法［J］．計算機集成制造系統(tǒng)，2010，16（4）：883－886．］

［7］YAN Sijie，ZHOU Yunfei，PENG Fangyu，et al．Optimization of allowance distribution on the workpieces with large sculptured surfaces in NC machining［J］．Journal of Huazhong University of Science and Technology：Nature Science，2002，30（10）：35－37（in Chinese）．［嚴思杰，周云飛，彭芳瑜，等．大型復雜曲面零件加工余量均布優(yōu)化問題研究［J］．華中科技大學學報：自然科學版，2002，30（10）：35－37．］

［8］LIU Bin，HUA Shungang，OU Zongying，et al．Plate prebending for complete fracture based on automatic registration of fractured bone model［J］．Journal of Optoelectronics Laser，2009，20（7）：977－982（in Chinese）．［劉 斌，華順剛，歐宗瑛，等．基于斷骨模型自動配準的完全性骨折鋼板預彎［J］．光電子·激光，2009，20（7）：977－982．］

［9］YAN Qingguang，LI Mingzhe，LI Dongcheng．Research on registration of 3Ddata in inspection of multi－point forming part［J］．China Mechanical Engineering，2003，14（19）：1648－1651（in Chinese）．［嚴慶光，李明哲，李東成．多點成形件檢測中三維數(shù)據(jù)配準方法的研究［J］．中國機械工程，2003，14（19）：1648－1651．］

［10］XU Faqiang，WAN Xiaoxia，ZHU Yuanhong．Color component prediction based on rotated principal component analysis［J］．Optics and Precision Engineering，2008，16（3）：518－523（in Chinese）．［許法強，萬曉霞，朱元泓．基于旋轉(zhuǎn)主成分分析的顏色組分預測［J］．光學精密工程，2008，16（3）：518－523．］

［11］JIANG Ming，ZHANG Guilin，HU Ruolan，et al．Research of an image matching method based on principal component a－nalysis［J］．Infrared and Laser Engineering，2000，29（4）：17－21（in Chinese）．［蔣 明，張桂林，胡若瀾，等．基于主成分分析的圖像匹配方法研究［J］．紅外與激光工程，2000，29（4）：17－21．］

［12］WU Feng，QIAN Zongcai，HANG Qiashi，et al．Principal axes algorithm based on image contour applied to medical image registration［J］．Journal of the Fourth Military Medical University，2001，22（6）：567－569（in Chinese）．［吳 鋒，錢宗才，杭洽時，等．基于輪廓的力矩主軸法在醫(yī)學圖像配準中的應用［J］．第四軍醫(yī)大學學報，2001，22（6）：567－569．］

［13］SHUI Wuyang，ZHOU Mingquan．An approach for model reconstruction based on multi－view scans registration［C］／／Proceedings of IEEE International Conference on Audio，Language and Image Processing．Washington，D．C．，USA：IEEE，2010：601－606．

［14］DAI Jinglan，CHEN Zhiyang，YE Xiuzi，et al．The application of ICP algorithm in point cloud alignment［J］．Journal of Image and Graphics，2007，12（3）：517－521（in Chinese）．［戴靜蘭，陳志楊，葉修梓，等．ICP算法在點云配準中的應用［J］．中國圖象圖形學報，2007，12（3）：517－521．］

［15］MODERSITZKI J．Numerical methods for image registration［M］．Oxford，UK：Oxford University Press，2003．

［16］SUN Jiaguang．Computer graphics［M］．2nd ed．Beijing：Tsinghua University Press，1995：344－348（in Chinese）．［孫家廣．計算機圖形學［M］．2版．北京：清華大學出版社，1995：344－348．］

［17］CIGNONI P，ROCCHINI C，SCOPIGNO R．Metro：measuring error on simplified surfaces［J］．Computer Graphics Forum，1998，17（2）：167－174．