姜曉通 郭保蘇 彭慶金 劉德利

1.常熟理工學院機械工程學院，常熟，2155002.燕山大學機械工程學院，秦皇島，066044 3.加拿大曼尼托巴大學機械工程學院，加拿大曼尼托巴，R3T 5V6

0 引言

在大部分三維打印中，模型的外部支撐結(jié)構(gòu)是不可避免的，是三維打印建模的重要研究內(nèi)容[1]。外部支撐結(jié)構(gòu)在打印過程中起到支撐模型懸空部分的作用，完成打印后需去除。外部支撐結(jié)構(gòu)的增加會增加打印模型的成本及時間，同時在去除支撐時，會影響支撐部分模型表面的質(zhì)量，有時甚至會破壞模型表面，因此如何在打印過程中減小支撐結(jié)構(gòu)的體積是三維打印建模的重要研究內(nèi)容之一。通常情況下，模型的支撐結(jié)構(gòu)體積與支撐本身的結(jié)構(gòu)和模型的打印方向兩個因素直接相關(guān)，因此，在支撐結(jié)構(gòu)相同的情況下，支撐的使用量主要取決于模型的打印方向。打印方向會在許多方面影響打印模型的質(zhì)量，如模型表面光滑度、支撐結(jié)構(gòu)材料的體積、模型強度及模型打印時間等。根據(jù)考慮影響打印模型質(zhì)量的因素數(shù)量，打印方向優(yōu)化可以分為單一目標屬性優(yōu)化和多目標屬性優(yōu)化。

單一目標屬性優(yōu)化是在確定模型的打印方向時只考慮影響模型打印質(zhì)量的單一屬性。KATTETHOTA等[2]在確定打印方向時考慮了模型表面的光滑度。GUPTA等[3]優(yōu)化了模型的打印時間，將單位球面分解為多個曲面多邊形，在每個曲面多邊形所定義的區(qū)間內(nèi)對打印時間進行優(yōu)化，計算最短打印時間對應的打印方向。大多數(shù)文獻計算支撐體積這一單一目標屬性來計算最優(yōu)方向。KHARDEKAR等[4]在計算最優(yōu)方向時采用GPU進行加速來提高算法的效率，相比CPU下的速度提高90%，且具有較好的精度。EZAIR等[5]采用與文獻[3]類似的方法，在GPU下求解最小打印體積所對應的方向。HU等[6]將模型分解為多個小的組件，計算每個組件不需要或只需要少量支撐的打印方向并分別進行打印，最后將小組件組裝成完整的模型。

多目標屬性優(yōu)化則是在確定打印方向時，同時考慮多個目標屬性。SANATI等[7]在確定最優(yōu)打印方向時，采用多目標遺傳算法來評估模型表面粗糙度、支撐體積及打印時間等屬性。BYUN等[8]同樣采用遺傳算法進行最優(yōu)打印方向的求解，在計算最優(yōu)打印方向時考慮打印時間及模型質(zhì)量兩個屬性。與單一目標屬性優(yōu)化相比，多目標屬性優(yōu)化更加耗時[9]。

現(xiàn)有的打印方向計算方法并沒有從理論上闡述最小支撐體積所對應方向的計算。本文將支撐體積目標屬性作為優(yōu)化對象，計算三維模型最小理論支撐體積所對應的打印方向。首先闡述了模型理論支撐體積的理論背景，針對離散網(wǎng)格模型的特點，提出了其理論支撐體積的計算方法，在計算過程中，利用光線追蹤算法處理凹模型。針對理論支撐體積函數(shù)的不可導性，提出了一種基于球面坐標采樣、無導數(shù)優(yōu)化算法進行最優(yōu)方向的計算方法。最后針對本文算法結(jié)果，從計算精度、時間及步長對結(jié)果的影響等方面對算法進行驗證分析。

1 理論支撐體積計算

