B樣條曲線的光順設(shè)計_浙江大學(xué)本科生畢業(yè)論文

1、本 科 生 畢 業(yè) 論 文（設(shè)計）文獻綜述和開題報告題目 B 樣條曲線的光順設(shè)計姓名與學(xué)號 指導(dǎo)教師 年級與專業(yè) 所在學(xué)院 浙江大學(xué)本科生畢業(yè)論文（設(shè)計）誠信承諾書1.本人鄭重地承諾所呈交的畢業(yè)論文（設(shè)計），是在指導(dǎo)教師的 指導(dǎo)下嚴格按照學(xué)校和學(xué)院有關(guān)規(guī)定完成的。2.本人在畢業(yè)論文（設(shè)計）中引用他人的觀點和參考資料均加 以注釋和說明。3. 本人承諾在畢業(yè)論文（設(shè)計）選題和研究內(nèi)容過程中沒有抄 襲他人研究成果和偽造相關(guān)數(shù)據(jù)等行為。4. 在畢業(yè)論文（設(shè)計）中對侵犯任何方面知識產(chǎn)權(quán)的行為，由 本人承擔相應(yīng)的法律責任。畢業(yè)論文（設(shè)計）作者簽名： 年 月 日一、 畢業(yè)論文二、 文獻綜述 三、 開題報告

2、四、 外文翻譯中文摘要在工業(yè)生產(chǎn)幾何設(shè)計中，人們大量應(yīng)用 B 樣條等數(shù)學(xué)工具來設(shè)計曲線。在 許多工業(yè)設(shè)計領(lǐng)域，比如造船業(yè)、汽車制造業(yè)等，設(shè)計師要求曲線要十分光順， 因為曲線的光順性會直接影響生成曲面的質(zhì)量。在本文中，我們先對光順的判 別準則進行討論，結(jié)合實際生產(chǎn)經(jīng)驗以得出更符合實際生產(chǎn)需要的光順判別準 則?；谒玫呐袆e準則，使用 L0/L1 范數(shù)優(yōu)化等理論給出目標函數(shù)，建立數(shù) 學(xué)模型。然后我們進行大量實驗，證明該模型的有效性。關(guān)鍵詞：B 樣條曲線光順L1 范數(shù)L0 范數(shù)壓縮感知 船體放樣AbstractB-spline are a modeling tool widely used in i

3、ndustrial geometric design. In many industrial design activities, curves need to be fair enough. The fairness of curves has a direct influence on the quality of the underlying surface. In this paper, we first discuss about the criterion of fairness, and then we present a reasonable criterion of fair

4、ness by considering the experience of actual production. According to the new criterion of fairness, we present a mathematical model by using the theories about L0/L1-norm optimization. We carry out large number of experiments which show that our solution is efficient.Key words: B-spline, fairness,

5、L1-norm, L0-norm, hull lofting目錄1 引言11.1 研究的目的和意義11.2 研究的問題和框架12 背景介紹與相關(guān)工作22.1 兩種通用光順判別準則22.1.1 光順判別準則 C1（N. Sapidis 等）22.1.2 光順判別準則 C2（蘇步青等）22.2整體能量優(yōu)化法32.3局部修改方法33 光順的判別準則（董光昌）33.1實例分析與光順判別準則 C333.1.1 典型的非光順曲線實例33.1.2 光順判別準則 C353.2函數(shù)樣條的光順算法簡介53.2.1 回彈法53.2.2 直尺卡樣法63.2.3 曲尺卡樣法73.2.4 一般曲線的光順設(shè)計算法83.3實

6、驗結(jié)果示例93.4算法的局限性94 基于 L1 范數(shù)優(yōu)化的光順算法104.1向量 c 與向量 e104.2數(shù)學(xué)模型的構(gòu)建104.2.1 光順判別準則 C3 與優(yōu)化目標114.2.2 曲線的拐點數(shù)與 c 的 L0 范數(shù)114.2.3 曲線的拐點數(shù)與 e 的 L0 范數(shù)144.2.4 曲線的拐點數(shù)與 e 的 L2 范數(shù)144.2.5 基于 L1 范數(shù)的優(yōu)化目標164.3實證結(jié)果分析164.3.1 c 的 L2 范數(shù)與 c 的 L1 范數(shù)164.3.2 e 的 L2 范數(shù)與 e 的 L1 范數(shù)204.4數(shù)學(xué)模型的優(yōu)化215 實驗結(jié)果216 結(jié)束語24參考文獻24附錄（matlab 相關(guān)代碼）261

7、引言在工業(yè)生產(chǎn)幾何設(shè)計中，“光順”是設(shè)計師們十分關(guān)心的概念。如在造船行業(yè)中， 船體若不夠光順，那么船在航行時會受到更大的阻力，也更容易被海水腐蝕，極大降 低船體壽命。隨著技術(shù)發(fā)展，B 樣條、NURBS 曲線/曲面等在生產(chǎn)設(shè)計中發(fā)揮了巨大 作用；人們也越來越關(guān)心如何對著這些數(shù)字曲線進行光順。1.1 研究的目的和意義在工業(yè)生產(chǎn)中，為了追求生成曲面或曲線具有更良好的物理特性或其他特性，往 往要求曲面或曲線具有光順性。例如，在造船業(yè)中，若船的水線、站線等足夠光順， 能夠減少船行駛遇到的阻力，且能夠延長船體使用年限。在實際的工業(yè)生產(chǎn)中，光順設(shè)計即將已有曲線變得更光順是十分重要的步 驟。而往往只有經(jīng)驗豐富

