三角函數(shù)逼近算法及其在光學(xué)條紋圖像分析中的應(yīng)用＊

周 珺，郭紅衛(wèi)

（上海大學(xué) 機電工程與自動化學(xué)院，上海 200072）

引 言

三角函數(shù)與反三角函數(shù)作為基本初等函數(shù)，在工程技術(shù)與科學(xué)研究中有廣泛的應(yīng)用，但在某些特定情況下則需通過逼近函數(shù)來計算其近似值。例如，在工業(yè)中常采用單片機、硬件計算設(shè)備或各種小型系統(tǒng)來實現(xiàn)過程的控制［1］，然而這些系統(tǒng)可能并不支持三角函數(shù)的庫函數(shù)計算，那么就需要尋找替代的逼近方法來計算其數(shù)值。另外，相對于一般算術(shù)運算，三角函數(shù)與反三角函數(shù)計算復(fù)雜度相對較高，計算耗時，難以滿足各種實時性要求。例如，在機器人動力學(xué)［2］方程的快速計算中，三角函數(shù)的計算也占很大的比例［3］。另外，在計算機三維建模與三維游戲中，三維模型的空間坐標變換也涉及到大量連續(xù)的三角函數(shù)運算［4］。GPU［5］作為一種可編程的圖形處理器，利用Vertex Shader計算反三角函數(shù)也很浪費時間。在這些情況下，采用三角函數(shù)的快速近似計算可以大大提高計算效率［6－7］，滿足各類實時性要求。

同樣，在光學(xué)檢測技術(shù)中，也涉及到大量的三角函數(shù)運算。例如，傅里葉變換位相分析方法［8－9］中，大量的復(fù)指數(shù)運算需通過正余弦函數(shù)的計算來實現(xiàn)；相移步距未知情況下的相移算法中，需使用反余弦函數(shù)來求解相對相移量［10］；各種位相分析方法一般需通過四象限反正切函數(shù)逐像素計算條紋圖像上的位相分布［11］；而在多視角（口徑）測量［12－13］中，則需利用坐標變換來實現(xiàn)各視角測量結(jié)果的配準與拼接，也要用到大量的三角函數(shù)計算。這些計算消耗大量時間，已成為快速測量或?qū)崟r測量的一個瓶頸問題，同時也使得測量數(shù)據(jù)的硬件或固件計算處理難以實現(xiàn)。

三角函數(shù)與反三角函數(shù)的快速近似計算可采用多種方法實現(xiàn)。其中，查表法首先對函數(shù)特定區(qū)間按特定分辨率進行采樣并計算其數(shù)值，并建表存儲在內(nèi)存中。使用時，直接訪問相應(yīng)地址即可獲得需要的三角函數(shù)值。查表法簡便易行，效率很高，但需要占用一定的內(nèi)存空間。特別是當精度要求提高時，所需內(nèi)存的大小也需相應(yīng)地增大。泰勒級數(shù)也常被用來計算函數(shù)的近似值。根據(jù)泰勒定理，一個無限可微函數(shù)可由其泰勒展開式無限逼近。但泰勒級數(shù)收斂較為緩慢，無法實現(xiàn)低階高精度的逼近。例如在區(qū)間［－1，1］上用泰勒級數(shù)計算反正切函數(shù)值時，若要精度達到10－3數(shù)量級，需將其展開至49階；若要精度達到10－4數(shù)量級，則需將其展開至499階。這一現(xiàn)象說明，當精度要求較高時，泰勒級數(shù)并不是一種實用有效的函數(shù)逼近方法。另一種方法是采用最小二乘多項式逼近三角函數(shù)，該多項式與目標函數(shù)之間的平方距離可達到極小值。由于最小二乘多項式系數(shù)的計算十分簡單，該方法較易于實現(xiàn)。采用這種方法，僅需5階和7階多項式即可分別使上述區(qū)間內(nèi)的反正切函數(shù)的求解精度達到10－3和10－4數(shù)量級，其逼近效率遠高于泰勒級數(shù)法。但是，這種最小二乘多項式并非是最優(yōu)的。換言之，存在階數(shù)更低的多項式可達到相同的精度。不同于上述方法，本文討論三角函數(shù)及反三角函數(shù)的最佳逼近方法，推導(dǎo)了基于無窮范數(shù)的最佳逼近多項式。

