擴散張量成像(DTI)是近年來快速發展的一種磁共振成像技術, 而擴散張量插值是DTI圖像處理過程中一個非常重要的步驟。傳統的譜四元數插值方法對插值點張量的方向進行了修正, 能夠保持張量的各向異性, 但沒有對張量的大小進行修正。本文在傳統譜四元數插值的基礎上提出了一種改進的譜四元數插值方法, 該方法首先對擴散張量進行譜分解, 并利用四元數表示張量的方向, 然后針對不同情況, 分別對張量的大小和方向進行修正, 最后求加權平均值得到插值點張量。通過仿真數據和真實數據將改進方法與譜四元數方法和對數-歐幾里德方法進行比較, 結果表明改進后的方法既能保持張量的部分各向異性指標(FA)和行列式的單調性, 又能保持張量的各向異性, 為擴散張量圖像處理提供了一種重要的插值方法。
引用本文： 徐永紅, 高上策, 郝小飛. 一種改進的擴散張量成像譜四元數插值方法. 生物醫學工程學雜志, 2016, 33(2): 362-367. doi: 10.7507/1001-5515.20160061 復制
引言
擴散張量磁共振成像(diffusion tensor magnetic resonance imaging, DTI)是近年來快速發展的一種磁共振成像技術。它能反映活體組織內水分子各方向的擴散運動，在神經生理、神經外科及腦部腫瘤研究中發揮了重要的作用。張量是圖像處理和分析的一種常用工具，張量本身也可以作為張量值圖像的像素值。目前張量值圖像在信息技術領域研究相對較多的就是DTI[1-2]。
由于現有技術限制，采集到的DTI圖像分辨率較低，為了得到高分辨率的數據，需要對擴散張量數據進行插值。張量值圖像插值是在給定的張量空間約束范圍內，填補和恢復低分辨率圖像中的張量值信息。近年來，一些研究人員致力于張量值圖像插值方法的研究，但目前這些方法都具有一定的局限性[3]。Zhukov等[4]提出張量圖像插值可以采用類似灰度圖像的雙線性插值方法，這種算法將標量圖像的雙線性插值方法直接使用到歐幾里德空間的矩陣中。Rodriguez-florido等[5]提出了局部結構張量的計算方法。為了克服歐幾里德空間中張量處理的缺陷，Pennec等[6]給張量空間賦予了一種仿射不變性度量，使張量空間變為黎曼流形，這樣可以使用黎曼空間中的方法方便地計算兩個張量間的測地線、張量的平均值，但計算較為復雜。Arsigny等[7]提出了一種計算更簡便的基于對數-歐幾里德插度量的張量值圖像插值方法，計算效率高，但不能保持張量各向異性。Fillard[8]提出了黎曼度量框架下張量值圖像線性插值方法，并推廣到多線性插值，該方法計算簡單，物理意義明確，但插值結果不能保持邊緣特性，會造成邊緣模糊。
Collard等[9]提出了一種基于黎曼度量的DTI各向異性保持方法即譜四元數插值方法，通過與對數-歐幾里德方法的比較，顯示該方法保持了張量值圖像的各向異性，降低了計算難度，提高了計算效率。但是，該方法只對插值點張量的方向進行修正，而沒有對張量的大小進行修正。
本文在文獻[9]方法的基礎上進行了改進。改進后的方法針對不同情況，分別利用簡化的部分各向異性(fractional anisotropy，FA)和相對各向異性對張量的大小和方向進行修正，然后求加權平均得到插值點張量。利用仿真數據和真實犬心DTI數據，將改進后的譜四元數(improved spectral quaternion, ISQ)方法與譜四元數(spectral quaternion, SQ)插值方法和對數-歐幾里德(Log-Euclidean, LE)插值方法進行比較，最后得出結論。
1 理論基礎
擴散張量成像是核磁共振成像的特殊形式，它的三維數據空間中每個體素為一個對稱正定的2階張量。擴散是一個三維過程，分子沿空間某一方向擴散的距離相等或不相等，可以將擴散的方式分為兩種：各向同性和各向異性。擴散張量將擴散運動表示為一個3×3的對稱正定矩陣[10]：
| $ D=\left\{ \begin{matrix} {{D}_{xx}}\ \ {{D}_{xy}}\ \ {{D}_{xz}} \\ {{D}_{xy}}\ \ {{D}_{yy}}\ \ {{D}_{yz}} \\ {{D}_{xz}}\ \ {{D}_{yz}}\ \ {{D}_{zz}} \\ \end{matrix} \right\} $ |
式中，矩陣的不同元素表達了不同方向上的擴散率，由擴散張量的矩陣形式分解可以得到三個實特征值，將其定義為λ1、λ2和λ3，每個特征值對應一個特征向量、和。
對稱正定矩陣譜分解為：
| $ S=U\Lambda {{U}^{T}} $ |
其中，Λ=diag(λ1, λ2, λ3)，并且，λ1≥λ2≥λ3 > 0，U是由特征向量組成的正交矩陣[11]。
譜四元數插值方法是在譜幾何的基礎上提出的，利用四元數來表示張量的方向，這樣降低了計算難度，提高了計算效率。具體算法見參考文獻[9]。譜四元數插值方法能夠很好地保持插值張量的各向異性, 但部分各向異性指標不穩定。當被插值張量FA指標差異較小時，插值點張量的FA值降低，這在臨床應用上很容易造成疾病誤診[12]。因此，本文在譜四元數插值的基礎上進行改進，即對方向進行修正的同時，針對被插值張量的不同情況，對張量譜分解后的特征值進行修正，然后求加權平均值得到最終結果。改進的譜四元數插值方法算法步驟如下：
(1)?給定兩個張量S1=U1Λ1U1T，S2=U2Λ2U2T。其中，Λi=diag(λi1, λi2, λi3)，并且，λi1≥λi2≥λi3 > 0，i=1, 2，U1、U2分別為由特征向量組成的正交矩陣。
(2)?分別用四元數q1、q2表示旋轉矩陣U1、U2。
(3)?計算插值點張量的特征值矩陣Λ(t)。
(4)?選擇夾角最小的兩個四元數進行插值，因為-q與q表示同一個旋轉矩陣。
(5)?計算插值點張量的方向q(t)。
(6)?將四元數q(t)轉換為旋轉矩陣U(t)。
(7)?插值點張量為：
| $ S\left(t \right)=U\left(t \right)\Lambda \left(t \right)U{{\left(t \right)}^{T}} $ |
其中，步驟(3)具體算法如下：
首先，引入過渡函數：
| $ f\left(x \right)=\frac{{{\left(\beta x \right)}^{4}}}{1+{{\left(\beta x \right)}^{4}}} $ |
其中，β為可調參數值。
①?當被插值張量的各項異性差值小于等于0.2時，保持張量的行列式更有利于插值結果。本文引入簡化的部分各項異性指標DA：
| $ DA=\frac{{{\left(\sum\limits_{i=1}^{3}{{{\lambda }_{i}}} \right)}^{2}}}{\sum\limits_{i=1}^{3}{\lambda _{i}^{2}}} $ |
插值權重分別為：
| $ w_{1}^{*}=\left(1-t \right)\frac{\gamma \left(D{{A}_{1}}, D{{A}_{t}} \right)}{\overline{\gamma }} $ |
| $ w_{2}^{*}=t\frac{\gamma \left(D{{A}_{t}}, D{{A}_{2}} \right)}{\overline{\gamma }} $ |
其中，γ(DA1, DA2)=f(min(DA1, DA2))，
此時，特征值的插值公式為：
| $ \Lambda \left(t \right)=\exp \left(w_{1}^{*}\log \left({{\Lambda }_{1}} \right)+w_{2}^{*}\log \left({{\Lambda }_{2}} \right) \right) $ |
②?當被插值張量的各項異性差值大于0.2時，保持張量的輪廓更有利于插值結果。
| $ {{D}_{1}}=\det \left({{S}_{1}} \right) $ |
| $ {{D}_{2}}=\det \left({{S}_{2}} \right) $ |
| $ D\left(t \right)=\left(1-t \right){{D}_{1}}+t{{D}_{2}} $ |
| $ f\left(t \right)=\frac{\log \left(D\left(t \right) \right)-\log \left({{D}_{1}} \right)}{\log \left({{D}_{2}} \right)-\log \left({{D}_{1}} \right)} $ |
此時, 特征值插值公式為：
| $ \Lambda \left(t \right)=\exp \left(\left(1-f\left(t \right) \right)\log \left({{\Lambda }_{1}} \right)+f\left(t \right)\log \left({{\Lambda }_{2}} \right) \right) $ |
其中，步驟(5)具體算法如下：
本文使用相對各向異性RA來修正插值點張量的方向，定義如下:
| $ RA=\frac{\sqrt{{{\left({{\lambda }_{1}}-\overline{\lambda } \right)}^{2}}+{{\left({{\lambda }_{2}}-\overline{\lambda } \right)}^{2}}+{{\left({{\lambda }_{3}}-\overline{\lambda } \right)}^{2}}}}{\sqrt{3\overline{\lambda }}} $ |
其中，，本文給出關于RA的插值公式：
| $ RA=\left(S\left(t; {{S}_{1}}, {{S}_{2}} \right) \right)=\left(1-t \right)RA\left({{S}_{1}} \right)+tRA\left({{S}_{2}} \right) $ |
下面，用RAt表示RA(S(t; S1, S2))。
運用四元數表示特征向量組成的旋轉矩陣，四元數線性插值為：
| $ q\left(t \right)=w_{3}^{*}\left(t \right){{q}_{1}}+w_{4}^{*}\left(t \right){{q}_{2}} $ |
其中，，。
并且，α(RA1, RA2)=f(min(RA1, RA2))，
2 實驗過程與分析
本文利用仿真數據和真實數據對該方法進行驗證，并且將ISQ方法與SQ方法和LE方法進行比較。為了充分考慮被插值張量的情況，本文構造的仿真數據包括各向異性和各向同性的大小不同的張量。首先確定張量的旋轉矩陣I，通過選擇不同的角度φ來確定方向，然后選擇不同的特征值來確定張量的大小。
針對真實數據，本文應用MedINRIA軟件提供的MATLAB函數工具包TensorInrIO，將原始擴散張量圖像數據經過MedINRIA軟件預處理后，保存為張量數據類型。通過TensorInrIO工具包導入到MATLAB中，在MATLAB2010下進行插值的研究，將研究結果通過TensorInrIO工具包保存為張量數據類型，應用MedINRIA軟件中的Tensor Viewer模塊進行張量的顯示。
首先，對兩個仿真張量進行插值。這兩個張量分別為S1和S2，S1的特征值為[5.3?2.5?0.2], S2的特征值為[6.6?2.6?1.1]，旋轉矩陣為
| $ I=\left[ \begin{matrix} \cos \varphi \ \-\sin \varphi \ \ 0 \\ \sin \varphi \ \ \cos \varphi \ \ 0 \\ 0\ \ \ \ \ \ \ \ \ 0\ \ \ \ \ \ \ \ 1 \\ \end{matrix} \right] $ |
其中，，結果如圖 1所示。
圖1
三種方法參數對比
Figure1.
Comparison of parameters among three methods
圖 1為三種方法插值結果的參數對比圖。由圖可知，LE方法的部分各向異性和相對各向異性分數降低，而ISQ的這兩個參數成線性關系；SQ方法的主特征向量角度也成線性，但其行列式有所降低，因此，ISQ方法比其他兩種方法更具優越性。
如圖 2所示為4個張量間的插值結果。四個已知張量分別為S1、S2、S3、S4，其中，，S2、S3、S4的旋轉矩陣均為，并且，。根據圖 2中張量的顏色變化可以看出，右半部分的張量有明顯差異，ISQ方法在插值過程中漸變趨勢更明顯。如圖 3所示為張量的部分各向異性等值曲線，由曲線可以看出，SQ方法的等值線間距逐漸變大，而ISQ方法中等值線間距更小，更均勻。由此可以得出，ISQ方法比SQ方法更能保持張量的各向異性。
圖2
SQ和ISQ方法在4個張量間的插值
其中顏色變化表示張量主特征向量的方向變化
Figure2. SQ and ISQ interpolation between four tensorsthe change of colors represented the main feature vector directions of tensors
圖3
SQ和ISQ方法部分各向異性等值線
Figure3.
Contour map of FA of SQ and ISQ method
由于技術有限，數據采集成本高，沒有條件采集實際的醫學診斷圖像，所以，本文利用權威數據庫的正常犬心臟DTI數據對本文方法進行驗證(http://cvrgrid.org/data/ex-vivo)。選取心臟第66個切片的部分數據，原始數據大小為256×256，為了能夠看清楚張量插值的細節變化，對原始數據進行減采樣，減采樣后數據大小為64×64，如
圖4
原始犬心DTI數據
Figure4.
Original DTI data of canine heart
選取圖 4中紅色矩形區域中的數據進行插值，數據大小為20×20，該區域的原始數據圖像如圖 5所示。分別利用LE、SQ和ISQ方法對該矩形區域中的數據進行插值，結果如圖 5所示。
圖5
三種方法插值結果對比圖
其中顏色變化表示張量主特征向量的方向變化
Figure5. Reults of the three interpolation methodsthe change of colors represented the main feature vector directions of tensors
其中，LE插值法在一些插值點比SQ方法更平滑，而SQ方法在一些插值點具有更好的輪廓，細節上更加清晰，但有些點的方向變化過于明顯。而ISQ插值方法介于兩者之間，既保證了輪廓的清晰，又能使局部的方向變化緩和，符合生物組織中水分子的擴散狀態。
為了更好地比較這三種算法，分別將插值后張量數據與原始數據進行均方誤差分析。主要比較張量的部分各向異性、平均擴散率(即張量特征值的均值)以及行列式的誤差，結果如表 1所示，可以看到SQ方法略優于LE方法，而ISQ插值法三個參數的均方誤差都是最小的，由此也可得出，ISQ方法更具優越性。
表1
三種插值方法的誤差分析
Table1.
Error analysis of three different interpolation methods
4 結論
本文提出了一種改進的譜四元數插值方法，該方法針對被插值張量的不同情況，對張量的大小進行了修正，同時對張量的方向進行修正，然后通過求加權平均值得到被插值點張量。最后，分別利用仿真數據和權威數據庫的真實犬心DTI數據將ISQ、SQ和LE方法進行比較，結果表明SQ插值法明顯優于LE插值法，而ISQ插值法優于SQ插值法。改進后的方法既能保持部分各向異性指標以及行列式的單調性，又能同時保持張量的各向異性，為擴散張量圖像處理提供了一種重要的插值方法。
文獻[9]中算法主要降低了擴散張量插值的計算復雜度，并且對插值方向算法進行了優化。本文主要是在插值大小上進行改進，同時對插值方向做了進一步修正。因此，本文的計算復雜度比文獻[9]稍高一些，但是提高了擴散張量插值的準確度。
引言
擴散張量磁共振成像(diffusion tensor magnetic resonance imaging, DTI)是近年來快速發展的一種磁共振成像技術。它能反映活體組織內水分子各方向的擴散運動，在神經生理、神經外科及腦部腫瘤研究中發揮了重要的作用。張量是圖像處理和分析的一種常用工具，張量本身也可以作為張量值圖像的像素值。目前張量值圖像在信息技術領域研究相對較多的就是DTI[1-2]。
由于現有技術限制，采集到的DTI圖像分辨率較低，為了得到高分辨率的數據，需要對擴散張量數據進行插值。張量值圖像插值是在給定的張量空間約束范圍內，填補和恢復低分辨率圖像中的張量值信息。近年來，一些研究人員致力于張量值圖像插值方法的研究，但目前這些方法都具有一定的局限性[3]。Zhukov等[4]提出張量圖像插值可以采用類似灰度圖像的雙線性插值方法，這種算法將標量圖像的雙線性插值方法直接使用到歐幾里德空間的矩陣中。Rodriguez-florido等[5]提出了局部結構張量的計算方法。為了克服歐幾里德空間中張量處理的缺陷，Pennec等[6]給張量空間賦予了一種仿射不變性度量，使張量空間變為黎曼流形，這樣可以使用黎曼空間中的方法方便地計算兩個張量間的測地線、張量的平均值，但計算較為復雜。Arsigny等[7]提出了一種計算更簡便的基于對數-歐幾里德插度量的張量值圖像插值方法，計算效率高，但不能保持張量各向異性。Fillard[8]提出了黎曼度量框架下張量值圖像線性插值方法，并推廣到多線性插值，該方法計算簡單，物理意義明確，但插值結果不能保持邊緣特性，會造成邊緣模糊。
Collard等[9]提出了一種基于黎曼度量的DTI各向異性保持方法即譜四元數插值方法，通過與對數-歐幾里德方法的比較，顯示該方法保持了張量值圖像的各向異性，降低了計算難度，提高了計算效率。但是，該方法只對插值點張量的方向進行修正，而沒有對張量的大小進行修正。
本文在文獻[9]方法的基礎上進行了改進。改進后的方法針對不同情況，分別利用簡化的部分各向異性(fractional anisotropy，FA)和相對各向異性對張量的大小和方向進行修正，然后求加權平均得到插值點張量。利用仿真數據和真實犬心DTI數據，將改進后的譜四元數(improved spectral quaternion, ISQ)方法與譜四元數(spectral quaternion, SQ)插值方法和對數-歐幾里德(Log-Euclidean, LE)插值方法進行比較，最后得出結論。
1 理論基礎
擴散張量成像是核磁共振成像的特殊形式，它的三維數據空間中每個體素為一個對稱正定的2階張量。擴散是一個三維過程，分子沿空間某一方向擴散的距離相等或不相等，可以將擴散的方式分為兩種：各向同性和各向異性。擴散張量將擴散運動表示為一個3×3的對稱正定矩陣[10]：
| $ D=\left\{ \begin{matrix} {{D}_{xx}}\ \ {{D}_{xy}}\ \ {{D}_{xz}} \\ {{D}_{xy}}\ \ {{D}_{yy}}\ \ {{D}_{yz}} \\ {{D}_{xz}}\ \ {{D}_{yz}}\ \ {{D}_{zz}} \\ \end{matrix} \right\} $ |
式中，矩陣的不同元素表達了不同方向上的擴散率，由擴散張量的矩陣形式分解可以得到三個實特征值，將其定義為λ1、λ2和λ3，每個特征值對應一個特征向量、和。
對稱正定矩陣譜分解為：
| $ S=U\Lambda {{U}^{T}} $ |
其中，Λ=diag(λ1, λ2, λ3)，并且，λ1≥λ2≥λ3 > 0，U是由特征向量組成的正交矩陣[11]。
譜四元數插值方法是在譜幾何的基礎上提出的，利用四元數來表示張量的方向，這樣降低了計算難度，提高了計算效率。具體算法見參考文獻[9]。譜四元數插值方法能夠很好地保持插值張量的各向異性, 但部分各向異性指標不穩定。當被插值張量FA指標差異較小時，插值點張量的FA值降低，這在臨床應用上很容易造成疾病誤診[12]。因此，本文在譜四元數插值的基礎上進行改進，即對方向進行修正的同時，針對被插值張量的不同情況，對張量譜分解后的特征值進行修正，然后求加權平均值得到最終結果。改進的譜四元數插值方法算法步驟如下：
(1)?給定兩個張量S1=U1Λ1U1T，S2=U2Λ2U2T。其中，Λi=diag(λi1, λi2, λi3)，并且，λi1≥λi2≥λi3 > 0，i=1, 2，U1、U2分別為由特征向量組成的正交矩陣。
(2)?分別用四元數q1、q2表示旋轉矩陣U1、U2。
(3)?計算插值點張量的特征值矩陣Λ(t)。
(4)?選擇夾角最小的兩個四元數進行插值，因為-q與q表示同一個旋轉矩陣。
(5)?計算插值點張量的方向q(t)。
(6)?將四元數q(t)轉換為旋轉矩陣U(t)。
(7)?插值點張量為：
| $ S\left(t \right)=U\left(t \right)\Lambda \left(t \right)U{{\left(t \right)}^{T}} $ |
其中，步驟(3)具體算法如下：
首先，引入過渡函數：
| $ f\left(x \right)=\frac{{{\left(\beta x \right)}^{4}}}{1+{{\left(\beta x \right)}^{4}}} $ |
其中，β為可調參數值。
①?當被插值張量的各項異性差值小于等于0.2時，保持張量的行列式更有利于插值結果。本文引入簡化的部分各項異性指標DA：
| $ DA=\frac{{{\left(\sum\limits_{i=1}^{3}{{{\lambda }_{i}}} \right)}^{2}}}{\sum\limits_{i=1}^{3}{\lambda _{i}^{2}}} $ |
插值權重分別為：
| $ w_{1}^{*}=\left(1-t \right)\frac{\gamma \left(D{{A}_{1}}, D{{A}_{t}} \right)}{\overline{\gamma }} $ |
| $ w_{2}^{*}=t\frac{\gamma \left(D{{A}_{t}}, D{{A}_{2}} \right)}{\overline{\gamma }} $ |
其中，γ(DA1, DA2)=f(min(DA1, DA2))，
此時，特征值的插值公式為：
| $ \Lambda \left(t \right)=\exp \left(w_{1}^{*}\log \left({{\Lambda }_{1}} \right)+w_{2}^{*}\log \left({{\Lambda }_{2}} \right) \right) $ |
②?當被插值張量的各項異性差值大于0.2時，保持張量的輪廓更有利于插值結果。
| $ {{D}_{1}}=\det \left({{S}_{1}} \right) $ |
| $ {{D}_{2}}=\det \left({{S}_{2}} \right) $ |
| $ D\left(t \right)=\left(1-t \right){{D}_{1}}+t{{D}_{2}} $ |
| $ f\left(t \right)=\frac{\log \left(D\left(t \right) \right)-\log \left({{D}_{1}} \right)}{\log \left({{D}_{2}} \right)-\log \left({{D}_{1}} \right)} $ |
此時, 特征值插值公式為：
| $ \Lambda \left(t \right)=\exp \left(\left(1-f\left(t \right) \right)\log \left({{\Lambda }_{1}} \right)+f\left(t \right)\log \left({{\Lambda }_{2}} \right) \right) $ |
其中，步驟(5)具體算法如下：
本文使用相對各向異性RA來修正插值點張量的方向，定義如下:
| $ RA=\frac{\sqrt{{{\left({{\lambda }_{1}}-\overline{\lambda } \right)}^{2}}+{{\left({{\lambda }_{2}}-\overline{\lambda } \right)}^{2}}+{{\left({{\lambda }_{3}}-\overline{\lambda } \right)}^{2}}}}{\sqrt{3\overline{\lambda }}} $ |
其中，，本文給出關于RA的插值公式：
| $ RA=\left(S\left(t; {{S}_{1}}, {{S}_{2}} \right) \right)=\left(1-t \right)RA\left({{S}_{1}} \right)+tRA\left({{S}_{2}} \right) $ |
下面，用RAt表示RA(S(t; S1, S2))。
運用四元數表示特征向量組成的旋轉矩陣，四元數線性插值為：
| $ q\left(t \right)=w_{3}^{*}\left(t \right){{q}_{1}}+w_{4}^{*}\left(t \right){{q}_{2}} $ |
其中，，。
并且，α(RA1, RA2)=f(min(RA1, RA2))，
2 實驗過程與分析
本文利用仿真數據和真實數據對該方法進行驗證，并且將ISQ方法與SQ方法和LE方法進行比較。為了充分考慮被插值張量的情況，本文構造的仿真數據包括各向異性和各向同性的大小不同的張量。首先確定張量的旋轉矩陣I，通過選擇不同的角度φ來確定方向，然后選擇不同的特征值來確定張量的大小。
針對真實數據，本文應用MedINRIA軟件提供的MATLAB函數工具包TensorInrIO，將原始擴散張量圖像數據經過MedINRIA軟件預處理后，保存為張量數據類型。通過TensorInrIO工具包導入到MATLAB中，在MATLAB2010下進行插值的研究，將研究結果通過TensorInrIO工具包保存為張量數據類型，應用MedINRIA軟件中的Tensor Viewer模塊進行張量的顯示。
首先，對兩個仿真張量進行插值。這兩個張量分別為S1和S2，S1的特征值為[5.3?2.5?0.2], S2的特征值為[6.6?2.6?1.1]，旋轉矩陣為
| $ I=\left[ \begin{matrix} \cos \varphi \ \-\sin \varphi \ \ 0 \\ \sin \varphi \ \ \cos \varphi \ \ 0 \\ 0\ \ \ \ \ \ \ \ \ 0\ \ \ \ \ \ \ \ 1 \\ \end{matrix} \right] $ |
其中，，結果如圖 1所示。
圖1
三種方法參數對比
Figure1.
Comparison of parameters among three methods
圖 1為三種方法插值結果的參數對比圖。由圖可知，LE方法的部分各向異性和相對各向異性分數降低，而ISQ的這兩個參數成線性關系；SQ方法的主特征向量角度也成線性，但其行列式有所降低，因此，ISQ方法比其他兩種方法更具優越性。
如圖 2所示為4個張量間的插值結果。四個已知張量分別為S1、S2、S3、S4，其中，，S2、S3、S4的旋轉矩陣均為，并且，。根據圖 2中張量的顏色變化可以看出，右半部分的張量有明顯差異，ISQ方法在插值過程中漸變趨勢更明顯。如圖 3所示為張量的部分各向異性等值曲線，由曲線可以看出，SQ方法的等值線間距逐漸變大，而ISQ方法中等值線間距更小，更均勻。由此可以得出，ISQ方法比SQ方法更能保持張量的各向異性。
圖2
SQ和ISQ方法在4個張量間的插值
其中顏色變化表示張量主特征向量的方向變化
Figure2. SQ and ISQ interpolation between four tensorsthe change of colors represented the main feature vector directions of tensors
圖3
SQ和ISQ方法部分各向異性等值線
Figure3.
Contour map of FA of SQ and ISQ method
由于技術有限，數據采集成本高，沒有條件采集實際的醫學診斷圖像，所以，本文利用權威數據庫的正常犬心臟DTI數據對本文方法進行驗證(http://cvrgrid.org/data/ex-vivo)。選取心臟第66個切片的部分數據，原始數據大小為256×256，為了能夠看清楚張量插值的細節變化，對原始數據進行減采樣，減采樣后數據大小為64×64，如
圖4
原始犬心DTI數據
Figure4.
Original DTI data of canine heart
選取圖 4中紅色矩形區域中的數據進行插值，數據大小為20×20，該區域的原始數據圖像如圖 5所示。分別利用LE、SQ和ISQ方法對該矩形區域中的數據進行插值，結果如圖 5所示。
圖5
三種方法插值結果對比圖
其中顏色變化表示張量主特征向量的方向變化
Figure5. Reults of the three interpolation methodsthe change of colors represented the main feature vector directions of tensors
其中，LE插值法在一些插值點比SQ方法更平滑，而SQ方法在一些插值點具有更好的輪廓，細節上更加清晰，但有些點的方向變化過于明顯。而ISQ插值方法介于兩者之間，既保證了輪廓的清晰，又能使局部的方向變化緩和，符合生物組織中水分子的擴散狀態。
為了更好地比較這三種算法，分別將插值后張量數據與原始數據進行均方誤差分析。主要比較張量的部分各向異性、平均擴散率(即張量特征值的均值)以及行列式的誤差，結果如表 1所示，可以看到SQ方法略優于LE方法，而ISQ插值法三個參數的均方誤差都是最小的，由此也可得出，ISQ方法更具優越性。
表1
三種插值方法的誤差分析
Table1.
Error analysis of three different interpolation methods
4 結論
本文提出了一種改進的譜四元數插值方法，該方法針對被插值張量的不同情況，對張量的大小進行了修正，同時對張量的方向進行修正，然后通過求加權平均值得到被插值點張量。最后，分別利用仿真數據和權威數據庫的真實犬心DTI數據將ISQ、SQ和LE方法進行比較，結果表明SQ插值法明顯優于LE插值法，而ISQ插值法優于SQ插值法。改進后的方法既能保持部分各向異性指標以及行列式的單調性，又能同時保持張量的各向異性，為擴散張量圖像處理提供了一種重要的插值方法。
文獻[9]中算法主要降低了擴散張量插值的計算復雜度，并且對插值方向算法進行了優化。本文主要是在插值大小上進行改進，同時對插值方向做了進一步修正。因此，本文的計算復雜度比文獻[9]稍高一些，但是提高了擴散張量插值的準確度。
