尤祖寰,董蘭芳

(中國科學技術(shù)大學 計算機科學與技術(shù)學院,合肥 230027)

布料仿真一直是計算機圖形學與虛擬現(xiàn)實領(lǐng)域研究的熱門課題.布料仿真技術(shù),旨在通過計算機構(gòu)建虛擬的布料模型,模擬逼真的褶皺及折痕,從而獲得高度的真實感與動態(tài)效果.近年來,隨著虛擬現(xiàn)實的日益火熱,人們更加注重真實感與代入感.同時伴隨著硬件性能的進步,各種大型虛擬現(xiàn)實應(yīng)用成為可能.因此布料仿真技術(shù)獲得了學者們的廣泛關(guān)注和進一步研究.

在進行布料仿真的過程中,布料的物理性質(zhì)通過其本構(gòu)模型得以體現(xiàn).本構(gòu)模型旨在盡可能真實地反映現(xiàn)實世界中布料所受應(yīng)力與應(yīng)變之間的關(guān)系.通常,將布料視為彈性材料,雖然能夠較好地仿真布料的形變特性,但容易造成過度拉伸或壓縮,即“超彈性”現(xiàn)象.應(yīng)變限制算法一般指對布料質(zhì)點應(yīng)變情況進行一定限制,在增強真實感的同時防止其發(fā)生“超彈性”現(xiàn)象.

另外,在布料仿真的一些實際應(yīng)用中,例如虛擬試衣場合,用戶可能對服裝進行局部拉伸操作,如果不加以限制,服裝有可能因所受應(yīng)力過大而發(fā)生非真實的過度形變.因此,需要對布料合理地施加應(yīng)變限制約束,增強布料仿真時的非線性塑性特征,防止“超彈性”現(xiàn)象的發(fā)生.

1 仿真流程概述

布料的基本仿真流程如圖1所示.

1)讀入模型:讀入布料模型與障礙物模型,并對讀取后發(fā)生穿透的部分進行簡單的分離處理;

2)添加物理屬性:根據(jù)布料的本構(gòu)模型及相應(yīng)力學性質(zhì)為布料模型添加物理屬性;

3)模型運動:根據(jù)運動時各節(jié)點受力情況計算各個節(jié)點的速度與位置變化,并更新其狀態(tài);

4)碰撞檢測與響應(yīng):根據(jù)模型包圍盒的距離、位置等信息來計算是否發(fā)生碰撞,并進行響應(yīng)處理.

圖1 布料仿真流程

一般情況下,只要實現(xiàn)了布料的基本仿真流程,即使不做任何優(yōu)化,也可以得到具有一定真實感的仿真布料,其仿真精度與讀入布料的面片數(shù)量、所使用的物理模型、碰撞檢測精度與響應(yīng)方法均密切相關(guān).

除此之外,在上述所有條件都已經(jīng)確定了的情況下,還有一些可進一步豐富布料細節(jié),增強布料真實感的仿真步驟,如網(wǎng)格重劃、應(yīng)變限制、塑性/脆性形變處理等等.

2 相關(guān)工作

2.1 網(wǎng)格重劃

網(wǎng)格的分辨率大小調(diào)節(jié)在布料仿真中相當重要,因為逼真的布料褶皺細節(jié)需要通過高分辨率的網(wǎng)格進行仿真,然而仿真高分辨率的網(wǎng)格又會大大降低計算效率.因此在布料仿真中,網(wǎng)格重劃無疑是兼顧仿真的實時性和真實感的有效方法.通過網(wǎng)格重劃,仿真系統(tǒng)可以自適應(yīng)地改變布料模型的局部分辨率.在仿真過程中,表現(xiàn)真實的褶皺細節(jié)時就需要用小而密的三角網(wǎng)格區(qū)域,而在布料較為平坦之處即可使用低分辨率的大塊三角面片.