設三維空間中的二維流形三角網(wǎng)格模型M=(C,E,V)∈R3，其中，C代表網(wǎng)格模型中點的拓撲關(guān)系，E代表網(wǎng)格模型中所有邊的集合，V={v0,v1,…,vn-1}(vi∈R3，i=0,1，…,n-1)是網(wǎng)格模型點的集合。定義S為模型的支撐結(jié)構(gòu)，V(S)為支撐結(jié)構(gòu)的體積。設+z為模型的打印方向，z方向上模型最低點的值為Zmin，即打印模型時的初始水平高度(打印平臺在z方向上的最小值)。在三維打印過程中，模型是“自上而下”逐層成形的過程，理論上模型懸空部分的所有網(wǎng)格面都需要添加支撐，該支撐為網(wǎng)格面到打印平臺的空間。然而，由于打印材料在成形過程中具有一定的延展性，模型的部分懸空網(wǎng)格面并不需要支撐，因此在實際打印過程中，支撐體積小于理論值。在計算模型最小支撐體積對應的打印方向時，在不確定支撐結(jié)構(gòu)的情況下不能以實際支撐體積作為計算依據(jù)。由此，在實際計算過程中，通常以理論支撐體積作為計算依據(jù)來計算實際支撐體積最小時對應的打印方向。基于此，本文在確定最小支撐體積對應的打印方向時同樣采用理論最小支撐體積。對于模型支撐S，其定義為

S:={(x,y,z)|z≥Zmin,?z′>z,
(x,y,z′)∈M且(x,y,z)?M}

(1)

式中，z′為投影點的z向值。

由式(1)可得，支撐S與模型M彼此互斥，即任何屬于模型M的部分都不屬于支撐S。

在計算模型給定方向的理論支撐體積時，由于網(wǎng)格模型并不是連續(xù)的曲面模型，其網(wǎng)格表面由許多離散的三角面片組成，故其求解并不是一個連續(xù)積分問題，而是一個離散積分問題。該離散積分問題求解的實質(zhì)是求和所有“懸空”三角面片對應的棱柱體的體積。根據(jù)凹凸性不通，模型可分為凹面模型和凸面模型，如圖1所示。在計算理論支撐體積時，它們的計算過程也不相同。圖1a為凸面模型的理論支撐體積，其值為所有“懸空”三角面片對應的三棱柱體積之和。圖1b為凹面模型的理論支撐體積，由于部分“懸空”三角面片對應的三棱柱并不在打印平臺上，而是位于模型上，因此在計算時需要計算三棱柱與模型交點的位置。

(a)凸面模型 (b)凹面模型圖1 模型理論支撐體積Fig.1 Theory support volume of models

(2)

若i=k-1，則i+1=0

圖2 多邊形面積計算Fig.2 Area calculation of polygon

與多邊形面積的求解相似，在計算模型的理論支撐體積時同樣采用離散積分的方法。假設T為模型的一個三角面片，3個頂點分別為vi=(xi,yi,zi)，i=0,1,2。三角面片T與打印平臺(打印平臺的z方向左邊為0)所組成三棱柱的體積計算公式為

V(T)=nz(z0+z1+z2)/6

(3)

n(T)=(nx,ny,nz),為三角面片T的非單位法矢，計算公式為n(T)=(v1-v0)×(v2-v0)。

對于凸面模型，其所有“懸空”三角面片對應的三棱柱的體積之和即為該模型的理論支撐體積，即

(4)

式中，n為模型中三角面片的個數(shù)。

對于凹面模型，在給定的打印方向，部分“懸空”三角面片對應的三棱柱的底部是落在網(wǎng)格模型上，而不是打印平臺,因此在計算凹面模型的理論支撐體積時，式(3)并不完全適用。此時需要計算哪些三角面片對應的三棱柱的底部落在打印平臺上，哪些不落在平臺上。本文利用光線追蹤算法判斷“懸空”三角面片這一屬性。