1 三角函數(shù)逼近算法

1.1 原理

假設(shè)y＝f（x）是區(qū)間［x1，x2］內(nèi)的連續(xù)函數(shù)，該區(qū)間其N 階逼近多項式可表示為：

pn（n＝0，1，…，N）為其系數(shù)。該方程涉及N（1＋N）／2次乘法和N 次加法運算。為簡化，式（1）可改寫為：

對于具體的三角函數(shù)，首先針對其某一特定區(qū)間［x1，x2］求解逼近多項式，然后再將其結(jié)果推廣到整個2π周期區(qū)間。根據(jù)數(shù)值分析理論，f＾（x）－f（x）在閉區(qū)間［x1，x2］內(nèi)至少有N＋1個不同的根（根的集合可能包含x1和x2），在開區(qū)間（x1，x2）內(nèi)至少有N 個交錯極值點，其極值的絕對值相等正負號交替。當將區(qū)間［x1，x2］上的逼近方程拓展到區(qū)間外時，不僅要保證分段逼近函數(shù)在整個函數(shù)定義區(qū)間內(nèi)盡可能是最優(yōu)的，而且要保證其在整個定義區(qū)間上的連續(xù)性。為此要求，在選定區(qū)間［x1，x2］的邊緣處，多項式與函數(shù)f（x）的關(guān)系滿足一定條件。假設(shè)x＝xi（i＝1或2）是一個邊緣點。如果函數(shù)f（x）關(guān)于x＝xi左右對稱則在x＝xi處有一個附加的交錯極值點，即：

其中，η為極值。關(guān)于極值的正負號，當i＝1（左邊緣點）時，m＝0，當i＝2（右邊緣點）時，m＝N＋1。如果函數(shù)f（x）在x＝xi處非左右對稱，則強制要求與f（x）精確相等，即：

第1步：任選N 個點ξk（k＝1，2，…，N），滿足x1＜ξ1＜ξ2＜…＜ξN ＜x2。

第2步：假設(shè)差值f＾（ξk）－f（ξk）（k＝1，2，…，N）具有相等的絕對值和交錯的正負號，即：

將式（6）與式（4）或式（5）聯(lián)立形成一個包含N＋2個方程式的線性方程組，該方程組包含pn（n＝0，1，…，N）和η共N＋2個未知數(shù)。

第3步：利用第二步計算所得系數(shù)pn構(gòu)建多項式。解方程：

這種方法可以構(gòu)建出一個函數(shù)在區(qū)間［x1，x2］內(nèi)的逼近多項式。利用該方法可以推導(dǎo)出基本的三角函數(shù)（正弦，正切）及其反函數(shù)在特定區(qū)間內(nèi)的逼近多項式，然后再將其結(jié)果拓展至整個2π周期或整個函數(shù)定義區(qū)間。

1.2 結(jié)果和討論

1.2.1 余弦函數(shù)逼近