8、的設(shè)計師或者工人才能將曲線光順好。設(shè)計師和工人們對已 有的光順準則并不滿意，認為只有靠經(jīng)驗才能解決光順的問題。我們希望結(jié)合實際生 產(chǎn)的經(jīng)驗，總結(jié)出更合理的光順判別準則，建立數(shù)學(xué)模型，最終能應(yīng)用到實際生產(chǎn)中 去。由于設(shè)計師和工人們對已有的光順準則并不認可，所以目前的光順方法并不信 任。因此，光順設(shè)計往往是依靠人力來完成。如果要光順一個較大的模型，比如船體 上的全部水線、站線、橫剖線，一般需要三個星期左右的時間。所以，我們希望利用 計算機進行輔助設(shè)計，找到一種合理的算法對曲線進行光順，減少設(shè)計師與工人的工 作量。B 樣條曲線被大量應(yīng)用于工業(yè)生產(chǎn)設(shè)計，我們將針對 B 樣條曲線的光順進行討論， 并找到

9、有效的辦法對 B 樣條曲線進行光順。論文開題的目的和意義，即研究出更符合實際生產(chǎn)需求的光順判別準則與基于該 準則的曲線光順設(shè)計算法，以解決實際生產(chǎn)中的問題，并將光順的概念數(shù)學(xué)規(guī)范化。1.2 研究的問題與框架在造船業(yè)中，放樣工人們會得到一系列的插值點，他們的工作是求得光順的曲線 來通過這些插值點。傳統(tǒng)的放樣中，一般是先在地板上畫出這些插值點（也叫型值點） 的位置，然后用細木條依次通過這些插值點。放樣工人們通過肉眼來判斷這些細木條 是否光順，如果不夠光順，在通過調(diào)整插值點的位置來達到光順的效果。我們研究的問題是：給定了一個點列 Pi(xi, yi) ， i = 1,.n ，我們需求的新的點列Pi(

10、xi, yi) ，其與 Pi 的距離不超過預(yù)先給定的容差值，使得我們用 B 樣條或者 NURBS等曲線對其進行插值，所得到的曲線是光順曲線。因此我們稱之為光順設(shè)計。 我們在第 2 章，會介紹該問題的背景以及相關(guān)的工作。給出了兩種通用的判別準則，以及基于這兩種判別準則的相關(guān)算法。在第 3 章，我們會介紹董光昌先生對光順的研究成果。我們的研究主要受到了董 先生的啟發(fā)，我們的算法是基于他給出的光順判別準則。在第 4 章，我們將依照董光昌先生給出的光順判別準則，使用 L0 范數(shù)與 L1 范數(shù) 對其建模，從而得到我們的光順算法。在第 5 章，則是實驗結(jié)果與分析。2 背景介紹與相關(guān)工作2.1 兩種通用光順

11、判別準則曲線“光順”的概念主要來源于實際的生產(chǎn)，并沒有一個明確的數(shù)學(xué)定義?，F(xiàn)在 普遍使用曲線光順的判別準則有兩種，分別由 G.Farin 等1，與蘇步青等人提出。2.1.1 光順的判別準則 C1（G.Farin 等）GFarin 等人認為，曲線的光順應(yīng)該滿足以下條件：（1）曲線 f ÎG2 ，即曲線至少有 2 階或 2 階以上的幾何連續(xù)性；（2） min ò k 2ds，在滿足一定的容差約束條件下。2.1.2 光順的判別準則 C2（蘇步青等）蘇步青等人認為，曲線的光順應(yīng)該滿足以下條件：（1） 曲線 f ÎG2 ，即曲線至少有 2 階或 2 階以上的幾何連續(xù)性；（2

12、） 曲線 f的拐點少；（3） 曲線 f的曲率變化均勻。 兩種不同的判別準則都被廣泛使用?；诓煌墓忭樑袆e準則，可以將現(xiàn)有的光順算法分為兩類：整體能量優(yōu)化法與局部調(diào)整法。2.2 整體能量優(yōu)化法整體能量優(yōu)化法是基于判別準則 C1 的能量優(yōu)化方法，主要有應(yīng)變能法、最小二 乘法和小波法。應(yīng)變能法即在滿足一定的容差約束條件下，使應(yīng)變能曲率平方的積分極小， 即判別準則 C1 的直接建模。最小二乘法則是在滿足一定的容差約束條件下，插值點（或控制點）處的曲率平 方之和求極小，即對判別準則 C1 的離散建模。小波法則是將曲線用小波進行分解，然后去除細節(jié)部分的小波。2.2 局部修改方法局部修改方法主要基于光順的

13、判別準則 C2。Kjellander2的工作給出了均勻參數(shù) 三次 B 樣條的光順方法。該方法的基本思想是曲線的幾何外形在大多數(shù)型值點處是光 順或比較光順的，只是在少數(shù)型值點處非光順，逐次找出這些非光順的點即“壞點”， 修改這些“壞點”，使曲線達到光順的要求。之后許多有許多人在發(fā)展了 Kjellander 的這套方法。G.Farin 在 Kjellander 方法的基礎(chǔ)上，給出了一種節(jié)點去除法。它是將一 部分壞點取出后，重新計算 B 樣條的控制頂點以實現(xiàn)曲線光順。N. Sapidis 與 G. Farin6結(jié)合以上兩種技巧，用曲率線圖的方式找到引起非光順的 “壞點”，每次對最“壞”的點進行處理，