在計算機圖形學中，光線追蹤算法是被用來生成逼真三維場景的主要算法之一。近年來，隨著可視化、虛擬現(xiàn)實的發(fā)展，該算法得到了越來越廣泛的應用。在計算幾何中，光線追蹤算法被廣泛應用于幾何求交,但該算法在求交過程中需要較大的運算開銷，運行效率較低。為了提高運算效率，在進行求交運算前，一般需要對模型進行空間劃分，構(gòu)建一定的空間層次結(jié)構(gòu)，如均勻網(wǎng)格劃分[10]、八叉樹[11]及KD-Tree[12](一種特殊的二叉空間分割等空間劃分方法)。這些空間層次結(jié)構(gòu)都能夠加快光線追蹤的速度，本文采用最常用的KD-tree來提高光線追蹤的求交速度。根據(jù)生成支撐求交的特殊性，在求交過程中只計算“懸空”三角面片在光線方向上距離模型最近的交點。設算法1為基于KD-Tree的光線追蹤求交算法(圖3)。根據(jù)算法1，當“懸空”三角面片對應的三棱柱的底面落在模型上時，式(3)變?yōu)?/p>

(5)

2 最小理論支撐體積所對應的打印方向計算

由于理論支撐體積函數(shù)并不光滑，故在求解最小理論支撐體積對應的方向時，并不能利用函數(shù)求導數(shù)的方法進行求解，通常需要利用無導數(shù)優(yōu)化算法[13-15]來求解最優(yōu)方向。此時，對于多極值函數(shù)，任意給定一個初值進行最優(yōu)值求解時，其求解結(jié)果往往只是一個局部最優(yōu)解，并不是全局最優(yōu)解。為得到函數(shù)的全局最優(yōu)解，在利用無導數(shù)優(yōu)化算法進行求解前，需要確定部分可信賴區(qū)間，然后計算各個可信賴區(qū)間的局部最優(yōu)值，進而得到全局最優(yōu)解。

本文利用球面坐標對打印方向進行采樣，以計算得到的理論支撐體積最小的前n個采樣方向所在的區(qū)間作為可信賴區(qū)間，在可信賴區(qū)間內(nèi)求得局部極值，再通過比較得到全局最優(yōu)解。在計算最小理論支撐體積對應的打印方向前，首先討論模型理論支撐體積函數(shù)V(S)的可導性。

2.1 理論支撐體積函數(shù)的可導性

為了簡化計算，本文以凸面模型為例討論其理論支撐體積函數(shù)的可導性。由于網(wǎng)格模型的理論支撐體積為模型“懸空”三角面片對應的棱柱體體積之和，因此這里只討論其中一個“懸空”三角面片T對應的三棱柱體積函數(shù)V(T)的可導性。當T繞著x軸旋轉(zhuǎn)角度θ時，式(3)中的zi變?yōu)閦icosθ+yisinθ(i=0,1,2)，nz變?yōu)閚zcosθ+nysinθ，則式(3)變?yōu)轶w積關(guān)于轉(zhuǎn)角θ的函數(shù)：

