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

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

2、四、 外文翻譯中文摘要在工業(yè)生產(chǎn)幾何設(shè)計(jì)中，人們大量應(yīng)用 B 樣條等數(shù)學(xué)工具來(lái)設(shè)計(jì)曲線。在 許多工業(yè)設(shè)計(jì)領(lǐng)域，比如造船業(yè)、汽車制造業(yè)等，設(shè)計(jì)師要求曲線要十分光順， 因?yàn)榍€的光順性會(huì)直接影響生成曲面的質(zhì)量。在本文中，我們先對(duì)光順的判 別準(zhǔn)則進(jìn)行討論，結(jié)合實(shí)際生產(chǎn)經(jīng)驗(yàn)以得出更符合實(shí)際生產(chǎn)需要的光順判別準(zhǔn) 則。基于所得的判別準(zhǔn)則，使用 L0/L1 范數(shù)優(yōu)化等理論給出目標(biāo)函數(shù)，建立數(shù) 學(xué)模型。然后我們進(jìn)行大量實(shí)驗(yàn)，證明該模型的有效性。關(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 研究的問(wèn)題和框架12 背景介紹與相關(guān)工作22.1 兩種通用光順判別準(zhǔn)則22.1.1 光順判別準(zhǔn)則 C1（N. Sapidis 等）22.1.2 光順判別準(zhǔn)則 C2（蘇步青等）22.2整體能量?jī)?yōu)化法32.3局部修改方法33 光順的判別準(zhǔn)則（董光昌）33.1實(shí)例分析與光順判別準(zhǔn)則 C333.1.1 典型的非光順曲線實(shí)例33.1.2 光順判別準(zhǔn)則 C353.2函數(shù)樣條的光順?biāo)惴ê?jiǎn)介53.2.1 回彈法53.2.2 直尺卡樣法63.2.3 曲尺卡樣法73.2.4 一般曲線的光順設(shè)計(jì)算法83.3實(shí)

6、驗(yàn)結(jié)果示例93.4算法的局限性94 基于 L1 范數(shù)優(yōu)化的光順?biāo)惴?04.1向量 c 與向量 e104.2數(shù)學(xué)模型的構(gòu)建104.2.1 光順判別準(zhǔn)則 C3 與優(yōu)化目標(biāo)114.2.2 曲線的拐點(diǎn)數(shù)與 c 的 L0 范數(shù)114.2.3 曲線的拐點(diǎn)數(shù)與 e 的 L0 范數(shù)144.2.4 曲線的拐點(diǎn)數(shù)與 e 的 L2 范數(shù)144.2.5 基于 L1 范數(shù)的優(yōu)化目標(biāo)164.3實(shí)證結(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 實(shí)驗(yàn)結(jié)果216 結(jié)束語(yǔ)24參考文獻(xiàn)24附錄（matlab 相關(guān)代碼）261

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

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

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

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

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

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

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

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

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

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

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

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

19、條的曲率與其二階導(dǎo)數(shù)十分接近，可以互相替代。 而三次函數(shù)樣條的二階導(dǎo)數(shù)可以用追趕法得到，從而得到插值點(diǎn)處的二階導(dǎo)數(shù) cici + 1 - ci（ i = 1.n），與 ei =xi + 1 - xi（ i = 2.n ） 。回彈法的過(guò)程如下：（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 代入后的重新計(jì)算（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）；直尺卡樣法的目的是減少多余的拐點(diǎn)，其分為三種應(yīng)對(duì)不同情況的過(guò)程。其過(guò)程 一如下：（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代入后的重新計(jì)算所得（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）；其過(guò)程二如下：（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代入后的重新計(jì)算所ci = 0得的值；求得最小的 d i +1，使 算所得的值，轉(zhuǎn)（4）；成立，其中 ci 為 yi +1 + d i +1代入后的重新計(jì)（4）如果| d i |<|d i + 1 |，則 yi = yi + d i，否則 yi + 1 = yi + 1 + d i + 1 ；然后在使用過(guò)程一中的方法，去除拐點(diǎn)；轉(zhuǎn)（5）；（5）如果

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

23、過(guò)程一和過(guò)程二。其過(guò)程一如下：（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代入后的重新計(jì)算所得（4） d = min(x ,d i) ， x 為預(yù)設(shè)的最大偏移量。更新（5）如果 i<n-1，i=i+1，轉(zhuǎn)（2）；否則結(jié)束。 其過(guò)程二如下：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代入后的重新計(jì)算所得的值；求得最小的 d i +1，使 計(jì)算所得的值，轉(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 ；然后在使用過(guò)程一中的方法，去除拐點(diǎn)；轉(zhuǎn)（5）；（5）如果 i<n-2，i=i+1，轉(zhuǎn)（2）；否則結(jié)束。3.2.4 一般曲線的光順設(shè)計(jì)算法在實(shí)際生產(chǎn)過(guò)程中，我們遇到

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

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

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

28、 c 與向量 e我們已知插值點(diǎn)（或型值點(diǎn)） Pi (xi , yi ) (i = 1,., n)，我們希望反求出這些型值點(diǎn)位置與其對(duì)應(yīng)的曲率向量間的關(guān)系。如果我事先給定了參數(shù) u1,.,un，它可以是均勻也可以是非均勻的。 那么，由參考文獻(xiàn)1，由三次 B 樣條的基本性質(zhì)，在每個(gè)插值T點(diǎn) Pi處的一階導(dǎo)數(shù)與二階導(dǎo)數(shù)都可由 P = x1, y1，., xn , yn 乘上一個(gè)矩陣得到，這個(gè)矩陣僅與參數(shù) u1,.,un 有關(guān)，與 Pi 所處的位置無(wú)關(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)化更為簡(jiǎn)單，一種方法是將二階導(dǎo)數(shù)當(dāng)作是曲率；另一種方法則是在優(yōu)化過(guò)程中 將一階導(dǎo)數(shù)作為常數(shù)來(lái)處理，因?yàn)楫?dāng)型值點(diǎn)的偏移量較小時(shí)，可以認(rèn)為一階導(dǎo)數(shù)不發(fā) 生太大的變化。令 c 為型值點(diǎn)處曲率所構(gòu)成的向量。當(dāng)我們將曲線的一階導(dǎo)數(shù)作為常數(shù)，曲率向 量 c 與 P 滿足線性關(guān)系，即可以得到一個(gè)矩陣 A，使得 c=AP。具體如何求，可以參 考附錄中所附的 matlab 代碼。而 ei = (ci+1 - ci ) / (ui+1 - ui ) ，當(dāng) u

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

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

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

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

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

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

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

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

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

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

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

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