14、或者修改其位置，或者刪除該點，或者在 附近增加輔助控制點，使其局部的拐點減少且曲率變得均勻。整體能量優(yōu)化法由于往往不考慮將曲率變均勻，或者去除拐點，它并不被工人們 所接受；而局部修改方法的缺點在于每次僅處理一點，運行效率慢，且沒有考慮整體， 因此，這些光順方法始終沒有解決工人們的實際需求。3 光順的判別準則（董光昌）3.1 實例分析與光順的判別準則 C33.1.1 典型的非光順曲線實例董光昌先生曾經(jīng)在船廠從事數(shù)學(xué)放樣的工作，他將所有非光順的曲線總結(jié)為三 類，并得到具有豐富光順經(jīng)驗的工人們的認可。第一類非光順曲線：曲線拐點較多。即使對沒有光順經(jīng)驗的的人來說，看到曲線 凹凸不平，也會認為該曲線并不

15、光順。根據(jù)拐點的定義，曲線通過拐點將改變曲線的 凹凸性，因此造成凹凸不平的原因是曲線的拐點較多。如 sin x ，它屬于 C¥ ，即無 窮階光滑，但是應(yīng)該不是光順曲線。圖 1 sin x第二類非光順曲線：曲率線的拐點較多。對于曲線 f= 0.5x2 + sin x，易得該曲線沒有拐點，其仍然令人直觀上覺得不夠“光順”，它像一條蛇一樣纏繞著 0.5x2 。這是因為其曲率線的拐點較多。曲率為當參數(shù)為弧長參數(shù)時的二階導(dǎo)數(shù)，對于一般參數(shù)， 曲率3k =f ¢¢（1）(1+ f ¢2 )2當 f ¢ 絕對值較小時，可以用一般二階導(dǎo)數(shù)近似曲率值。圖 2 紅

16、色線表示 f= 0.5x2 + sin x ,藍色線為 0.5x2第三類非光順曲線：曲率的變化幅度大。如圖中的圓弧曲線，在 x=0 左邊，其為 半徑為 1 的圓，在 x=0 右邊，為半徑為 5 的圓。由于其曲率變化幅度較大，讓感覺左 邊部分較“鼓”，右邊部分較“癟”，不夠光順。圖 3 鼓癟段實例3.1.2 光順的判別準則 C3根據(jù)實例所示的非光順曲線，為避免這些非光順情形，董光昌先生給出了新的光 順判別準則。判別準則 C1 與判別準則 C2 都將曲線限定了至少 2 階的幾何連續(xù)性， 這樣在工業(yè)設(shè)計中常用的圓弧曲線就無法滿足其要求。我們可以僅僅要求曲線有 1+ 階的幾何連續(xù)性， >0。光順

17、的判別準則 C3：在滿足一定的約束條件小，（1） 曲線 f ÎG1+d ， d > 0 ；（2） 曲線 f的拐點少；（3） 曲線 f的曲率線的拐點少；（4） 曲率 f 的曲率線變化幅度小。其中，我們稱曲線 f 的拐點為第一振動數(shù)，曲線 f 的曲率線的拐點為第二振動數(shù)。 第一振動數(shù)比第二振動數(shù)與曲率線的變化幅度要重要得多，因為第一振動數(shù)說明了曲 線本身發(fā)生了質(zhì)變。因此，減少曲線本身的拐點數(shù)應(yīng)該是我們的第一目標。在下文中，我們用向量 c 來表示插值點（或控制頂點）處的曲率構(gòu)成的向量，用 向量 e 來表示曲率線的曲率構(gòu)成的向量。如果 c 或者 e 不存在，那么用有限差分等方 法對其近

18、似替代。3.2 函數(shù)樣條的光順方法簡介基于以上的判別準則，董光昌先生等人以函數(shù)樣條為主要對象，得出一種十分有 效的光順算法。這個算法又由回彈法、直尺卡樣法與曲尺卡樣法構(gòu)成9-12。我們的 主要工作受到了董先生這套方法的啟發(fā)，因此，我們在此簡單地介紹回彈法、直尺卡樣法和曲尺卡樣法。3.2.1 回彈法由公式（1）可知，當曲線的一階導(dǎo)數(shù)絕對值較小時，其曲率可以近似為其二階導(dǎo)數(shù)。令 Pi(xi, yi)，i = 1,.n表示插值點，向量 PiPi + 1 （ i = 2,.n -1） 分別與向量P0 P1構(gòu)成拐角ai，a = maxai稱作這一系列點的最大拐角。當最大拐角不超過 120度時，三次函數(shù)樣

19、條的曲率與其二階導(dǎo)數(shù)十分接近，可以互相替代。 而三次函數(shù)樣條的二階導(dǎo)數(shù)可以用追趕法得到，從而得到插值點處的二階導(dǎo)數(shù) cici + 1 - ci（ i = 1.n），與 ei =xi + 1 - xi（ i = 2.n ） ?；貜椃ǖ倪^程如下：（1）初始化 i = 1 ，轉(zhuǎn)（2）；（2）如果 eiei + 1 < 0，轉(zhuǎn)（3）；否則轉(zhuǎn)（5）；（3）求得最小的d i所得的值，轉(zhuǎn)（4）；，使 eiei + 1=0 成立，其中 ei 為 yi + d i 代入后的重新計算（4） d = min(x, 0.5d i) ， x 為預(yù)設(shè)的最大偏移量。更新（5）如果 i<n-1，i=i+1，轉(zhuǎn)（2