網(wǎng)格重劃大大提高了仿真效率,同時又能最大程度地保證了應(yīng)有的布料細節(jié).目前已經(jīng)存在幾種不同的方法,通過使用自適應(yīng)地網(wǎng)格重劃來進行布料仿真.Hutchinson等人是最早提出網(wǎng)格重劃算法的,其通過曲率作為質(zhì)點-彈簧模型的網(wǎng)格重劃標準[1].Villard等人將類似的方法推廣到了常規(guī)四邊形網(wǎng)格的質(zhì)點彈簧模型中[2].二者當時并未考慮布料褶皺的形成.Simnett等人在前者基礎(chǔ)上增加了壓縮和碰撞作為額外細化標準,進一步地提高了細化效率[3].Lee等人則是結(jié)合質(zhì)點彈簧模型與Loop細分的方法對三角網(wǎng)格模型進行網(wǎng)格重劃,通過預(yù)先計算細分步驟,得到多分辨率的層次結(jié)構(gòu),從而自適應(yīng)地降低線性方程組的維數(shù),更好地進行隱式積分求解[4].Narain等人則通過動態(tài)張量場重劃網(wǎng)格,使用邊的分割,翻轉(zhuǎn)以及塌陷這三種基本操作,已經(jīng)能夠較好地展現(xiàn)出布料的褶皺效果[5].

總而言之,網(wǎng)格重劃主要根據(jù)節(jié)點法向、曲率等參數(shù),對碰撞響應(yīng)結(jié)束之后的網(wǎng)格進行重劃分,以滿足展現(xiàn)布料微觀形變的網(wǎng)格精度要求.

2.2 應(yīng)變限制

應(yīng)變限制是對布料應(yīng)變量的約束,旨在防止過度拉伸或壓縮的情況發(fā)生.應(yīng)變限制算法目前可分為直接約束方法與額外約束方法.

早期的研究者們普遍通過使用直接約束方法對應(yīng)變進行限制.Provot最早提出將應(yīng)變限制作為構(gòu)建剛性彈簧的方法,對每個環(huán)節(jié)施加最大最小允許應(yīng)變限制約束[6].Hauth等人通過增大布料剛度,從而盡可能地減少過度拉伸[7],但由于布料剛度的增大,布料數(shù)值求解的精度將無法保證.

綜上所述,當布料受力較小時,采用直接約束處理的方法限制應(yīng)變可以得到較為精確的結(jié)果,而當布料受力較大時,直接約束方法的效率將大幅降低.因此,有學者提出采用額外約束的方法對應(yīng)變進行限制.Thomaszewski等人提出基于連續(xù)介質(zhì)的方法,將應(yīng)變限制應(yīng)用到有限元方法之中,獨立地約束應(yīng)變張量的三個部分[8].Narain等人提出上述方法均不適用于各向異性的高度不均勻網(wǎng)格,通過將應(yīng)變限制視為一個不等式約束非線性規(guī)劃問題,對三角網(wǎng)格進行各向同性的應(yīng)變限制約束[5].Ma等人在Narain等人的工作基礎(chǔ)上,直接使用布料平面坐標與世界坐標的映射關(guān)系提取出經(jīng)緯方向的應(yīng)變,以此近似奇異值分解而得到的兩個主方向的應(yīng)變,得到了各向異性的應(yīng)變限制約束[9].

3 應(yīng)變限制

3.1 應(yīng)變限制定義

假設(shè)布料厚度及其交織模式相比其彈性行為對布料形變的影響要小得多,因此本文將布料視為一個二維連續(xù)體.將布料離散為三角形網(wǎng)格并進行分析,設(shè)三角形頂點i的平面坐標為(ui,vi),世界坐標為xi,如圖2所示.

圖2 三角形網(wǎng)格示意圖

不妨設(shè)w(u,v)為將平面坐標映射到世界坐標系的一個線性函數(shù),這里我們并不關(guān)心w(u,v)的具體形式.考慮其偏導(dǎo)數(shù),其大小決定了在對應(yīng)方向上頂點的拉壓程度.根據(jù)上述假設(shè),給定任意一個三角形元素,關(guān)于其三個頂點i,j,k,定義Δx1=xj-xi,Δx2=xk-xi,Δu1=uj-ui,Δu2=uk-ui,Δv1=vj-vi,Δv2=vk-vi,則可得到Δx1=wuΔu1-wvΔv1,Δx2=wuΔu2-wvΔv2,即:

式(1)可簡寫為:

其中,X=[xixjxk],Δ為 3×2 的矩陣差分算子,β-1為該三角形有限元基礎(chǔ)矩陣.利用矩陣乘法結(jié)合律,我們可以得到:

一般情況下,平面坐標(ui,vi)是不變的,僅在網(wǎng)格重劃之后發(fā)生改變.因此,可以將視為世界坐標xi的函數(shù).