算法 1:光線追蹤求交輸入: (1)pNode:與網(wǎng)格模型M相關(guān)聯(lián)的Kd-Tree(2)vOrigin:光線追蹤源點(3)vDirection:光線方向輸出:(1)fNearest:浮點數(shù),初始值為FLT_MAX(2)vHit:光線與模型的交點(3)nHitTri:交點所在的三角面片 執(zhí)行:bool rayIntersection(pNode, vOrigin, vDirection, fNearest,vHit, nHitTri){ Point pHit; bool bHit = isIntersection (pNode, vOrigin, bHit); //判斷光線與pNode的包圍盒是否相交 If (isLeaf(pNode))If(bHit) Triangles=GetTriangleID (pNode); //得到節(jié)點pNode內(nèi)的三角面片 bool bHitTri = IsInsectionTri (Triangle,vOrigin, vDirection, pHit); //判斷光線與三角面片是否相交,若相交則得到交點pHit If(bHitTri&&dis(vOrigin, pHit)pLeft, vOrigin, vDirec-tion, fNearest,vHit, nHitTri); bool bRight =rayIntersection (pNode->pRight, vOrigin, vDi-rection, fNearest,vHit, nHitTri); If(bLeft | |bRight)return true; return false;}

圖3光線追蹤求交算法

Fig.3Raytracingintersectionalgorithm

(6)

式(6)變形后為

(7)

為了進一步簡化計算，這里假設T的初始位置與水平面平行，則式(7)可以簡化為

(8)

(a)立方體模型 (b)RockArm模型圖4 測試模型Fig.4 Test models

2.2 基于無導數(shù)優(yōu)化算法的最優(yōu)方向計算

2.2.1球面坐標采樣

在利用無導數(shù)優(yōu)化算法求解最優(yōu)打印方向前，首先需要確定求解的部分可信賴區(qū)間。由于本文算法的求解目標為打印方向，因此采用球面坐標的方法對空間方向進行采樣，計算每個采樣點對應打印方向的理論支撐體積，選擇體積最小的前n個方向?qū)膮^(qū)間作為無導數(shù)優(yōu)化算法求解的可信賴區(qū)間。

(a)立方體模型

(b)RockArm模型圖5 模型理論支撐體積曲線圖Fig.5 Curve graph of models’ theory support volume

圖6 球面坐標示意圖Fig.6 Schematic diagram of spheroidal coordinates

圖7 球面坐標采樣結(jié)果Fig.7 Results of spheroidal coordinate sampling

球面坐標以坐標原點為參考點，由方位角、仰角及距離構(gòu)成。如圖6所示，假設p(x,y,z)為三維空間內(nèi)一點，則點p的空間位置可用3個有次序的數(shù)(r,θ,φ)來確定，其中,r為原點o與點p間的距離，θ為從正z軸來看自x軸按逆時針方向轉(zhuǎn)到op′所轉(zhuǎn)過的角，φ為有向線段op與z軸正向的夾角，這里p′為點p在xoy面上的投影。r、θ、φ的取值范圍分別為：r∈[0,+∞)，φ∈[0,2π]，θ∈[0,2π]。球面坐標采樣的實質(zhì)為對θ及φ方向進行球面坐標均勻采樣，并利用一定的轉(zhuǎn)化關(guān)系將其轉(zhuǎn)化成直角坐標系內(nèi)的空間點。以圖6中的點p為例，該轉(zhuǎn)化關(guān)系為：xp=rsinθcosφ，yp=rsinθsinφ，zp=rcosθ。這樣，在三維空間內(nèi)生成一系列的采樣點pi(xi,yi,zi)，則其采樣方向opi=(xi,yi,zi)。圖7為利用球面坐標生成的采樣點，在生成采樣點過程中，在θ及φ方向上進行步長為6°的采樣，共得到1 800個采樣方向，然后計算每個采樣方向?qū)睦碚撝误w積。

2.2.2基于無導數(shù)優(yōu)化算法的方向優(yōu)化

在對打印方向進行采樣并計算對應的理論支撐體積后，將最小的幾個理論支撐體積對應的采樣打印方向作為初始值，利用無導數(shù)優(yōu)化算法在選定的采樣方向?qū)目尚刨噮^(qū)間內(nèi)計算局部最優(yōu)方向，則計算得到的所有局部最優(yōu)方向的最小值為全局最優(yōu)方向。在利用無導數(shù)優(yōu)化算法計算局部最優(yōu)解時，求解的可信賴區(qū)間是保證算法能夠收斂的基礎。當區(qū)間過大時，區(qū)間內(nèi)部可能會有兩個甚至更多的極值，利用無導數(shù)優(yōu)化算法得到的局部最優(yōu)解不可信。該區(qū)間的大小取決于采樣步長，為了保證計算區(qū)間可靠，計算收斂，采樣時步長不宜過大。