20、）；否則結(jié)束。3.2.2 直尺卡樣法yi = yi + d，轉(zhuǎn)（5）；直尺卡樣法的目的是減少多余的拐點，其分為三種應(yīng)對不同情況的過程。其過程 一如下：（1）初始化 i = 2 ，轉(zhuǎn)（2）；ci = 0（2）如果 ci-1ci <= 0且cici + 1 <= 0，轉(zhuǎn)（3）；否則轉(zhuǎn)（5）；（3）求得最小的 d i，使 的值，轉(zhuǎn)（4）；成立，其中 ci 為 yi + d i代入后的重新計算所得（4） d = min(x ,d i) ， x 為預(yù)設(shè)的最大偏移量。更新（5）如果 i<n-1，i=i+1，轉(zhuǎn)（2）；否則結(jié)束。yi = yi + d，轉(zhuǎn)（5）；其過程二如下：（1）初始化

21、i = 2，轉(zhuǎn)（2）；（2）如果 ci-1ci < 0、cici + 1 > 0且ci + 1ci + 2 < 0，轉(zhuǎn)（3）；否則轉(zhuǎn)（5）；ci +1 = 0（3）求得最小的 d i，使成立，其中 ci 為 yi + d i代入后的重新計算所ci = 0得的值；求得最小的 d i +1，使 算所得的值，轉(zhuǎn)（4）；成立，其中 ci 為 yi +1 + d i +1代入后的重新計（4）如果| d i |<|d i + 1 |，則 yi = yi + d i，否則 yi + 1 = yi + 1 + d i + 1 ；然后在使用過程一中的方法，去除拐點；轉(zhuǎn)（5）；（5）如果

22、i<n-2，i=i+1，轉(zhuǎn)（2）；否則結(jié)束。其過程三的主要面對的情況是 ci - 1ci < 0, cici +1 > 0,., cici + k - 1 > 0, cici + k < 0 ，k>2，一般認為 k<=5。首先，利用回彈法修改插值點使 ei+1 = 0,., ei + k - 2=0 。這樣，相當于將中間的插值點 Pi + 1,., Pi + k - 2去除。3.2.3 曲尺卡樣法刪除，于是可以用過程二的方法將拐點給我們將直尺卡樣法中的 c 都用 e 來替代，就是曲尺卡樣法。即曲尺卡樣法就是 對曲率線的直尺卡樣法，而曲尺卡樣法一般只執(zhí)行

23、過程一和過程二。其過程一如下：（1）初始化 i = 2 ，轉(zhuǎn)（2）；ei = 0（2）如果 ei-1ei <= 0且eiei + 1 <= 0，轉(zhuǎn)（3）；否則轉(zhuǎn)（5）；（3）求得最小的 d i，使 的值，轉(zhuǎn)（4）；成立，其中 ei 為 yi + d i代入后的重新計算所得（4） d = min(x ,d i) ， x 為預(yù)設(shè)的最大偏移量。更新（5）如果 i<n-1，i=i+1，轉(zhuǎn)（2）；否則結(jié)束。 其過程二如下：yi = yi + d，轉(zhuǎn)（5）；（1）初始化 i = 2 ，轉(zhuǎn)（2）；ei = 0（2）如果 ei-1ei < 0、eiei + 1 > 0且ei +

24、1ei + 2 < 0，轉(zhuǎn)（3）；否則轉(zhuǎn)（5）；ei +1 = 0（3）求得最小的 d i，使成立，其中 ei 為 yi + d i代入后的重新計算所得的值；求得最小的 d i +1，使 計算所得的值，轉(zhuǎn)（4）；成立，其中 ei +1 為 yi +1 + d i +1代入后的重新（4）如果| d i |<|d i + 1 |，則 yi = yi + d i，否則 yi + 1 = yi + 1 + d i + 1 ；然后在使用過程一中的方法，去除拐點；轉(zhuǎn)（5）；（5）如果 i<n-2，i=i+1，轉(zhuǎn)（2）；否則結(jié)束。3.2.4 一般曲線的光順設(shè)計算法在實際生產(chǎn)過程中，我們遇到

25、的曲線很多都不是函數(shù)型曲線。盡管三次函數(shù)樣條 還有二次導(dǎo)數(shù)平方的積分最小這樣良好的性質(zhì)，它面對最大拐角大于 120 度的插值點 序列時，上述的光順方法就顯得也無能為力了。因此，我們需要對插值點進行分段，逐段進行光順，然后將結(jié)果拼接在一起。為 了避免拼接后兩段之間有明顯的差異，處理的方法是相鄰的兩段間共用 4 到 5 個插值 點。比如，1 到 10 為一段，下一段為 7 到 15，這樣 7、8、9、10 四個點是共用的插 值點。分段光順后，被共用的點可能對應(yīng)有有兩個不同的新位置，這兩個新位置的中 間點的位置為該共用點的新位置，這個步驟稱之為混合。圖 4 該線分成了兩段（紅與綠），黑色部分為共用點

26、其算法的主要過程如下：（1）將 Pi(xi, yi) 分段，每相鄰兩段有 4 到 5 個共用點，且每一段滿足該段的插值 點列最大拐角不超過 120 度；（2）對于每一段，分別對其依次使用：回彈法、直尺卡樣法、曲尺卡樣法進行 光順。（3）每一段得到的光順后的結(jié)果，在共用點處進行混合。圖 5對圖 4 所示曲線的光順結(jié)果的曲率線圖，藍色為原曲率線，紅色為光順后曲率線。左圖為分段后沒有共用點所得的結(jié)果，右圖為分段共用 4 個點并混合后所得的結(jié)果，右圖的拐點顯然比左圖少。3.3 實驗結(jié)果示例圖 6 鞋印，左圖是原曲線，右圖是曲率線。紅色為光順后，藍色為光順前。圖 7圖 83.4 算法的局限性如 3.2.