對式(4)進行奇異值分解:

其中,U的列向量即為兩個應(yīng)變主方向,Σ中的兩個對角線元素η1、η2即為其主應(yīng)變率.通過對η1、η2進行限制,即可控制布料模型的應(yīng)變情況.

3.2 約束非線性規(guī)劃[3]

傳統(tǒng)的求解應(yīng)變限制問題的方法是通過雅可比迭代或高斯賽德爾迭代進行求解的,但當方程組較大時收斂慢是不可避免的一大問題.另外針對高度不規(guī)則與各向異性網(wǎng)格時,這種求解方法很難保證其穩(wěn)定性.因此,將應(yīng)變限制問題轉(zhuǎn)化為約束非線性規(guī)劃問題,可以更為有效地處理上述問題,同時提高計算效率.定義最優(yōu)化模型如下:

其中,mn為頂點n的質(zhì)量(即其所在三角形的平均值),為未受約束時所處的位置(即頂點的當前世界坐標).關(guān)于兩個主應(yīng)變率ηj(j=1,2)的約束條件c(x)為:

其中,a為材質(zhì)平面空間下三角形的面積,tj,min為壓縮應(yīng)變閾值,tj,max為拉伸應(yīng)變閾值.又ηj關(guān)于xn的偏導(dǎo)數(shù)為:

則xi方向上的應(yīng)變梯度為:

3.3 基于經(jīng)緯方向的應(yīng)變限制[6]

Ma等人借鑒了如下拉壓能量狀態(tài)函數(shù)[7]

通過將3.2節(jié)中由SVD得到的兩個主應(yīng)變率η1、η2分別用‖wu(x)‖、‖wv(x)‖進行近似,則可以得到關(guān)于經(jīng)緯方向的最優(yōu)化約束條件

又定義剪切向量為ws=wu+wv,類似地,則可以增加兩個關(guān)于剪切向量的約束條件

該方法一方面直接使用經(jīng)緯方向的應(yīng)變率近似SVD得到的主方向應(yīng)變率,從而避免了SVD的開銷,提高了計算效率,另一方面通過定義剪切向量保留了剪切變形,限制了拉伸應(yīng)變,從而使得布料能夠不表現(xiàn)出各向同性的行為.

通過該方法施加的應(yīng)變限制約束,一方面在形變后經(jīng)緯方向無法保證正交,另一方面在每次網(wǎng)格重劃后,三角形網(wǎng)格的經(jīng)緯方向都將發(fā)生變化,如圖3所示.這樣得到的應(yīng)變限制雖然大體上是各向異性的,但是局部約束方向并不一定總能沿著所希望的方向.為此,本文提出一種方法,即通過內(nèi)積進行映射,使得布料每個三角形網(wǎng)格的局部應(yīng)變限制在網(wǎng)格重劃之后仍然能夠和全局應(yīng)變限制保持一致,從而提高了布料仿真的準確性和真實感.

4 內(nèi)積倒數(shù)映射

考慮到布料的針織顯然是具有方向性的,因此經(jīng)向與緯向的應(yīng)變一般并不相同,然而在每一幀對網(wǎng)格進行重劃分之后,重劃分的網(wǎng)格中每個三角形的主方向均發(fā)生了變化.Ma等人的工作雖然實現(xiàn)了各向異性的應(yīng)變限制,但無法較好地處理上述情景,因此針對網(wǎng)格中每個三角面片的主方向不固定的特點,需要定義一種新的應(yīng)變限制來處理這些問題.為此我們提出了內(nèi)積倒數(shù)映射的概念.

為了滿足需求,考慮滿足如下性質(zhì)的一個映射關(guān)系:

1)若三角形面片的主方向與布料全局應(yīng)變限制方向相同,則約束值保持不變;

基于上述性質(zhì),定義映射函數(shù):

其中,ei表示布料的全局應(yīng)變方向,eu表示三角面片內(nèi)的局部應(yīng)變單位方向,εi表示ei方向上的全局應(yīng)變,εu,i表示由ei方向上的全局應(yīng)變εi映射到三角面片eu方向上的局部應(yīng)變.

下面確定求解約束非線性規(guī)劃問題中所需的約束條件c(x)及其關(guān)于xn的梯度向量?xnc,分別就關(guān)于SVD主應(yīng)變方向與關(guān)于布料經(jīng)緯方向的應(yīng)變限制約束進行討論.