正弦和余弦函數(shù)是兩種基本的三角函數(shù)。因為sin（x）＝cos（π／2－x），正弦函數(shù)sin（x）可以由通過余弦函數(shù)cos（x）求得，所以只需針對余弦函數(shù)進行討論。余弦函數(shù)是一個定義在（－∞，＋∞）內(nèi)的周期函數(shù)，利用上述方法可以在一個基本周期［0，2π］內(nèi)求解其逼近多項式。但是這種分段方式無法獲得效率較高的逼近方程。例如，采用4次逼近多項式，最大誤差達到0.035；若要使逼近精度達到10－3數(shù)量級，多項式次數(shù)至少為6次。這意味著執(zhí)行較多的乘法運算，要消耗較多的計算時間。因此，函數(shù)逼近區(qū)間局限在［0，π／2］范圍內(nèi)，最后通過三角恒等式得到整個周期上的分段逼近多項式。采用第2.1節(jié)所述的方法，對［0，π／2］區(qū)間內(nèi)的余弦函數(shù)進行多項式逼近。由于余弦函數(shù)是偶函數(shù)，在區(qū)間左邊界x＝0處，采用式（4）作為邊界條件；而在右邊界x＝π／2處，以式（5）作為邊界條件。計算所得的3至5次多項式系數(shù)及最大誤差如表1所示，其中pn為式（1）所定義的多項式的n次項系數(shù)。

表1 區(qū)間［0，π／2］內(nèi)余弦函數(shù)逼近多項式系數(shù)和最大誤差Tab.1 Coefficients and maximum errors of approximation polynomials for the cosine function in［0，π／2］

整周期［0，2π］區(qū)間內(nèi)的余弦函數(shù)則可以利用上述逼近多項式，以分段逼近多項式計算得到，即

1.2.2 正切函數(shù)逼近

正切函數(shù)tan（x）是一個周期函數(shù)。在周期（－π／2，π／2）之內(nèi)，當x 趨近于±π／2 時，該函數(shù)取值趨近于±∞。在整個周期上就無法獲得正切函數(shù)單一的最佳逼近多項式。但是可以選取一個相對較小的區(qū)間計算其逼近多項式。例如，表2所示是采用1.1節(jié)方法求得的正切函數(shù)在［0，π／2］之間的逼近多項式的系數(shù)及最大誤差。為保證連續(xù)性，計算采用的邊界條件由式（5）確定。

圖1 余弦函數(shù)3次（短劃線）、4次（點線）和5次（實線）分段逼近多項式誤差Fig.1 The approximation errors of cos（x）with polynomial degree being 3（dashed），4（dot）and 5（solid lines）

表2 區(qū)間［0，π／4］內(nèi)正切函數(shù)逼近多項式系數(shù)與最大誤差Tab.2 Coefficients and maximum errors of approximation polynomials for the tangent function in［0，π／4］

利用上述［0，π／4］區(qū)間內(nèi)的逼近多項式，有多種方法可以實現(xiàn)整個周期區(qū)間內(nèi)的正切函數(shù)的近似計算。第一種方法就是利用三角恒等式，可得：

式（9）表示正切函數(shù)在一個整周期（－π／2，π／2）上的分段逼近函數(shù)。其誤差分布如圖2所示。從中可見，根據(jù)式（9）所得正切函數(shù)近似值在一個周期內(nèi)仍然是連續(xù)分布的。在［－π／4，π／4］區(qū)間上最大逼近誤差被最小化，其逼近函數(shù)是最佳逼近。但是，在［－π／4，π／4］區(qū)間之外，誤差分布是不均勻的，這種逼近不再是最佳逼近。特別是在靠近周期邊緣處，當接近π／2時，誤差增大較快。這是由于正切函數(shù)在此處趨近于∞的性質(zhì)決定的。

第二種獲取整周期函數(shù)的方法是利用前述正弦函數(shù)sin（x）和余弦函數(shù)cos（x）的逼近多項式，通過三角關(guān)系式tan（x）＝sin（x）／cos（x）進行計算。模擬實驗證明在接近±π／2處，這種方法精度略高于第一種方法，但是其誤差仍然不是均勻分布，這種方法也不是最優(yōu)的，并且需要執(zhí)行較多的乘法運算。實際上，正切函數(shù)在（－π／2，π／2）周期內(nèi)不存在一個最佳逼近多項式或分段的最佳逼近多項式。在具體應(yīng)用中，所要求解的正切函數(shù)角度范圍往往被限制在一個較小的已知區(qū)間內(nèi)，此時就可以利用1.1節(jié)所述方法直接求解該區(qū)間內(nèi)的最佳逼近多項式。

