TWI555514B

TWI555514B - 用於實施具遞迴最佳化之高效率電腦斷層掃描之系統、工作站及方法

Info

Publication number: TWI555514B
Application number: TW100125177A
Authority: TW
Inventors: 沃爾夫藍Ｒ杰利奇
Original assignee: 沃爾夫藍Ｒ杰利奇
Priority date: 2010-07-16
Filing date: 2011-07-15
Publication date: 2016-11-01
Also published as: WO2012009535A1; EP2593906A1; CN103052959A; EP2593906B1; TW201204327A; EP2593906A4; US20100278413A1; US8660330B2; CN103052959B

Description

用於實施具遞迴最佳化之高效率電腦斷層掃描之系統、工作站及方法

本發明係相關於影像重建之領域。

斷層掃描(tomography)，即，在一n維空間中的一區域內的密度之計算(估計)，其係依據該區域(以像素來代表；通常0<m<n)的m維投射(projection)，通常落入兩種主要類別之內：濾波反投射(filtered back-projection，FBP)，與疊代重建(iterative reconstruction，IR)，後者係習知為代數重建技術(Algebraic Reconstruction Technique，ART)的一種變化技術。FBP與習知之IR技術皆各有其缺陷，特別是在每一次投射常會包含了由2000個像素(單維)至2000×2000個像素(二維)，其常有多達10⁶-10¹⁰個影像元素(即，立體像素(voxel))需要進行計算的現代組構之中為然。

與FBP相關的主要缺陷包括，其需要有龐大次數的投射才能達到有限的量化精確度。其投射次數通常數以百計，但其投射並未被有效率地加以運用。例如，其計算得出之密度估計值可能出現負值，但一般皆知物理密度值不會是負數。其有限的量化精確度有時可以利用後續的疊代微調加以改進。與FBP相關的其他缺陷包括，其無法依據立體像素之位置而改變資料權重，如Wood,S.L.,Morf,M.,A fast implementation of a minimum variance estimator for computerized tomography image reconstruction,IEEE,Trans. on Biomed. Eng.,Vol. BME-28,No. 2,pp. 56-68,1981(此後稱此文件為Wood文件)中所討論者，以及其無法有效地考量其各種限制。

由於FBP方法通常是依賴快速傅立葉轉換(fast Fourier transformation)，故其強項包括其計算速度，中央切片理論，以及其在適當及預先計算濾波上的可能利用。

FBP的變化技術以及習知的IR技術皆使用矩陣運算，其因此可能適用於其體積或圖像元素(立體像素或二維像素)在數千個程度的數量範圍內的「小型」問題。不過，對於典型的斷層掃描影像重建組構而言，在可見的將來，其計算上的負擔是無法以此等以矩陣為基礎的技術來承擔的。這是由於其電腦運算次數的規模是以N的冪次來計算的。例如，在3-D情況之中，其運算次數規模可以高過N⁵(其中N是為一個要予重建立方體，即所要檢視區域，其一側的像素之數量)，如同Wood文件以及2001年12月8日之美國專利6,332,035所討論者。

IR技術的好處包括其可以顧慮到其限定條件的能力，特別是其可以確保其對密度的估算不會出現負值的能力。非負值密度估算值的確保，可以導致影像對比上的顯著改進，以及較佳的量化結果。

習知IR技術之相關缺陷包括，其需要反覆地解反運算以及前向投射運算才能獲得後續疊代的更正項。其對於許多此種反覆出現的成對運算，通常會使IR技術的執行變慢，並且要依賴特定的，不穩定的執行方式。

由於此些習知IR技術使用了各種的(疊代)最佳化方法，其朝向解答的收斂速度需依賴某特定技術在赫士矩陣(Hessian matrix)上的執行效果如何，而此矩陣是為所要最佳化的目的函數(亦即，相對於影像元素的性能標準)的二次微分矩陣。其最常見的技術是為IR的變化技術。不過，由於影像重建所使用的赫士矩陣之大尺度，其構造(以及其導數矩陣之構造)通常是被忽略或僅只粗略地加以估算。此外，由於赫士矩陣的特徵值(eigen-value)之廣域分佈，目前的最佳化技術在一定次數的疊代(通常為數十至數百次)之後，皆傾向於顯現不出更進一步的進展，其亦只能處理大的特徵值。

此類演算法的多格柵變化可能有所幫助，但由於赫士矩陣牽涉到較細的格柵，故其最終仍會失敗。多格柵解析度在此係指隨著疊代的進行所使用的漸進較細的解析度。

當一次投射包含了由2,000個像素(一維)至65,000 x 65,000個像素(二維)時，現代斷層掃描的組構通常需要計算10⁶-10¹⁴個影像元素(即，立體像素)。近似法(approximation techniques)，諸如線性化(linearization)與次級問題維度的減低，可能有助於使高效率電腦斷層掃描(High Efficiency Computed Tomography，HECT)變為可行。不過，以現代斷層掃描組構所進行的影像重建，由於典型例如非線性估算的緣故，通常是相當複雜的。

本發明可提供一系統，其包含：一或更多個發射器，以發射一激發能量進入被觀察的一物體內；一或更多個偵側器以產生投射空間資料，以將由該一或多個偵測器反應於被發射入被觀察物體內之激發能量而所接收到的一能量加以編碼；一控制器，其控制該一或更多個發射器發射激發能量及該一或更多個偵側器產生該投射空間資料；與一影像重建器，具有至少一處理器以利用以下方式而接收該投射空間資料並處理該投射空間資料：計算表現該投射空間資料與預測投射資料間之差異特徵之一第一量，其中投射空間資料與預測投射資料間之差異係與一投射增益相關聯；於一資料儲存裝置內記錄下將該差異加以編碼之一第一資料；依據來自先前疊代所記錄的第一資料而計算經修改過的第一量；計算對應於至少一脈衝反應之一第二量，每一脈衝反應係對應於物體空間中一位置處具單位強度的一參考立體像素；利用修改過第一量與第二量之一反函數而計算一更新值，其中更新值係與一立體像素增益相關聯；與利用更新值而重建代表被觀察物體之一物體空間影像。

本發明可提供一工作站，其包含：一或更多個處理器；與一或更多個資料儲存裝置以儲存可由一或更多個處理器所執行的軟體，該軟體包含：可接收代表由一物體空間至一投射空間之一物體的輸入資料的軟體碼；可計算表現該投射空間資料與預測投射資料間之差異特徵之一第一量之軟體碼，其中投射空間資料與預測投射資料間之差異係與一投射增益相關聯；計算對應於至少一脈衝反應之一第二量，每一脈衝反應係對應於物體空間中一位置處具單位強度的一參考立體像素；可利用第一量與第二量之反函數而計算一更新值之軟體碼；可於一或更多個資料儲存裝置，或另一或其他更多個資料儲存裝置內記錄下將該更新值加以編碼之一第二資料之軟體碼；依據來自先前疊代所記錄的第二資料而修改該更新值之軟體碼；與可利用被修改過更新值而重建代表被觀察物體之一物體空間影像之軟體碼，其中該被修改過之更新值係與一立體像素增益相關聯。