4.1 約束條件的確定

首先由式(13),我們可以得到三角面片應(yīng)變閾值為

其中,εmin為壓縮應(yīng)變閾值,εmax為拉伸應(yīng)變閾值.結(jié)合式(7)和式(14)可得

4.2 梯度向量的確定

要求?xnc,則需要先確定?xnηj與?xneu.由式(9)可得:

為確定?xneu,首先需要確定eu與ei的關(guān)系,基于式(5),我們可以得到如下關(guān)系式:

其中,V的列向量即是全局應(yīng)變方向到局部應(yīng)變主方向的旋轉(zhuǎn)矩陣.

血栓形成可以導(dǎo)致嚴重心腦血管方面疾病，其中肺栓塞、大面積腦梗死可以直接致死和致殘，此類患者往往預(yù)后很差。靜脈血栓方面疾病具有一定的發(fā)病率，在住院患者中尤其危重病患者亦有較高的發(fā)病率，可增加死亡率[1] 。深靜脈血栓形成（deep venous thromboembolism，DVT）是最常見的靜脈血栓疾病，在大手術(shù)后及危重病患者人群中，其發(fā)病率較高，往往可能導(dǎo)致較嚴重的并發(fā)癥。有研究顯示，可導(dǎo)致肺動脈高壓[2] 等并發(fā)癥。甚至出現(xiàn)意外死亡和嚴重傷殘，這對該病的治療及預(yù)防提出了巨大挑戰(zhàn)。下面詳細闡述國內(nèi)外研究進展情況。

因此,從ei到eu的旋轉(zhuǎn)映射為:

借鑒Papadopoulo等人[10,11]的工作,由SVD得到的矩陣V對原矩陣元素wij的偏導(dǎo)數(shù)為:

由鏈式法則,可以得到:

5 改進乘子法

經(jīng)過上一節(jié)的詳細推導(dǎo),我們得到了本文提出的自適應(yīng)各向異性應(yīng)變限制算法的目標函數(shù)、約束條件及其關(guān)于自變量的梯度向量.下面則需要對算法得到的約束非線性規(guī)劃問題進行求解.

求解約束非線性規(guī)劃問題,首先想到的是罰函數(shù)法和拉格朗日乘子法.考慮到二種求解算法各有各的優(yōu)點,因此本文將二者相結(jié)合,提出通過改進乘子法來處理該最優(yōu)化規(guī)劃問題.結(jié)合后的改進方法相當于在目標函數(shù)中同時引入了二次懲罰項與線性逼近項.迭代過程中仍是從小到大地依次增大懲罰因子序列,這樣在算法過程早期當懲罰因子很小時,拉格朗日乘子的更新可以很大,即通過引入二次懲罰項使得算法的收斂速度更快,同時又因為拉格朗日乘子的加速更新,一般可以在懲罰因子趨于無窮之前就收斂到最優(yōu)解,從而避免了之前提到的目標函數(shù)海塞矩陣隨懲罰因子增大而失效的情況.

下面具體介紹求解方法,結(jié)合式(6)得到如下約束非線性規(guī)劃問題:

記可行域D={x|x∈Rn,hi(x)≤0,i=1,2,...,m}.若存在不等式約束gj(x)≤0,可加入松弛變量將其轉(zhuǎn)換為等式約束,即:

下面將罰函數(shù)法和拉格朗日乘子法的思想分別引入,即通過拉格朗日乘子和懲罰函數(shù)的組合,將約束條件加入到目標函數(shù)當中,構(gòu)造函數(shù):

其中,拉格朗日乘子λ∈Rm,懲罰因子M∈R.整個過程分兩步:

1)對于固定的λ和M,通過最小化LA更新x.

2)使用前文提及的方法更新拉格朗日乘子:

與懲罰因子M.

對于足夠大的M,迭代一定會收斂,且收斂時λ→λ?,x→x?.

6 實驗結(jié)果

本文的實驗仿真是在Intel(R)Core(TM)i3-2120 CPU @ 3.30 GHz,4 GB內(nèi)存的機器下進行的,采用C++語言進行編寫,并通過OpenGL對所有模型進行渲染與顯示.