無導數(shù)優(yōu)化算法可以分為有限差分方法[13]、基于模型的方法[14]、模式搜索方法[15]等3種方法。模式搜索方法是一種直接的局部搜索方法，在搜索過程中直接使用函數(shù)值，并不需要評估函數(shù)的導數(shù)。該算法的基本思想從幾何意義上來說，是尋找具有較小函數(shù)值的“山谷”，通過迭代的方法使搜索沿著“山谷”的走向來逼近極小值。該算法從給定的起始點開始，通過局部移動尋找下降的方向來得到局部的最優(yōu)值，其最經(jīng)典的算法為坐標搜索。坐標搜索是一種無約束的優(yōu)化求解算法，求解目標函數(shù)f(x)的最小值。在計算時首先選定一個初始值及步長，然后通過迭代的方法得到一個局部最小值。本文在計算最優(yōu)打印方向時，同樣是找到最優(yōu)打印方向所對應的球面坐標，初始化值為上面所計算的最小的幾個理論支撐體積所對應的方向，步長閾值為1°。在每個作為初始值的方向上，利用坐標搜索算法得到局部的最小理論支撐體積。這些局部最小理論支撐體積的最小值所對應的打印方向即為最優(yōu)的打印方向。

3 實驗結(jié)果分析

本文算法開發(fā)環(huán)境為VS2008，編程語言為C++，基于OpenGL圖形庫進行模型顯示，測試用電腦的處理器為Inter Core 10 ms i5-2500，主頻為3.30GHz，內(nèi)存為4GB。本節(jié)主要從算法的精度及計算時間兩個方面來驗證本文算法的有效性與實用性，同時分析采樣步長對計算結(jié)果的影響。在計算過程中，為了加快運算速度，采用OpenMP技術(shù)進行并行運算，以提高運算效率。

3.1 理論支撐體積的計算精度與計算時間

圖8 理論支撐體積的計算誤差Fig.8 Calculation error of theory support volume

圖9 理論支撐體積的計算時間Fig.9 Calculation time of theory support volume

圖9所示為理論支撐體積的計算時間與模型三角面片個數(shù)之間的關(guān)系。測試所用的模型為圖4b的RockArm模型，且在測試過程中模型的網(wǎng)格數(shù)不斷增加。對于每個測試模型，在測試過程中的采樣步長為10°，共648個采樣方向。模型在每個采樣方向上都進行理論支撐體積的計算，然后得到單次計算的平均時間。由圖9中的曲線可得到，單次理論支撐體積的計算時間與模型三角面片個數(shù)成線性關(guān)系，隨三角面片個數(shù)的增加而增長。當三角面片個數(shù)為70 000時，計算時間為10 ms左右，完全能夠滿足時間要求。

3.2 采樣步長對計算結(jié)果的影響

在計算過程中，采樣步長會對計算時間及最終的優(yōu)化結(jié)果產(chǎn)生較大影響。這里利用不同的采樣步長來測試其對計算結(jié)果的影響。

圖10為采樣步長對采樣時間及優(yōu)化時間的影響曲線，測試模型為圖4b的RockArm模型，測試所用的采樣步長分別為2°、3°、4°、5°、6°、9°、10°、12°、15°、18°、20°、22°、25°、30°、36°、45°、55°、60°、80°、90°、180°。在利用無導數(shù)優(yōu)化算法進行優(yōu)化計算時，計算理論支撐體積最小的前100個采樣方向，在這100個采樣方向確定的區(qū)間內(nèi)進行優(yōu)化計算。由圖10可以得到，對于采樣時間，當采樣步長增大時，采樣方向減小，在單次計算時間一定的情況下總的采樣時間縮短。對于優(yōu)化時間，當采樣步長較小時，在優(yōu)化方向數(shù)量確定的情況下，由于優(yōu)化區(qū)間較小，單次優(yōu)化時間也會縮短，因此總的優(yōu)化時間也較短。隨著采樣步長的增大，優(yōu)化區(qū)間增大，優(yōu)化時間也會相應增長。當采樣步長為25°時，采樣方向的數(shù)量剛好為100，此時理論上的優(yōu)化時間應該最長。在實際采樣過程中，由于模型的初始位置不同，用于局部優(yōu)化的各初始方向也不同，單次優(yōu)化時間也會所有變化，因此最長優(yōu)化時間對應的采樣步長約為25°，但會有所波動。當采樣步長大于25°時，采樣方向的個數(shù)小于100，隨著采樣步長增大，采樣方向的個數(shù)也隨之減少，因而優(yōu)化時間也隨著減短。