1.2.3 反余弦函數(shù)逼近

反正弦和反余弦函數(shù)是兩種基本的反三角函數(shù)。因為arcsin（x）＝π／2－arccos（x）成立，反正弦函數(shù)arcsin（x）的近似值可以通過反余弦函數(shù)arccos（x）的逼近函數(shù)求得。反余弦函數(shù)是一個定義在區(qū)間［－1，1］內(nèi)的函數(shù)。如果在整個區(qū)間內(nèi)求解單一的逼近多項式，也無法獲得效率較高的逼近方程。例如，采用3次逼近多項式，最大誤差達到0.14；若要使逼近精度達到10－4數(shù)量級，多項式次數(shù)至少為11次。這是因為反余弦函數(shù)在x趨近于±1時，導(dǎo)數(shù)趨近于±∞，這就要求采用非常高階的多項式才能逼近該曲線。采用高階多項式意味著需執(zhí)行更多的乘法運算，消耗更多的時間。因此，將函數(shù)逼近的區(qū)間局限在范圍內(nèi)，最后通過三角恒等式得到整個定義域上的分段逼近多項式。

采用1.1節(jié)所述的方法，以式（5）為邊界條件，分別用1、3、5次多項式對區(qū)間內(nèi)的余弦函數(shù)進行多項式逼近，所得的多項式系數(shù)及最大誤差如表3所示。

圖2 正切函數(shù)的3次（短劃線）和5次（實線）分段逼近多項式誤差Fig.2 The approximation errors of tan（x）with polynomial degrees being 3（dashed）and 5（solid lines）

表3 區(qū)間內(nèi)反余弦函數(shù)逼近多項式的系數(shù)和最大誤差Tab.3 Coefficients and maximum errors of approximation polynomials for the arccosine function in

表3 區(qū)間內(nèi)反余弦函數(shù)逼近多項式的系數(shù)和最大誤差Tab.3 Coefficients and maximum errors of approximation polynomials for the arccosine function in

多項式的 多項式系數(shù) 最大逼近次數(shù)p0 p1 p2 p3 p4 p5誤差／rad 1 π／2 －1.110 720 0.033 00 3 π／2 －1.019 450 0.126 854 －0.361 939 0.000 87 5 π／2 －1.002 146 0.036 749 －0.364 830 0.432 948 －0.420 863 0.000 038

可以看出式（10）所得到的逼近多項式在各分段邊界處是連續(xù)的，在各個分段區(qū)間內(nèi)的最大誤差是一致的。與前述整個定義區(qū)間上單一的逼近多項式相比，其多項式次數(shù)較少，計算量要小得多。

1.2.4 反正切函數(shù)逼近

反正切函數(shù)arctan（x）是定義在（－∞，＋∞）上的奇對稱函數(shù)，無法獲得整周期內(nèi)的單一逼近多項式。為此可以選取一個小區(qū)間計算其逼近多項式，再把其拓展到整個區(qū)間。表4 所示是采用1.1 節(jié)方法，以式（5）為邊界條件，求得的正切函數(shù)在［0，1］之間的逼近多項式的系數(shù)及最大誤差。從表中可見，由于反正切函數(shù)的奇對稱性質(zhì)，這些多項式的0次項（常數(shù)項）為0。

圖3 反余弦函數(shù)的3次（短劃線）與5次（實線）分段逼近多項式誤差Fig 3 The approximation errors of arccos（x）with polynomial degrees being 3（dashed）and 5（solid lines）

表4 區(qū)間［0，1］內(nèi)反正切函數(shù)逼近多項式的系數(shù)和最大誤差Tab.4 Coefficients and maximum errors of approximation polynomials for the arctangent function in［0，1］

得到區(qū)間［0，1］內(nèi)的逼近多項式后，整個定義域區(qū)間（－∞，＋∞）內(nèi)的反余弦函數(shù)可以利用分段逼近多項式計算得到，即：