為了驗證算法,本文設(shè)計了兩個對比實驗,即分別模擬兩端懸掛布料以及仿真服裝,通過觀察每組實驗對象的形變情況,進行算法比較與分析.

6.1 兩端懸掛布料實驗

實驗分為三組,模擬兩端懸掛布料在穩(wěn)定狀態(tài)下的形變情況,分別實現(xiàn)了基于SVD主應(yīng)變方向的應(yīng)變限制、基于布料經(jīng)緯方向的應(yīng)變限制以及本文提出的基于內(nèi)積倒數(shù)映射的應(yīng)變限制三種方法,具體實驗構(gòu)成參數(shù)如表1所示.其中,圖3(a)(b)(c)垂直方向閾值為[0.99,1.01],水平方向閾值為[0.95,1.05],顯然此時垂直方向要比水平方向具有更高的剛性;圖3(d)、圖3(e)、圖3(f)則與圖3(a)、圖3(b)、圖3(c)相反,垂直方向閾值為[0.95,1.05],水平方向閾值為[0.99,1.01],此時水平方向較垂直方向更具有剛性.

圖3(a)、圖3(d)施加的是基于SVD主應(yīng)變方向的各向同性應(yīng)變限制,即使規(guī)定的兩個主應(yīng)變方向下的限制不同,展現(xiàn)出來的整體布料形態(tài)仍十分相似,僅因為閾值不同而使得細微褶皺上有所區(qū)別,即閾值越小,布料應(yīng)變越小,其褶皺也越少.

表1 布料實驗參數(shù)

圖3 布料對比實驗結(jié)果圖

圖3(b)、圖3(e)展示的是基于經(jīng)緯方向的應(yīng)變限制,因為是各向異性的應(yīng)變限制,所以當經(jīng)向與緯向的應(yīng)變限制約束閾值不同時,影響的不僅僅是細節(jié)的褶皺部分,其展現(xiàn)出來的整體布料形態(tài)也不相同,即布料應(yīng)變方向上的形變量受到其閾值的限制.

圖3(c)、圖3(f)施加的是本文所定義的應(yīng)變限制,與圖3(b)、圖3(e)展示的整體布料形態(tài)比較相似,因為二者均受到各向異性的應(yīng)變限制.但仔細觀察兩種方法的實驗結(jié)果可以發(fā)現(xiàn),在相同限制閾值的情況下,圖3c、圖3(f)相對于圖3(b)、圖3(e)而言展現(xiàn)了更加豐富的褶皺,以及更大的形變狀態(tài).究其原因在于本文所定義的應(yīng)變限制考慮到整體閾值到局部的映射,因此圖3(c)、圖3(f)在垂直方向與水平方向上的效果差異要比前面兩種定義更加明顯.

6.2 仿真服裝實驗

將應(yīng)變限制運用到服裝仿真中,通過三組實驗進一步地驗證了本文方法較基于經(jīng)緯方向的方法限制效果更加準確,實驗參數(shù)如表2所示.如圖4展示了通過三種不同的應(yīng)變限制約束方法下仿真人偶在靜止狀態(tài)下其身著T恤的形變情況.其中,前兩組水平、垂直方向閾值均為±5%,不同之處在于第一組施加的是基于布料經(jīng)緯方向的應(yīng)變限制,第二組施加的是本文提出的基于內(nèi)積倒數(shù)映射的應(yīng)變限制;為進一步地控制實驗參數(shù)變量,第三組施加基于布料經(jīng)緯方向的應(yīng)變限制,并將其水平、垂直方向閾值設(shè)置為±3%.

表2 服裝實驗參數(shù)

顯然地,細微褶皺應(yīng)變的變化程度較為劇烈,而平坦或較大范圍褶皺區(qū)域的應(yīng)變變化程度則較緩慢.因此在限制了T恤的應(yīng)變后,最終應(yīng)該得到較為平整且僅保留較大范圍褶皺的結(jié)果.可以看到,第一組使用本文方法,在應(yīng)變限制閾值為±5%的情況下得到了預(yù)期結(jié)果,即細微褶皺被如期限制,如圖4(a)所示.

而第二組使用基于經(jīng)緯方向的應(yīng)變限制方法,在相同的應(yīng)變限制閾值(±5%)下并未取得預(yù)期結(jié)果.如圖4(b)所示,其T恤的胸口、領(lǐng)口處仍存在理應(yīng)受到限制的細微褶皺.