27、4 所述，對于一般曲線，這套方法需要先分段再混合。這樣，在相鄰兩段的拼接處，原來的光順結(jié)果可能會被破壞，從而得不到更好的結(jié)果。圖 9 左圖是圖 4 的例子加了噪音后的原曲線，右圖是曲率線圖。藍色表示光順前的曲率線，紅色則表示按 3.2.4 所述算法光順后所得的曲率線。即使有分段和混合，中間部分的光順性反而更差了。該算法比起第 2 章所述的整體方法，更多地考慮了局部性；而比起局部方法，又 是一段一段地去光順，考慮了部分地整體性。因此，它可以說是“大局部”方法。我 們希望得到類似的方法，既考慮局部，又不忽略全局整體性，從而得到更好且更魯棒 的光順效果。4 基于 L1 范數(shù)優(yōu)化的光順算法4.1 向量

28、 c 與向量 e我們已知插值點（或型值點） Pi (xi , yi ) (i = 1,., n)，我們希望反求出這些型值點位置與其對應(yīng)的曲率向量間的關(guān)系。如果我事先給定了參數(shù) u1,.,un，它可以是均勻也可以是非均勻的。 那么，由參考文獻1，由三次 B 樣條的基本性質(zhì)，在每個插值T點 Pi處的一階導(dǎo)數(shù)與二階導(dǎo)數(shù)都可由 P = x1, y1，., xn , yn 乘上一個矩陣得到，這個矩陣僅與參數(shù) u1,.,un 有關(guān)，與 Pi 所處的位置無關(guān)。又由參數(shù)平面曲線 s(t) = (x(t), y(t) ，其曲率為k = x¢(t) y ¢(t) - x ¢(t) y

29、¢(t)(x¢2 (t)+y¢2 (t)3/ 2盡管一階導(dǎo)數(shù)與二階導(dǎo)數(shù)與 P 是線性關(guān)系，但是曲率與它卻是非線性關(guān)系。為了 使優(yōu)化更為簡單，一種方法是將二階導(dǎo)數(shù)當作是曲率；另一種方法則是在優(yōu)化過程中 將一階導(dǎo)數(shù)作為常數(shù)來處理，因為當型值點的偏移量較小時，可以認為一階導(dǎo)數(shù)不發(fā) 生太大的變化。令 c 為型值點處曲率所構(gòu)成的向量。當我們將曲線的一階導(dǎo)數(shù)作為常數(shù)，曲率向 量 c 與 P 滿足線性關(guān)系，即可以得到一個矩陣 A，使得 c=AP。具體如何求，可以參 考附錄中所附的 matlab 代碼。而 ei = (ci+1 - ci ) / (ui+1 - ui ) ，當 u

30、1,.,un 為事先給定的常數(shù)時，它與曲率向量 c 存在線性關(guān)系。因此，存在矩陣 B，使得 e=BP 成立。具體如何求，可以參考附錄中所 附的 matlab 代碼。4.2 數(shù)學(xué)模型的構(gòu)建4.2.1 光順判別準則 C3 與優(yōu)化目標如第 3 章中所述，光順判別準則 C3 有三個優(yōu)化目標：（1）曲線的拐點數(shù)；（2） 曲率線的拐點數(shù)；（3）曲率線的變化幅度。而第三章中，董光昌先生等提出的曲線光順方法，回彈法對應(yīng)減少曲率線變化的 幅度，直尺卡樣法對應(yīng)減少曲線的拐點數(shù)，曲尺卡樣法對應(yīng)減少曲率線的變化幅度。 我們先考慮曲率線的變化幅度，這是三個優(yōu)化目標中相對容易量化的目標。向量e 中的每一個值表示該插值點處

31、曲率的變化幅度，因此，減少向量 e 的 L1 范數(shù)或 L2范數(shù)，即可減少曲率線的變化幅度；最小化| e |1化幅度?；蛘遼 e |2 即可對應(yīng)最小化曲率的變再考慮前兩個優(yōu)化目標，它們都是希望減少拐點數(shù)。曲線的拐點數(shù)，即第一振動數(shù)，我們用 Nc來表示；曲率線的拐點數(shù)，即第二振動數(shù)，我們用 Ne來表示。由第3 章所述，第一振動數(shù)代表的是曲線光順性發(fā)生“質(zhì)變”的次數(shù)，它是最重要的優(yōu)化 目標。因此，我們首要考慮如何優(yōu)化 Nc 。而第一振動數(shù)與第二振動數(shù)是相似的，只 要我們知道如何優(yōu)化 Nc ，總可以用類似的方法去優(yōu)化 Ne ，就如直尺卡樣法和曲尺卡 樣法一樣。4.2.2 曲線的拐點數(shù)與 c 的 L0

32、范數(shù)曲率向量 c 的 L0 范數(shù)，即曲率向量 c 中非零元的個數(shù)。圖 10 曲率線圖 c=1,2,-1,1,-1,3圖 11 曲率線圖 c=1,2,0,1,-1,3圖 10 中的拐點數(shù)為 4，而其 c 的 L0 范數(shù)為 6。圖 11 中拐點數(shù)為 2，其 c 的 L0 范 數(shù)為 5。易有不等式Nc<=| c |0 - 1因為只有當向量 c 中相鄰的兩個非零元的符號不同，才會產(chǎn)生一個拐點，那么當 c 的 L0 范數(shù)固定，最多只會產(chǎn)生非零元個數(shù)減一個拐點。從這點來看，c 的 L0 范數(shù) 為拐點數(shù)的上界，優(yōu)化 c 的 L0 范數(shù)，能在一定程度上減少拐點數(shù)。如果我們將向量 c 中某個非零元 ci

33、變成 0，而其他值保持同號：（1）當 ci 的兩個相鄰的非零元都與它異號時，它變成零，拐點數(shù)會減少 2。如將 圖 10 中第 3 個元素-1 變成 0 就變成了圖 11，而圖 11 的拐點數(shù)比圖 10 的拐點數(shù)少了 2。（2）當 ci 得相鄰非零元有一個與它同號，它變成零，拐點數(shù)不會改變。 我們知道優(yōu)化 c 的 L0 范數(shù)，可以近似轉(zhuǎn)化為優(yōu)化 c 的 L1 范數(shù)。當我們?nèi)?yōu)化 c的 L1 范數(shù)時，一般來說，c 中的每個元素在優(yōu)化后不會與原來異號，或者是同號， 或者是變成零，即確保不會增加拐點數(shù)。根據(jù)上述討論，優(yōu)化 c 的 L1 范數(shù)，應(yīng)該能 夠有效地減少曲線的拐點數(shù)。下面圖示是實驗論證。圖 1