本發明可提供由一或更多個處理器執行儲存於一或更多個資料儲存裝置中之軟體碼之一方法，該方法包含：由一投射影像獲取系統，該一或多個資料儲存裝置，另一或其他更多個資料儲存裝置，由該一或多個處理器所執行之一轉體模擬，或由另一或其他更多個處理器所執行之軟體模擬之至少其中之一之處接收代表由一物體空間至一投射空間之一物體投射之輸入資料；計算表現該輸入資料在一解析度格柵上與預測投射資料在解析度格柵上其間之差異特徵之一第一量，其中輸入資料與預測投射資料間之差異係與一投射增益相關聯；於一或多個資料儲存裝置，或另一或其他更多個資料儲存裝置內記錄下將該差異加以編碼之一第一資料；依據來自先前疊代所記錄的第一資料而計算經修改過的第一量；計算對應於至少一脈衝反應之一第二量，每一脈衝反應係對應於物體空間中一位置處具單位強度的一參考立體像素；利用第一量或修改過之第一量其中之一，與第二量之反函數而計算一更新值，其中該更新值係與一立體像素增益相關聯；於該一或多個資料儲存裝置，或該另一或其他更多個資料儲存裝置內記錄下將該更新值加以編碼之一第二資料；依據來自先前疊代所記錄的第二資料而修改更新值；利用被修改過更新值而重建代表被觀察物體之一影像，與將該被重建之影像輸出至一輸出裝置，其中該輸出裝置係至少為該一或更多個資料儲存裝置，該另一或其他更多個資料儲存裝置，一顯示裝置，或一列印裝置其中之一；其中接收輸入資料，計算第一量，記錄第一資料，計算修改之第一量，計算第二量，記錄第二量，修改更新值，重建影像，與輸出經重建之影像係由該一或更多個處理器依據儲存於該一或更多個儲存裝置中之軟體碼而執行。

以下將詳細說明本發明說明性質之各實施例。為清楚起見，以下說明性質實施例之討論將採用特定名詞。不過，本發明不應受限於此些特定的名詞用語上。習於本技藝者皆可理解，其他對等的部件亦皆可以適用，且其他方法亦可在不偏離於本發明之廣泛觀念的情況下被發展出來。本說明書所援引之各參考資料亦各皆各自包容於本發明之習知技藝範圍內。

圖1顯示依據本發明某些實施例之系統圖。投射像素資料104包含有可將一物體之內部構造加以編碼的資訊，其可利用一影像重建器105而轉換為可在一輸出裝置107上觀看的一個重建影像106。投射像素資料104可來自於一次實驗性質的擷取101，或在一部電腦上的一次模擬103，或者包含有來自於先前較早實驗或模擬時所記錄下來的投射像素資料104之一資料儲存裝置102。實驗性擷取101，資料儲存裝置102，以及模擬103皆可由遠端直接或經由網路，例如區域網路(LAN)，廣域網路(WAN)，網際網路等，而對影像重建器105提供投射像素資料104。雖然圖中只繪示投射像素資料104，但通常該資料可以是為將任何形式的激發能量發送進入被觀察物體所獲得的結果，其因此即可能包括來自於，例如，電子顯微術(electron microscopy)，磁共振(magnetic resonance，MR)，正子放射斷層造影術(positron emission tomography，PET)，單一正子放射電腦斷層造影術(single positron emission computed tomography，SPECT)，超音波，螢光，多光子顯微術(multi-photon microscopy，MPM)，光同調斷層掃描(optical coherence tomography，OCT)等，以及電腦斷層掃描(computed tomography，CT)。資料儲存裝置102可為，例如，電腦記憶，CD-ROM，DVD，磁光碟(magneto-optical，MO)，硬碟，軟碟，zip碟，快閃碟(flash-drive)，等等。

圖2顯示依據本發明某些實施例之一成像系統。此成像系統可包括：一控制器205；一介面206(例如圖形使用者介面，鍵盤，搖桿，等)，以供與一研究者211(例如，一名人類使用者，一部電腦等)之信號通訊，以及供與控制器205進行同步；發射器201(一或多部)，供反應於來自於控制器205的一控制信號而產生並發射激發能量202；以及偵測器207(一或多部)，其被組構來產生投射像素資料104。投射像素資料104可為數位形態並被儲存於一儲存裝置之中。成像系統亦可包括：一平移或旋轉台204，其被組構可在其上承接一物體203，並可被操作而相對於發射器201(一或多部)及偵測器207(一或多部)進行移動或轉動；以及一部影像重建器105，209，其係電性地互相耦接，或者經由選擇性的資料傳輸器208，諸如網際網路，其他網路，或資料匯流排而耦接至偵測器207(一或多部)與控制器205。資料匯流排可為，例如，可在一部電腦的電腦組件之間，或在電腦之間傳輸資料的電纜線的一個次級系統。雖然在相連結的組件之間所流動的資料流是為單向的，但在連結的組件之間所進行的通訊亦可以是雙向的。

激發能量202可為放射性形態的能量，諸如X射線能量，電磁(EM)能量，光能量，紅外線(IR)能量，粒子能量(例如電子，中子，原子射束)，振動能量，諸如超音波等。激發能量202可被發射至一物體203上，其可為，例如，一假體(phantom)，一名人類病人，一樣本，或其任何之組合。激發能量202可由對應能量形態的發射器201(一或多部)進行放射。激發能量202可傳遞通過物體203且其一部份可被適當的偵測器207(一或多部)所接收。偵測器207(一或多部)可將其接收到的能量轉換成為可量測的電信號，並再進一步以數位格式將被量測的電信號轉換成為投射像素資料104。

控制器205可為連接發射器201(一或多部)與偵測器207(一或多部)的一個電路，並可對發射器201(一或多部)與偵測器207(一或多部)送出控制信號，以使激發能量202的傳輸與偵測器207(一或多部)的操作同步。此電路可為類比，數位，或混合信號的形式。控制器205亦可為具有一或多個處理器與一或多個記憶的一部電腦，以將控制信號發送至發射器201(一或多部)與偵測器207(一或多部)。

影像重建器105，209可對控制器205進行反應，並可接收投射像素資料104以經由依據本發明某些實施例之方法而重建物體203之影像，以利用其高計算效率而產生出物體的高真實度影像。影像重建器105，209可為具有儲存了依據本發明之說明性質實施例而進行操作的軟體碼的一或多個資料儲存裝置的一部電腦。此電腦可包括有一或多個處理器以及一或多個記憶體，以便讀取並執行儲存在一或多個資料儲存裝置內的軟體碼。在另一說明性質之實施例中，影像重建器105，209可包括有一或多個程式儲存及執行裝置，其可依據本發明之說明性質實施例而操作。影像重建器105，209可產出重建影像106。

一輸出裝置107，210可以接收重建影像106。輸出裝置107，210可為一視覺顯示裝置或為一資料儲存裝置。視覺顯示裝置可為，例如，一顯示器或一列印裝置。顯示器之實例可包括，例如，陰極射線管(CRT)，發光二極體(LED)顯示器，液晶顯示器(LCD)，數位光學投影(digital light projection，DLP)監視器，真空螢光顯示器(vacuum florescent display，VFDs)，表面傳導電子放射顯示器(surface-conduction electron-emitter display，SED)，場放射顯示器(field emission display，FEDs)，矽上液晶顯示器(liquid crystal on silicon，LCOS)等。列印裝置之實例可包括，例如，使用碳粉之印表機，液態噴墨印表機(liquid ink-jet printers)，固態墨印表機(solid ink printers)，染料昇華印表機(dye-sublimation printers)，以及諸如熱感印表機與紫外線印表機(ultraviolet(UV) printers)等的無墨印表機。

成像系統更可包括有一研究者211。研究者211可為人類或執行程式之電腦以接收來自於輸出裝置210或影像重建器209的重建影像106，以便接著施用一演算法(例如，預先規劃的程式，人工智慧，機器學習，等)而抽取出有關於物體203的診斷資訊，或者進行發射器201(一或多部)，偵測器207(一或多部)，台204，影像重建器209等的控制參數之微調。在某些實施例中，介面206或輸出裝置210可能並非必要。在某些實施例中，研究者211，控制器205，以及影像重建器105，209可以建置在相同的一部電腦上或不同的電腦上。