圖4顯示在使用3次與5 次逼近時誤差的比較。可以看出式（11）所得到的逼近多項式在各分段邊界處是連續(xù)的，在各個分段區(qū)間內(nèi)的最大誤差是一致的。

2 應(yīng) 用

2.1 條紋分析算法

圖4 反正切函數(shù)的3次（短劃線）與5次（實線）分段逼近多項式誤差Fig 4 The approximation errors arctan（x）with polynomial degrees being 3（dashed）and 5（solid lines）

光學(xué)三維測量技術(shù)常以二維條紋圖像的形式提供其原始測量數(shù)據(jù)［14－15］。其中，第三維信息就包含于條紋圖像的位相之中。因此，分析條紋的位相信息是實現(xiàn)光學(xué)測量的關(guān)鍵步驟，其中涉及大量的三角函數(shù)計算。下面以相移步距未知情況下的相移算法為例來說明本文方法在條紋圖像分析中的應(yīng)用。當相移步距未知并沿圖像為不均勻分布時，第k 幅相移條紋圖可表示為

其中，（x，y）是圖像像素坐標；Ik、a和b 分別表示光強、背景光強與調(diào)制度；位相φ 與相對相移量α 均為未知量。條紋分析的第一步是針對每一像素估計其相對相移量：

再利用最小二乘相移算法求解位相。解方程組：

則各像素的位相為

上述過程對每一像素點來說總共涉及2 K 次的正弦或余弦計算，一次反余弦計算和一次反正切計算。這些運算都可以利用本文方法計算，以簡化計算。

2.2 實驗

上述相移步距未知情況下的相移算法，適用于匯聚光相移干涉測量［16］、陰影莫爾形貌測量［17］、鏡像莫爾形貌測量［18］等。在這些方法中，位相與相移量均依賴于被測表面的深度或梯度等形貌信息，因而未知且呈非均勻分布。圖5所示為一幅相移鏡像莫爾條紋圖。總相移步數(shù)為K＝64。圖6比較了條紋分析的結(jié)果。圖6（a）和6（b）為通過式（13），分別利用標準反余弦函數(shù)及其逼近算法計算所得相對相移量。其中，逼近多項式階數(shù)為3，系數(shù)見表3。圖6（c）顯示兩者的差值，說明逼近方法誤差小至10－4數(shù)量級。圖6（d）和6（e）為通過式（14）與式（15），分別利用標準三角函數(shù)及其逼近算法計算所得包裹位相。其中，逼近多項式階數(shù)為3，余弦函數(shù)逼近多項式系數(shù)見表1，反正切函數(shù)逼近多項式系數(shù)見表4。圖6（f）為兩者之差，從中可見，采用3階多項式逼近算法，導(dǎo)致的最大位相計算誤差僅為0.027rad。采用更高階多項式，可以進一步提高精度。對整幅圖像進行處理，采用標準函數(shù)算法與逼近算法所用時間分別為45.70s和39.15s。后者略快于前者。除計算復(fù)雜性外，計算速度還取決于計算機性能，程序代碼的優(yōu)化等等。在某些情況下，函數(shù)變量局限于一個較小的已知區(qū)間時，可以減少比較運算的次數(shù)，大大提高計算的效率。

圖5 莫爾條紋圖Fig 5 A Moiréfringe pattern

圖6 條紋圖像分析結(jié)果Fig 6 Fringe pattern analysis results

3 結(jié) 論