34、2 左邊為原曲線，右邊為曲率圖。藍色為優(yōu)化 c 的 L1 范數(shù)前，紅色為優(yōu)化后。圖 13 左邊為原曲線，右邊為曲率圖?？梢钥闯?，優(yōu)化后拐點明顯減少。圖 14 左邊為原曲線，右邊為曲率圖。通過實驗，我們認為優(yōu)化 L0 范數(shù)，可以有效地去除多余的拐點。4.2.3 曲線的拐點數(shù)與 e 的 L0 范數(shù)從上小節(jié)的分析，我們知道減少曲率向量 c 的 L0 范數(shù)，可以有效地減少曲線的 拐點數(shù)。而向量 e 可以看成是曲率線的曲率向量，減少向量 e 的 L0 范數(shù)，可以有效 減少曲率線的拐點數(shù)。當曲線的拐點數(shù)很多時，向量 c 中元素符號經(jīng)常發(fā)生改變，曲率線應(yīng)該是沿著橫 軸上下振動，此時曲率線本身的拐點也很多。反

35、之，當曲率線的拐點較少時，相對應(yīng) 的曲率線的振動也會較少，那么曲率線穿越橫軸的振動也很可能相應(yīng)的較少，如此， 曲線的拐點數(shù)就會比較少。因此，減少 e 的 L0 范數(shù)，也可能會減少曲線的拐點數(shù)。4.2.4 曲線的拐點數(shù)與 e 的 L2 范數(shù)由前文可知，e 的 L2 范數(shù)應(yīng)該直接與曲率線的變化幅度相關(guān)。在本節(jié)中，我們認 為減少 e 的曲率的變化幅度，也能在一定程度上減少曲線的拐點數(shù)。我們在本節(jié)之所以考慮 e 的 L2 范數(shù)，而不是 e 的 L1 范數(shù)，是為了不受 e 的稀疏 性（L0 范數(shù)）影響。圖 15 中所示的曲率線圖有兩個拐點，在中間點處的曲率變動幅度較大。而在未 被光順的曲線中，在局部類似

36、這樣的曲率線圖有很多。當我們減小中間點的變化幅度， 只要減少達到一定的量，如圖 16 所示，就很可能將拐點減少。圖 15 曲率線圖。曲率向量 c=2,-1,2圖 16 曲率線圖。紅色曲率向量為 c=1.8,0.1,1.8我們進行以下實驗：在同一容差條件下，即型值點的最大偏移量不超過一個固定 值，對比優(yōu)化 c 的 L1 范數(shù)與 e 的 L2 范數(shù)的結(jié)果。圖 17 實驗對象的原曲線圖圖 18：最大偏移量為 0.001 時的曲率線圖。藍色為原曲率，綠色為 c 的 L1 范數(shù)優(yōu)化后結(jié)果，紅色為 e 的 L2 范數(shù)優(yōu)化后的結(jié)果。圖 19：最大偏移量為 0.0028 時的曲率線圖?？梢钥闯觯t線 的拐點數(shù)

37、已經(jīng)減少 2，而綠線與藍線拐點數(shù)相同。實驗說明，在有些情況下，優(yōu)化 e 的 L2 范數(shù)會比優(yōu)化 c 的 L1 范數(shù)更快地減少拐 點數(shù)。4.2.5 基于 L1 范數(shù)的優(yōu)化目標綜合以上所述，光順的判別準則 C3 可以分別對應(yīng)三個優(yōu)化目標：| c |0| e |2 。、| e |0 和我們知道，關(guān)于 L0 范數(shù)的優(yōu)化可以轉(zhuǎn)化為 L1 范數(shù)的優(yōu)化，因此，目標函數(shù) | c |0可以用| c |1 來替代以進行優(yōu)化。而| e |0 可以用| e |1 來替代進行優(yōu)化。而 e 的 L1 范數(shù)代表了曲率線的變化幅度， 因此，我們可以不用再去優(yōu)化 | e |2 。光順判別準則 C3 的后兩項目標，我們可以用 |