本發明之某些實施例可以提供包括有一或多個處理器的工作站，其被組構來以相類似於影像重建器105，209的方式而重建影像。工作站可以接收來自於成像系統，資料儲存裝置，或電腦等之至少其中之一的輸入資料。輸入資料可經由一資料匯流排，一電纜，一有線網路，一無線網路等而接收。工作站更可包括有一輸出裝置107，210以接收重建影像。輸出裝置可為一資料儲存裝置，一顯示裝置，一列印裝置等。資料儲存裝置，顯示裝置與列印裝置的實例係如前述。

圖3顯示依據本發明一說明性質實施例之流程圖。

在方塊301中，可選擇一適當的限定(調節)函數f(x)，諸如：

-a<f(x)<b　式(1)

其中-a<0<b且就x=0的某些鄰區而言，依據下式，函數f(x)係近似於x

|f(x)-x|<ε,就適當的較小ε>0.　式(2)

函數f可微調超過之數值，特別是在早期的疊代之中。例如，其一種函數是為atan(x)。限定函數可被釋鬆，且其可能沒有嚴謹的截斷值a及b。另一個限定函數的例子可為h(x)=c‧x+(1-c)‧atan(x)),0<c<1。Box,G.E.P.,Jenkins,GM.,Time Series Analysis: Forecasting and Control,San Francisco,CA: Holden-Day,1976之中包含了可使用做函數f的非線性資料轉換的一種選擇。

在方塊302中，m維(m-D)像素的輸入投射資料104，就每一給定的投射方向而言，可以被定置於一系列的變動解析度的解析度格柵上(例如，每一個接續的解析度格柵皆有比其前者為少的像素，但同時則維持其影像的範圍)。例如，系列中的一個解析度格柵可以具有該系列中之前者一半的解析度。其最粗的尺寸可為，例如，單一m-D像素。

應注意通常1<m<n且在所有的投射之中m可能全皆相等。不過，一般而言，m可能隨著每一投射而變動，且m>=0(例如，就某些吵雜的估算問題而言)。

同樣亦在方塊302中，一啟始組的n維(n-D)立體像素可能被置於n-D空間之中，以使其投射涵蓋該列之解析度格柵上的可得m-D投射像素。例如，此等啟始立體像素的適當組的數值可能為對應於一平均投射密度的一個數值；且就每一立體像素j而言，可能會設定v_Lj=log(v_j)。除了對數以外的其他單調函數皆可能可以應用。例如，啟始單一n-D立體像素密度v_i可能被設定等於平均投射密度。

在方塊303中，目前格柵解析度的n-D立體像素與m-D像素兩者皆可選擇性地利用適當核心(kernel)而予以平滑化，其中m-D核大約是為n-D核的投射。此操作可對應於含雜訊卡爾曼濾波器(Kalman filter)的預測性演算法(不含更新，其中v係由v_L而得)。其一實例函數係為，例如，v=exp(v_L)。其他的函數，包括不連續性，皆可被用來計算額外的密度限制。一說明性質的平滑化核可為具有西格碼寬度(sigma-width)的高斯形狀(Gaussian shape)，例如，約在0.4與2立體像素的間隔之間。不過，如同以下配合方塊307所說明的，其改變是有可能的。

若有需要，立體像素v可利用針對某些數值的對映而予以限制。例如，在某些區域中，其可能得以推知其立體像素密度為零。同樣地，為了避免極小立體像素密度情況下的數值下溢(numerical underflow)，v_Lj可被限定在一個極小可忽略值上。一說明性質的限定數值可代表最大立體像素密度值的約為10^-15的密度。

在方塊304中，預測投射m-D像素p_p可依據n-D空間中該組n-D立體像素v而進行計算。預測投射m-D像素p_p可能是在該系列解析度格柵的一個解析度格柵上。若其格柵解析度已由一先前解析度格柵上有所增加，則預測投射與立體像素估算值v便可能由粗/先前解析度格柵上利用內插而獲得。

在方塊305中，描繪預測投射資料與解析度格柵上輸入投射資料104之間之差異的第一量便可進行計算。

若輸入投射資料104的投射像素值p_i大於零，則第一量便可能為，例如，

u_i=f(log(L_i))　式(3)

其中L_i是為輸入投射資料104的測得p_i與預測投射資料的預測像素值p_pi之間的誤差比，其係以，例如，下式表示

L_i=p_i/p_pi　式(4)

L_i的其他合適關係式同樣亦可使用，例如，L_i會隨p_i增加並隨p_pi減小的該些函數。此些函數可以補捉到輸入投射資料104與預測投射資料之間誤差的大小。

上列之對數函數可為錯誤比例的對數，但其亦可為輸入相對於預測投射值的一個更為一般性的函數(例如，應用在傳輸電子顯微術(transmission electron microscopy)之遺失楔塊問題(missing wedge problem)中的函數)。

若輸入投射資料104的投射像素值p_i不大於零，則第一量便可能為，例如，

u_i=-a.　式(5)

其中a表示一個限制。

在方塊306中，對應於n-D空間中具單位強度之一參考立體像素的至少一脈衝反應R_p之量是可加以計算的。依據目前之解析度格柵，脈衝反應R_p可為所有p投射與反投射的每一所需之參考n-D立體像素而計算出來。脈衝反應R_p可隨參考n-D立體像素之位置而定。脈衝反應R_p與其反數可為m-D投射或為n-D立體像素而計算出來，此係依較佳之操作程序而定，而此程序則依所要精確度而定。此係類比於FBP。參考n-D立體像素的一個有限子集可為脈衝反應R_p之評估而選定。在某些應用之中，在重建空間內中心處可以選定最少單一個的n-D立體像素。

脈衝反應可以數種方式而為對應之參考n-D立體像素計算出來，而此些脈衝反應的反數則可用於與第一量的後續結合，如同以下針對方塊307所說明的。

在一說明性質之實施例中，在n-D重建/物體空間中一區域內的一特定位置上的脈衝反應R_p可為當考量單位強度之一參考立體像素時的一線性系統之脈衝反應R_p。

當在一特定位置上的此單一參考n-D立體像素被投射入所有p投射方向上時，即可以產生出p個單一像素投射(其不必然全皆具有相同的維度m)。當此些p個投射皆為重建所需的反投射時，便可以產生出原始立體像素的一個模糊版本。此模糊版本可以原始立體像素線性相關聯，但亦可能分散於整個n-D空間中。此模糊函數R_P可能為此特定位置上重建系統之脈衝反應(或點散佈函數(point spread function)，PSF)。

就一有限重建體積而言，此脈衝反應R_P可能隨參考n-D立體素的位置而定。例如，當投射的數量足夠少且其格柵亦夠細，則脈衝反應R_P便可能為蜘蛛形(例如，其可能顯現出「蜘蛛腳」擴張狀或由「蜘蛛」中心位置條狀發散而出)。其相鄰之立體像素可能具有類似的「蜘蛛腳」形狀的延伸，但亦可能因其有限重建區域可能在其邊緣處截段其延伸而具有不同的長度。此外，在使用p個錐形射束投射的諸如錐形射束之幾何形狀之中，「蜘蛛腳」的方向及相互角度可能會隨位置而改變，因「蜘蛛腳」可能會指向p錐的固定尖端。