42、. Farin 的局部修改法。這兩種算法分 別屬于整體方法和局部方法，也是目前光順?lè)椒ㄖ凶罹叽硇缘膬煞N算法。我們進(jìn)行了大量實(shí)驗(yàn)，證明了我們的算法的在光順效果上的優(yōu)越性。圖 23 左圖是曲線圖；右圖是三種方法光順后的曲率線圖。藍(lán)色為原曲率，紅色為我們的光順?lè)椒ǎ?黑色為 N. Sapidis 與 G. Farin 的局部修改法，綠色為最小二乘法（求 c 得 L2 范數(shù)極?。Ｇ拾魧?duì)比圖如下：圖 24 曲率棒對(duì)比圖。左上為原圖，右上為 N. Sapidis 與 G. Farin 的局部修改法的結(jié)果，左下為我們 的方法的結(jié)果，右下則是最小二乘法的結(jié)果。圖 25 左邊是原圖；右邊曲率線對(duì)比圖。顏色

43、對(duì)應(yīng)與圖 23 相同。圖 26 曲率棒對(duì)比圖。左上為原圖，右上為 N. Sapidis 與 G. Farin 的局部修改法的結(jié)果，左下為我們 的方法的結(jié)果，右下則是最小二乘法的結(jié)果。實(shí)驗(yàn)證明，我們的光順?lè)椒ū绕鹌渌麅煞N方法的光順效果更好。 我們的算法基于光順判別準(zhǔn)則 C3，而董光昌先生的算法也是基于光順判別準(zhǔn)則C3。我們認(rèn)為在董先生的算法局限性的場(chǎng)合我們的算法依然能有效地光順，以下是實(shí) 驗(yàn)結(jié)果。圖 27 曲率線圖。該圖的原曲線時(shí)圖 9 中的曲線。藍(lán)色為原曲率，綠色為董先生方法光順后的曲率， 紅色為我們的方法光順后所得的曲率。6 結(jié)束語(yǔ)我們的光順?lè)椒ㄖ饕菍?duì)向量 c 與向量 e 的 L1 范數(shù)進(jìn)

44、行多目標(biāo)優(yōu)化。通過(guò)討論， 我們了解了影響曲線光順性的主要原因，并針對(duì)這些原因建立了數(shù)學(xué)模型。我們認(rèn)為通過(guò)這篇論文所述的算法，曲線的光順問(wèn)題已經(jīng)能夠得到較為良好的解 決。我們希望能將這套方法往曲面上推廣，從而直接解決曲面光順的問(wèn)題。參考文獻(xiàn)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é)報(bào) A 輯(中文版)，

48、2003(04)。11 董光昌，張但，劉志誠(chéng)，馬利莊，Curve fairingJ, Progress in Natural Science, 1997(05)。12 劉志誠(chéng)，董光昌，樣條曲線光順概念及指標(biāo)回彈法的應(yīng)用J, 數(shù)學(xué)的實(shí)踐與 認(rèn)識(shí)，1985(03)。13 董光昌，梁友棟，何援軍，樣條曲線擬合與雙圓弧逼近J，應(yīng)用數(shù)學(xué)學(xué)報(bào)，1978(04)。附錄（相關(guān)代碼）主函數(shù)：function nx,ny = l1_fairing( x,y,tolerance ) n = length(x);% 設(shè)定多目標(biāo)優(yōu)化的迭代步驟kk = 8;%設(shè)定步長(zhǎng)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 為控制點(diǎn)的位置% t 為節(jié)點(diǎn)位置，即對(duì)應(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 合并為一個(gè)長(zhǎng)向量,矩陣對(duì)應(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;% 多目標(biāo)迭代優(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