38、 e |1 來表示。由此，我們確定了曲線光順問題的兩個優(yōu)化目標函數(shù)：| c |1 、| e |1 。4.3 實證結(jié)果分析4.3.1 c 的 L2 范數(shù)與 c 的 L1 范數(shù)如第 2 章所述，如果我們基于光順判別準則 C1，那么可以通過優(yōu)化 c 的 L2 范數(shù) 來達到光順的效果。而在光順判別準則 C3 中，減少拐點數(shù)才是我們的光順的首要目標。以下的實驗圖，是在同一容差條件下，分別優(yōu)化 c 的 L2 范數(shù)與 c 的 L1 范數(shù)的對 比結(jié)果。圖 20（a）原曲線圖圖 20(b) 優(yōu)化 L2 范數(shù)后曲率線與原曲率線對比圖。藍色為原曲率，綠色為優(yōu)化后。圖 20(c) 優(yōu)化 L1 范數(shù)后曲率線與原曲率線對

39、比圖。藍色為原曲率，紅色為優(yōu)化后。如圖 20 所示，圖 20（c）中的紅線的拐點數(shù)明顯比圖 20（b）中的綠線拐點數(shù)要 少。從減少拐點數(shù)這個角度來說，優(yōu)化 c 的 L1 范數(shù)比優(yōu)化 c 的 L2 范數(shù)更加有效。圖 21(a) 原曲線圖 21(b)c 的 L2 范數(shù)優(yōu)化后曲率對比圖圖 21（c）c 的 L1 范數(shù)優(yōu)化后曲率對比圖綜合圖 21 的例子，可以明顯看出，在同一容差條件下（每個點的偏移量不超過(max( y) - min( y) /1000 ），優(yōu)化 c 的 L1 范數(shù)比起優(yōu)化 c 的 L2 范數(shù)更能夠減少第一振動數(shù)，即曲線的拐點數(shù)。4.3.2 e 的 L2 范數(shù)與 e 的 L1 范數(shù)圖

40、 22 曲率線圖。藍色為原曲率，綠色為 e 的 L2 范數(shù)優(yōu)化，紅色為 e 的 L1 范數(shù)優(yōu)化。圖 22 的原曲線為圖 17 中的曲線，圖 22 設(shè)定偏移量不超過 0.0027 的結(jié)果?？梢?看出，優(yōu)化 e 的 L1 范數(shù)比優(yōu)化 e 的 L2 范數(shù)可能更有效地減少拐點。圖中藍、綠線穿 過了橫軸 2 次，而紅線沒有穿過橫軸。綜合以上所述，我們討論了在光順問題上優(yōu)化 L1 范數(shù)相比 L2 范數(shù)的優(yōu)越性。因 為 L1 范數(shù)還與 L0 范數(shù)有直接相關(guān)的關(guān)系，而 L0 范數(shù)是離散的量，它的變化相比起 L2 范數(shù)的變化是一種“質(zhì)變”，而我們想要減少的拐點，也是曲線發(fā)生了“質(zhì)變”的 點，因此，L1 范數(shù)會比

41、 L2 范數(shù)更加有效。4.4 數(shù)學(xué)模型的優(yōu)化綜合以上的模型的討論，我們可以建立多目標優(yōu)化的數(shù)學(xué)模型：其中矩陣 A，B 為 4.1 中所討論的矩陣，c0、e0分別為初始值，所求的 d為偏移量， 為給定的偏移量最大范圍，一般不能太大，因為根據(jù) 4.1，我們在優(yōu)化過程中 將一階導(dǎo)數(shù)當成常數(shù)來處理，如果型值點位置變化較大，讓一階導(dǎo)數(shù)也變化較大，會 明顯影響到結(jié)果。我們采用網(wǎng)上提供的 cvx matlab 凸優(yōu)化包來對其進行優(yōu)化。具體的代碼會在附錄 中給出。5 實驗結(jié)果本章我們將比較我們的算法與其他算法的結(jié)果。這里我們參照的算法主要是最小 二乘法（曲率的 L2 范數(shù)求極小）、N. Sapidis 與 G

42、. Farin 的局部修改法。這兩種算法分 別屬于整體方法和局部方法，也是目前光順方法中最具代表性的兩種算法。我們進行了大量實驗，證明了我們的算法的在光順效果上的優(yōu)越性。圖 23 左圖是曲線圖；右圖是三種方法光順后的曲率線圖。藍色為原曲率，紅色為我們的光順方法， 黑色為 N. Sapidis 與 G. Farin 的局部修改法，綠色為最小二乘法（求 c 得 L2 范數(shù)極?。?。曲率棒對比圖如下：圖 24 曲率棒對比圖。左上為原圖，右上為 N. Sapidis 與 G. Farin 的局部修改法的結(jié)果，左下為我們 的方法的結(jié)果，右下則是最小二乘法的結(jié)果。圖 25 左邊是原圖；右邊曲率線對比圖。顏色

43、對應(yīng)與圖 23 相同。圖 26 曲率棒對比圖。左上為原圖，右上為 N. Sapidis 與 G. Farin 的局部修改法的結(jié)果，左下為我們 的方法的結(jié)果，右下則是最小二乘法的結(jié)果。實驗證明，我們的光順方法比起其他兩種方法的光順效果更好。 我們的算法基于光順判別準則 C3，而董光昌先生的算法也是基于光順判別準則C3。我們認為在董先生的算法局限性的場合我們的算法依然能有效地光順，以下是實 驗結(jié)果。圖 27 曲率線圖。該圖的原曲線時圖 9 中的曲線。藍色為原曲率，綠色為董先生方法光順后的曲率， 紅色為我們的方法光順后所得的曲率。6 結(jié)束語我們的光順方法主要是對向量 c 與向量 e 的 L1 范數(shù)進

44、行多目標優(yōu)化。通過討論， 我們了解了影響曲線光順性的主要原因，并針對這些原因建立了數(shù)學(xué)模型。我們認為通過這篇論文所述的算法，曲線的光順問題已經(jīng)能夠得到較為良好的解 決。我們希望能將這套方法往曲面上推廣，從而直接解決曲面光順的問題。參考文獻1 G.Farin, Curves and Surfaces for CAGD:A Practical Guide. Academic Press, 5thedition, 2002.2 J.Kjellander, Smoothing of cubic parametric splinesJ. Computer Aided Design, 15(3):175-

45、179, 1983.3 H. Hagen, S.Hahmann and G.Schulze, Automatic smoothing with geometric surface patchesJ. Computer Aided Geometric Design, 4(3):231-235, 1987.4 J.Poliakoff, An improved algorithm for automatic fairing of nonuniform parametric cubic splinesJ. Computer-Aided Design, 28(1):59-66,1996.5 X.Ye,

46、Curvature continuous interpolation of curve meshesJ. Computer Aided Geometric Design, 14:169-190, 1997.6 N. Sapidis and G. Farin, Automatic fairing algorithm for b-spline curvesJ.Computer-Aided Design, 22(2):121-129, 1990.7 W. Renz, Interactive smoothing of digitized point dataJ, Computer Aided Desi

47、gn,14:267-269, 1982.8 J. Hoschek, Smoothing of curves and surfaceJ. Computer Aided Geometric Design, 15:175-179, 1983.9 C.Zhang, P.Zhang and F.Cheng, Fairing spline curves and surfaces by minimizing energyJ. Computer Aided Design, 33(13):919-923, 2001.10 董光昌，吳明華，洪安祥，樣條曲線光順的數(shù)學(xué)模型分析J，高校應(yīng)用數(shù)學(xué)學(xué)報 A 輯(中文版)，

48、2003(04)。11 董光昌，張但，劉志誠，馬利莊，Curve fairingJ, Progress in Natural Science, 1997(05)。12 劉志誠，董光昌，樣條曲線光順概念及指標回彈法的應(yīng)用J, 數(shù)學(xué)的實踐與 認識，1985(03)。13 董光昌，梁友棟，何援軍，樣條曲線擬合與雙圓弧逼近J，應(yīng)用數(shù)學(xué)學(xué)報，1978(04)。附錄（相關(guān)代碼）主函數(shù)：function nx,ny = l1_fairing( x,y,tolerance ) n = length(x);% 設(shè)定多目標優(yōu)化的迭代步驟kk = 8;%設(shè)定步長step = tolerance / kk / 2;%

49、 獲得參數(shù) Qx,Qy: Qx = 0;x;0; Qy = 0;y;0;% A : Px = A * Qx; Py = A * Qy;其中 Px、Py 為控制點的位置% t 為節(jié)點位置，即對應(yīng)文中的 u% N 為 A 的逆矩陣% N1 一階導(dǎo)數(shù)求導(dǎo)矩陣，即 Px' = N1 * Px; Py' = N1 * Py;% N2 二階導(dǎo)數(shù)求導(dǎo)矩陣 即 Px'' = N2 * Px; Py'' = N2 * Py;% DF 為差分矩陣，即 e = DF * c; A,Qx,Qy,t,N =Cubic_BSpline_M(x,y); N1,N2,DF =

50、Df_Mat(t);% 將 Qx，Qy 合并為一個長向量,矩陣對應(yīng)修改T = zeros(size(A); A = A,T;T,A;A = sparse(A); Q = Qx'Qy' P = A * Q;% c = matC * P;matC = GetCurvatureMatrix(P,N1,N2);% c = matC * Q; matC = matC * A;% e = matE * A; matE = DF * matC; n = length(P);% 初始值c0 = matC * Q; e0 = matE * Q;% 多目標迭代優(yōu)化for i = 1 : kkcvx

51、_beginvariable dy(n); minimize(norm(matC * dy + c0,1) subject tonorm(matE * dy + e0,1) <= norm(e0,1); norm(dy ,inf) <= step;cvx_end hold on;Q = Q + dy; P = A * Q;matC = getCurvatureMatrix(P,N1,N2); matC = matC * A;matE = DF * matC; c0 = matC * Q;e0 = matE * Q; cvx_beginvariable dy(n); minimize

52、(norm(matE * dy + e0,1) subject tonorm(dy ,inf) <= step;endcvx_endQ = Q + dy; P = A * Q;matC = getCurvatureMatrix(P,N1,N2); matC = matC * A;matE = DF * matC; c0 = matC * Q;e0 = matE * Q;end相關(guān)函數(shù)：function A,Qx,Qy ,t,N = Cubic_BSpline_M( x,y ) n=length(x);l=0;for i=1:n-1l=l+norm(x(i+1)-x(i),y(i+1)-y

53、(i);end t(1)=0;t(2)=0;t(3)=0; s=0;for i=4:n+2s=s+norm(x(i-2)-x(i-3),y(i-2)-y(i-3); t(i)=s/l;end t(n+3)=t(n+2);t(n+4)=t(n+2);h1=t(4)-t(3); h2=t(5)-t(4); h3=t(6)-t(5);f1=h1(-2); f2=-h1(-2)-h1(-1)*(h1+h2)(-1)-(h1+h2)(-2);f3=h1(-1)*(h1+h2)(-1)+(h1+h2)(-2)+(h1+h2)(-1)*(h1+h2+h3)(-1); f4=-(h1+h2)(-1)*(h1+h2+h3)(-1);Qx(1)=0;Qy(1)=0;h1=t(n+2)-t(n+1); h2=t(n+1)-t(n); h3=t(n)-t(n-1);g1=-h1(-2); g2=h1(-2)+h1(-1)*(h1+h2)(-1)+(h1+h2)(-2);g3