不過，為了達成近似逼近之目的，此等脈衝反應函數的緩慢變動特點，可能被利用來在n-D空間中讓立體像素之足夠小的鄰近區域內可以分享相同的脈衝反應函數。「蜘蛛腳」的近似程度上的幾何考量，可以決定分享了相同近似脈衝反應且其符合良好之立體像素區之大小。其在「蜘蛛腳」的尾端點之處的例如數個(比如，數十個)百分點之相對於其內容的整體不符合度，可能相當可以被接受，特別是當其接續的疊代足以更正此些錯誤時為然。相較之下，FBP法通常並不考量此等截斷效應，而諸如Wood所教示的矩陣建構則會加以考量。

同樣亦在方塊306中，一給定立體像素的脈衝反應R_P接著即可與來自一雜訊變化的部份結合以獲得R，其即為第二量。輸入投射資料即可予以量測以獲得雜訊變化。此結合可以增加「蜘蛛」中心之處的值。R及其反數接著即可為每一m-D投射或為每一n-D立體像素而計算出來，此係隨操作的較佳序列而定，而此較佳操作則隨所需之精確度而定。其較小的增加，例如當所假設之雜訊位準為低時，則可以顯著地增加對R反數之精確計算的需求(例如，對應於卡爾曼/威納濾波器(Kalman/Wiener filter)增益)。不過，相當方便地，在多解析度格柵上的計算則可將殘餘投射誤差中之信號部份維持在中度/低度位準上，如此即得以支援R^-1的有效率計算。此信號部份可為進行重建時殘餘投射誤差中所未被考量到之構造資訊。

在方塊307中，計算所得之第二量接著即可計算反數以獲得P，其，例如，可表示為P=R^-1(以取代赫士反數之使用)。此步驟可在傅立葉域中進行。本發明之說明性質實施例在傅立葉域中可將此反函數歸類為卡爾曼或威納增益，如同以下所將討論者。

測得之輸入投射資料Y之一觀察模型可為，例如，Y=H X+W，其中X為可代表物體密度，H可代表系統轉換函數，而Y則代表輸入投射資料104。

在矩陣之表示之中，X可依可得之資料Y而估算出來。資料Y可包括來自於雜訊之量測所貢獻者。X可具有以下之統計性質：

X₀=E[X],and　式(6)

V_X=E[(X-X₀)‧(X-X₀)^t]　式(7)

而投射雜訊W統計則可表示為：

E[W]=0,　式(8)

V=E[W‧W^t],及　式(9)

E[W‧X^t]=0　式(10)

最小變動(卡爾曼)增益，K，如同前述，在例如Sage之文件(即，Sage,A.P.,J.L. Melsa,Estimation Theory with Applications to Communications and Control,New York: McGraw-Hill,1971,page 262)中可表示為

K=V_XH^t[H V_X H^t+V]^-1=V_X H_t R_Y ¹　式(11)

或為

K={H^t[I+V(H V_X H^t)^-1]}^-1　式(12)，

利用複數頻譜表示，一物體之功率-頻譜可為S_XX，其觀察轉換濾波器可為H，投射雜訊頻譜功率密度可為S_VV，則其威納(例如，見Wiener,Norbert,Extrapolation,Interpolation,and Smoothing of Stationary Time Series. New York: Wiley,1949)增益即可表示為：

W=H* S_XX/(H*S_XXH+S_VV)　式(13)