在三角函數(shù)及反三角函數(shù)的特定區(qū)間內(nèi)推導(dǎo)了基于無窮范數(shù)的最佳逼近多項式，得到了各個三角函數(shù)和反三角函數(shù)分段的多項式。利用三角函數(shù)的關(guān)系式將這些三角函數(shù)推廣到整周期區(qū)間，或利用反三角函數(shù)的關(guān)系式將這些反三角函數(shù)推廣到整個定義區(qū)間。文中也給出了一定精度內(nèi)的各函數(shù)逼近多項式的系數(shù)及最大誤差。將上述結(jié)果應(yīng)用于光學(xué)條紋圖像的分析，通過實驗證明了方法的有效性。所述三角函數(shù)逼近算法，由于只是涉及一些簡單的邏輯和算術(shù)運算，非常適合于光學(xué)測量信息處理中硬件與固件計算的實現(xiàn)。另外，在一些具有實時性要求的應(yīng)用中，可以替代庫函數(shù)運算，提高計算速度。此外，這些結(jié)果具有普遍性，也可用于其它科學(xué)研究及工程應(yīng)用目的。

［1］ 胡天亮，張洪波，劉日良.基于STEP－NC的數(shù)控車削系統(tǒng)的加工工藝參數(shù)優(yōu)化［J］.工藝與檢測，2011（8）：128－131.

［2］ 魏振忠，張 博，張廣軍.雙機器人系統(tǒng)的快速手眼標定方法［J］.光學(xué) 精密工程，2011，19（8）：1895－1902.

［3］ 于華東，方 濱，周東輝.一種機器人動力學(xué)方程快速計算的三角函數(shù)發(fā)生器［J］.機器人，2002，24（6）：7－9.

［4］ BELLIS S J，MARNANE W P.A CORDIC arctangent FPGA implementation for ahigh－speed 3D－camerasystems［J］.Lecture Notes in Computer Science，2000，1896：485－494.

［5］ 許 楨.關(guān)于CPU＋GPU 異構(gòu)計算的研究與分析［J］.科技信息，2010（17）：613.

［6］ 李加文，陳宗雨，李從心.基于函數(shù)逼近的三角函數(shù)加減速方法［J］.機床與液壓，2006（3）：66－67.

［7］ 郭新貴，李從心.一種新型柔性加減速算法［J］.上海交通大學(xué)學(xué)報，2003，37（2）：205－207.

［8］ TAKEDA M，MUTOH K.Fourier transform profilometry for the automatic measurement of 3－D object shapes［J］.Appl Opt，1983，22（24）：3977－3982.

［9］ TAKEDA M，INA H，KOBAYASHI S.Fourier transform method of fringe－pattern analysis for computer－based topography and interferometry［J］.J Opt Soc Am，1982，72（1）：156－160.

［10］ GUO H W，CHEN M Y.Least－squares algorithm for phase－stepping interferometry with an unknown relative step［J］.Appl Opt，2005，44（23）：4854－4859.

［11］ GUO H W，LIU G Q.Approximations for the arctangent function in efficient fringe pattern analysis［J］.Opt Express，2007，15（6）：3053－3066.

［12］ GUO H W，CHEN M Y.Multiview connection technique for 360－deg three－dimensional measurement［J］.Opt Eng，2003，42（4）：900－905.

［13］ HE H T，CHEN M Y，GUO H W，et al.Novel multiview connection method based on virtual cylinder for 3－D surface measurement［J］.Opt Eng，2005，44（8）：083605.

［14］ 陶 濤，郭紅衛(wèi)，何海濤.鏡面反射面形光學(xué)三維測量技術(shù)綜述［J］.光學(xué)儀器，2005，27（2）：90－95.

［15］ 蔣尊興，郭紅衛(wèi)，陳明儀.自適應(yīng)空間載波相移算法［J］.光學(xué)儀器，2006，28（6）：54－58.

［16］ ROBINSON D W，REID G.Interferogram analysis：digital fringe pattern measurement［M］.Bristol：IOP，1993：94－140.

［17］ LADAK H M，DECRAEMER W F，DIRCKX J J J，et al.Systematic errors in small deformations measured by use of shadow－Moire topography［J］.Appl Opt，2000，39（19）：3266－3275.

［18］ GUO H W.Phase－shifting deflectometric Moirétopography［J］.Structural Longevity，2011，5（1）：

39－48.