圖10 不同采樣步長的測試時間Fig.10 Test time of different sampling step

圖11所示為在不同采樣步長下，采樣得到的最小理論支撐體積與優(yōu)化后得到的最小理論支撐體積。從圖11中可以得到，當測試過程中的采樣步長足夠小時(采樣步長的數(shù)值為2°、3°、6°時)，經(jīng)采樣得到的最小理論支撐體積即為優(yōu)化后的最小理論支撐體積，不需要后續(xù)的優(yōu)化即可得到最優(yōu)打印方向；當采樣步長較大時需要利用優(yōu)化算法才可得到最優(yōu)的打印方向；當采樣步長超過一定的值時，即使利用優(yōu)化算法也得不到最優(yōu)打印方向。

圖11 不同采樣步長的計算結(jié)果Fig.11 Calculation results of different sampling step

表1所示為不同網(wǎng)格模型的測試對比結(jié)果。在測試過程中，采樣步長為5°，共計算36×72=2 592個方向的理論支撐體積。在優(yōu)化過程中，選取理論支撐體積最小的前100個方向進行優(yōu)化計算，分別計算模型三角面片個數(shù)、網(wǎng)格模型體積V(M)、最優(yōu)采樣方向理論支撐體積、最優(yōu)方向理論支撐體積、采樣時間及優(yōu)化時間。從表1中的計算結(jié)果可以得到，計算時間隨著模型三角面片個數(shù)的增加而增長。當模型網(wǎng)格數(shù)量較小時，本文算法效率非常高，即使模型三角面片個數(shù)為200 000，計算時間也短于180 s。當模型三角面片個數(shù)太多時，可先對模型進行基于特征的簡化，減少模型三角面片個數(shù)，然后再進行打印方向的優(yōu)化。為了進一步驗證本文算法有效性，將本文算法與Autodesk公司的Meshmixer軟件進行對比測試。在進行對比測試時，除了一個模型外，其余所有模型的最終打印方向與本文算法的結(jié)果基本一致，如圖12所示。在得到各自的最優(yōu)打印方向后，利用Meshmixer軟件提供的算法生成樹狀支撐結(jié)構(gòu)，分別計算兩個不同最優(yōu)打印方向下的實際支撐體積。從結(jié)果可以看出，對于該測試模型，本文算法結(jié)果更優(yōu)。

表1 不同模型測試結(jié)果

(b)本文算法計算結(jié)果圖12 結(jié)果比較Fig.12 Comparison of results

4 結(jié)論

(1)本文提出了一種基于最小理論支撐體積的模型三維打印方向的計算方法，該方法能夠快速準確地計算最小理論支撐體積所對應的模型打印方向，為模型三維打印方向提供了輔助的優(yōu)化方案。

(2)闡述了模型理論支撐體積的理論背景，給出了離散網(wǎng)格模型理論支撐體積的計算方法，同時針對凹模型的理論支撐體積的特殊性，利用光線追蹤的方法來計算凹模型的理論支撐體積。

(3)針對理論支撐體積函數(shù)不可導的特點，利用無導數(shù)優(yōu)化算法來計算最小理論支撐體積所對應的方向，該方法可采用并行運算以加快計算速度。

(4)實驗表明本文算法在計算精度和計算效率上都優(yōu)于已有的方法。

本文算法在計算方向時只考慮支撐體積這一單一因素，下一步研究將進一步考慮三維打印的表面粗糙度、打印時間等多方面因素，增強算法的適用性。