=[H(U+S_VV/(H* S_XX H))^-1　式(14)

頻譜公式就斷層掃描用途而言在計算上可能是可行的。在斷層掃描重建上，物體密度的不確定性(例如，S_XX或V_X可能相對為大)可能獨佔了投射雜訊(例如，S_XX或V可能相對為小)，而尋求威納增益W(或卡爾曼增益K)的主要任務即可能變成尋求式(13)(或等效式(11))的有用的近似。

如此，輸入投射資料104即可被測量而得到V以及V_X，如上式所示。H V_X H^t與H* S_XX H亦可計算，例如，依據所量測得之V_X而計算。方塊306中所提及，一雜訊變動之貢獻可能是指因，例如，V，H V_X H^t，或H* S_XX H所引起之貢獻。雖然卡爾曼增益之式(12)以及威納濾波器增益之式(14)皆可直接進行計算，但其進一步之簡化仍是可能的。例如，就雜訊W量測之小交叉關聯(例如，可能是以對角為主)以及投射物體密度的中度交叉關聯(例如，H V_X H^t可能是以對角為主，特別是在疊代程序時)而言，式(12)或(14)的複合表示便可利用為個別(或其群組)投射進行區段式增益計算而得以簡化。例如，區段式處理可以利用，例如，在立體像素至立體像素(n-D處理)或像素至像素(m-D近似處理)基礎上的修改過局部脈衝反應之反數，而進一步予以簡化。局部脈衝反應之修改可以，例如，式(12)中所見之雜訊對信號構造之元素V(H V_X H^t)^-1，或簡單僅利用式(14)中所見之S_VV/(H* S_XX H)而給定。此修改可在立體像素至立體像素的基礎上施行。此修改可利用增加「蜘蛛」中心區(例如，接近脈衝反應函數之峰值之處，對應於(接近)式(12)中之矩陣V(H V_X H^t)^-1或式(14)中之S_VV/(H* S_XX H)的對角元素)的值而進行逼近。在W為白雜訊且投射物體之不確定性為白值且具有估計均方根(RMS)雜訊信號比為r的情況之中，蜘蛛中心即可以，例如，增加一個(1+r²)的係數。在此修改之後，卡爾曼增益與威納濾波器增益即可利用數值反數，根據，例如，式(12)及(14)而予以計算。

注意到，R_P，如同前述，係與上述矩陣(或轉換濾波器)H相關，如式(11)及(13)中所顯示的。處理程序可以利用將卡爾曼增益或威納增益中之單一步驟打斷成為兩個步驟的向量操作而加以改變。此兩步驟可為，首先反投射，例如，投射資料，並再以參考立體像素K的一n-D濾波器P_K=R_K ^-1進行濾波；某些相鄰的環繞立體像素K可以利用其結果之濾波器輸出。此兩步驟亦可為，首先，例如，以p個m-D濾波器P_iK(0<=i<p)對投射資料進行預濾波，其中每一P_iK可由經雜訊對信號比加以修改的n-D脈衝反應(PSF)R_PK之投射的反數而獲得，並接著再將濾波器輸出反投射；再一次地，某些相鄰的環繞立體像素K(或其投射)可以利用其結果之濾波器輸出(更新)。為了覆蓋整個的n-D空間，可能需要使用數個立體像素位置及其鄰近區域。如此，顯示於式(11)及(13)中的極大矩陣的反數便可以分解為遠為較小數量的R_K組的中等尺度組的反數，其中每一R_K只會在物體空間中緩慢地變動，以容許其重覆的使用。再者，參考立體像素的每一R_K ^-1便可利用傅立葉技巧而極有效率地進行計算。

總結之，用來計算增益P的近似法之兩要點可能對高統計效率有所助益。

首先，當低頻時，卡爾曼/威納濾波器增益(現在以增益P加以取代)可利用在脈衝反應函數的原點增加一個小常數而進行估算，此可以涵蓋到在足夠寬廣的頻率範圍內的雜訊信號比。靠近脈衝反應函數中心的一個小鄰近區可利用加入一個「狹窄」的高斯鐘(Gaussian bell)而非僅只是原點上的一個常數而予修改，以便，例如，考量到「粉紅」雜訊量測。

其次，當高頻時，對應於增益P的濾波器便可能被高頻信號雜訊比的特性所傾斜/平滑化。為了增進數值精確度，此傾斜/平滑化函數可在目前的解析度格柵上於影像密度及/或估算之影像密度上執行。被估算影像密度之低通濾波器，其在其系列之解析度格柵上的疊代計算期間的每一階段，皆可輔助將高頻影像部份傾斜下來，以減低計算上的負擔。

在方塊308中，利用第一量與P可以計算一更新值。例如，P可以與第一量結合以獲得一給定立體像素(或像素)影像元素的一個原始更新值i。例如，P可以與第一量迴旋結合(convolve)以獲得原始更新值。當P在n-D重建空間內緩慢地變動時，利用將P做為給定立體像素鄰近的一個合適但足夠小的一個近似值，便可以節省計算時間。利用與先前所討論之式(1)及(2)中之函數f相類似而建構的一個限定函數g，接著即可用來供應予原始更新值I以便獲得一個修改過的更新值i_g=g(i)。此修改可以保護避免發生過度的更正。在本發明說明性質實施例發展的期間，許多直接利用原始更新值的實驗皆失敗了。本案發明人接著即構思並發展出限制函數以及其使用法則以限制過度的更新值之出現。修改過之更新即可為後續處理的更新值。

在方塊304-308中，對錐形射束幾何構形重建(例如，應用於血管放射攝影術，angiography中者)，其中投射之數量p可能相對為小(p<<N，其中N可代表沿著通過重建體積內的一代表性投射路徑的立體像素的數量)者而言，額外的簡化考量有其助益。在此種情況下，個別的「蜘蛛腳」可能不同，特別是在多重重建格柵的高解析度最終階段期間為然。在此高解析度n-D空間中為許多小區域計算脈衝反應R_P可能是耗費時間的事。做為使用n-D脈衝反應的一個近似，可以使用維度m_i(0<=i<p)的p個低維度投射脈衝反應。其對應之修改過反數或卡爾曼/威納增益接著即可疊代應用到投射的個別誤差影像上。此些m-D脈衝反應可為上述n-D脈衝反應於對應之m-D投射空間上之投射，其可能為蜘蛛形者。其性質可能即如前述，亦即，在投射空間內緩慢地變動。不過，計算p個投射的脈衝反應，此些計算得脈衝反應之反數，以及後續的迴旋計算(此時應用於m-D投射上)時的計算負擔，可能遠比其在n-D空間，特別是當p<<N時為少。計算所得之脈衝反應可依與R類似方式進行修改。計算所得之脈衝反應反數亦可進行頻帶濾波以降低斷層掃描重建時的雜訊。斷層掃描反數的此種技巧可以類似於其投射在進行重建之前先行濾波的反投射濾波處理。此技術可能受限於錐形射束幾何形狀，但亦可以應用於其他的幾何形狀之中，諸如，平形射束的幾何形狀。

在方塊309中，修改過之更新值i_g可被加至該組n-D立體像素中之對應影像密度值上以更新該組n-D立體像素v。此新增可能會受到前述式(1)及(2)類似的一個限制函數的限制。此新增可以在整組的n-D立體像素的一部份上進行。為了精確度的緣故，立體像素值v可比例調節以使其投射能與原始影像資料最佳地相符。

在方塊310中，若解析度格柵是最細的且若計算所得之預測值與輸入投射資料之間的差異(例如，殘餘值)(統計上)足夠小，則代表密度估算值v，v_L，或兩者(並且，若有必要，估算之預測值與誤差)的更新組的n-D立體像素，便可輸出至一輸出裝置，如方塊311中顯示。平滑化可以在輸出之前進行。否則，解析度格柵便可改變為具有比目前解析度格柵為大尺度的一個解析度格柵(或者，n-D平滑核心的嚴謹度/範圍便能減低)且在方塊312中，其流程可以進至方塊303。

在方塊311結束時即可獲得物體106的一個高真實度重建影像。此重建影像106可以非侵入及非破壞方式顯現物體203的內部資訊，此在診斷醫療影像與非破壞性評估上具有廣泛範圍的用途。重建影像106可在一部視覺裝置上觀看或儲存於電腦記憶或儲存媒體之中，如同以上所討論的輸出裝置107所可以做到的情形。

來自於Trachtenberg,S.,Dorward,L.M.,Speransky,V.V.,Jaffe,H.,Andrews,S.B.,Leapman,R.D.,"Structure of the Cytoskeleton of Spiroplasma melliferum BC3 and Its Interactions with the Cell Membrane," J. Mol. Biol.,378,pp. 776-787,2008之投射像素資料104被用來展現經由本發明一說明性質實施例所可獲得的改進。圖4A顯示依據投射像素資料104之濾波反投射(FBP)所重建的一個重建影像106。圖4B顯示依據本發明一說明性質實施例之重建影像106。圖4A與4B之間的比較顯示，以諸如細胞間對細胞外空間，細胞膜，核糖體，脊髓液之細節等各種密度的較尖銳對比作為量測標準，圖4B顯然可以獲得遠為更優良的影像清析度。圖4B中的影像可更為使用者所需要。兩種重建皆在約+/- 70度的範圍內使用了約166個雙傾斜傳輸電子顯微術(double-tilt transmission electron microscopy，TEM)投射。

圖5顯示一模擬假體(phantom)之覆面樣本投射。此模擬假體包含有在尺寸為256×256×64個立體像素的一3-D盒體內隨機位置上的約100個尺寸為6×6×6立體像素的立方體。該些立體像素具有理置於0.85單位的環繞密度中的5單位密度。電腦模擬103獲得了具有添加之白雜訊的總共49個投射(單傾斜，+/- 72度)以產生投射像素資料104。

依據統計模式驗證的習知黃金標準，重建量可以最為嚴謹地利用殘餘之檢查，亦即，真實影像(例如，圖5中之推知之模擬假體)與重建影像106之間的差異而評估出來。高真實度的重建可以導致隨機並包含有來自於3-D物體之構造(例如，圖5中之模擬假體)的極小影響的殘餘。圖6A顯示依據同時疊代重建技術(simultaneous iterative reconstruction technique，SIRT)，而圖6B則顯示依據本發明一說明性質實施例之方法之重建誤差。圖6A中相當醒目的殘餘方形圖案(15次疊代之後)以及圖6B中所出現的主要白雜訊(在一次高解析度疊代之後)更進一步地顯示了本發明在影像真實度與計算效率上的超越性。SIRT的疊代次數為此例之最佳，此乃因SIRT的更多或更少次的疊代皆僅只會導致影像品質的衰退。

本發明亦可應用於投射之數量為小數量的情況，例如，血管放射攝影術之錐形射束幾何形狀的情況之中。圖7A顯示在碘對比染料(例如，一般常見心導管血管放射攝影術所採用者)被注入由塑膠管所製作一跳動液力冠狀假體內時之X射線投影。其所使用之投射像素資料104，由於注入動作的暫態性質只由其底下假體的6個數位減除投射所構成。圖7B顯示依據本發明一說明性質實施例之方法所獲得之假體之重建影像106，其精確度代表跳動液力假體的幾何造形。由於投射之數量通常很小且/或疊代之次數以目前之電腦資源而言變得不切實際，此種重見影像是無法利用濾波反投射(FBP)或同時疊代重建技術(SIRT)，甚或最大熵法(maximum entropy，MENT)而獲得的。

如同上述，本發明某些說明性質之實施例結合了計算/處理技術並獲得在目標區域內所需n-D立體像素的快速且精確的計算/估算結果，同時並容許使用任意的m-D投射(及方向)。應注意的是，m並不需固定，而是可隨投射指標而定。此些技術之組合效率可被應用於廣泛投射數量範圍內。可得投射(量測)資訊的有效使用可使此方法當投射的數量極大時或受限時(小數量或不完整)，諸如心血管攝影時，皆可有效。

本發明之非線性反數技術可結合：(i)經由將n-D影像估算分解為重建區域內的兩非線性相關域而將一非線性卡爾曼濾波器的某些部份線性化(近似貝爾新估算，approximating Bayesian estimation)；(ii)一多格柵計算，其可與「階梯」(“stair-cases”)相比擬，並利用使用減低的影像平滑化程度而進一步分割解析度；與(iii)當有需要時，線性系統之反脈衝反應(或轉換函數)的有效率計算，該系統係由利用反赫士而推估的影像雜訊之估算加以修改，並為物體空間中的選定區域而獲得近似卡爾曼/威納濾波器增益。

特別地，非線性卡爾曼濾波器已在Sage,Jarisch W. R. & Detwiler J. S.,"Statistical Modeling of Fetal Heart Rate Variability," IEEE Trans. On Biomed Eng.,Vol. BME-27,No. 10,pp. 582-589,Oct. 1980中有所討論，且卡爾曼濾波器在影像重建中的應用亦已在Jarisch W.,Determination of Material Properties by Limited Scan X-ray Tomography,Technical Report USAF AFWAL-TR-81-4050,NTIS,HC A07/MF A01,Defense Technical Information Center,Cameron Station,Alexandria,Virginia 22314,1981中有所討論。其他的方法，諸如FBP，FFT與威納濾波器，亦可取代反數的線性部份而應用。在重建區內的兩非線性相關域之其中任何一個皆可代表卡爾曼系統狀態。

利用上述結合，其每一次疊代皆可達成錯誤的有效縮減，類似於赫士在最佳化中的有效應用。其所達成之效率可接近於四次收斂(quadratic convergence)。此種作法可以跳過其他方法所會遇到的高維度與不良收斂速度的嚴重問題，諸如Kolehmainen,V.,Brandt,S.S.,Structure-From-Motion Without Correspondence From Tomographic Projections by Bayesian Inversion Theory,IEEE Trans. Medical Imaging,ITMID4,Vol.26,No. 2,pp. 238-248,Feb. 2007中所討論者。

由另一角度觀之，本發明之一說明性質實施例可被視為是IR與(空間變數/非靜止/分段化)FBP的一種結合，在此稱為S-FBP。依近似之需要，可能會有一或多個分段會與S-FBP相關聯。通常，FBP在重建之區域內不會是靜止。不過，區域性質的緩慢改變可能會因將區域的重建煞停進入重疊區/段內，其在此些區之間有加權轉換而容許近似。S-FBP的速度，效率，以及彈性，經由以下所討論數種作法的組合而皆可達成。

例如，S-FBP部份可被應用在非線性轉換資料上，而非原始投射資料。

例如，多格柵計算可支援線性化並使每次疊代的操作次數保持極低。例如，在最終計算最細解析度格柵之前的總操作次數可低於S-FBP在最細解析度格柵上的次數。

例如，立體像素密度可以兩非線性相關域(以一映對加以連結)來代表：其一代表經由S-FBP而更新(其在n維空間中為可變，依所需之近似技術而定)，而另一則代表計算投射預測。S-FBP部份的濾波器加權可依立體像素位置而定，但將立體像素結合成群在一起以供應用於相同的權重上可以增進計算的效率。

例如，利用計算來自線性化之系統脈衝反應(或轉換函數)反數的S-FBP部份的濾波器加權，則任意投射維度與方向便皆可使用。線性化之系統脈衝反應(或轉換函數)可能已與由輸入投射資料之量測值所決定的一雜訊變動結合。

例如，適當選定的限制(微調)函數可以保護避免計算進行期間的數值奇點(numerical singularity)與不穩定性。

例如，適當選擇的平滑化函數，在每一次疊代進行時依規則而施加在立體像素與像素值上，可有助於數值運算的穩定性及精確度。立體像素與像素值的此種平滑化在功能上可比擬於產生更新值之前的卡爾曼預測演算法。平滑化的此疊代應用可類似於卡爾曼濾波器中的狀態預測演算法。

由於最細解析度格柵計算之前疊代的精確度可能已屬足夠，最細解析度格柵的疊代可能在計算上並不必要，因此即得以省下誤差更正項的計算。在此情況下，重建速度可能變得與S-FBP相當。此可能是因最細解析度格柵上的計算可以達成與FBP相當的操作計算次數(在>1的係數之內)，而在較粗解析度格柵上的所有先前較低解析度疊代的結合操作可能少於最細解析度格柵上的計算，特別是較高維度的影像重建。

由於各種原因，上述的非線性更新技巧可能無法在單一次疊代中執行一次完美的更正。相反的，在接續的疊代之中，可能會有些區域其更新太小而其殘餘則在接續的幾次疊代中持續存在，或者，其中可能有些區域的更新太大而其殘餘則顯現出一種限制下的循環(例如，在正負號之間震盪)並持續數次疊代─一種較高增益情況下的常見回授迴路效應。來自於先前疊代的此些殘餘可能便得以有利地被利用來達成更為有效率的更新。

特別地，在每一影像點上的增益(不論是投射空間中給定像素的投射增益或物體空間中一立體像素的立體像素增益)可以做為先前更正/殘餘r_h之歷史函數而加以評估。例如，目前的更正/殘餘可標示為r₀(在二維中的一個給定位置ij或三維中的ijk)且前一疊代者則為r₁。當r₀與r₁具有相同的正負號時，在一位置(例如，當m=2時在位置ij)上的增益因素通常會增加；類似地，當r₀與r₁具有相反的正負號(且未為其他因素所驅動)時，其增益便會減小。可能影響增益修改之選擇的其他因素包括，例如，與更正r_h相關聯的不確定性，來自其他投射之估算驅動力的不確定性(例如，當來自其他投射而施加於殘餘上的增益因素獨佔了重建的更正時)，先前增益的有限精確度，投射影像的雜訊位準，及/或一穩定解的估計收斂程度。

如此，局部投射與立體像素增益便可以依先前疊代中的歷史圖像而進行適應性調整。圖8A顯示依據本發明一說明性質實施例之一開展流程圖。圖8A之開展流程圖除了額外的方塊以外，其與圖3相類似。

在方塊801中，為預測投射資料與輸入投射資料間之差異編碼的第一資料即可被記錄下來(做為投射殘餘)，例如，記錄於一資料儲存裝置之中，諸如儲存裝置102或另一裝置。其差異可由方塊305中獲得。

同樣亦在方塊801中，疊代t時的第一量r_t(來自方塊305)可以一增益矩陣G _t進行加權，以便依據先前疊代所記錄的第一資料而獲得一個更新第一量。

就一給定低(例如，4×4或更低)格柵解析度的所有像素i而言，一個適當的初始增益矩陣G ₀通常會與一單一固定值(G_0i=G₀>0)相關聯。增益矩陣G _t的初始固定元素值G_t接著即可利用對第t^th次疊代時增益矩陣G _t的對應區施以矩陣g _t之增益改變因素g_t(例如，利用一元素接一元素地在對應的位置上相乘擴大)而增加或減小。其增益改變因素g _t係隨觀察到的穩定性/更新的震盪/遞迴(recursions)而定。

一個聰明選定的初始值G₀可能小到足以安全地避免估算物體密度的震盪但又大到足以達成估算物體密度的有效更正。某些數值嚐試可能會引導到一給定重建設定的合適初始值。一最佳化值的選擇在此點上可能不重要，因為最佳值終將會在後續的疊代中被找出來，如同以下所將討論的。

在首先前兩疊代之後，當獲得二或更多接續組的投射殘餘時，每一個別投射像素的增益矩陣G _t的調整即可開始。

調整增益矩陣G _t的一個目的可以是避免密度估算之震盪(例如交替變動)，其為太高增益的一個特性。解決交替變動投射殘餘(或其等效之物體密度震盪，如同以下配合方塊803所討論者)的一種可能調整實例顯示於式15中：

其中g_t是為估算之增益改變因素，而r _s與r _t則表示來自於一給定像素位置之接續疊代的殘餘s與t(s=t-1)。

式15是為增益改變因素g_t的一個實例。不過，式15因出現在殘餘r_s與r_t的雜訊通常可能導致由一系列增益改變因素g_t所獲得的增益矩陣G _t之過度變動，而可能不被直接使用。

為了避免增益矩陣G _t之元素因，例如，被r_s-r_t的差值之近乎零值所除而成為過大的數值，有數種預防作法可以執行，例如：

(i)　殘餘中雜訊的一個局部量測可利用，例如，估算投射殘餘中的局部雜訊功率(例如，依所表示的一參考區域的標準差的量測)而獲得。

(ii)　當r_s與r_t具相反正負號但|r_s|與|r_t|皆小於sig_ref，或當r_s與r_t具相反正負號且|r_s-r_t|<2‧sig_ref時，即可能被定為一單位。

(iii)　在(ii)並不適用的情況下，當式15的分母中r_s或r_t的絕對值低於sig_ref以下時，r_s或r_t的可隨需要而直接由sign(r_x)‧sig_ref加以取代，以使|1-g_t|具最小值，其中x{r,t}。

(iv)　當r_t的大小持續地大於r_s而r_t與r_s具相同正負號時，在特定位置上的潛在性控制喪失情形即可能發生。局部增益因素g_t的增加可能使控制喪失的問題和緩。在該位置上局部增益改變(增加)因素的一種可能選擇可能為：

其中s=t-1,q=s-1,且

(v)　一局部增益改變因素g_t利用一最大與一最小限定值而可能在範圍上被限制：

g_min<g_t<g_max　式(18)

其中

g_max=C‧n_pro　式(19)

g_min=1/g_max,　式(20)

且其中n_pro是為所使用之投射數。C的合適值通常在I以下，例如0.25。G_t的元素值可以類似的方式加以限制，例如，在一投射的所有位置i與j上，max{G_ti}/min{G_tj}<0.5*n_pro，其中n_pro是為可獲得投射的數量。

(vi)　局部增益改變因素g_t中的雜訊可利用將局部增益改變因素在一個小鄰近區域內加以平均化(平滑化)。就2-D投射資料而言，一3×3平均化陣列可能傾向於在疊代重建程序的收斂性上獲致顯著的改進。

在方塊801中，如同前述，當被記錄之第一資料接續其先前疊代而具有相同的正負號時，被更新之第一量可能傾向於增加投射增益，而當被記錄之第一資料接續其先前疊代而具有相反的正負號時，被更新之第一量則傾向於減小其投射增益。

如此，隨著來自於先前疊代的投射殘餘而變的一適應性增益矩陣即已說明如上。此適應性增益矩陣亦可隨像素在投射空間內的位置而變。

在方塊802中，一更新值可依類似於先前方塊308所說明方式，利用P，來自方塊307的反函數，以及來自方塊801的更新第一量或來自方塊305的第一量而被計算出來。

在方塊803中，將更新值編碼的一第二資料可被記錄於資料儲存裝置之中，諸如儲存裝置102或其他。更新值可依據先前疊代所記錄的第二資料而進行更新。適應性微調更新值的一個增益矩陣G _t可由更新值的序列u_t計算出來，可能以及投射殘餘的次組，其可以類似於先前方塊801所描述之方式而由先前之疊代中計算出來。與方塊801相比，利用更直接地減低估算錯誤，對更新值應用增益矩陣G _t即可以適應性地估算出物體密度估計值。在方塊803中，其更新值可與一立體像素增益相關聯。當所記錄的第二資料具有與其接續先前疊代中相同的正負號時，其修訂過之更新值可能傾向於增加立體像素增益，且當所記錄的第二資料具有與其接續先前疊代中相反的正負號時，其修訂過之更新值可能傾向於降低立體像素增益。

在方塊804中，來自803的修改過更新值或來自308的更新值可被加至該組n-D立體像素v中之對應物體密度之值，以便更新該組n-D立體像素v。此值之加入可能會被類似於先前所討論式(1)與(2)的一限定函數所限制。此外，所估算之物體密度可能受限於某推定之已知值，諸如，零/可忽略之密度。再者，此值之加入可能是在整組n-D立體像素的一部份上進行。為了精確度之故，立體像素值v可比例調整以使其投射可以最佳地與原始投射資料匹配。

如此即說明了依據來自先前疊代之物體估算誤差的一種適應性增益矩陣。此適應性增益矩陣可以考量到重建之幾何因素以及資料相關因素。幾何因素可為，例如，一物體已知之中空區域(零密度)，其不需要有增益，且其估計物體密度被重置為零。

圖9A至9D分別顯示10，50，130，與170度軸傾斜角時之輸入投射資料之實例。輸入投射資料可由圖8A之方塊302中獲得。就圖8A中由方塊303至方塊312的每一疊代循環而言，其每一投射角度皆會獲得一個預測投射矩陣。顯示於圖9A至9D之輸入投射資料係來自於傳輸電子顯微術/掃描傳輸電子顯微術(TEM/STEM)，其係由卡爾斯如技術學院以一度角度漸進增加的投射角提供。

圖10A至10D分別顯示10度軸傾斜角時1，2，5，及6次疊代之預測投射矩陣，其係由圖8A之方塊304中獲得。一組大約140次的投射被用來進行重建。其預測投射資料係對應於10度的選定單軸傾斜角。

圖11A至11D分別顯示10度軸傾斜角時1，2，5，及6次疊代之原始輸入投射資料與預測投射矩陣之間之殘餘實例(於最高解析度格柵)，其係由圖8A之方塊305中獲得。

圖12A至12D分別顯示10度軸傾斜角時1，2，5，及6次疊代之結果增益矩陣G _t，其係由方塊801中獲得。具有持續殘餘值的區域可利用增加的投射增益加以處理，如同其中，例如，某些中心珠之增加的亮度所顯示的，而其他顯現了限制循環的區域(亦即，在接續的先前疊代中出現了交替的正負號之震盪情形)，則以減低的增益進行處理。其顯示係為對數比例尺。每一影像左側的參考尺度表示其每一灰階位準上10%的改變。典型的動態範圍增益(亦即，像素或立體像素的高至低比例)係在10至100的係數範圍，此表示其計算上的對應加速情形。

上述更新第一量的適應性增益矩陣與更新值，兩者皆是在圖8A之方塊303至310的回授循環內，供計算物體密度之估計值。傳統上，一個「良好」更正增益的選擇可被應用在疊代的誤差量測(即，殘餘)上，其可能與一塑形矩陣結合，以做為進行有效率計算所需的一個重要設計參數，如同，例如，Pan,T.-S.,Yagle,A.E.,Acceleration and Filtering in the Generalized Landweber Iteration Using a Variable Shaping Matrix,IEEE Trans. Med. Imaging,ITMID-4,Vol.12,pp.278-292,June 1993文件中所討論者。一個不適當地小的更正增益係數可能導致朝向一個解的極為緩慢的收斂情形，而一個過大的更正增益係數則可能造成限制循迴或疊代的震盪，並有可能造成解的發散。限制循迴與震盪可能對影像的重建造成摧毀性的效果。一個「良好」更正增益係數的有規則選擇因此即可能顯著地增進疊代演算法的執行性能。

如同本案發明人所發現的，以一個中等程度的高更正增益係數所進行的殘餘分析，其結果顯示限制循迴/震盪通常可以被限定在重建影像(或對應之投射影像)的殘餘的某些區域內。此外，投射的限制循迴/震盪可在對應的重建影像區中找到，此表示重建物體密度的某些部份可能會受限於限制循迴或震盪。根據發明人的此一新穎發現，在投射空間及/或物體空間內，可以應用中等區域性增益的空間性改變，以使有觀察到限制循迴/震盪的區域，其投射增益或立體像素增益得以被減低，而其收斂緩慢的區域，其投射增益或立體像素增益則可以被增加。

此外，一增益矩陣的空間構造可被利用來將增益值最佳化。如圖9A至9D與12A至12D中所可見，增益值之圖形傾向於追隨輸入投射資料的圖形。若要應用此些系統性的局部關係並減低一吵雜增益矩陣中不確定性之數量，智慧性的平滑化可被利用來獲得更為精確的增益估計值。例如，當使用平面形投射(維度m=2)時，利用將一回歸模型安裝在環繞著投射空間或者物體空間中的一元素的一個適當的，通常小的，鄰近區域上，則所估算的增益矩陣g _t便可被平滑化。可供安裝在覆蓋著相關區域上的數個像素g _ti的一種可能回歸模型可為：

其中t標示疊代t，I可代表一相關區域內的一小組的元素，α_k標示一組K回歸係數的一元素，p可代表投射像素值，而ε則表示不確定性。例如，值為2的K可能即足以利用適當的速度執行其演算法。此外，g _ti的一個估算，諸如最小平方估算，可以被利用來獲得局部估算值。遠較為平滑的可被用來修改顯示於方塊801及繪示於圖11A至11D的m維第一量。相類似的一種作法可被應用於與物體密度估算的n維更新─顯示於方塊803中─相關的增益改變係數，依據此種作法，具有持續性正負號的更新的區域內的增益可被加大，而在其限定循環/變動正負號的該些區域內的增益則可被減小。如方塊303中所顯示的，可為m維投射p_i...j，或計算得之n維物體密度p^obj _i...k，或參考值的其他組合施行一次平滑化程序。

最後，由於來自於方塊801與803中先前疊代的第一資料與/或第二資料需予維持之故，前述程序的單純直接施行可能極度需求記憶體的使用。為了降低改進估算時對於記憶體的需求，可以使用降低精確度的作法。例如，其殘餘與改正可能僅只是原始影像資料或估算物體密度的尺度的一小部份。因此，第一資料或第二資料，以及方塊305中的第一量或方塊306中的第二量，皆可以低於原始輸入投射資料或物體密度代表值的精確度的壓縮格式而以數位形式儲存。此外，具中等精確度(例如，數個位元)的增益因素的估算可能即已足夠。總結言之，減低精確度的資料表示可將對記憶體的需求維持在中等程度，並且只會在計算效率上產生極少的不利衝擊。

圖8B顯示依據本發明另一說明性質實施例之另一開展流程圖。除了圖8A的方塊306，307與802係以圖8B的方塊805加以取代之外，圖8B之開展流程圖係類似於圖8A。在方塊805中，更新值可利用適當的斷層掃描重建演算法，諸如使用習知之中央切塊理論(central-slice-theorem，CST)，濾波反投射演算法等的演算法，依據來自方塊805的修改過之第一量而加以計算。對於高品質的重建，特別是當使用錐形束幾何時，物體空間中一個相關區域的一代表性立體像素的局部投射幾何(例如，由交叉之投射束所給定者)，可被應用在計算第一量，並應用做為選定斷層掃描重建演算法的一個輸入。在錐形束重建的情況之中，可能會需要選擇一組代表性的立體像素。就為了該組代表性立體像素所執行的第一量之每一全體積斷層掃描重建而言，代表性立體像素的間隔，可能會被交叉之投射束的特徵，如同前述，例如，方塊306，變為重要時的距離所驅動。接著，環繞著代表性立體像素的個別區域之輸出可適當地權重其輸出而被接合。輸出對應於一特定重建區的加權可能掩蔽了所有其他的加權，例如，在參考立體像素上逼近於一。接著，重建之接合輸出可被一增益矩陣G _t '所修改，並被用來修改該組的n-D立體像素。雖然不同的增益矩陣G _t與G _t '可分別在方塊805的輸入及輸出處設計，但有數種原因使得單一增益矩陣較為有利。例如，使用數個矩陣可能會有通過方塊312而回到方塊303的回授循環的增加計算不穩定性之風險。

本發明之說明性質實施例，其相對於圖3，8A及8B所顯示的方法所討論者，可以軟體碼的形式而提供，並儲存於一資料儲存裝置中，諸如，CD-ROM，DVD，磁光碟(magneto-optical，MO)，硬碟，軟碟，zip碟，快閃碟(flash-drive)等。其所儲存的軟體碼可由具有一或更多個處理器，一或更多個記憶裝置，諸如隨機存取記憶(RAM)裝置，動態RAM(DRAM)裝置，快閃記憶裝置，以及靜態RAM(SRAM)裝置等的電腦讀取並執行，以便實施前述相對於圖3，8A及8B所討論的說明性質方法。

本發明之說明性質實施例可以提供一或多個程式儲存及執行裝置，例如，特定用途集積電路(application specific integrated circuits，ASIC)，可程式閘陣列(field programmable gate arrays，FPGA)，複雜可程式邏輯元件(complex programmable logic devices，CPLD)，特定用途指令集處理器(application specific instruction-set processors，ASIP)等，以供實施前述相對於圖3，8A及8B所討論的說明性質方法。

前述例子與實施例係為非限定性質之實例。

本發明已經利用說明係質實施例說明如上，習於本技藝者皆可理解，前述說明可予變化及修改而不脫離於本發明之精神範疇。本發明因此即應依據後附申請專利範圍乙節所界定。

101‧‧‧投射影像獲取

102‧‧‧儲存裝置

103‧‧‧模擬

104‧‧‧投射像素資料編碼物體資訊

105‧‧‧影像重建器

106‧‧‧物體之重建影像

107‧‧‧輸出裝置

201‧‧‧發射器

203‧‧‧物體

204‧‧‧表

205‧‧‧控制器

206‧‧‧介面

207‧‧‧偵測器

208‧‧‧資料傳輸器

209‧‧‧影像重建器

210‧‧‧輸出裝置

211‧‧‧研究者

本發明之說明性質實施例配合所附圖式予以詳細說明。圖式之中：圖1顯示依據本發明一說明性質實施例之系統圖。

圖2顯示依據本發明一說明性質實施例之一成像系統。

圖3顯示依據本發明一說明性質實施例之流程圖。

圖4A及4B分別顯示依據一習知濾波反投射(filtered back-projection，FBP)技術所重建與依據本發明一說明性質實施例之方法所重建之細胞影像。

圖5顯示一模擬假體(phantom)之覆面樣本投射。

圖6A及6B分別顯示依據同時疊代重建技術(simultaneous iterative reconstruction technique，SIRT)與依據本發明一說明性質實施例之方法之重建誤差。

圖7A顯示在碘對比染料被注入一跳動液力冠狀假體內時之X射線投影。

圖7B顯示依據本發明一說明性質實施例之方法所獲得之假體之重建影像。

圖8A顯示依據本發明一說明性質實施例之一開展流程圖。

圖8B顯示依據本發明另一說明性質實施例之另一開展流程圖。

圖9A至9D分別顯示10，50，130，與170度軸傾斜角時之輸入投射資料。

圖10A至10D分別顯示10度軸傾斜角時1，2，5，及6次疊代之預測投射矩陣。

圖11A至11D分別顯示10度軸傾斜角時1，2，5，及6次疊代之原始輸入投射資料與預測投射矩陣之間之殘餘。

圖12A至12D分別顯示10度軸傾斜角時1，2，5，及6次疊代之結果投射增益。