第三組表明若在施加基于經(jīng)緯方向的應(yīng)變限制的情況下希望得到上述效果,則需要將應(yīng)變限制閾值進一步縮小.如圖4(c)所示,在閾值限定為±3%的情況下,基于經(jīng)緯方向的應(yīng)變限制方法才能較好地保證網(wǎng)格中每個三角形的應(yīng)變都得到了如期的限制.

圖4 服裝對比實驗結(jié)果圖

綜上所述,同為自適應(yīng)各向異性應(yīng)變限制方法,本文所定義的內(nèi)積倒數(shù)映射充分考慮到了整體閾值到局部閾值的映射,使得布料在預(yù)期方向上的應(yīng)變限制效果差異要比基于經(jīng)緯方向的方法更加準確.

7 結(jié)語

本文在基于網(wǎng)格重劃的布料仿真流程中,對布料模型的應(yīng)變限制添加了一種新的約束,解決了網(wǎng)格重劃后各個三角形形變主軸的局部應(yīng)變限制與宏觀上規(guī)定的全局應(yīng)變限制不一致的情況,并使用改進乘子法對約束非線性規(guī)劃問題進行求解.同時,利用OpenGL圖形庫進行顯示與交互,成功地完成了該布料仿真實驗,取得了較其他工作更為科學準確的仿真結(jié)果.

在本文的實驗案例中,內(nèi)積倒數(shù)映射的作用還沒有得到完全地展現(xiàn).在之后的工作中,可以針對具體的布料材質(zhì)可以由力學實驗測得其相應(yīng)的極限應(yīng)變,例如布料橫紋(即布料經(jīng)向)、直紋(布料緯向)與斜紋(布料經(jīng)緯45度方向)的拉伸剛度就具有很大的差異,并通過材料力學的相關(guān)理論確定其不同方向的應(yīng)變限制閾值,從而取代人為設(shè)定的閾值.再結(jié)合本文所提出的基于內(nèi)積倒數(shù)映射的應(yīng)變限制,可以精確地仿真出各向異性的布料材質(zhì),從而進一步地提高布料仿真實驗的準確性與真實感.

1Hutchinson D,Preston M,Hewitt T.Adaptive refinement for mass/spring simulations.Boulic R,Hégron G.Computer Animation and Simulation’96.Vienna,Austria:Springer,1996.31-45.

2Villard J,Borouchaki H.Adaptive meshing for cloth animation.Engineering with Computers,2005,20(4):333-341.[doi:10.1007/s00366-005-0302-1]

3Simnett TJR,Laycock SD,Day AM.An edge-based approach to adaptively refining a mesh for cloth deformation.Proceedings of 2009 EG UK Theory and Practice of Computer Graphics.Cardiff University,United Kingdom.2009.77-84.

4Lee Y,Yoon SE,Oh S,et al.Multi-resolution cloth simulation.Computer Graphics Forum,2010,29(7):2225-2232.[doi:10.1111/cgf.2010.29.issue-7]

5Narain R,Samii A,O’Brien JF.Adaptive anisotropic remeshing for cloth simulation.ACM Transactions on Graphics (TOG),2012,31(6):152.

6Provot X.Deformation constraints in a mass-spring model to describe rigid cloth behaviour.Graphics interface.Canadian Information Processing Society.Toronto,Canada.1995.147-147.

7Hauth M,Etzmu? O,Stra?er W.Analysis of numerical methods for the simulation of deformable models.The Visual Computer,2003,19(7-8):581-600.[doi:10.1007/s00371-003-0206-2]

8Thomaszewski B,Pabst S,Stra?er W.Continuum-based strain limiting.Computer Graphics Forum,2009,28(2):569-576.[doi:10.1111/cgf.2009.28.issue-2]

9Ma GH,Ye JT,Li JT,et al.Anisotropic strain limiting for quadrilateral and triangular cloth meshes.Computer Graphics Forum,2015,35(1):89-99.

10Baraff D,Witkin A.Large steps in cloth simulation.Proceedings of the 25th Annual Conference on Computer Graphics and Interactive Techniques.New York,NY,USA.1998.43-54.

11Papadopoulo T,Lourakis MIA.Estimating the Jacobian of the singular value decomposition:Theory and applications.Proceedings of European Conference on Computer Vision.Berlin Heidelberg,Germany.2000.554-570.