快速分層x線斷層造影方法和裝置的制作方法

            文檔序號:6496127閱讀:708來源:國知局
            專利名稱:快速分層x線斷層造影方法和裝置的制作方法
            技術領域
            本發明涉及X線斷層造影,尤其涉及用于從投射創建像素圖像的方法和裝置。
            背景技術
            斷層圖像重建是一種眾所周知的技術,幾乎是所有診斷顯像醫療器械的基礎,包括計算機X線斷層造影(CT)、正電子斷層掃描儀(PET)、單光子斷層掃描儀(SPET)、以及用于磁振造影(MRI)的某些獲得方法。它還發現在關于非破壞性評價(NDE)、安全掃描的制造業、合成孔徑雷達(SAR)、射電天文學、地球物理學和其它區域中的應用。
            斷層圖像數據有若干種主要格式(i)平行光束,其中沿多組平行線執行線積分;(ii)發散光束,其中沿多組發散為扇形或圓錐形的線執行線積分;(iii)曲線,其中沿多組諸如圓、橢圓或其它閉合或開曲線的曲線執行積分。斷層圖像重建的一個問題是從其線積分投射集中重建2D或3D圖像。斷層圖像重建的另一個問題是從其面積分投射集(即其在曲面族上的積分)中重建3D圖像。例如,3D Radon變換涉及圖像在距離原點各個方位和距離的2D平面族上的積分。斷層圖像重建的一些問題、以及一些重建方法在標準參考書中描述,諸如F.Natterer的“TheMathematics of Computerized Tomography”(計算機X線斷層造影的數學),1986年John Wiley出版社在Chichester出版;F.Natterer和F.Wubbeling的“MathematicalMethods in Image Reconstruction”(圖像重建中的數學方法),2001年在費城工業與應用數學學會發表;A.C.Kak和M.Soaney的“Principles of ComputerizedTomographic Imaging”(計算機X線斷層造影的原理),1988年IEEE出版社在紐約出版;以及S.R.Deans的“The Radon Transform and Some of Its Applications”(“Radon變換及其一些應用”),1983年Wiley出版社在紐約出版。
            選擇用于斷層圖像重建的方法是過濾反投射(FBP)或回旋反投射(CBP),使用不加權的(平行光束或Radon變換情形中)或加權的(大多數其它情形中)反投射步驟。該步驟是技術中的計算瓶頸,其中計算需求規模為對2D的N×N像素圖像為N3,對3D的N×N×N體素圖像為至少N4。因而,從N到2N倍增圖像的分辨率導致計算量增加約8倍(或者3D中為16倍)。當計算機已變得快得多時,同時出現了能夠實時收集更大量數據的新技術(例如使用多行檢測器的心臟成像、介入成像)、以及3D采集幾何結構的盛行,對快速重建技術有越來越多的需求。快速重建可加快成像過程、降低專用圖像重建計算機的成本,或兩者皆可。
            反投射的二元操作是再投射(reprojection),它是計算電子存儲圖像的投射的過程。該過程也在斷層圖像重建中起了基本作用。反投射和再投射的組合也可用來構建對螺旋型錐面光束幾何結構中長對象問題的快速重建算法,這對人體的實際3D成像是關鍵的。此外,在各種應用中,使用迭代重建算法是有利甚至必要的,其中對單個圖像的重構執行反投射和再投射步驟若干次。加快反投射和再投射步驟將確定這些迭代方法的經濟可行性。這些年來已提出了若干種加快重建的方法。例如,Brandt等人的美國專利5,778,038描述了使用多層分解、每個階段產生覆蓋整個視野且分辨率遞增的圖像的2D平行光束斷層造影方法。Nillson等人的美國專利6,151,377揭示了其它分層反投射方法。盡管這些系統都具有價值,但仍然需要產生更準確圖像、提供準確度和速度之間的更大靈活性的方法和裝置。
            因此,本發明的一個目的是提供用于計算機X線斷層造影(CT)掃描的新穎和經改進的方法和裝置。
            另一個目的是提供產生更準確圖像、并提供準確度和速度之間的更大靈活性的、用于CT掃描的新穎和經改進的方法和裝置。

            發明內容
            這些目的由本發明實現或超越。通過反投射選定投射來產生中間圖像、并執行數字圖像坐標變換和/或對中間圖像重新采樣,像素圖像得以根據投射創建。數字圖像坐標變換被選擇用來說明中間圖像的組成投射的視角、及其傅立葉特征,從而中間圖像可通過稀疏樣本來準確地表示。結果的中間圖像被聚合成子集,且該過程以遞歸方式重復,直到已處理和聚合了足夠的投射和中間圖像,以形成像素圖像。
            數字圖像坐標變換可包括旋轉、切變、拉伸、收縮等。重新采樣可包括上采樣、下采樣等。
            通過執行數字圖像坐標變換和/或重新采樣和/或抽選,并重新投射最后的中間圖像,投射可得以從像素圖像創建。


            通過參閱附圖和本發明一實施例的以下描述,本發明的上述和其它特征以及獲得它們的方式將變得顯而易見,且本發明本身將得到最佳的理解,在附圖中圖1是用于本發明的裝置的框圖;圖2A、2B和2C是用于本發明一些實施例中的采樣圖(sampling pattern)的示圖;圖3A、3B和3C是用于本發明一些實施例中的其它采樣圖;圖4是示出一種已知反投射方法的示圖;圖5A是示出本發明一實施例的算法的示圖;圖5B是示出在圖5A的實施例中產生中間圖像的方式的示圖;圖6是示出用于產生圖5A的中間圖像的傅立葉特征的示圖;圖7A、7B和7C是示出當坐標變換是數字圖像旋轉時用于圖5A中示出的反投射算法的中間圖像的傅立葉支持的示圖;圖8是示出用于本發明的另一實施例中的算法的示圖;圖9是示出圖8算法中光譜支持的演變的示圖,其中框(1……9)對應于圖8中的相應各點;圖10A是描述用于本發明實施例的一種算法的示圖;圖10B示出用于圖10A的實施例中的圖像坐標變換;圖11A是示出切變比例反投射的示圖,而圖11B是示出分層切變比例反投射的示圖;圖12A和12B是示出圖像在中間圖像的光譜支持上的切變效果的示圖;圖13是示出本發明另一個實施例的算法的示圖;圖14是找出最優切變因子的算法;圖15示出用于本發明另一實施例的算法;圖16示出用于本發明又一實施例的算法;圖17是示出具有圓形掃描軌跡的通用扇形光束幾何結構的示圖;圖18示出用于本發明再一實施例的算法;圖19示出圖18中所述的算法中中間圖像的加權;
            圖20A、20B和20C示出圖18的第二分層結構層的采樣點;圖21A和21B是用于圖18算法中的采樣圖的示圖;圖22是示出使用圖20A-22C中所示方法獲得的原始交點的示圖;圖23A、23B、23C和23D示出用于圖18中使用的上采樣的旋轉采樣點;圖24示出在由圖18的算法產生的中間圖像的一點上的本地光譜支持;圖25是在圖18的算法中使用的不統一采樣圖的示圖;圖26A和26B示出用于再次采樣和坐標變換的采樣圖;圖27A和27B示出可用于本發明的選擇性采樣方案;圖28包括用于本發明的發散光束的兩個示圖;圖29A是示出錐面光束的示圖,而圖29B是示出再次采樣的示圖;圖30是用于再次采樣投射的算法的示圖;圖31是用于在本發明一實施例中再次采樣的另一種算法;圖32是用于快速分層再投射的算法的示圖;圖33是用于快速分層再投射的另一種算法的示圖;圖34是示出使用本發明的實驗結果的示圖;圖35包括用本發明產生的樣本圖像;圖36A是使用常規算法的重建圖像的顯示,而圖36B顯示使用本發明的快速算法獲得的結果;以及圖37A和37B是與常規算法相類似的算法的點分布函數的示圖,而圖37C顯示本發明快速算法的點分布函數。
            具體實施例方式
            符號和字體以下數學符號和字體系統將用來提高明確度。
            空間域中的函數由小寫字母表示(例如f(x)),而其傅立葉變換由大寫字母表示(F(ω))。
            兩個變量函數的索引取決于方便性可不同地表示。函數f的以下三種表示是等同的f(x1,x2),f(x),以及f(x1x2).]]>連續域和離散域函數分別通過其索引中括號的樣式來區分f(x1,x2)是兩個連續變量的函數(即, ),而f[m1,m2]是f(x)的采樣版本,因此是2-D數組。
            線性算子及其相應矩陣通過字體樣式來區分。假設 是矩陣,則A是其相關聯的線性算子。例如,如果A是坐標變換,則g(x)=(Af)(x)=f(Ax)。
            相同的算子有時在框圖內外會不同地表示。在框圖外它可表示為A(α),在框圖內它表示為Aa。
            硬件概述本發明在包括CT掃描儀的各種成像裝置中都有應用。典型的成像裝置10(圖11)包括掃描儀12,它從諸如頭的對象中采集數據,并向投射預處理器16發送與例如具有不同光束幾何形狀的線積分投射14相對應的原始數據。該投射預處理器16對數據應用各種轉換、歸一化、和校正、以及加權和可以是漂移變化的濾波。投射預處理器16的輸出是預處理后投射的集合,下文中簡稱為投射,也稱為sinogram1(竇腔X線造影)18,它可提供給竇腔X線造影更新處理器20。該竇腔X線造影更新處理器20使用來自sinogram234的信息,可能更改所輸入的sinogram118,例如校正包括射線束硬化的各種人工因素、或作為多步驟或迭代重建過程的一部分。
            竇腔X線造影更新處理器20的輸出是sinogram322,它被輸入給快速反投射處理器24。該快速反投射處理器24通常是經編程和/或接線用以執行這里所述算法的任何適當類型的計算機或專用硬件,或其任意組合。
            該快速反投射處理器24的輸出是電子圖像126,它被輸入給圖像調整處理器28。該圖像調整處理器28執行電子圖像的必要后處理,可能包括非源自腦中的電波圖像、或用于在多步驟或迭代重建過程中作進一步處理的圖像的標識和提取。
            在需要時,圖像調整處理器28可產生提供給快速再投射處理器32的電子圖像230。該快速再投射處理器32通常是經編程和/或接線用以執行這里所述算法的任何適當類型的計算機或專用硬件,或其任意組合。需要時該處理器能共享反投射處理器24采用的同一計算機和硬件的運用。
            快速再投射處理器32的輸出是反饋回竇腔X線造影更新處理器20的sinogram234。反投射/再投射過程可繼續,直到獲得適當結果。盡管再投射并不總是需要的,但在許多情形中都有幫助。
            當電子圖像126適當時,圖像調整處理器28產生提供給存儲/分析/顯示設備38的電子圖像336。預期電子圖像336可被存儲在計算機存儲器中,和/或為例如異常或危險物質對其作電子分析,和/或顯示,和/或用一些可查看形式打印。
            本發明的反投射和再投射方法的概述本發明的反投射方法使用各種技術來創建由像素(圖片元素)和/或體素(3D圖片元素)(下文中統稱為像素)組成的圖像,以下將用一般方法進行介紹。
            本說明書使用常常在多維信號處理中使用的術語和過程,例如1983年Prentice-Hall出版社在Englewood Cliffs出版的D.Dudgeon和R.Mersereau的“Multidimensional Digital Signal Processing”(多維數字信號處理)中所述。本發明說明書中的一些術語用于以下上下文中。術語“采樣圖”指其位置相對于坐標系統限定的空間點集。采樣圖的示例如圖2A-2C以及3A-3C所示。笛卡爾采樣圖指由兩組相互垂直的平行線的交點形成的點集。術語連續圖像指在坐標系統上定義的函數,例如f(x,y)和f(x,y,z)分別是2D和3D函數。數字圖像是采樣圖上連續圖像的各個值的陣列。更廣泛地,連續圖像可由數字陣列表示,這些數字用作參照諸如樣條的一些基集的序列擴展中的系數。零階樣條的笛卡爾積產生用于將數字圖像顯示為持續圖像的熟悉的方形像素形式。在下文中,該數字陣列也被稱為數字圖像。所有存儲在數字計算設備中的圖像必須是數字的。為簡短起見,數字圖像和連續圖像兩者在下文中都經常被簡稱為圖像,其意思可從上下文中推斷。使用該術語,像素圖像是與格子采樣圖相對應的數字圖像,即均勻分隔的周期性圖案、通常并非必然是笛卡爾。
            如果一種采樣圖產生比連續圖像少的樣本總量,則它被稱為比另一種采樣圖稀疏。通常較為稀疏的采樣圖將具有較低的采樣密度。過采樣指使用比以期望準確度表示連續圖像所必需的樣本更多的樣本。相應的數字圖像將稱為被過采樣。給定與一采樣圖的連續圖像相對應的數字圖像,用期望準確度產生與不同采樣圖上的相同連續圖像相對應的新數字圖像的過程稱為數字圖像重新采樣。上采樣和下采樣分別是對較密集或較稀疏采樣圖重新采樣的特定情形。此外,上采樣或下采樣1倍不涉及數字圖像中的變化,被視為是上采樣或下采樣的一種形式。數字圖像附加指在相對于基的擴展情形中,相對于同一采樣圖或同一基定義的數字圖像的逐點附加。低維數字過濾指多維陣列僅在其多個維的子集上的數字過濾,例如2D矩形數字圖像的每個列的單獨1D過濾。
            連續圖像的坐標變換指諸如旋轉、切變、拉伸、或收縮的運算。要定義一數字坐標變換,考慮通過坐標變換關聯的兩個連續圖像、以及表示相對于公共采樣圖的連續圖像的相應數字圖像。從另一個數字圖像產生一個數字圖像的過程稱為數字圖像坐標變換。這可通過數字過濾來實現,即通過對數字圖像的離散指數運算。各特定示例包括(但不限于)數字圖像旋轉、數字圖像切變、數字圖像拉伸、數字圖像收縮、以及這些運算的組合。用于執行數字圖像坐標變換的各種方法是眾所周知的,例如,如IEEE Transactions Image Processing 1995年第4卷1371-1381頁的M.Unser、P.Thevenaz、和L.Yaroslavsky的“Convolution-based interpolation for fast,high-quality rotation ofimages”(用于圖像的快速、高質量旋轉的基于回旋內插)中所述。
            一些數字圖像坐標變換在圖2A-2C和3A-3C中示出。圖2A示出連續圖像(矩形)的輪廓、以及表示連續圖像的數字圖像的采樣圖。較大點上連續圖像的值包括在數字圖像中。圖2-B和2-C分別示出對同一采樣圖的旋轉后連續圖像和切變后連續圖像,其中較大點示出包括在圖2-A數字圖像的旋轉后/切變后數字版本中的值。
            圖3A也示出一連續圖像以及定義數字圖像的采樣圖。圖3B示出在x和y維上拉伸某常數倍后的連續圖像。該數字上拉伸后的圖像由較大點上拉伸后連續圖像的值來定義。圖3C示出與圖3A相同的連續圖像,但其中采樣圖要更密一定拉伸倍。圖3C中的對應的數字圖像將與圖3B中的數字圖像相同。更一般地,對于具有一定規則度的采樣圖,數字圖像拉伸或收縮對數字圖像上采樣或下采樣是等同的。
            注意,某些坐標變換的應用不會使數字圖像改變,諸如旋轉0度、切變參數為零地切變、或拉伸或收縮1倍,因此在過程的描述中可包括在內或略去,而不會改變結果。
            連續圖像的傅立葉變換使人能通過采樣理論確定采樣圖的屬性,從而相應數字圖像將連續圖像表示到期望準確度,如例如1983年Prentice-Hall出版社在Englewood Cliffs出版的D.Dudgeon和R.Mersereau的“Multidimensional DigitalSignal Processing”(多維數字信號處理)中所述。類似地,數字圖像的離散時間傅立葉變換(DTFT)使人能確定什么數字圖像坐標變換將產生將變換后的連續圖像表示到期望準確度的數字圖像。連續圖像的傅立葉變換的相關屬性、以及數字圖像的DTFT將統稱為傅立葉特征。
            加權反投射運算需要對每個投射加權取決于反投射所產生像素的位置的權重。取決于例如獲取投射的源的位置,不同的權重可用于不同的投射,如1988年IEEE出版社在紐約出版的A.C.Kak和M.Slaney的“Principles of ComputerizedTomographic Imaging”(計算機X線斷層造影的成像原理)中所述。作為一特定情形,加權倍數可等于1,它對應于無加權但的確是加權倍數。除非特別指出,經加權和未經加權的反投射將統稱為反投射。
            有了這樣的背景技術信息,將描述本發明的若干實施例。
            反投射處理器28(圖1)被編程為執行用來實踐本發明的算法。這些算法將詳細地進行討論,但首先將以更普通的術語進行描述。反投射過程中的步驟用標號在框圖中示出。
            圖4示出基于旋轉的反投射,作為通過步驟50中以0度反投射各個投射q1...qp而形成的旋轉后中間圖像的和。反投射產生在52進行坐標變換的圖像,并在54聚合以生成圖像 該結構單獨地等同于常規反投射,并在運算計數中不提供縮減。然而,它用作引入本發明的一些快速分層反投射方法的進階。
            圖5A示出用于從多個投射q1...qp中創建像素圖像 的一種分層反投射方法。每個投射qm都在100反投射以產生多個中間圖像I1,1...I1,p這是分層反投射的零層或準備層。在102對選定中間圖像執行數字圖像坐標變換。變換后中間圖像的子集(在此為成對)在104聚合以產生聚合中間圖像I2,I...I2,p/2。這是分層反投射的第一層。第一層的聚合中間圖像用作輸入到分層反投射的下一層的新的中間圖像。在106將數字圖像坐標變換應用于選定中間圖像、以及在108聚合選定中間圖像以產生新的中間圖像的過程繼續,直到已處理完所有中間圖像、并在116聚合成最終圖像 在需要時,可組合跨一些層和這些層內的運算。例如,對于一些投射、以及來自通過在兩個或多個(對于圖5A中的實施例正好為2)選定投射qp的選定視角上的反投射112產生的集I2,1...I2,p/2的相應中間圖像,第0和1層可組合,如圖5B所示。或者,一些初始中間圖像可通過不涉及顯式反投射的等同過程生成,諸如直接傅立葉重建方法,如2001年在費城召開的工業和應用數學學會上發表的F.Natterer和F.Wubbeling的“Mathematical Methods and Image Reconstruction”(數學方法和圖像重建)中所述。
            各個數字圖像坐標變換的參數被選擇,用以說明中間圖像的構成投射的視角,并說明中間圖像的傅立葉特征,從而中間圖像的聚合可由稀疏樣本準確地表示,如以下將解釋的。這些傅立葉特征集中在中間圖像的必要光譜支持上,即傅立葉(頻率)域中的區域,其中中間圖像的傅立葉變換遠不同于零。采樣理論教導連續圖像的光譜支持確定產生表示基礎連續圖像的、以及從中能可靠地重建連續圖像的數字圖像的采樣圖的特性。特別地,圖6示出分層算法中中間圖像 的典型楔形光譜支持。
            圖7(A)和7(B)示出在坐標變換被選為數字圖像旋轉時,圖5(A)所示的二進制分層反投射算法中中間圖像的傅立葉支持。圖7(A)示出虛擬中間圖像I的傅立葉域支持。虛擬圖像 由包括在相應圖像Il,m中的多個投射的反投射組成。圖7(B)示出通過選擇坐標變換的參數來說明中間圖像的構成投射的視角以及中間圖像的傅立葉特征,中間圖像的縱向帶寬(虛線矩形的高度)可最小化,從而允許稀疏采樣并降低計算需求。圖7(C)示出Il,m、Il-1,2m和Il-1,2m-1的空間域采樣方案。采樣點在空間域中虛線的交點處。
            圖8示出本發明的另一個實施例,它執行中間圖像的三重聚合,其中P=2×3L個投射。為了具體說明,示出P=18個投射的情形。象圖5(A)中一樣,投射qθ1,qθ2,…,qθ18在100反投射以產生多個中間圖像I1,1d…I1,18d。這是分層反投射的第零層。
            各投射以3個為一組,且每個組中的選定投射(圖8的3個中的2個)在102上進行數字圖像的坐標變換。變換后中間圖像I1d的子集(圖8中的三元組)在104上聚合,以生成聚合中間圖像I2,1d…I2,6d。這是分層反投射的第一層(1=1)。
            在分層反投射的第二層上(1=2),選定聚合中間圖像I2,md在106進行由沿y坐標拉伸或上采樣組成的坐標變換,而在108對選定的聚合中間圖像進行另一次坐標變換K2,m。在此,變換后的中間圖像在110還是聚合成三個一組,以生成中間圖像I3,1d,I3,2d。來自層2的選定中間圖像在112進行數字圖像坐標變換,并在114聚合以生成下一層的中間圖像I4,1d。取決于投射的數量,該過程重復到在116生成圖像 所必需的程度。在112最后一層上由K標注的數字圖像坐標變換是一數字圖像旋轉,且在前面各層102和108上的坐標變換也可選為數字圖像旋轉。在此,各個數字圖像坐標變換的參數被選為說明中間圖像的構成投射的視角、和中間圖像的傅立葉特征,從而中間圖像的聚合可由稀疏樣本來準確地表示,如以下將解釋的。
            圖9示出圖8的三重反投射算法中光譜支持的進展。編號(1),...,(9)對應于圖8算法框圖中的對應點。所示光譜是數字圖像Il,m的離散時間傅立葉變換(DTFT)。
            圖10(A)描述了本發明另一個實施例的算法,該算法使用三重基于兩切變的分層反投射。圖10(A)的實施例與圖8的實施例相似,且來自圖8的標號也在適用時使用。如在圖8中,示出P=2×3L個投射,且L=2的情形。圖10(a)的實施例與圖8的不同之處有兩個方面。首先,沿x坐標的附加上采樣步驟被包括在101第一層的坐標變換中。其次,在102和108由K標注的數字圖像坐標變換的至少一部分由兩個數字圖像切變運算序列組成,第一個(120)沿y坐標(122),第二個沿x坐標,如圖10(b)所示。
            圖11A示出切變縮放(shear-scale)反投射,而圖11B示出等效的分層的切變縮放反投射。在圖11A中,多個投射q1,…,q4在140反投射,并在142進行切變縮放坐標變換。最后的中間圖像在144聚合。
            在圖11B中,多個投射q1,…,q4在146反投射,在148進行切變縮放坐標變換。并在150聚合成子集。通過聚合產生的中間圖像在152再次進行切變縮放變換,且最后的中間圖像在154聚合。該過程以遞歸方式繼續,直到所有的投射和中間圖像都已得到處理并聚合,用以形成圖像 在此,數字切變變換和中間圖像采樣的參數被選為說明視角和傅立葉特征,以便于降低采樣需求和所需計算量。
            圖12示出圖像切變對中間圖像的光譜支持的效果。圖12A示出當投射應出現在最終圖像中時其某個子集的光譜支持。圖12(B)示出由相同投射組成的中間圖像的光譜支持,其中坐標變換參數被選為使縱向最高角頻率最小化。
            圖13示出本發明又一實施例—三重分層切變縮放反投射的算法。投射qθ1,…,qθ18在160反投射并在162進行切變縮放數字圖像坐標變換。結果圖像的子集在164聚合,且結果中間圖像在166和168進行上采樣,并在168進行數字切變坐標變換。這些圖像的子集在170聚合,且選定結果圖像在172進行其它坐標變換。該過程繼續直到所有投射和中間圖像已處理并在174聚合,生成圖像 在172所示的最終坐標變換K-π/2僅涉及當像素圖像 的采樣圖是笛卡爾時像素的重新排列。
            圖14是使用圖12所示中間圖像的傅立葉屬性的屬性,尋找先前所述實施例的坐標變換參數的一種算法。
            圖15示出本發明的再一實施例—過采樣三重兩切變分層反投射的算法。圖15的實施例與圖10A的實施例相似,且來自圖10A的標號在適用處使用。然而,除了圖10A的各個步驟之外,圖15的實施例包括最后步驟之前一個步驟中的下采樣步驟109,它在圖15中為第二層。在此,各個數字圖像坐標變換的參數被選為說明中間圖像的構成投射的視角和中間圖像的傅立葉特征,從而中間圖像的聚合可由稀疏樣本來準確表示。然而,一定程度的過采樣用來改進隨后處理的準確度,如將作解釋的。如果需要,則為了改進計算效率,下采樣步驟109可與第二個x軸切變組合,使過程108和109一起被圖15(B)中所示的過程所替換,其中第二個x軸切變包括產生切變縮放變換的兩切變數字圖像坐標變換108(如圖10B所示)。
            圖16示出本發明的又一實施例—過采樣三重分層的切變縮放反投射的算法。圖16的實施例與圖13的實施例相似,且來自圖13的標號在適用處使用。然而,圖16包括中間層中上采樣161的附加步驟,以及在最后層之前一層中的下采樣169,其中中間圖像被分別上采樣和下采樣。與圖15的實施例相似,在圖16的實施例中,當沿x坐標的下采樣步驟169↓xUl,m在x切變步驟168 Sxαl,m之后時,為了計算效率兩者可組合成單個數字圖像切變縮放。
            用于快速分層反投射算法的非笛卡爾采樣方案圖8、10和13中所示的分層反投射的不同實施例中的中間圖像具有便于有效非笛卡爾采樣的特殊光譜支持。特別地,第1層中的基礎連續域圖像Il,m在傅立葉空間中占據一楔形,如例如在圖6、7、9和12中可見。多維采樣理論告訴我們,對于具有這種光譜支持的圖像,在笛卡爾柵格上的采樣不如在適當非笛卡爾柵格上的采樣有效。非笛卡爾采樣可通過將基帶頻譜的副本更緊地封裝于傅立葉平面而改進采樣效率。對于2D采樣的解釋參見[?]。特別地,梅花形采樣方案減少中間圖像的大小,因此使該算法的計算成本降低幾乎2倍。對周期性的非笛卡爾采樣圖的數字圖像坐標變換可使用一維移位不變濾波器來有效地執行。類似地,所有先前描述的用于選擇數字圖像變換的參數的方法也在這種采樣圖情形中應用。因此,前述所有實施例擴展到使用周期性非笛卡爾采樣圖的各個實施例。
            已經描述的本發明各個實施例的變體可應用于各種3D投射形式的3D反投射,例如正電子斷層掃描儀中產生的X線變換;以及在磁振造影中產生的3D Radon變換。
            三維(3D)X線變換是3D中沿平行線集在各個方位上的積分的集合。每個3D X線投射是可由線定向的3D角度表征的兩維函數。3D X線數據的分層反投射的框圖與前述的相似,諸如圖8、10、13、15或16。該情形中的中間圖像是三維的,而非前述示例中的二維。每個中間圖像在各構成投射的3D角度的平均值方向上最稀疏的3D采樣圖上采樣。隨著該算法的展開,沿該稀疏采樣(緩慢)方向上樣本的密度增大,以適應該方向上的漸增帶寬,如3D X線投射的傅立葉分析所述。因此在該算法的每個階段,在聚合圖像之前每個圖像都必須沿該緩慢方向上采樣。3D實施例中可用的額外維數還提供可用于算法中各個步驟的更多坐標變換,諸如3D中的旋轉。像2D情形中一樣,這些數字圖像坐標變換的參數可選為說明構成視角和中間圖像的傅立葉特征。這些數字圖像坐標變換可被分解成一維運算序列,諸如如前所述的切變和切變縮放。像2D情形中一樣,可迫使三維的任何子集中的過采樣,以改進圖像質量。
            3D Radon變換投射是一維函數——沿平行2D平面集的以平面離原點的距離為參數的積分的集合。投射視角是與平面集垂直的矢量的3D方位角的視角。3DRadon投射的分層反投射的框圖象圖8、10、13、15或16中一樣。在第一層中多個投射被Radon反投射到3D圖像域上。這些圖像是沿兩維恒定的,且因此僅需要在垂直于該兩個恒定方向的方向上采樣。當多個組的3D Radon投射組合時,聚合圖像的帶寬可取決于這些構成投射的視角在一個或兩個維上增加。因此,有必要在坐標變換中間圖像并添加到其它中間圖像之前,在可能兩個維中上采樣該中間圖像。像2D情形中一樣,坐標變換可單獨執行、可與上采樣運算組合、并且過采樣可在中間圖像上實施。
            在本發明的各個實施例中,數字圖像坐標變換、以及下采樣或上采樣運算可通過一系列較低(一)維線性數字過濾運算來執行。此外,當所使用的采樣圖具有一些周期性時,這些數字濾波器可以是移位不變的,如將要更詳細描述的。為了計算效率,數字濾波器可使用遞歸濾波器結構、或使用眾所周知的FFT來實現。一種確定數字濾波器的較佳值的方法是使用樣條內插理論,如IEEE Transactions onPattern Analysis and Machine Intelligence 1991年第13卷第277-285頁的M.Unser、A.Aldroubi和M.Eden的“Fast B-Spline transforms for continuous imagerepresentation and interpolation”(用于連續圖像表示和內插的快速B樣條變換);IEEE Transactions Signal Processing 1993年第41卷第821-832頁的M.Unser、A.Aldroubi和M.Eden的“B-spline signal processingPart I-theory”(B樣條信號處理部分I-理論);以及IEEE Transactions Signal Processing 1993年第41卷第834-848頁的M.Unser、A.Aldroubi和M.Eden的“B-spline signal processingPart II-efficientdesign and applications”(B樣條信號處理部分II-有效設計和應用)中所述。
            本發明的分層反投射方法可應用于廣泛范圍的X線斷層成像幾何結構,諸如具有任意源軌跡的包括扇形光束和錐面光束的發散光束幾何結構。具有圓形掃描軌跡的通用扇形光束幾何結構在圖17中示出。射線源繞半徑為R的圖像在半徑為D的圓形掃描軌跡上移動。源角度β上的扇形光束投射對應于沿以扇形角度γ為參數的射線集的線積分度量。TST是沿源射線到圖像盤的最近邊緣的距離,TEND是沿源射線到圖像盤的最遠邊緣的距離。
            圖18示出用于本發明一實施例的一種算法——分層的三重基于旋轉反投射,可應用于扇形光束加權反投射。圖18的實施例與圖8的實施例相似,且圖8的標號在適用時使用。圖18的實施例與圖8所示實施例的不同之處在于若干方面。第一,示出P=4×3L個投射,且L=2的情形。第二,在第0層上,初始中間圖像Il,md通過在99的標示為W的加權反投射從扇形光束投射中生成。第三,在最后階段,聚合4個而不是2個中間圖像。第四,在分層結構的大多數較早層中使用的采樣圖最好選為非笛卡爾,如將要解釋的。在此,如果笛卡爾采樣圖用于倒數第二層中的中間圖像I3,md,則最后的數字圖像坐標變換同樣僅需要重新排列圖像像素。此外,像圖5(B)中一樣,分層結構的第0層和第1層中的運算可組合,從而中間圖像I2,md可通過其三個構成投射的直接加權扇形光束的反投射、或通過其它手段生成。
            在選定數量的層中,通過將附加加權步驟包括在上采樣步驟和數字圖像旋轉(圖18中的步驟106和108)的級聯前后(如圖19所示),來更改圖18的實施例是有利的。中間圖像I2,m分別在步驟106和108前后在180和182通過在空間上更改權重來加權。如將要解釋的,該加權可用來減少中間圖像的采樣需求,從而減少計算量。
            中間圖像的采樣方案影響算法的性能。期望采樣方案是使用最少樣本但不丟失圖像信息的方案。在此,各個數字圖像坐標變換和重新采樣運算的參數可選為說明中間圖像的構成投射的視角、以及中間圖像的傅立葉特征,從而中間圖像的聚合可通過稀疏樣本來準確地表示,如將要解釋的。如將要解釋的,用于選擇這些參數的選擇性方法基于特殊射線或曲線集的交點。
            圖20A-20C示出在源角度以間隔Δβ均勻隔開的情形中,中間圖像通過圖18所示三重分層結構算法的各個層次的進度。示出了算法中組成中間圖像的多個投射的扇形。在零層上,每個中間圖像I1,md由具有方位為β=0的扇形的單個投射組成,如圖20(A)所示。在遞歸的第一層之后,每個中間圖像由三個投射組成,其中扇形由圖20(B)所示。在第二層之后,每個圖像由9個反投射扇形組成,如圖20(C)所示。
            算法中用于對具有如圖20(C)所示構成扇形的中間圖像I3,m選擇坐標變換和采樣/重新采樣圖的基于交點方法如圖21A和21B所示。I3,m的采樣圖由位于中央構成扇形的點組成,與圖20(A)中所示的一致。如圖21(A)和21(B)所示,I3,m的采樣點由中央扇形的一半與中央扇形兩側上的極值構成扇形的交點所確定。
            通過應用兩個約束條件來更改結果采樣圖,以改進反投射的準確度并減少計算需求是有利的。首先,沿射線的樣本密度受到限制,不超過最終圖像的期望采樣密度。其次,采樣圖被迫使在每條射線上包含每一側的圖像盤外的至少一個樣本。圖22中的曲線顯示采樣點沿中央扇形的扇形角度為γr的特定射線上的位置。采樣點位于距離γ’的固定間隔上,其中γ’是圖21(a)、(b)和(c)所示極值扇形之一的扇形角度。落在圖22中連續曲線上的各個點是使用圖21A和21B中所示的方法獲得的原始交點。虛曲線上的各個點是使用以上兩個約束條件更改的。圖23D示出使用該更改交點方法獲得的采樣圖的一個示例。所示扇形是已為之生成該采樣圖的中間圖像的中央扇形。
            當在本發明前述實施例的情形中時,將重新采樣和數字圖像旋轉運算分解成一維運算系列是有利的。圖18中標為↑U的框表示將對中央扇形采樣的中間圖像上采樣成同一扇形上采樣點的較小集。這可通過沿扇形的每條射線的單獨上采樣來實現。對于圖18中標有↑U和K的上采樣和旋轉的級聯,聯合分解成1D運算是方便的,如圖23(A)-23(D)所示。在4個畫面的每一個中,兩個虛線圓表示半徑為D的圖像的邊界和半徑為R的(更大)圓上的源軌跡。采樣點由小圓圈表示。具有圖23(A)所示的中央扇形(Fana)和采樣圖的數字中間圖像要重新采樣、并旋轉到圖23(D)所示中央扇形(Fand)的角度,具有圖23(d)所示的最終采樣圖。這用以下兩個步驟來完成(i)分別沿著Fana的每條射線將圖像上采樣成圖23(b)中所示的采樣圖,它由Fana和Fand的交點定義。這些點因此將位于圖23(c)所示的Fand上;(ii)分別沿著Fand的每條射線將圖像上采樣成圖23(d)中所示的采樣圖。
            這些步驟使用1D上采樣運算來實現重新采樣和旋轉的組合運算。
            在對周期性采樣圖采樣的中間圖像情形中,為本發明各實施例描述的傅立葉分析技術被擴展成設計分散光束情形中的空間變化采樣方案。這些技術足夠普遍以應用于兩維和三維中任意軌跡上的投射和反投射。對于非周期性直線或曲線系統上的反投射情形,局部光譜支持概念代替光譜支持概念。這在圖24中示為由圖24左側示出的單個扇區的反投射所產生的中間圖像。如以下將解釋的,在所示點上以r和θ為參數的連續中間圖像的局部光譜是圖24(b)右側所示傅立葉域中的線段。由多個扇形光束投射的反投射組成的中間圖像的點上的局部光譜支持如圖25所示。在左側,示出了點的位置、以及用于中間圖像的構成投射的視角范圍。右側的蝶形領結狀區域是相應的局部光譜支持。
            該局部光譜支持分析用來確定中間圖像的局部采樣需求。結果在空間上不均勻的采樣圖由圖26A和26B中的虛線弧表示,用于兩個典型的中間扇形反投射圖像。圖像邊界由虛線圓表示,且僅有極少的點需要在該邊界外取得。該局部傅立葉采樣方法可直接應用于在直線、曲線或任意維的平面上尋找任意投射幾何結構的采樣方案。
            確定用于重新采樣的采樣圖、以及用于分層扇形光束反投射的坐標變換的基于傅立葉的方法可進一步地擴展到位于除中間圖像的中央扇形之外的直線或曲線上的采樣圖,如下所述。該有用集的示例如圖27(A)和27(B)所示。
            普通發散光束算法所述用于扇形光束反投射的方法直接擴展到其它發散光束幾何結構,包括現代診斷放射醫學中最重要內容之一的3D錐形光束,如下所述。類似于扇形光束幾何結構,其中發散射線的源在圍繞成像對象的軌跡(即圓)上移動,在普通發散光束的幾何結構中發散射線的源在圍繞成像對象的3D空間中的軌跡(并非必須為圓)上移動。檢測器表面通過成像對象測量線積分。
            用于普通發散光束幾何結構的分層反投射算法的一個實施例可再次通過類似于圖18的框圖描述,但用以下方法更改。首先,在第0層,發散光束投射在99用與名義“零”源位置相對應的適當發散光束的單視圈反投射W1來零反投射,從而產生初始中間圖像。其次,因為源的軌跡并非必須是圓,中間圖像的構成發散光束不僅可像在扇形光束幾何結構一樣旋轉,而且還能彼此平移。相應地選擇18中的坐標變換K。再次,取決于源軌跡和位置中的對稱性的出現,可以有或沒有“自由”的坐標變換,諸如替代扇形光束算法中π/2旋轉的像素重新排列。像扇形光束情形一樣,初始中間圖像由該算法分層處理。與扇形光束情形相類似地,在位置和方位上彼此接近的中間圖像被聚合,以便于聚合后的中間圖像可稀疏地采樣。通過研究旋轉和相對彼此平移的構成加權反投射發散光束的結構,中間圖像的3D采樣圖得以確定。如圖28所示的一種采樣圖,是樣本沿與中間圖像的中央構成光束相對應的發散光束的射線分布的。沿每條射線的樣本間隔被選擇為確保該中間圖像的所有構成發散光束都足以沿該射線采樣。或者,更普通的方法是使用局部傅立葉方法來尋找采樣圖,如前所述的扇形光束情形。知道了中間圖像如何由其構成投射組成,每個中間圖像的局部3D傅立葉結構得以確定。中間圖像的每一點上的3D局部采樣矩陣函數被發現與局部傅立葉支持的采樣需求相匹配,如扇形光束情形中所述。該矩陣函數然后在圖像域上(可能數字地)積分,以確定樣本的位置。
            組合旋轉和平移地上采樣成新的采樣光束的可分開方法類似于扇形光束情形實現。它將3D坐標變換和重新采樣運算減少成1D重新采樣運算序列。如圖29(a)所示,每個發散光束被視為按照方位角分布的一組垂直坡面扇形光束、與在不同仰角上分布的傾斜坡面扇形光束集的交點。可分開坐標變換的步驟如下(如圖29(b)中所示)1.原始發散光束被重新采樣成原始射線與新發散光束的垂直平面的交點集。這些點因此位于最終采樣點所共享的平面上。
            2.來自扇形光束情形中可分開重新采樣的步驟1和2對每個平面單獨執行,以重新采樣成最終點集。
            使用適當的有效采樣方案,用于發散光束的快速分層反投射算法可期望實現具有期望準確度的較大加速。如在扇形光束情形中可能使用在連續層中更改的偽光束(例如使用移到距離原點更遠的頂點位置)。
            對本領域技術人員顯而易見的是,本發明諸方法并不限于在此所述的成像幾何結構或特定實施例的示例。這些方法同樣可應用于具有其它幾何結構的反投射的一般問題。
            圖30示出基于重新采樣的反投射,作為通過源位置βp上的各個投射的一般化反投射所形成的上采樣后中間圖像的和。它與基于旋轉的反投射4相似,但不同之處在于兩個方面。第一,在第一步驟中(可能是處理后的)第p個投射在源位置或方位βp進行加權反投射184,而不是4情形中的零反投射。例如在扇形光束投射的情形中,每個投射是在其源角度方位上反投射的。其次,在188通過累加聚合之前,每個初始中間圖像在186進行上采樣運算。這是必須的,因為這些P個單投射中間圖像中的每一個的采樣圖被選為有效和稀疏的圖,所以通常對每個投射是不同的,且常常是非笛卡爾的。在中間圖像聚合之前,它們需要被重新采樣成普通和更密集的柵格。該結構本身并不提供運算計數的減少。然而,它用作引入本發明的一些快速分層反投射方法的進階。
            圖31示出本發明的另一個實施例——用于從多個投射q1...qp中創建像素圖像 的基于重新采樣的分層反投射。圖31是圖30的二進制分層版本。首先,與非分層情形中一樣,每個投射在其各個方位上反投射184成適于該方位的采樣圖,以產生中間數字圖像。這是分層結構的第0層。在分層結構的第一層中,這些中間數字圖像在186被上采樣成圖像選定對共用的一種更密集采樣圖。該上采樣通常將是非笛卡爾采樣圖,但可通過一維重新采樣運算序列執行,如圖23或圖29(b)所示。選定的結果上采樣圖像在190成對地聚合,以產生新的中間圖像。在第二層中,新的中間圖像在192再次上采樣并在194聚合,以產生新的中間圖像。該過程繼續直到所有中間圖像和投射都已得到處理,以在最后聚合步驟198之后產生最終圖像 像前述實施例一樣,運算可在層內或跨層組合。
            分層結構中每個階段上的采樣圖和圖31所示實施例中的上采樣運算的參數可通過任一前述方法選擇。例如,在扇形光束投射的情形中,對給定中間圖像的一種可能采樣圖將位于兩個扇形的交點上第一個扇形朝向中央構成源位置;第二個扇形朝向極值構成源位置。或者上采樣步驟的采樣圖和參數可基于傅立葉或局部傅立葉屬性、以及包括在中間圖像中的投射視角來確定。特別地,對于非周期性的采樣圖,所述用于扇形光束的基于旋轉算法的局部傅立葉方法可用來尋找采樣圖知道投射如何組合以形成中間圖像導致局部光譜支持的確定,該局部光譜支持用來計算在圖像域上積分時產生中間圖像的采樣圖的局部采樣矩陣函數。
            圖32是快速分層重新投射的框圖。重新投射是從給定電子圖像計算X線斷層投射的過程。重新投射算法通過將流程圖換位的過程應用于反投射算法的任何框圖來發現,其中在加權運算中可能有一些改變。在該過程中,運算被其共軛或二元運算所替代。因此重新投射的框圖顯現為與反投射的相應框圖的反向版本類似。圖32中所述的重新投射過程是從圖8所示反投射算法中獲取的一種這樣的重新投射實施例。
            在第一層中,輸入圖像f的副本保留在示圖的頂部分支中作為頂部中間圖像,并且在底部分支圖像f在200旋轉-π/2,以產生底部中間圖像。在第二層中,在示圖的上半部分,未旋轉的中間圖像在202進行三次獨立的數字圖像坐標變換,其中一部分不會使圖像改變,從而產生三個不同的頂部中間圖像。類似的過程應用于底部分支,從而產生三個底部中間圖像。頂部和底部中間圖像中的每一個(一共6個)然后在204進行抽選過程(低通濾波然后是下采樣),以產生新的中間圖像。在圖32所示實施例的實例中,僅有2·3L(其中L=2)個視角,所以第三層是遞歸分層結構中的最后一層。在第三和最后一層上,中間圖像進行單獨的坐標變換(206),其中一部分不會使圖像改變,從而產生18個中間圖像。不是遞歸一部分的最后一個步驟是不同的每個中間圖像在208進行零度的重新投射。零度重新投射等同于累加垂直像素行以產生一維信號。這些一維信號(518)是該算法產生的輸出投射。
            算法中數字圖像坐標變換的參數根據中間圖像的傅立葉特征知識來選擇。這些參數僅相關于相應反投射算法的參數。易于看到,因為重新投射框圖是反投射框圖的流程圖換位,重新投射框圖的每個分支在反投射框圖中都有相應的分支,且相應分支中的坐標變換是彼此的數學共軛。在與圖10雙切變反投射算法相對應的該重新投射算法的版本中,第二層的坐標變換(202)是x切變加y切變,而最后一層的坐標變換(206)是x切變加y切變,以及x中的部分抽選(這三個運算可減為切變縮放)。這些切變參數是用于反投射算法的相應參數的負。x中抽選的參數與反投射算法的第一層中x上采樣的相同。在本算法的切變縮放版本中(對應于圖13所示的切變縮放反投射算法),第二層的坐標變換(202)是x中的切變(且每個切變的參數是圖13中相應參數的負)。最后一層中的坐標變換(206)是切變縮放。在重新投射算法中使用的切變縮放是在反投射中使用的相應切變縮放的數學共軛。抽選因子的參數也與反投射中相應的上采樣因子相同。
            與反投射算法相似,算法中的中間圖像的過采樣可通過在算法開始時首先上采樣圖像并在算法結束時下采樣相同因數來實施。此外,這些運算可在層內或跨層組合。
            圖33是基于抽選的加權重新投射的框圖。它是圖31的流程圖換位。它顯示將圖像重新投射成18個不同源角度上的投射集。開始時,給定圖像沿兩條平行路徑進行處理。在每條路徑中圖像進行三次平行重新采樣(210),成為三個不同的更稀疏采樣圖。使用一維抽選(低通濾波、隨后是重新采樣成更稀疏柵格),該重新采樣能以可分開方法執行。濾波器的參數根據中間圖像的傅立葉特征和期望投射來確定。期望投射的局部傅立葉分析用于投射不排列于平行直線上的情形中。在算法的最終層中,每個中間圖像再次進行三次平行重新采樣(212),成為更稀疏的采樣圖。最后對圖像執行加權投射以產生投射pθ。這涉及要產生一維投射的圖像像素的加權和。
            實現和實驗結果本發明的較佳實施例用C編程語言實現,并用數字實驗在具有384MB RAM的Sun Ultra 5工作站上測試。該測試圖像是Shepp-Logan頭部模型(head phantom)(用于X線斷層造影算法的數字估算中的標準測試圖像)。通過改變該算法的參數,可在準確度和計算成本之間找到平衡點。準確度指重建圖像的質量。盡管視覺質量并不易于計量,所以我們測量重建圖像和原始圖像之間的誤差,從中可用數字計算Radon變換。所使用誤差的測量是相對均方根誤差(RRMSE)。從f[m2,m1]的X線斷層投射重建N×N圖像 的RRMSE如下進行計算RRMSE=1N2Σm2=1NΣm1=1N|f^[m1,m2]-f[m1,m2]|2maxm1,m2f[m1,m2]---(1)]]>對于平行光束數據,測試圖像大小為256×256,視角數量為486,且Shepp-Logan濾波器(理想V型濾波器乘以正弦函數)用來濾波各個投射。圖30顯示在過采樣參數γ為0.75~1.0之間的各個值上兩種算法的RRMSE誤差對運行時間。雙切變算法由圓圈表示,而切變縮放算法由方塊表示。每種算法使用兩種類型的濾波器——稱為MOMS 16的三階(虛線)和五階(實線)樣條濾波器——來運行。對于該算法的每種特性,當γ減小時算法的誤差降低而運行時間增加。未與任何其它點相連的曲線點是由相連點表示的算法的非過采樣版本。作為比較,使用線性內插的常規算法的運行時間為14秒,且其輸出圖像的RRMSE誤差為0.0486(比在此顯示的快速算法差)。
            來自用于平行光束數據的算法輸出的一些樣本圖像,在圖35中顯示。從左到右按列地,它們是來自常規反投射、雙切變、和切變縮放算法的輸出圖像。γ=0.82的過采樣應用于兩種快速算法。下面一行圖像詳細地顯示了上面一行相應圖像的一部分。
            用于扇形光束幾何結構的快速分層算法在512×512 2D Shepp-Logan模型上成功測試。所考慮的獲得幾何結構具有源半徑D=1.06×N=544,972個源角度,以及1025個等角檢測器。快速算法的常規和過采樣版本使用各種內插方法來實現。結果的重建圖像與使用線性內插的常規扇形光束算法作比較。在所有實驗中,投射都用Shepp-Logan濾波器濾波。
            圖36(A)和36(B)顯示重建圖像,圖36(A)來自常規算法,而圖36(B)來自快速算法。快速算法的結果可與常規算法的作比較。快速算法的點分布函數可與常規的作比較。圖37(B)顯示常規算法的PSF,圖37(C)顯示使用線性內插和γ=0.4過采樣的快速算法的PSF。圖37(c)通過常規算法、非過采樣的和γ=0.4的過采樣快速算法的psf的x軸比較各樣條。PSF的相似性確認快速算法中相當的圖像質量。
            所使用的采樣方案是基于扇形交點的方法,無需下文中提議的增強。除了先前提及的本方案缺點之外,為了簡單實現,算法準確運算所不需要的許多采樣點用于實驗中測試的實施例中。盡管這些低效性,對于N=512和D=1.41N,常規算法中(數據相關)乘法運算與快速(線性)算法的比例為6.4,且加法運算的比例為3.0。注意,在該實驗中使用的具有極接近原點的源(D=2.06R)的幾何結構,因為接近源頂點的樣本的高密度所以對于該未增強實現特別具有挑戰性。在更實用的系統中源更遠一些,減少了這些效應。此外,使用前述選擇性采樣方案將導致高得多的加速因數。
            用于本發明的反投射和重新投射算法的詳細描述反投射是用來從已處理投射中創建圖像的過程。重新投射是逆過程,用來從給定圖像計算投射。兩種運算都用于根據投射的圖像重建,如圖1所示。首先將描述常規反投射,隨后描述本發明的反投射和重新投射方法。將首先描述對2D圖像的平行光束投射,然后描述更一般的2D和3D幾何結構。
            用來從其投射估算圖像的典型和較佳方法是濾波反投射(FBP)算法。FBP首先涉及用所謂V型或更改后V型濾波器來濾波投射,然后反投射這些過濾后的投射。附加預處理步驟也可在反投射之前應用于該投射。在視角θp(p=1,...,P)上取得的處理后投射是投射索引,取決于對變量t或θ的相關性是否需要顯式地示出,將被交替地表示為q(t,p)或qp或qθr。為簡單起見,處理后的投射在下文中將稱為投射。
            從一組P投射(在由長度為P的矢量 的分量指定的角度上)的FBP重建 被示為f^(x1,x2)=(Bθ→q)(x1,x2)---(2)]]>其中定義0.1反投射運算符由以下定義(Bθ→q)(x1,x2)=ΔΣp=1Pq(x1cosθp+x2sinθp,p)---(3)]]>本發明的基于O(N2logP)坐標變換的分層反投射算法基于反投射按照坐標變換的分層分解。若干較佳實施例的細節用坐標變換的不同選擇來描述,以最一般的坐標變換結束。
            特別地,本發明的基于O(N2logP)旋轉的分層反投射算法基于反投射按照圖像旋轉的分解。
            將圖像 映射成其旋轉版本(K(θ)f) 的旋轉θ運算符由以下定義定義0.2(κ(θ)f)(x→)=Δf(Kθx→)]]>其中定義0.3在平面上旋轉角度θ的矩陣為Kθ=cosθ-sinθsinθcosθ]]>1D函數q(t)在單個角度θ上的反投射被定義為(Bθq)(x→)=Δπq(x1cosθ+x2sinθ).]]>特別地,零角度反投射(對于θ=0)為(B0q)(x→)=q(x1cos0+x2sin0)=q(x1)=q(10x→)---(4)]]>定義0.4單個θp反投射的中間圖像是Iθp(x→)=Δ(BθPqp)(x→)]]>如圖4可見,方程式(3)中的反投射可重寫如下。在此qp指其視角為θp的投射。
            f^(x→)=Σp=1P(κ(-θp)B0qp)(x→)---(5)]]>
            分層反投射的一個實施例基于若干連續旋轉的累積結果仍然是旋轉的事實。特別地,κ(θ1)κ(θ2)...κ(θN)=κ(θ1+θ2+...θN)(6)從方程式(6)可得出,對于某整數L的p=2L,圖4的框圖可被重新安排成圖5所示的分層樹結構。第1層的第m個分支中的中間圖像被示為Il,m。在初始層中,I1,m=ΔB0qmm=1,2,...,P]]>(7)中間圖像通過反投射個別投射來產生,并且在后續層中,即,對于l=1,2,3,...,log2PIl+1,m=Δκ(δl,2m-1)Il,2m-1+κ(δl,2m)Il,2m]]>(8)中間圖像進行坐標變換(在該實施例中為旋轉)并通過累加聚合。對于任意組的投射角θi,i=1,...,P,中間旋轉角δl,m。可選為確保圖4和5中的結構等效,從而圖5所示的分層反投射算法產生期望圖像 實際上,因為有比視角(約束)多得多的中間旋轉角(自由參數),所以選擇δl,m時有許多自由度。此外,對實踐中使用的數字圖像使用數字圖像坐標變換,其準確度和計算成本取決于采樣圖和實現的選擇。這些各種自由度或參數在本發明中用于降低分層反投射算法的計算需求,如以下所述。
            定義0.5中間圖像Il,m的構成投射的指數集Nl,m是有從圖5中輸入到中間圖像Il,m的路徑的投射的指數集。例如,如圖5中可見,N3,P/4={P-3,P-2,P-1,P}。使用方程式(8)或在檢查框圖之后,便于顯示Nl,m={2l-1(m-1)+k∶k=1,2,...,2l-1} (9)因而Nl,m列示組成中間圖像Il,m的投射。如果與由圖5所示的方法處理相反,這些由集Nl,m指示的相同投射在其相應視角上一起直接反投射,則這將產生虛擬中間圖像 定義0.6I~l,m=Σp∈Nl,mIθp]]>在檢查圖5中框圖之后,可看到中間圖像中的投射之間的相對角度在最終的重建圖像 中得到保留。換言之,對于算法中的所有中間圖像,即l,m,Il,m=κ(αl,m)I‾l,m]]>對于某αl,m∈(-π ,π](10)分層算法的中間旋角δl,m完全根據αl,m如下確定。從Nl,m的定義或方程式(9)中可簡便地得出Nl+1,m=Nl,2m-1UlNl,2m,且從而通過定義I得出
            I~l+1,m=I‾l,2m-1+I~l,2m---(11)]]>方程式(8)、(10)和(11)隱含δl,2m-1=αl+1,m-αl,2m-1和δl,2m=αl+1,m-αl,2m(12)中間旋轉角(δl,m)和因此相對角αl,m被選擇成用以減少中間圖像的帶寬。較小帶寬的圖像可由少數樣本表示,并因此降低整個算法的計算成本。通過理解傅立葉域中的算法、并跟蹤中間圖像的光譜支持通過分層結構的進展,這些中間圖像的帶寬得以確定。
            傅立葉域中的X線斷層造影投射數字圖像坐標變換的參數被選為說明選定投射的視角θp,如前所述。參數還為它們對中間圖像的傅立葉特征的效應進行選擇。這些考慮也被用來確定選擇將哪些中間圖像聚合在一起。
            解釋傅立葉域中的反投射算法的關鍵是使投射P0f的一維傅立葉變換(F1)與圖像f的兩維傅立葉變換(F2)相關的投射-片(projection-slice)定理。投射-片定理4稱(F1P0f)(ω)=(F2f)(ωcosθ,ωsinθ)圖4的反投射算法可在傅立葉域內解釋。零反投射的運算符產生其光譜支持受限于水平頻率軸ω1的圖像(B0q)(x1,x2)=q(x1)(F2B0q)(ω1,ω2)=2πQ(ω1)δ(ω2)(13)眾所周知,在空間域中旋轉一圖像導致其傅立葉變換也旋轉相同角度。
            g(x→)=(B0q)(Kθx→)⇔G(ω→)=(F2B0q)(Kθω→)---(14)]]>可得出方程式(5)中的反投射運算符的功能是將每個投射的光譜分量旋轉到重建圖像的適當角度,并將它們累加在一起。
            假設具有單側帶寬(t變量中)等于π的投射,分層算法中典型的虛擬連續中間圖像Il,m因此在由構成投射確定的角度范圍內具有楔形的光譜支持,如圖6所示。特別地,設Il,m為圖5(A)中所示算法的框圖中的中間圖像,其中Nl,m={b,b+l,...,e},即該圖像的構成投射的視角為{θp∶p=b,b+l,...,e}。因為方程式(10)中的關系,Il,m的光譜支持僅是Il,m的具有旋轉角αl,m的旋轉版本。可得出,圖6中楔形占據的角度的范圍為[γ-β,γ+β]=[θb-αl,m,θe-αl,m]。分別在第一和第二坐標系中Il,m的帶寬Ωs和Ωf在楔形的外圍矩形上示出,如圖6中的虛線所示。因而,αl,m的選擇提供了Il,m帶寬的控制方法。
            采樣理論指示這種連續中間圖像Il,m如何由笛卡爾采樣圖上的樣本表示。特別地,第一和第二坐標系中具有帶寬Ωs和Ωf的連續圖像可在與坐標軸對齊的笛卡爾圖上采樣,其中第一和第二坐標系中采樣間隔為Δs和Δf,從而Δs<π/Ωs和Δf<π/Ωf。為了使表示連續圖像所需的樣本數量最小化,應當最小化與外圍矩形的面積成比例的數量1/(ΔsΔf)>π2ΩsΩf。使該區域最小化、并因此使對中間圖像Il,m的采樣需求最小化的旋轉角為αl,m*=(θb+θe)/2---(15)]]>選擇用于坐標變換的參數方程式(15)用于方程式(12)中,以在圖5所示的基于旋轉的分層反投射實施例中確定用于數字圖像坐標旋轉的最優旋轉參數δl,m。指數b和e被確定為分別對應于集Nl,m中的最小值和最大值。有了這種選擇,中間圖像Il,m=K(αl,m*)I‾l,m]]>及其在第一坐標系中的帶寬(慢帶寬)和第二坐標系中的帶寬(快帶寬)分別為Ωs(Il,m)=πsinθe-θb2Ωf(Il,m)=π---(16)]]>為了使帶寬進一步最小化,應當最小化差值θe-θb,這通過選擇在每個步驟上聚合從其視角之間具有最小的最大角距的投射產生的中間圖像來實現。
            除了旋轉角和對要聚合哪些中間圖像的選擇以外,有必要選擇算法中用于各個中間圖像的數字圖像坐標變換的適當采樣圖。這些使用采樣要求進行選擇,而采樣要求又使用中間圖像的光譜支持和帶寬來確定。
            特別地,對于通過單個投射的反投射形成的初始中間圖像,如圖5中的Il,m=B0qθm,慢帶寬Ωs(Il,m)=0,且沿x2坐標的單個樣本充足。慢帶寬將更大,且沿x2坐標的更多樣本將為如圖3-b通過反投射一個以上投射形成的初始中間圖像所需。
            兩個虛擬中間圖像之和(I~l,m=I‾l-1,2m-1+I~l-1,2m)]]>及其相應中間圖像如圖7(A)-7(C)在空間域和頻域中示出。圖7(A)和7(B)示出I和I的傅立葉域支持。圖7(C)示出Il,m、Il-1.2m、Il-1,2m-1的空間域采樣方案。各采樣點在空間域中的水平線和垂直線的交點上。聚合后的中間圖像Il,m具有比其每個分量更大的慢帶寬,需要更密集的采樣圖。相反,在分層算法較早層上的中間圖像需要更稀疏的采樣圖,從而節約計算。
            計算成本和節約計算成本是易于估算的。假設在
            ]>=N2sin(π(2l-1-1)/(2P))]]><N2(π2l-1/(2P))=O(N22l/P).]]>假設旋轉S個像素的數字圖像的成本是O(S)且size(Il,md)=O(N22l/P),]]>該算法的復雜度是Σl=1logPP2l-1O(N22lP)=Σl=1logPO(N2)=O(N2logP)---(17)]]>對于典型的P=O(N),這是O(N2logN),對該圖像大小成本要比常規反投射的O(N3)在比例上有利得多。
            經改進的三重基于旋轉的分層反投射容易看出采樣圖像旋轉某特定角度——0,±π/2和π——是計算上廉價的運算,因為它們不涉及內插而最多只是現有像素的重新排列。本發明算法的計算效率可通過將這些免運算(free operation)結合其中來改進。二進制算法(圖5)可更改成在每個階段累加三個而非兩個圖像。每個這種三元組的中央圖像旋轉0弧度-一免運算。為了使用-π/2的自由旋轉,它被包括在最后階段中構成圖像之一旋轉-π/2并累加到另一未旋轉圖像上。雙重和三重階段的這個特定組合導致一組大小為2×3L的視角(L為某整數)。盡管算法可適用于任意組投射和任意個投射,該算法的描述和分析可通過假設其數量正好為2×3L且均勻分布而簡化如下θl=-π4(1-3-L)+Δ0(i-1)]]>其中i=1,2,...,2·3L且Δθ=π/(2·3L)(18)因此,視角可被分成兩組,每組3L個,一組以θ=0為中心而另一組以θ=π/2為中心。該三重算法的框圖以L=2的特定情形在圖8中示出,即一組P=2×32=18個投射的特定情形。數字中間圖像Il,md是基礎的連續中間圖像Il,m的采樣版本。在該框圖中,標為B0的框表示零角度反投射運算符。標為Kδ的框表示數字圖像旋轉運算符。如下所述,用于數字圖像和用于數字圖像旋轉的采樣周期和旋轉角度被選為使計算成本最低。
            選定的聚合中間數字圖像可伸展以產生新的中間圖像。圖8中標有↑yUl,m的框表示沿第二個(慢)坐標系上采樣,這對算法進展時容納中間圖像的漸增帶寬是有必要的。對于(18)和(22)中視角的配置,第1層上中間圖像的帶寬為Ωs(Il,m)=πsin(Δθ(3l-1-1)/2)和Ωf(Il,m)=π。在算法標為(1),...(9)的各關鍵點上的標準化光譜支持如圖9所示。隨著算法的進展(當1增加時),慢帶寬增大,而所需的慢采樣周期Δs<π/Ωs(Il,m)=sin((3l-1-1)/2)-1減小。各個上采樣塊中的上采樣因數Ul,m可選擇為滿足這些要求。
            小節0.0.0.1描述涉及兩個切變的各個離散圖像的可分旋轉如何結合到反投射算法中。為了改進反投射后圖像的質量,中間圖像的系統過采樣被引入。小節0.0.0.1描述了經更改包括了該過采樣的算法。
            用于坐標變換的中間旋轉角可如下選擇。對于一般L,分層結構包括(L+1)層。在第零層,每個投射q0被零反投射(B0)以產生相應圖像B0qθ。然后各圖像被分成三組,并組合成1/3多的圖像。在層1上的分組由以下關系式定義。在后續層中,l=1,...,L。
            Il+1,m=κ(δl,3m 2)Il,3m-2+κ(δl,3m-1)Il,3m-1+κ(δl,3m)Il,3m(19)和 IL+2,1=κ(δL+1,1)IL+1,1+κ(δL+1,2)IL+1,2(20)對于(18)中的視角配置,最優中間旋轉角是δl,3m..1=0,l=1,2,...,L、δL+1,1=-π/2和δL+1,2=-π/2。在檢查圖10之后,容易看出構成圖像Il,m的投射的指數Nl,m如下Nl,m={3l-1(m-1)+1,3l-1(m-1)+2,...,3(l-1)m}l=1,2,...,L+1 (21)通過方程式(21),b=3l-1(m-1)+1和e=3l-1(m-1)+3l-1,且使用方程式(18)和(15),αl,m*=-π4(1-3-L)+Δθ(3l-1(m-1/2)-1/2)]]>。與二元的相似,可得出對于l=1,...,L, 對于最后階段我們使用自由旋轉π/2和0弧度。
            0.0.0.1坐標變換可基于切變,如下所述。中間圖像的數字圖像旋轉可由圖10(B)所示的兩個數字圖像切變序列代替,其中x和y坐標中的切變定義如下定義0.7 x切變和y切變由以下定義(Sx(α)f)(x→)=f(Sxαx→)]]>(Sy(α)f)(x→)=f(Syαx→)]]>其中x切變和y切變矩陣分別為Sxα=Iα01]]>和Syα=10α1.]]>該選擇性實施例從旋轉矩陣的雙切變因數分解導出Kθ=SytanθSx-sinθcosθSc,]]>其中Sc=cosθ001/cosθ.]]>該分解暗指在數字圖像的情形中,使用完美的限帶寬內插和帶寬充分受限的數字圖像作為雙切變數字圖像坐標變換的輸入,該雙切變圖像僅是圖像旋轉θ,且有效地在x下采樣并在y上采樣1/cosθ因數。為了校正采樣圖中的這些變化,在每次執行雙切變運算時都需要非整數的重新采樣。相反,通過不同地采樣中間圖像,僅需要對最終圖像解決累計效果。
            方便的選擇是通過在x上采樣初始數字中間圖像并在y上下采樣,在一開始分層算法時就校正重新采樣。當如圖10所示實施例中通過單視圈反投射產生初始中間圖像時,y上下采樣是自由的,因為它僅涉及初始反投射到具有不同y密度的采樣圖。因此未在圖10中框圖明顯示出。第一階段中在x上采樣避免了在后面階段中的偽信號(aliasing)問題。對于初始圖像中的每一個,對最終圖像的累計有效下采樣/上采樣因素1/ПKk=1cosθk在對到最終圖像的路徑的所有雙切變變換上計算,而逆運算應用于所述初始圖像。可得出初始圖像Il,pd,p=1,...,P中的每一個用不同的x和y采樣密度進行重新采樣。
            因此10(A)所示的結果算法包括首先重新采樣所有圖像適當因數B0qθp,p=1,...,P(取決于θp),然后執行雙切變數字圖像坐標變換,如圖10(B)所示。該雙切變考慮每個中間圖像的變化采樣圖來執行。
            切變縮放算法基于切變縮放的分層反投射算法以反投射的定義為基礎,該反投射是被縮放(拉伸/收縮)和切變的單投射圖像之和。
            定義0.8 x切變縮放運算符 由以下定義(C(σ→)f)(x→)=f(Cσ→x→)]]>其中x切變縮放矩陣Cσ→=σ1σ201.]]>所以,方程式(3)可重寫為f^(x→)=Σp=1P[C(cosθp,sinθp)B0qp](x→)---(23)]]>以下是描述可便于顯示的累積切變和切變縮放的效果的一些結果Πi=1nSx(αi)=Sx(Σi=1nαi)---(24)]]>C(σ→′)C(σ→′′)=C(σ1′σ1′′,σ1′′+σ1′σ2′′)---(25)]]>與方程式(24)相似,方程式(25)可歸納地應用,以證實一個n切變縮放序列也是切變縮放。使用該屬性,方程式(23)可用分層結構算法結構來重寫。圖11(A)是方程式(23)的框圖表示(P=4),而圖11(B)是其分層等效體。
            兩種結構的等同性要求建立給定投射角θp和未知切變縮放參數σl,m之間相等的系統,其解總是存在。
            使用整數縮放因數而不是非整數縮放因數在計算上是有利的(根據運算或存儲空間)。例如,使用x切變縮放 (其中整數切變參數σ1=1.0)是有利的,它實際上是純粹的x切變Sx(σ2)。將算法設計成在層l=2,3,...,L上僅使用純粹的x切變并避免x切變縮放是可能的。
            與用一系列切變縮放替換方程式(23)中的每個切變縮放相反,使用方程式(24)和(25)每個切變縮放都用單個切變縮放加一系列切變來替換C(cosθp,sinθp)=Sx(αn)Sx(αn-1)...Sx(α1)C(cosθp,σ)]]>s.t.σ+cosθpΣi=1nαi=sinθp---(26)]]>具有自由參數α1,α2,...,αn。
            結果分層配置在圖13中顯示。方便的選擇是在層l=1上執行x切變縮放 并在后續層l=2,3,...,L上x切變。所以第一層中基礎的連續域圖像的表達式如下I2,m=Σi=02C(σ→3m-i)B0q3m-i---(27)]]>對于l=2,3,...,L,中間圖像Il+1,m由前一層中的三個中間圖像({Il,3m-1}l=02)的x切變版本組成Il+1,m=Sx(αl,3m..2)Il,3m-2+Il,3m-1+Sx(αl,3m)Il,3m(28)并且為了使用0和-π/2弧度的自由旋轉,在最后一層中IL+2,1=IL+1,1+κ(-π/2)IL+1,2圖13中所示算法中坐標變換的參數(用于各種切變縮放和切變)尚未確定。實際上,僅有第一層的(非整數)縮放因數σ1根據相關聯投射的角度被完全確定。中間切變因素({αl,m}l=2,3...L)被選為使中間圖像的帶寬最小,以便于可最稀疏地采樣。給定任一組{αl,m},在層l=1上σ2(切變因數)總是可選為確保由圖13的整個分層算法產生的反投射后圖像 等于反投射表達式(23)。特別地,分層算法初始層中的切變縮放參數取決于組{αl,m},如下所示
            σ2p=sinθp-cosθpΣl=2Lαl,μ(p,l):p≤P2sin(θp-π2)-cos(θp-π2)Σl=2Lαl,μ(p,l):p>P2]]>其中 描述構中從層1中第p個投到根的路徑。
            通過檢查圖13中框圖的上半部分p≤P/2,然后從最后一層反向序列切變,可以看到在該分層結構的上半部分,基礎連續中間圖像Il,m通過切變與其相關聯的虛擬中間圖像 相關,即Il,m=Sx(βl,m)I‾l,m---(29)]]>與方程式(12)的推導相似,可示出中間切變因數取決于{βl,m},如下所示αl,3m-i=βl+1,m-βl,3m-ii=0,1,2. (30)選擇切變系數αl,m的自由度提供選擇β參數中的自由度,它可用于使中間投射的帶寬最小化。這將降低其采樣要求,并改進算法的計算效率。
            設中間圖像Il,3m-1的光譜支持由Wl,3m-i表示,且其y帶寬(即慢帶寬)由Ωy(Wl,3m-i)表示。然后使Ωy(Wl,3m-i)最小化的最優βl,3m-i對i=0,1,2為βl,3m-i*=argminβ∈R{Ωy(Sy(β)W‾l,3m-i)}---(31)]]>如果各視角均勻分布,則可發現βl,3m-1*≈βl+1,m*.]]>然后一種有利的選擇是βl,3m-1*=Δβl+1,m*,]]>結果αl,3m-1*=0,]]>使所需切變運算消除約1/3。
            組{βl,m}通過從分層結構中最后一層到第一層反向地跟蹤算法中中間圖像的光譜支持來確定。根據片投射定理,顯然Il,m的光譜支持是W~l,m={(γcosθp,γsinθp):p∈Nl,mγ∈[-π,π]}---(32)]]>為了采樣要求,光譜支持Wl,m由其標示為 的端點來表征E~l,m={(cosθp,sinθp):p∈Nl,m}---(33)]]>并如圖12所示。該集 的上下帶寬邊界分別由 和 表示。所以如果Nl,m={b,b+1,...,e},則E~l,m+=(cosθe,sinθe)]]>且E~l,m-=(cosθb,sinθb).]]>由于三重組合規則,這些光譜支持端點的上、中、下三組被標示為E~l+1,m0,E~l+1,m1,E~l+1,m2,]]>即E~l+1,mi=E~l,3m-i.]]>并從(30)和有關仿射變換的傅立葉理論得出El,m=Sy(βl,m)E~l,m.]]>然后最優切變因數為αl,3m-i*=-(El+1,mi)y++(El+1,mi)y-(El+1,mi)x++(El+1,mi)x-i=0,2---(34)]]>
            圖14中示出的算法FINDALPHAS尋找當在切變縮放反投射算法中從層L到層2向下遍歷該分層結構時的最優切變因數αl,m*。該算法的輸入是在第(L+1)層開始時圖像的光譜支持的端點集EL+1,m,m=1,2。特別地,在均勻分布角度的情形中,EL+1,1=EL+1,2={(cosθp,sinθp)p=1,2,...,P/2}。
            尋找上采樣因數諸算法和本發明不同實施例的y上采樣因數{Ul,m}確定慢方向中中間圖像的采樣密度。這些上采樣因數可選為符合中間圖像的采樣要求,同時限制為整數值。限制為整數有計算上的優點,并出于同樣原因可擴展成具有小值分母的有理數。
            選擇這些因數的問題在切變縮放算法中比雙切變情形略簡單些。在前一情形中,查看圖13容易看出中間圖像Il,m的慢采樣間隔是Πl′=lLUl′,μ(l′,l,m),]]>其中 述了從Il,m到算法最后節點的路徑。采樣理論要求慢方向采樣間隔如下Δsl,m=Πl′=lLUl′,μ(l′,l,m)π/Ωy(Wl,m)---(35)]]>其中Ωy(Wl,m)是Il,m的y帶寬。假設一組上采樣因數為u=Δ{Ul,m∈Z:l=2,3,...,L;]]>m=1,2,...,6.3L-1}整個算法的計算成本可示為j(u)=Σl=2LΣm=16.3L-1(cl,mΠl′=lLUl′,μ(l′,l,m))+c---(36)]]>其中常數cl,m表示應用于中間圖像Il,md的數字坐標變換的相對計算成本。這些常數通過所使用數字濾波器的階、以及坐標變換的特定實現來確定,如將要在坐標變換的有效實現中所述。滿足約束方程式(35)的使J(U)最小化的最佳集U可通過動態編程或綜合性搜索來解出。對于一給定視角組,上采樣因數的最佳集可預先計算并存儲在查詢表中。
            在雙切變旋轉算法的情形中,問題因為使用雙切變旋轉導致中間圖像的實際分數上下采樣的事實而略為復雜化。特別地,算法相鄰層中的慢采樣周期Δl,m相關如下 其中 所以對慢方向采樣間隔的限制變成如下Δsl,m=Σl′=lL(Ul′,(l′,l,m)κl′,(l′,l,m))π/Ωy(Wl,m)---(37)]]>
            其中Ωy(Wl,m)是Il,m的y帶寬。方程式(36)給出的計算成本現在可以方程式(37)中的約束為條件最小化。
            算法的過采樣版本過采樣被結合在本算法中,以改進重建的準確度。一種這樣做的方法是在算法中的所有中間圖像上統一在算法中執行內插時,制作運算的圖像至少過采樣某預定常數γ。特別地,進行數字坐標變換或重新采樣的每個中間圖像的Nyquist(尼奎斯特)速率與采樣頻率(慢方向和快方向上)之比最好應小于γ。該過采樣最好在盡量少地更改本發明算法時結合。
            慢方向上的過采樣在先前所述的本發明算法中,慢方向上的采樣頻率由上采樣因數Ul,m控制。這證實是在算法過采樣版本中保持過采樣條件的有用工具,因此導致對本發明算法的結構沒有改變。以1/γ倍上采樣通過簡單地更改慢方向采樣間隔上的約束γ倍來實現,即需要Δsl,m≤γπ/Ωs(Il,m),l=2,3,...,L.]]>快方向上的過采樣在快方向上,計算上廉價的整數上采樣并不用來控制過采樣,所以算法被更改為包括分數上采樣。在雙切變分層反投射算法中,這種慢方向上的分數上采樣因數{Ul,m∶m=1,2,...,2.3L}已包括在l=1層中。因此我們簡單地增大上采樣因數Ul,m以在執行了上一數字坐標變換之后結合該過采樣并下采樣該圖像,以將該圖像返回給所需采樣方案(其中Δf=Δs=1.0)。這種對雙切變算法的更改因此僅涉及分數重新采樣的一個附加層。該x坐標重新采樣與第L層中的數字x切變組合,用于改進計算效率。該算法的框圖如圖15所示。在切變縮放分層反投射算法中,x上采樣運算與x切變縮放(在算法開始時為 )組合,以避免額外層重新采樣。容易看出,x上采樣僅僅是x切變縮放的特定情形,其中切變因數設定為0。該兩種運算的組合仍然是x切變縮放C(σ1,σ2)C(1/U,O)=C(σ1/U,σ2/U)。
            算法結束時的下采樣與第L層中的x切變(SαL,m|m=1,...,6)相組合,以使它們成為有效的數字圖像的切變縮放。可示出x切變α加上x下采樣U是有效的x切變縮放C(U,O)Sx(α)=C(U,1+Uα)。
            這些上采樣和下采樣因數的確切值根據參數γ和中間圖像的光譜結構確定。在快方向中,過采樣條件是快方向(y坐標)采樣間隔滿足Δfl,m≤γπ/Ωf(Il,m),l=2,3,...,L.]]>所以所必需的附加上采樣為
            η=max2≤l≤LΩf(Il,m)γπ---(38)]]>因此,過采樣算法的上采樣因數 根據非過采樣的雙切變算法的Ul,m更改如下U~1,m=U1,m,U~L+1,m=η---(39)]]>換言之,第一層中的上采樣因數被更改成確保算法中的所有中間圖像都根據參數γ來過采樣。因此,在最后旋轉之后,圖像具有快速采樣間隔1/η;因此,在第L層,它們必須被下采樣η以向單元采樣方案返回它們。
            在此必須注意,因為輸入投射被假設為確切地以Nyquist速率采樣,快方向中的過采樣條件將不被第一層圖像B0q0p,p=1,...,P滿足。
            該三重過采樣的基于切變縮放的算法框圖如圖16所示。
            相似運算的壓縮序列在同一單個坐標上有兩個運算序列時,它們可相組合以改進計算成本和重新采樣準確度。以下是除了前述的組合x切變與x上采樣/下采樣、或x切變縮放與x-上采樣/下采樣的另一個示例。
            圖10的非過采樣的雙切變算法包括第L階段中6個中間圖像的4個的y切變加上x切變。我們通過將兩個運算序列——該階段的x下采樣加上x切變——壓縮成一個,來把圖像的最終下采樣結合成圖15的過采樣版本。這使內插級聯的長度與非過采樣情形中未有改變。
            任意基數/分支系數的分層對于視角組不是形式2*3L的情形,所有這些算法都可簡便地更改。盡管較佳實施例用于包括聚合中間圖像的三元組的分層算法的所有分支,或使用±π/2或0的旋轉,但可在每個階段上組合任意數量的中間圖像。
            假設在算法的特定層上有任意數量(假設為M)的中間圖像,它們中的3×M/3可組合成三個組(其中x是小于或等于x的最大整數)。如果剩余數量的中間圖像是兩個,則它們可聚合成對以在下一層上產生圖像。如果僅剩下單個中間圖像,則它可不加更改地傳遞給分層算法的下一層。
            分層算法的分支系數(聚合在分層算法節點上的分支數量)可變化,不僅是為了容納任意數量的視角,而且是為了減少分層算法的深度,從而改進圖像質量。在該情形中,不成對或不以三元組、而是以更大的組來聚合圖像是有用的。
            坐標變換的參數的先前規定可簡便地擴展成具有任意分支系數的節點。例如,在基于旋轉的算法中,其中Il+1,m=Σi=1Mκ(δl,mi)Il,mi,]]>程式(10)和(15)中所述的關系仍然成立,因此旋轉角度由δl,mi=αl+1,m-αl,m,指定。在基于切變縮放算法的情形中,其中Il+1,m=Σi=1MSx(αl,mi)Il,mi,]]>方程式(29)仍然成立。將概念El,ml擴展成指正在聚合的M組投射的第i組,可獲得等于方程式(34)右側的αl,m*的方程式。
            0.1基于其它圖像變換的分層算法對于任何矩陣Aθ,反投射方程式(3)可寫成形式f^(x→)=1PΣp=1P[B0qp](AθPx→)---(40)]]>只要Aθ的第一行為[cosθsinθ]。對Aθ的剩余兩個條目選擇任意值的自由度允許在設計分層算法中使用的坐標變換時的靈活性。矩陣Aθ可被表示為與θ相關的某些參數向量δl的因式Aθ=Πl=1LA(δl),]]>但由于Aθ兩個底部條目的自由度而沒有完全確定。該因式分解可用來得出反投射方程式(40)的分層分解,諸如圖5所示的相應框圖,其中坐標變換步驟由表示定義為(κ(δl,m)f)(x→)=f(A(δl,m)x→)]]>的圖像坐標變換運算符的κ(δl,m)標示。
            顯然,在此所述的本發明特定實施例是矩陣Aθ及其因式分解、及相關聯坐標變換的較普遍選擇的特定情形。這些坐標變換對中間圖像的傅立葉光譜的效果與已經描述的情形相似地加以分析,因為由矩陣A在空間域中的仿射變換的效果是矩陣A-T(A的逆轉置矩陣)在頻率域中的仿射變換。因此類似的考慮用來選擇變換中的自由參數,其目標是降低計算要求。因而,用于本發明的分層反投射算法中的一類數字圖像坐標變換除了所述特定實施例之外還包括許多其它變換。此外,因為矩陣A(δl,m)可被因式分解成三角矩陣的積,所以坐標變換可按需作為單軸坐標變換的級聯來執行。
            0.2數字圖像坐標變換和重新采樣的有效和準確實現本發明的分層反投射和重新投射算法的準確度和速度取決于各個數字圖像坐標變換和重新采樣運算的實現細節。經改進的準確度一般需要通常增加計算成本的高階濾波器或內插。高階濾波或內插的成本可通過將運算分解成較低階運算來降低。如果所使用的濾波器是移位不變的,則可獲得計算和/或存儲器需要的其它降低。高階有限脈沖響應濾波器可通過使用低階遞歸濾波器、或通過使用快速傅立葉變換(FFT),來有效地實現。
            特別地,如果對連續圖像假設可分的表示基礎,則數字圖像的數字y切變可通過數字圖像陣列中各個豎直行的部分延遲來獲得。類似地,數字x切變可表達為某時一行的部分延遲。從而,1D信號的部分延遲可使用移位不變濾波器來完成。3D圖像和切變運算的類似分解是眾所周知的。
            重新采樣數字圖像通常也可使用常常是移位不變的低維運算來執行。例如,沿一坐標上采樣整數倍U的圖像可被分解成U個不同的有效計算的部分延遲。這實際上是數字信號處理中眾所周知的所謂多相分解。沿一坐標的有理但非整數重新采樣可被分解成每個都有效執行的整數下采樣和上采樣的級聯。
            更普通的數字圖像重新采樣也可被分解成低維運算。考慮數字2D圖像f從位于標示為CF1的曲線族上樣本S1的采樣圖,重新采樣成具有位于不同曲線族CF2上的點S2的另一個圖,從而產生數字圖像h。如果兩個曲線族足夠密度地相交,則參照圖23所述的用于從一個扇形重新采樣成另一旋轉后扇形的方法可用于一般曲線。否則可引入以所需密度與CF1和CF2相交的第三曲線族CF3。數字圖像可從CF1重新采樣到其與CF3的交點,然后到CF3與CF2的交點,最后到CF2,成為所需采樣圖。
            通過例如考慮面而非曲線,該過程歸一化為3D,以首先減少對表面的重新采樣之一的過程,然后使用2D過程進一步減少對表面的重新采樣以對曲線重新采樣。
            數字坐標變換和重新采樣常常可組合,從而改進計算效率。例如,用于圖13中切變縮放算法的數字x切變縮放運算可被分解成對各條水平線的重新采樣運算。
            1發散光束的快速分層反投射算法1.1扇形光束投射和反投射考慮在圍繞對象的圓上放置的等角分隔檢測器的情形。兩維圖像f(x,y)的源角度為β的扇形光束X線斷層造影投射由(Rβf)(γ)標示,并被定義為沿以原點為圓心半徑為D的圓上的源位置為中心、參數為γ的扇形射線的線積分組。函數 被假設為在半徑為R的盤外為零值。
            如圖17所示的扇形光束射線V(β,γ,T)定義如下V(β,γ,T)=D-sinβcosβ+Tsin(β+γ)-cos(β+γ)---(41)]]>在源角度為β扇形角度為γ的扇形光束投射為(Rβf)(γ)=∫T=0∞f(V(β,γ,T))dT.]]>因為f在半徑為R的盤外為0,所以僅需在盤內即TST和TEND之間執行積分。
            在使用扇形光束幾何結構的計算X線斷層造影中,投射在離散源角度集{βp∶p=1,2,...,P}上可用,且在每個扇形內射線的角度由{γJ∶j=1,2,...,J}索引。在等角扇形光束幾何結構的情形中,檢測器均勻分布在以該源為中心的圓弧上,所以扇形角度是均勻間隔的。在等間隔檢測器的情形中,檢測器均勻分布在與從源位置到原點的線垂直的線上。在此所述的用于扇形光束幾何形狀的快速反投射算法假設為等角分布,其中J為奇數,且γ3=Δγ·(j-(J+1)/2)。但是,諸算法可簡便地擴展成其它扇形角幾何結構。
            根據一組P個扇形光束投射{Rβi(γ):i=1,2,...,P}]]>的重建算法,可表達為每個扇形光束投射加上加權反投射(5)的縮放和濾波f^(x→)=Σp=1PW(T(x→,βp))qβp(γ(x→,βp))---(42)]]>其中W(T)是適當的加權函數, 和 是分別是沿射線的源和源位置β的圖像點 之間的距離,以及該射線的扇形角度,而qβ(γ)是加權和濾波后的投射Qβ(γ)=R′β(γ)*g(γ),R′β(γ)=Rβ(γ)·D·cosγ其中g(γ)是適當的濾波器。
            定義1.1單個源角度β上函數q的加權反投射定義為(Wβq)(x→)=ΔW(T(x→,β))q(γ(x→,β))---(43)]]>并且加權反投射運算符Wβ→(β→∈
            ]>定義為(Wβ→{qp}p=1P)(x→)=ΔΣp=1P(Wβpqp)(x→)---(44)]]>因此f^(x→)=(Wβ→{qp}p=1P)(x→).]]>容易示出Wβ=Κ(-β)W0(45)其中Κ(β)表示β運算符的旋轉,因此(Wβ→{qp}p=1P)(x→)=Σp=1P(κ(-βp)W0qp)(x→)---(46)]]>因而,像平行光束情形中一樣,P個扇形光束投射的加權后投射可被表達為加權后零反投射圖像之和,如圖4可見,其中B0=W0。實際上,方程式(43)和(44)中定義的反投射運算符的接近相似性更普遍地應用于從一般形式的投射重建兩維或更多維的函數,其中函數 和 適當定義。本發明的諸方法擴展到這些具有適當定義的坐標變換K的其它應用。
            1.2快速分層反投射用于扇形光束幾何結構的快速分層反投射算法與平行光束情形的相似。它使用了由投射角為B彼此接近的投射所形成的中間圖像可稀疏采樣的事實,來組合三重分層結構中的扇形光束投射。框圖如圖18所示。相似的算法可對雙重或其它(可能混合)基數結構得出。
            支配圖18中基礎連續圖像的組合的方程式如下Il,m=ΔW0qm---(47)Il+1,m=Δκ(δl,3m-2)Il,3m-2+Il,3m-1+κ(δl,3m)Il,3m,l=1,2,...,L---(48)IL+1,1=ΔIL,1+κ(-π/2)IL,2+κ(-π)IL,3+κ(-3π/2)IL,4---(49)]]>在隨后所述的實施例中,假設源角度在
            ]>其中△β=π/(2·3L)然后中問旋轉角使用方程式(22)像平行光束情形一樣進行選擇,其中Δβ替換Δθ。
            具有圖18框圖的算法的三層中的中間圖像的構成扇形如圖20所示。
            1.2.1扇形光束采樣方案扇形光束算法的一個實施例使用類似于平行光束情形導出的算法中的中間圖像的采樣方案,這接下來描述。可選方案和改進在后面描述。
            單個扇形光束反投射圖像(由單個扇形光束的加權反投射產生的圖像W0q)具有可進行稀疏采樣的結構。可從方程式(43)導出(Wβq)(V(β,γ,T))=W(T)q(γ)。因此在方位為β=0的扇形中由γ標示的射線上的特定T=T1上的單個樣本足以沿整條射線指定圖像的值當T是從源到射線上點的距離時,它僅與W(T)成比例。
            算法中的中間圖像是這種單扇形光束的反投射圖像的若干旋轉后版本之和,每一個都從構成投射中生成。我們將在各構成投射所占角度間隔中央的角度上的扇形指為中央構成扇形。這可對應于實際投射——例如在三重算法的情形中,或對應于虛擬投射——例如在雙重算法的情形中。例如,圖20(a)中的扇形是圖20(b)和20(c)的中央構成扇形。
            回想在平行光束算法中,中間圖像在與中央成分對齊的笛卡爾采樣圖上采樣。在基于旋轉的平行光束情形中,中間圖像沿平行的豎直射線采樣。樣本沿這些平行射線中的每一條的間隔被選為確保圖像的其它構成投射足以沿該射線采樣(根據Nyquist采樣標準)。該必要間隔實際上等于豎直射線的交點的間隔,其中一組旋轉后平行射線對應于極值構成投射,即視角中距離中央投射最遠的投射。
            在平行光束情形的直接模擬中,在此中間圖像沿中央構成扇形的射線采樣。沿中央扇形的射線的采樣點是中央扇形與極值構成扇形的交點。這導致樣本未沿每條射線均勻分布。考慮兩個扇形V(β,γj,·)|j=1J]]>和V(β′,γj′,·)|j=1J.]]>假設等角扇形光束幾何結構具有奇數J個檢測器,γi=Δγ·(j-(J+1)/2)。與β’扇形相交的β扇形的第γ條射線上的點(對應于扇形角γr),由方程式V(β,γr,T)=V(β′,γj′,T′)|j=1J,]]>確定,其解為Tβ,β′,γr(γj′)=D[sin(β-β′-γj′)+sin(γj′)]sin(β-β′+γr-γj′)]]>函數Tβ,β′,γr(γ′)描述扇形V(β′,·,·)如何沿扇形V(β,·,·)的第r條射線變化。它傳送有關沿扇形射線的變化采樣速率的信息。斜率越大,則樣本越稀疏。任何γ’上的本地采樣速率與(dT/dγ′)-1成正比。典型的T(γ’)由圖22中的虛線所示。
            對Tβ,β′,γr(γ′)作兩次更改就得到新的采樣函數 1.中間圖像要在半徑R的盤內采樣,其中邊界為盤外每條射線上至少一個樣本(具有T<TST的一個樣本和具有T>TEND的另一個)。定義γEND′=Tβ,β′,γr-1(TEND):]]>γ’的值對應于TEND。具有相鄰源角度的扇形可因為T>TEND而缺少交點。為了校正之,對T>TEND, 的斜率固定為Tβ,β′,γr(γ′END)的斜率。
            2.最終圖像使用單元采樣間隔在笛卡爾圖上采樣。這是中間圖像需要采樣的最小間隔。獲得沿第γ條射線的單元采樣速率的點γ′通過解方程式d/(dγ)Tβ,β,γr(γ)==1/Δγ來得出γ。為了保持單元采樣間隔的約束, γ′≥γ′的斜率被設置為1/Δγ。
            圖22示出典型β和γ的Tβ,β′,γr(γ)和 它還顯示是集{T~β,β′,γr(jΔγ):j∈Z}]]>的射線上實際采樣點的T值。整數j被選擇為在圖像邊界的每一側上都有至少一個樣本的邊界。
            顯然,由在此概述的原理規定的樣本點位置根據扇形的幾何形狀、構成投射集Nl,m的方程式(21)、中間圖像的選定旋轉角δl,m(對于等距離分布視角情形,在方程式(22)中給出)計算,并從 中選擇。
            用于扇形光束采樣的可分旋轉和上內插如在概述中所述,上采樣運算、以及與旋轉運算組合的上采樣可被分解成在計算上有效的一維重新采樣運算。
            1.2.2基于本地傅立葉結構的采樣方案假設圖像f(x→):Rn→R,]]>我們想要找出一采樣函數t(x→):Rn→Rn,]]>使 具有較小的基本帶寬,因此可在點集{t(m→):m→∈Zn}]]>上采樣。我們發現該采樣函數如下。如將要簡述的,知道了 如何由其構成投射組成,我們發現描述函數f應如何在點 上局部采樣的矩陣函數υ(x→)=υ11(x→)υ12(x→)υ21(x→)υ22(x→).]]>然后我們在圖像域上組合該矩陣函數,以得到采樣函數 在我們的算法中,我們知道中間圖像f是這樣的形式對于某些角δp∈[βmin,βmax],f=Σpκ(δp)W0qp.]]>略去W(T)的加權,并假設投射qp實際上以Nyquist速率采樣,我們知道圖像f中每個點上的光譜支持的局部結構。因為每個帶寬受限投射qp的扇形反投射,從單個投射進行扇形反投射的圖像在扇形的輻條方向上具有可略去的空間帶寬,而在垂直方向上其帶寬與距扇形頂點的距離成反比,如圖24所示。當一組這些扇形被旋轉并累加在一起,該結果圖像的局部光譜支持是各個扇形的光譜支持的合并。例如當構成扇形具有βmin和βmax之間的源角度時,如圖25所示,圖像域中點 上的局部光譜支持是方位在θmin和θmax之間具有半徑、以及角度θ上徑向帶寬ΩR/rθ的彎曲楔形。此處Ω是在圖像平面中央經反投射的投射的空間帶寬,假設以Nyquist速率采樣各投射。在等角情形中,Ω是投射的樣本數量除以弧通過經反投射的投射所覆蓋圖像平面的原點的長度。數學上光譜支持是{(cosθωθ,sinθωθ):θ∈[θmin,θmax]]]>和ωθ∈[-Ωθ,Ωθ]}]]>(50)其中θmax=tanRsinβmax-x1Rcosβmax+x2]]>θmin=tanRsinβmin-x1Rcosβmin+x2]]>Ωθ=RΩ/rθrθ=(x1-Rsinθ)2+(x2-Rcosθ)2]]>給定圖像中點 上的局部兩維傅立葉結構的該認知,該點上的矩陣函數 是這樣的采樣矩陣如果該點上的光譜支持在整個圖像上都是統一的(諸如平行光束情形),則會有效采樣圖像。在感興趣的中間圖像中,這產生兩個不同的小帶寬和大帶寬采樣方向。我們將第一個坐標固定為小帶寬方向。
            我們通過(s100s21α10110α21)T]]>來近似該采樣矩陣函數 此時α2=-Ωminsinθmin+ΩmaxsinθmaxΩmincosθmin+Ωmaxcosθmax]]>和α1=Ωmincosθmin-ΩmaxcosθmaxΩmaxsinθmin+Ωminsinθmin]]>S1=(101α10110α21cosθmidsinθmidΩmid)-1]]>和S2=(0110α21cosmaxsinmaxΩmax)-1---(51)]]>此處Ωmin=ΔΩθmin]]>且Ωmax=ΔΩθmax.]]>角度θmid指與源角度(βmin+βmax)/2相對應的投射角度,且Ωmid=ΔΩθmid.]]>通過光譜支持上幾何變量選擇的這些參數確保對采樣矩陣變換的光譜支持受限于[-π,π]×[-π,π]。
            整個圖像的采樣函數通過在整個圖像上求該采樣矩陣函數的積分來尋找——解微分方程集dtidnj=υij(t(n→)),i,j=1,2.]]>這可得出用數字表示的解。即使未使用該確切指定的圖,局部傅立葉支持分析將估計任何采樣圖的有效性。
            一些中間圖像的結果采樣圖如圖26(A)和26(B)所示。對于扇形光束情形,該局部的基于傅立葉方法產生與因前述基于交點方法所導致的相似的采樣圖。從一個采樣圖重新采樣成另一個的前述可分方法也可用于這些圖。該局部傅立葉采樣方法可直接應用,以尋找用于直線上、任意維上曲線或平面的任意投射幾何結構的采樣方案。在平行光束幾何結構的情形中,它減少成先前所述的采樣圖。
            1.2.3可選采樣方案樣本位于其頂點在源軌跡上的扇形上的采樣方案(稱為“采樣扇形”)具有一些缺點。對每條射線單獨選擇采樣點,導致在兩維意義上并非必然最優的方案。盡管算法中的最終圖像在均勻矩形網格上用單位間隔采樣,但是使用上述方案的中間圖像比圖像特定區域內(例如接近扇形頂點)所必需地更密集地采樣。為了調整這個,上述方案可被更改成使該過采樣區域中的采樣稀疏。兩種這樣的可能性在圖27中示出。這兩種可能性結合這樣的事實圖像在第一層中對扇形有效采樣,而在最后一層則需要在矩形網格上采樣。這些方案嘗試結合從對扇形的采樣到網格的逐步轉換。
            在圖27(A)中,樣本在其頂點位于距離原點比源半徑更遠的偽扇形上選擇。在連續各層中,中間圖像從更大范圍的源角度的扇形中形成,且偽扇形頂點離原點的距離增大。在頂點在無窮遠處的最終層上;即,射線是矩形網格的平行線。在圖27(B)中,樣本不位于扇形上,而是位于距離源位置越近越不發散的光束上。在連續各層上,光束變得越來越平行但越來越不發散。取兩維(或者對于3D圖像,三維)視點的這些或其它采樣方案將有助于更快速的算法。
            1.2.4其它優化·在具有全軌跡的扇形光束情形中,源角度分布在
            小。不同之處在于反投射中使用的加權額外的加權因數(例如Parker加權)被引入以說明短期掃描。類似地,諸算法歸一化為源角度的、甚至不允許完整圖像重建的多個源角度小組的不均勻分布。此外,算法中運算的順序可用本領域技術人員所眾所周知的各種方法更改(例如流水線技術)。例如,運算可被安排為單個投射一可用處理就開始的順序。
            權利要求
            1.一種用于從投射(q1...qp)創建像素圖像 的方法(圖5A),包括以下步驟(a)從選定投射(q1...qp)產生(100)中間圖像(I1,P);(b)對選定中間圖像(I1,P)執行數字圖像坐標變換(102),坐標變換的參數被選為說明從中產生中間圖像的投射的視角、并說明所述中間圖像的傅立葉特征;(c)聚合在步驟(b)中生成的變換后中間圖像的子集(104),以生成聚合中間圖像(I2,P/2);以及(d)用遞歸方式重復步驟(b)、(c),直到已處理和聚合所有投射和中間圖像以形成像素圖像 其中坐標變換參數被選為中間圖像的聚合(104)可由稀疏樣本用所需準確度表示。
            2.如權利要求1所述的方法,其特征在于,所述聚合(104、108)通過累加數字圖像來執行。
            3.如權利要求1所述的方法,其特征在于,至少部分中間圖像(In,m)通過反投射選定投射(q1...qp)在步驟(a)生成。
            4.如權利要求1所述的方法,其特征在于,至少部分中間圖像(In,m)的每一個都通過在步驟(a)反投射兩個或更多個選定投射(q1...qp)而形成。
            5.如權利要求1所述的方法,其特征在于,至少部分聚合中間圖像(In,m)的每一個都通過在步驟(d)聚合三個或更多個選定變換后中間圖像而形成。
            6.如權利要求1所述的方法,其特征在于,數字圖像坐標變換使用數字濾波來執行。
            7.如權利要求1所述的方法,其特征在于,選定坐標變換包括數字圖像旋轉。
            8.如權利要求1所述的方法,其特征在于,選定坐標變換包括數字圖像切變(圖10B 120、122)或切變縮放。
            9.如權利要求1所述的方法(圖15),其特征在于,選定坐標變換包括數字圖像的上采樣(101、106)和/或下采樣(109)。
            10.如權利要求8所述的方法,其特征在于,所述數字圖像切變通過一維線性數字濾波器來執行。
            11.如權利要求9所述的方法,其特征在于,所述數字圖像上采樣和/或下采樣通過一維線性數字濾波器來執行。
            12.如權利要求10所述的方法,其特征在于,至少部分所述數字濾波器是移位不變的。
            13.如權利要求11所述的方法,其特征在于,至少部分所述數字濾波器是移位不變的。
            14.如權利要求12或13所述的方法,其特征在于,至少部分所述數字濾波器是遞歸的。
            15.如權利要求12或13所述的方法,其特征在于,至少部分所述數字濾波器使用快速傅立葉變換(FFT)來實現。
            16.如權利要求1所述的方法,其特征在于,選定過采樣應用于選定中間圖像、和/或變換后中間圖像、和/或聚合中間圖像。
            17.如權利要求1所述的方法,其特征在于,使用了非笛卡爾采樣圖。
            18.如權利要求1所述的方法,其特征在于,多個選定坐標變換可在分層算法的一個層內或跨相鄰層組合。
            19.一種用于沿線、曲線或表面的集合從投射(q1...qp)創建像素圖像 的方法(圖31),包括以下步驟(a)產生(184)中間圖像(I1,m);(b)對選定中間圖像執行數字圖像重新采樣(186),所選定樣本的位置被選為說明選定投射的視角、并說明所述中間圖像的傅立葉特征;(c)聚合(190)重新采樣后中間圖像的選定子集,以生成聚合中間圖像(Iz,m);以及(d)用遞歸方式重復步驟(b)、(c),在遞歸的每一層上增加中間圖像的樣本密度,直到已處理和聚合所有投射和中間圖像以形成所述像素圖像;其中采樣方案被選為重新采樣后的中間圖像的聚合可由稀疏樣本用所需準確度表示。
            20.如權利要求19所述的方法,其特征在于,至少部分中間圖像通過選定投射的加權反投射(圖19 180、182)在步驟(a)生成。
            21.如權利要求19所述的方法,其特征在于,至少部分中間圖像的每一個都通過在步驟(a)加權反投射兩個或更多個選定投射(圖19 180、182)而形成。
            22.如權利要求19所述的方法,其特征在于,至少部分聚合中間圖像的每一個都通過在步驟(d)聚合三個或更多個選定變換后中間圖像而形成。
            23.如權利要求19所述的方法,其特征在于,所述中間圖像具有位于線、曲線或表面的族上的樣本。
            24.如權利要求19所述的方法,其特征在于,所述數字圖像重新采樣通過使用位于線、曲線或平面族的交點處的中間采樣方案,由一系列低維數字濾波運算來執行。
            25.如權利要求19所述的方法,其特征在于,選定過采樣度應用于所選定的重新采樣后的中間圖像、以及聚合后的中間圖像。
            26.如權利要求19所述的方法,其特征在于,所述聚合通過累加數字圖像來執行。
            27.如權利要求19所述的方法,其特征在于,所述重新采樣和聚合可跨連續層組合。
            28.如權利要求30所述的方法,其特征在于,至少一個中間圖像在所包括的重新采樣步驟(b)前后加權。
            29.如權利要求19所述的方法,其特征在于,采樣密度的改變通過數字濾波來完成。
            30.一種用于從投射(q1...qp)創建像素圖像 的方法(圖18),包括以下步驟(a)產生(99)多個中間圖像(I1,m),其中至少一個中間圖像對應于非笛卡爾和/或非周期性采樣圖;(b)對選定中間圖像執行數字圖像上采樣或下采樣(106);(c)對上采樣后/下采樣后的中間圖像執行數字圖像坐標變換;(d)聚合(110)在步驟(c)中生成的變換后中間圖像的子集,以生成聚合中間圖像;以及(e)用遞歸方式重復步驟(b)、(c)和(d),直到已處理和聚合所有投射和中間圖像以形成像素圖像;其中數字圖像坐標變換的至少之一用非笛卡爾和/或非周期性的采樣圖來執行,且坐標變換參數被選為中間圖像的聚合可由稀疏樣本用所需準確度表示。
            31.如權利要求30所述的方法,其特征在于,所述聚合通過累加數字圖像來執行。
            32.如權利要求30所述的方法,其特征在于,至少部分中間圖像通過反投射選定投射在步驟(a)生成。
            33.如權利要求30所述的方法,其特征在于,至少部分中間圖像的每一個都通過在步驟(a)反投射兩個或更多個選定投射而形成。
            34.如權利要求30所述的方法,其特征在于,至少部分聚合中間圖像的每一個都通過在步驟(c)聚合三個或更多個選定變換后中間圖像而形成。
            35.如權利要求30所述的方法,其特征在于,數字圖像坐標變換使用數字濾波來執行。
            36.如權利要求30所述的方法,其特征在于,選定坐標變換包括數字圖像旋轉。
            37.如權利要求30所述的方法,其特征在于,所述數字圖像重新采樣的上采樣和/或下采樣通過一維線性數字濾波器來執行。
            38.如權利要求37所述的方法,其特征在于,至少部分所述數字濾波器是移位不變的。
            39.如權利要求37所述的方法,其特征在于,至少部分所述數字濾波器是遞歸的。
            40.如權利要求38所述的方法,其特征在于,至少部分所述數字濾波器使用快速傅立葉變換(FFT)來實現。
            41.如權利要求30所述的方法,其特征在于,選定過采樣應用于選定中間圖像、和/或變換后中間圖像、和/或聚合中間圖像。
            全文摘要
            通過反投射(100)選定投射以生成中間圖像(I
            文檔編號G06T15/08GK1867926SQ200480029785
            公開日2006年11月22日 申請日期2004年9月9日 優先權日2003年9月9日
            發明者A·K·喬治, Y·布瑞斯勒 申請人:伊利諾斯州立大學托管會
            網友詢問留言 已有0條留言
            • 還沒有人留言評論。精彩留言會獲得點贊!
            1
            婷婷六月激情在线综合激情,亚洲国产大片,久久中文字幕综合婷婷,精品久久久久久中文字幕,亚洲一区二区三区高清不卡,99国产精品热久久久久久夜夜嗨 ,欧美日韩亚洲综合在线一区二区,99国产精品电影,伊人精品线视天天综合,精品伊人久久久大香线蕉欧美
            亚洲精品1区 国产成人一级 91精品国产欧美一区二区 亚洲精品乱码久久久久久下载 国产精品久久久久久久伊一 九色国产 国产精品九九视频 伊人久久成人爱综合网 欧美日韩亚洲区久久综合 欧美日本一道免费一区三区 夜夜爽一区二区三区精品 欧美日韩高清一区二区三区 国产成人av在线 国产精品对白交换绿帽视频 国产视频亚洲 国产在线欧美精品 国产精品综合网 国产日韩精品欧美一区色 国产日韩精品欧美一区喷 欧美日韩在线观看区一二 国产区精品 欧美视频日韩视频 中文字幕天天躁日日躁狠狠躁97 视频一二三区 欧美高清在线精品一区二区不卡 国产精品揄拍一区二区久久 99久久综合狠狠综合久久aⅴ 亚洲乱码视频在线观看 日韩在线第二页 亚洲精品无码专区在线播放 成人亚洲网站www在线观看 欧美三级一区二区 99久久精品免费看国产高清 91麻豆国产在线观看 最新日韩欧美不卡一二三区 成人在线观看不卡 日韩国产在线 在线亚洲精品 亚洲午夜久久久久中文字幕 国产精品成人久久久久久久 精品国产一区二区在线观看 欧美精品国产一区二区三区 中文在线播放 亚洲第一页在线视频 国产午夜精品福利久久 九色国产 精品国产九九 国产永久视频 久久精品人人做人人综合试看 国产一区二区三区免费观看 亚洲精品国产电影 9999热视频 国产精品资源在线 麻豆久久婷婷国产综合五月 国产精品免费一级在线观看 亚洲国产一区二区三区青草影视 中文在线播放 国产成人综合在线 国产在线观看色 国产亚洲三级 国产片一区二区三区 久久99精品久久久久久牛牛影视 亚洲欧美日韩国产 四虎永久免费网站 国产一毛片 国产精品视频在 九九热在线精品 99精品福利视频 色婷婷色99国产综合精品 97成人精品视频在线播放 精品久久久久久中文字幕 亚洲欧美一区二区三区孕妇 亚洲欧美成人网 日韩高清在线二区 国产尤物在线观看 在线不卡一区二区 91网站在线看 韩国精品福利一区二区 欧美日韩国产成人精品 99热精品久久 国产精品免费视频一区 高清视频一区 精品九九久久 欧美日韩在线观看免费 91欧美激情一区二区三区成人 99福利视频 亚洲国产精品91 久热国产在线 精品久久久久久中文字幕女 国产精品久久久久久久久99热 成人自拍视频网 国产精品视频久久久久久 久久影院国产 国产玖玖在线观看 99精品在线免费 亚洲欧美一区二区三区导航 久久久久久久综合 国产欧美日韩精品高清二区综合区 国产精品视频自拍 亚洲一级片免费 久久久久久九九 国产欧美自拍视频 视频一区二区在线观看 欧美日韩一区二区三区久久 中文在线亚洲 伊人热人久久中文字幕 日韩欧美亚洲国产一区二区三区 欧美亚洲国产成人高清在线 欧美日韩国产码高清综合人成 国产性大片免费播放网站 亚洲午夜综合网 91精品久久一区二区三区 国产无套在线播放 国产精品视频网站 国产成人亚洲精品老王 91在线网站 国产视频97 欧美黑人欧美精品刺激 国产一区二区三区免费在线视频 久久久国产精品免费看 99re6久精品国产首页 久久精品91 国产成人一级 国产成人精品曰本亚洲 日本福利在线观看 伊人成综合网 久久综合一本 国产综合久久久久久 久久精品成人免费看 久久福利 91精品国产91久久久久久麻豆 亚洲精品成人在线 亚洲伊人久久精品 欧美日本二区 国产永久视频 国产一区二 一区二区福利 国产一毛片 亚洲精品1区 毛片一区二区三区 伊人久久大香线蕉综合影 国产欧美在线观看一区 亚洲国产欧洲综合997久久 国产一区二区免费视频 国产91精品对白露脸全集观看 久久亚洲国产伦理 欧美成人伊人久久综合网 亚洲性久久久影院 久久99国产精一区二区三区! 91精品国产欧美一区二区 欧美日韩亚洲区久久综合 日韩精品一二三区 久久久夜色精品国产噜噜 国产在线精品福利91香蕉 久久久久久久亚洲精品 97se色综合一区二区二区 91国语精品自产拍在线观看性色 91久久国产综合精品女同我 日韩中文字幕a 国产成人亚洲日本精品 久久国产精品-国产精品 久久国产经典视频 久久国产精品伦理 亚洲第一页在线视频 国产精品久久久久三级 日韩毛片网 久久免费高清视频 麻豆国产在线观看一区二区 91麻豆国产福利在线观看 国产成人精品男人的天堂538 一区二区三区中文字幕 免费在线视频一区 欧美日韩国产成人精品 国产综合网站 国产资源免费观看 亚洲精品亚洲人成在线播放 精品久久久久久中文字幕专区 亚洲人成人毛片无遮挡 国产一起色一起爱 国产香蕉精品视频在 九九热免费观看 日韩亚洲欧美一区 九九热精品在线观看 精品久久久久久中文字幕专区 亚洲欧美自拍偷拍 国产精品每日更新 久久久久国产一级毛片高清板 久久天天躁狠狠躁夜夜中文字幕 久久精品片 日韩在线毛片 国产成人精品本亚洲 国产成人精品一区二区三区 九九热在线观看 国产r级在线观看 国产欧美日韩精品高清二区综合区 韩国电影一区二区 国产精品毛片va一区二区三区 五月婷婷伊人网 久久一区二区三区免费 一本色道久久综合狠狠躁篇 亚洲综合色站 国产尤物在线观看 亚洲一区亚洲二区 免费在线视频一区 欧洲精品视频在线观看 日韩中文字幕a 中文字幕日本在线mv视频精品 91精品在线免费视频 精品国产免费人成在线观看 精品a级片 中文字幕日本在线mv视频精品 日韩在线精品视频 婷婷丁香色 91精品国产高清久久久久 国产成人精品日本亚洲直接 五月综合视频 欧美日韩在线亚洲国产人 精液呈暗黄色 亚洲乱码一区 久久精品中文字幕不卡一二区 亚洲天堂精品在线 激情婷婷综合 国产免费久久精品久久久 国产精品亚洲二区在线 久久免费播放视频 五月婷婷丁香综合 在线亚洲欧美日韩 久久免费精品高清麻豆 精品久久久久久中文字幕 亚洲一区网站 国产精品福利社 日韩中文字幕免费 亚洲综合丝袜 91精品在线播放 国产精品18 亚洲日日夜夜 伊人久久大香线蕉综合影 亚洲精品中文字幕乱码影院 亚洲一区二区黄色 亚洲第一页在线视频 一区二区在线观看视频 国产成人福利精品视频 亚洲高清二区 国内成人免费视频 精品亚洲性xxx久久久 国产精品合集一区二区三区 97av免费视频 国产一起色一起爱 国产区久久 国产资源免费观看 99精品视频免费 国产成人一级 国产精品九九免费视频 欧美91精品久久久久网免费 99热国产免费 久久精品色 98精品国产综合久久 久久精品播放 中文字幕视频免费 国产欧美日韩一区二区三区在线 精品久久蜜桃 国产小视频精品 一本色道久久综合狠狠躁篇 91在线免费观看 亚洲精品区 伊人成综合网 伊人热人久久中文字幕 伊人黄色片 99国产精品热久久久久久夜夜嗨 久久免费精品视频 亚洲一区二区三区高清不卡 久久久久国产一级毛片高清板 国产片一区二区三区 久久狠狠干 99久久婷婷国产综合精品电影 国产99区 国产精品成人久久久久 久久狠狠干 青青国产在线观看 亚洲高清国产拍精品影院 国产精品一区二区av 九九热在线免费视频 伊人久久国产 国产精品久久久久久久久久一区 在线观看免费视频一区 国产精品自在在线午夜区app 国产精品综合色区在线观看 国产毛片久久久久久国产毛片 97国产免费全部免费观看 国产精品每日更新 国产尤物视频在线 九九视频这里只有精品99 一本一道久久a久久精品综合 久久综合给会久久狠狠狠 国产成人精品男人的天堂538 欧美一区二区高清 毛片一区二区三区 国产欧美日韩在线观看一区二区三区 在线国产二区 欧美不卡网 91在线精品中文字幕 在线国产福利 国内精品91久久久久 91亚洲福利 日韩欧美国产中文字幕 91久久精品国产性色也91久久 亚洲性久久久影院 欧美精品1区 国产热re99久久6国产精品 九九热免费观看 国产精品欧美日韩 久久久久国产一级毛片高清板 久久国产经典视频 日韩欧美亚洲国产一区二区三区 欧美亚洲综合另类在线观看 国产精品自在在线午夜区app 97中文字幕在线观看 视频一二三区 精品国产一区在线观看 国产欧美日韩在线一区二区不卡 欧美一区二三区 伊人成人在线观看 国内精品91久久久久 97在线亚洲 国产在线不卡一区 久久久全免费全集一级全黄片 国产精品v欧美精品∨日韩 亚洲毛片网站 在线不卡一区二区 99re热在线视频 久久激情网 国产毛片一区二区三区精品 久久亚洲综合色 中文字幕视频免费 国产视频亚洲 婷婷伊人久久 国产一区二区免费播放 久久99国产精品成人欧美 99国产在线视频 国产成人免费视频精品一区二区 国产不卡一区二区三区免费视 国产码欧美日韩高清综合一区 久久精品国产主播一区二区 国产一区电影 久久精品国产夜色 国产精品国产三级国产 日韩一区二区三区在线 久久97久久97精品免视看 久久国产免费一区二区三区 伊人久久大香线蕉综合电影网 99re6久精品国产首页 久久激情网 亚洲成人高清在线 国产精品网址 国产成人精品男人的天堂538 香蕉国产综合久久猫咪 国产专区中文字幕 91麻豆精品国产高清在线 久久国产经典视频 国产精品成人va在线观看 国产精品爱啪在线线免费观看 日本精品久久久久久久久免费 亚洲综合一区二区三区 久久五月网 精品国产网红福利在线观看 久久综合亚洲伊人色 亚洲国产精品久久久久久网站 在线日韩国产 99国产精品热久久久久久夜夜嗨 国产综合精品在线 国产区福利 精品亚洲综合久久中文字幕 国产制服丝袜在线 毛片在线播放网站 在线观看免费视频一区 国产精品久久久精品三级 亚洲国产电影在线观看 最新日韩欧美不卡一二三区 狠狠综合久久综合鬼色 日本精品1在线区 国产日韩一区二区三区在线播放 欧美日韩精品在线播放 亚洲欧美日韩国产一区二区三区精品 久久综合久久网 婷婷六月激情在线综合激情 亚洲乱码一区 国产专区91 97av视频在线观看 精品久久久久久中文字幕 久久五月视频 国产成人福利精品视频 国产精品网址 中文字幕视频在线 精品一区二区三区免费视频 伊人手机在线视频 亚洲精品中文字幕乱码 国产在线视频www色 色噜噜国产精品视频一区二区 精品亚洲成a人在线观看 国产香蕉尹人综合在线 成人免费一区二区三区在线观看 国产不卡一区二区三区免费视 欧美精品久久天天躁 国产专区中文字幕 久久精品国产免费中文 久久精品国产免费一区 久久无码精品一区二区三区 国产欧美另类久久久精品免费 欧美精品久久天天躁 亚洲精品在线视频 国产视频91在线 91精品福利一区二区三区野战 日韩中文字幕免费 国产精品99一区二区三区 欧美成人高清性色生活 国产精品系列在线观看 亚洲国产福利精品一区二区 国产成人在线小视频 国产精品久久久久免费 99re热在线视频 久久久久久久综合 一区二区国产在线播放 成人国产在线视频 亚洲精品乱码久久久久 欧美日韩一区二区综合 精品久久久久免费极品大片 中文字幕视频二区 激情粉嫩精品国产尤物 国产成人精品一区二区视频 久久精品中文字幕首页 亚洲高清在线 国产精品亚洲一区二区三区 伊人久久艹 中文在线亚洲 国产精品一区二区在线播放 国产精品九九免费视频 亚洲二区在线播放 亚洲狠狠婷婷综合久久久久网站 亚洲欧美日韩网站 日韩成人精品 亚洲国产一区二区三区青草影视 91精品国产福利在线观看 国产精品久久久久久久久99热 国产一区二区精品尤物 久碰香蕉精品视频在线观看 亚洲日日夜夜 在线不卡一区二区 国产午夜亚洲精品 九九热在线视频观看这里只有精品 伊人手机在线视频 91免费国产精品 日韩欧美中字 91精品国产91久久久久 国产全黄三级播放 视频一区二区三区免费观看 国产开裆丝袜高跟在线观看 国产成人欧美 激情综合丝袜美女一区二区 国产成人亚洲综合无 欧美精品一区二区三区免费观看 欧美亚洲国产日韩 日韩亚州 国产欧美日韩精品高清二区综合区 亚洲午夜国产片在线观看 精品久久久久久中文字幕 欧美精品1区 久久伊人久久亚洲综合 亚洲欧美日韩精品 国产成人精品久久亚洲高清不卡 久久福利影视 国产精品99精品久久免费 久久久久免费精品视频 国产日产亚洲精品 亚洲国产午夜电影在线入口 精品无码一区在线观看 午夜国产精品视频 亚洲一级片免费 伊人久久大香线蕉综合影 国产精品久久影院 久碰香蕉精品视频在线观看 www.欧美精品 在线小视频国产 亚洲国产天堂久久综合图区 欧美一区二区三区不卡 日韩美女福利视频 九九精品免视频国产成人 不卡国产00高中生在线视频 亚洲第一页在线视频 欧美日韩在线播放成人 99re视频这里只有精品 国产精品91在线 精品乱码一区二区三区在线 国产区久久 91麻豆精品国产自产在线观看一区 日韩精品成人在线 九九热在线观看 国产精品久久不卡日韩美女 欧美一区二区三区综合色视频 欧美精品免费一区欧美久久优播 国产精品网址 国产专区中文字幕 国产精品欧美亚洲韩国日本久久 日韩美香港a一级毛片 久久精品123 欧美一区二区三区免费看 99r在线视频 亚洲精品国产字幕久久vr 国产综合激情在线亚洲第一页 91免费国产精品 日韩免费小视频 亚洲国产精品综合一区在线 国产亚洲第一伦理第一区 在线亚洲精品 国产精品一区二区制服丝袜 国产在线成人精品 九九精品免视频国产成人 亚洲国产网 欧美日韩亚洲一区二区三区在线观看 在线亚洲精品 欧美一区二区三区高清视频 国产成人精品男人的天堂538 欧美日韩在线观看区一二 亚洲欧美一区二区久久 久久精品中文字幕首页 日本高清www午夜视频 久久精品国产免费 久久999精品 亚洲国产精品欧美综合 88国产精品视频一区二区三区 91久久偷偷做嫩草影院免费看 国产精品夜色视频一区二区 欧美日韩导航 国产成人啪精品午夜在线播放 一区二区视频在线免费观看 99久久精品国产自免费 精液呈暗黄色 久久99国产精品 日本精品久久久久久久久免费 精品国产97在线观看 99re视频这里只有精品 国产视频91在线 999av视频 亚洲美女视频一区二区三区 久久97久久97精品免视看 亚洲国产成人久久三区 99久久亚洲国产高清观看 日韩毛片在线视频 综合激情在线 91福利一区二区在线观看 一区二区视频在线免费观看 激情粉嫩精品国产尤物 国产成人精品曰本亚洲78 国产成人精品本亚洲 国产精品成人免费视频 国产成人啪精品视频免费软件 久久精品国产亚洲妲己影院 国产精品成人久久久久久久 久久大香线蕉综合爱 欧美一区二区三区高清视频 99热国产免费 在线观看欧美国产 91精品视频在线播放 国产精品福利社 欧美精品一区二区三区免费观看 国产一区二区免费视频 国产午夜精品一区二区 精品视频在线观看97 91精品福利久久久 国产一区福利 国产综合激情在线亚洲第一页 国产精品久久久久久久久久久不卡 九色国产 在线日韩国产 黄网在线观看 亚洲一区小说区中文字幕 中文字幕丝袜 日本二区在线观看 日本国产一区在线观看 欧美日韩一区二区三区久久 欧美精品亚洲精品日韩专 国产日产亚洲精品 久久综合九色综合欧美播 亚洲国产欧美无圣光一区 欧美视频区 亚洲乱码视频在线观看 久久无码精品一区二区三区 九九热精品免费视频 久久99精品久久久久久牛牛影视 国产精品成久久久久三级 国产一区福利 午夜国产精品视频 日本二区在线观看 99久久网站 国产亚洲天堂 精品国产一区二区三区不卡 亚洲国产日韩在线一区 国产成人综合在线观看网站 久久免费高清视频 欧美在线导航 午夜精品久久久久久99热7777 欧美久久综合网 国产小视频精品 国产尤物在线观看 亚洲国产精品综合一区在线 欧美一区二区三区不卡视频 欧美黑人欧美精品刺激 日本福利在线观看 久久国产偷 国产手机精品一区二区 国产热re99久久6国产精品 国产高清啪啪 欧美亚洲国产成人高清在线 国产在线第三页 亚洲综合一区二区三区 99r在线视频 99精品久久久久久久婷婷 国产精品乱码免费一区二区 国产在线精品福利91香蕉 国产尤物视频在线 五月婷婷亚洲 中文字幕久久综合伊人 亚洲精品一级毛片 99国产精品电影 在线视频第一页 久久99国产精品成人欧美 国产白白视频在线观看2 成人精品一区二区www 亚洲成人网在线观看 麻豆91在线视频 色综合合久久天天综合绕视看 久久精品国产免费高清 国产不卡一区二区三区免费视 欧美国产中文 99精品欧美 九九在线精品 国产中文字幕在线免费观看 国产一区中文字幕在线观看 国产成人一级 国产精品一区二区制服丝袜 国产一起色一起爱 亚洲精品成人在线 亚洲欧美精品在线 国产欧美自拍视频 99精品久久久久久久婷婷 久99视频 国产热re99久久6国产精品 视频一区亚洲 国产精品视频分类 国产精品成在线观看 99re6久精品国产首页 亚洲在成人网在线看 亚洲国产日韩在线一区 久久国产三级 日韩国产欧美 欧美在线一区二区三区 国产精品美女一级在线观看 成人午夜免费福利视频 亚洲天堂精品在线 91精品国产手机 欧美日韩视频在线播放 狠狠综合久久综合鬼色 九一色视频 青青视频国产 亚洲欧美自拍一区 中文字幕天天躁日日躁狠狠躁97 日韩免费大片 996热视频 伊人成综合网 亚洲天堂欧美 日韩精品亚洲人成在线观看 久久综合给会久久狠狠狠 日韩精品亚洲人成在线观看 日韩国产欧美 亚洲成aⅴ人片在线影院八 亚洲精品1区 99久久精品免费 国产精品高清在线观看 国产精品久久久免费视频 在线亚洲欧美日韩 91在线看视频 国产精品96久久久久久久 欧美日韩国产成人精品 91在线亚洲 热久久亚洲 国产精品美女免费视频观看 日韩在线毛片 亚洲永久免费视频 九九免费在线视频 亚洲一区网站 日本高清二区视频久二区 精品国产美女福利在线 伊人久久艹 国产精品久久久久三级 欧美成人精品第一区二区三区 99久久精品国产自免费 在线观看日韩一区 国产中文字幕一区 成人免费午夜视频 欧美日韩另类在线 久久99国产精品成人欧美 色婷婷中文网 久久天天躁夜夜躁狠狠躁2020 欧美成人伊人久久综合网 国产精品福利资源在线 国产伦精品一区二区三区高清 国产精品亚洲综合色区韩国 亚洲一区欧美日韩 色综合视频 国语自产精品视频在线区 国产高清a 成人国内精品久久久久影 国产在线精品香蕉综合网一区 国产不卡在线看 国产成人精品精品欧美 国产欧美日韩综合精品一区二区三区 韩国电影一区二区 国产在线视频www色 91中文字幕在线一区 国产人成午夜免视频网站 亚洲综合一区二区三区 色综合视频一区二区观看 久久五月网 九九热精品在线观看 国产一区二区三区国产精品 99久热re在线精品996热视频 亚洲国产网 在线视频亚洲一区 日韩字幕一中文在线综合 国产高清一级毛片在线不卡 精品国产色在线 国产高清视频一区二区 精品日本久久久久久久久久 亚洲国产午夜精品乱码 成人免费国产gav视频在线 日韩欧美一区二区在线观看 欧美曰批人成在线观看 韩国电影一区二区 99re这里只有精品6 日韩精品一区二区三区视频 99re6久精品国产首页 亚洲欧美一区二区三区导航 欧美色图一区二区三区 午夜精品视频在线观看 欧美激情在线观看一区二区三区 亚洲热在线 成人国产精品一区二区网站 亚洲一级毛片在线播放 亚洲一区小说区中文字幕 亚洲午夜久久久久影院 国产自产v一区二区三区c 国产精品视频免费 久久调教视频 国产成人91激情在线播放 国产精品欧美亚洲韩国日本久久 久久亚洲日本不卡一区二区 91中文字幕网 成人国产在线视频 国产视频91在线 欧美成人精品第一区二区三区 国产精品福利在线 久久综合九色综合精品 欧美一区二区三区精品 久久国产综合尤物免费观看 久久99青青久久99久久 日韩精品免费 久久国产精品999 91亚洲视频在线观看 国产精品igao视频 色综合区 在线亚洲欧国产精品专区 国产一区二区三区在线观看视频 亚洲精品成人在线 一区二区国产在线播放 中文在线亚洲 亚洲精品第一国产综合野 国产一区二区精品久久 一区二区三区四区精品视频 99热精品久久 中文字幕视频二区 国产成人精品男人的天堂538 99精品影视 美女福利视频一区二区 久久午夜夜伦伦鲁鲁片 综合久久久久久久综合网 国产精品国产欧美综合一区 国产99视频在线观看 国产亚洲女在线精品 婷婷影院在线综合免费视频 国产亚洲3p一区二区三区 91成人爽a毛片一区二区 亚洲一区二区高清 国产欧美亚洲精品第二区首页 欧美日韩导航 亚洲高清二区 欧美激情观看一区二区久久 日韩毛片在线播放 亚洲欧美日韩高清中文在线 亚洲日本在线播放 国产精品一区二区制服丝袜 精品国产一区二区三区不卡 国产不卡在线看 国产欧美网站 四虎永久在线观看视频精品 国产黄色片在线观看 夜夜综合 一本色道久久综合狠狠躁篇 欧美亚洲综合另类在线观看 国产91在线看 伊人久久国产 欧美一区二区在线观看免费网站 国产精品久久久久三级 久久福利 日韩中文字幕a 亚洲午夜久久久久影院 91在线高清视频 国产亚洲一区二区三区啪 久久人精品 国产精品亚洲午夜一区二区三区 综合久久久久久 久久伊人一区二区三区四区 国产综合久久久久久 日韩一区精品视频在线看 国产精品日韩欧美制服 日本精品1在线区 99re视频 无码av免费一区二区三区试看 国产视频1区 日韩欧美中文字幕一区 日本高清中文字幕一区二区三区a 亚洲国产欧美无圣光一区 国产在线视频一区二区三区 欧美国产第一页 在线亚洲欧美日韩 日韩中文字幕第一页 在线不卡一区二区 伊人久久青青 国产精品一区二区在线播放 www.五月婷婷 麻豆久久婷婷国产综合五月 亚洲精品区 久久国产欧美另类久久久 99在线视频免费 伊人久久中文字幕久久cm 久久精品成人免费看 久久这里只有精品首页 88国产精品视频一区二区三区 中文字幕日本在线mv视频精品 国产在线精品成人一区二区三区 伊人精品线视天天综合 亚洲一区二区黄色 国产尤物视频在线 亚洲精品99久久久久中文字幕 国产一区二区三区免费观看 伊人久久大香线蕉综合电影网 国产成人精品区在线观看 日本精品一区二区三区视频 日韩高清在线二区 久久免费播放视频 一区二区成人国产精品 国产精品免费精品自在线观看 亚洲精品视频二区 麻豆国产精品有码在线观看 精品日本一区二区 亚洲欧洲久久 久久中文字幕综合婷婷 中文字幕视频在线 国产成人精品综合在线观看 91精品国产91久久久久福利 精液呈暗黄色 香蕉国产综合久久猫咪 国产专区精品 亚洲精品无码不卡 国产永久视频 亚洲成a人片在线播放观看国产 一区二区国产在线播放 亚洲一区二区黄色 欧美日韩在线观看视频 亚洲精品另类 久久国产综合尤物免费观看 国产一区二区三区国产精品 高清视频一区 国产精品igao视频 国产精品资源在线 久久综合精品国产一区二区三区 www.五月婷婷 精品色综合 99热国产免费 麻豆福利影院 亚洲伊人久久大香线蕉苏妲己 久久电影院久久国产 久久精品伊人 在线日韩理论午夜中文电影 亚洲国产欧洲综合997久久 伊人国产精品 久草国产精品 欧美一区精品二区三区 亚洲成人高清在线 91免费国产精品 日韩精品福利在线 国产一线在线观看 国产不卡在线看 久久99青青久久99久久 亚洲精品亚洲人成在线播放 99久久免费看国产精品 国产日本在线观看 青草国产在线视频 麻豆久久婷婷国产综合五月 国产中文字幕一区 91久久精品国产性色也91久久 国产一区a 国产欧美日韩成人 国产亚洲女在线精品 一区二区美女 中文字幕在线2021一区 在线小视频国产 久久这里只有精品首页 国产在线第三页 欧美日韩中文字幕 在线亚洲+欧美+日本专区 精品国产一区二区三区不卡 久久这里精品 欧美在线va在线播放 精液呈暗黄色 91精品国产手机 91在线免费播放 欧美视频亚洲色图 欧美国产日韩精品 日韩高清不卡在线 精品视频免费观看 欧美日韩一区二区三区四区 国产欧美亚洲精品第二区首页 亚洲韩精品欧美一区二区三区 国产精品视频免费 在线精品小视频 久久午夜夜伦伦鲁鲁片 国产无套在线播放 久热这里只精品99re8久 欧美久久久久 久久香蕉国产线看观看精品蕉 国产成人精品男人的天堂538 亚洲人成网站色7799在线观看 日韩在线第二页 一本色道久久综合狠狠躁篇 国产一区二区三区不卡在线观看 亚洲乱码在线 在线观看欧美国产 久久福利青草精品资源站免费 国产玖玖在线观看 在线亚洲精品 亚洲成aⅴ人在线观看 精品91在线 欧美一区二三区 日韩中文字幕视频在线 日本成人一区二区 日韩免费专区 国内精品在线观看视频 久久国产综合尤物免费观看 国产精品系列在线观看 一本一道久久a久久精品综合 亚洲免费播放 久久精品国产免费 久久人精品 亚洲毛片网站 亚洲成a人一区二区三区 韩国福利一区二区三区高清视频 亚洲精品天堂在线 一区二区三区中文字幕 亚洲国产色婷婷精品综合在线观看 亚洲国产成人久久笫一页 999国产视频 国产精品香港三级在线电影 欧美日韩一区二区三区四区 日韩国产欧美 国产精品99一区二区三区 午夜国产精品理论片久久影院 亚洲精品中文字幕麻豆 亚洲国产高清视频 久久免费手机视频 日韩a在线观看 五月婷婷亚洲 亚洲精品中文字幕麻豆 中文字幕丝袜 www国产精品 亚洲天堂精品在线 亚洲乱码一区 国产日韩欧美三级 久久999精品 伊人热人久久中文字幕 久热国产在线视频 国产欧美日韩在线观看一区二区三区 国产一二三区在线 日韩国产欧美 91精品国产91久久久久 亚洲一区小说区中文字幕 精品一区二区免费视频 国产精品视频免费 国产精品亚洲综合色区韩国 亚洲国产精品成人午夜在线观看 欧美国产日韩精品 中文字幕精品一区二区精品