欠采樣的ct圖像重建方法
【技術領域】
[0001] 本發明涉及數字圖像處理技術領域，特別涉及一種對欠采樣CT掃描獲得的原始 投影數據進行數據處理的圖像重建方法。
【背景技術】
[0002] 現有技術中具有代表性的CT圖像重建算法是一種基于濾波反投影的解析重建方 法。此算法以中心切片定理為依據，即二維圖像在某一方向上的投影數據的一維傅里葉變 換函數的值與該圖像的二維傅里葉變換函數在頻域平面上沿同一方向過原點的直線上的 值相等。因此，濾波反投影算法通過各個方向上的投影數據得到一系列一維傅里葉變換函 數，將這些函數按照對應方向加權疊加并進行濾波變換后得到待重建CT圖像的二維傅里 葉變換函數，對其做二維傅里葉反變換得到重建結果。
[0003] 為達到精確重建的目的，上述解析重建方法的采樣頻率必須大于信號最高頻率的 2倍，（奈奎斯特抽樣定理表明若頻帶寬度有限，要從抽樣信號中無失真地恢復原信號，抽 樣頻率應大于2倍信號最高頻率。當抽樣頻率小于2倍頻譜最高頻率時，信號的頻譜有混 疊，解析重建方法從CT圖像頻率域還原到空間域時，會引入因為頻譜混疊產生的偽影，降 低重建圖像的質量。）這使得投影數據的信息量遠遠大于重建圖像本身，采樣數據冗余度過 高，高采樣率的要求與降低CT輻射劑量的需求相悖。而如果降低CT掃描的采樣率，當不滿 足奈奎斯特采樣定律的條件時，重建圖像會存在比較明顯的偽影，影響CT圖像的質量；另 一方面，解析重建方法對投影數據的噪聲也很敏感。
【發明內容】
[0004] 針對現有技術存在的不足，本發明的目的在于通過構建一種全新的數據處理方 法，使得可以通過欠采樣CT掃描獲得的原始數據重建出高質量的CT圖像。
[0005] 解析重建和迭代重建是CT圖像重建的兩種基本方法。與解析重建算法相比，迭代 重建算法具有可在數據不完全和低信噪比（低劑量引起的噪聲影響）條件下高質量重建CT 圖像的優點。針對濾波反投影算法無法應對欠采樣投影數據的缺點，本案提出一種結合加 權懲罰因子的基于字典學習理論的迭代重建算法。該算法首先通過已經訓練構建好的自適 應參數選擇模型確定迭代算法目標函數中正則化參數的值，然后建立基于字典學習理論的 迭代重建模型，在目標函數的正則約束項上添加加權懲罰因子，形成加權字典學習迭代重 建模型，該模型的每次迭代過程交替更新加權懲罰因子、自適應字典、圖像塊的稀疏表示以 及重建圖像，最終達到迭代停止條件。
[0006] 為實現上述目的，本發明通過以下技術方案實現：
[0007] -種欠采樣的CT圖像重建方法，包括以下步驟：
[0008] 步驟1)加載待處理的原始CT掃描數據；
[0009] 步驟2)基于字典學習的迭代重建算法，假設正則化參數λ為正無窮大，得到此條 件下的初步重建圖像
[0010] 步驟3)計算出所述初步重建圖像的前向投影數據與所述原始CT掃描數據的 相對誤差δ λ
[0011] 步驟4)構建一個正則化參數λ的自適應選擇模型，該自適應選擇模型通過以所 述相對誤差S 為自變量來獲取與所述原始CT掃描數據相匹配的正則化參數λ的值；
[0012] 步驟5)構建一個基于圖像塊全變分的加權懲罰因子ws;
[0013] 步驟6)將所述加權懲罰因子Ws與字典學習的約束項相結合，構建一個基于加權 字典學習迭代重建模型的目標函數一；
[0014] 步驟7)對所述目標函數一最小化求解，得到圖像μ ;
[0015] 步驟8)根據圖像μ更新加權懲罰因子ws;
[0016] 步驟9)重復執行步驟7)和步驟8)，直至滿足迭代終止條件，輸出最終的CT重建 圖像μ。
[0017] 優選的是，所述的欠采樣的CT圖像重建方法，其中，步驟2)中的初步重建圖像 KL30通過以下方法獲得：
[0018] 步驟a):構建一個基于字典學習迭代重建模型的目標函數二，該目標函數二的表 達式為：
[0020] 其中，A =丨&} e RM:，其為CT成像的系統矩陣；/ = (He Ru，其為原始投 影數據；街qxpH/ll/lU，其為統計權重；μ = %/%，...，//、.:/，其為待求圖像排成的列向量； 21(叫/2)([鄭];-為數據保真項；Es: = !<}eR^，其為從圖像中提取NqXNq維小圖像塊 的算子；I> = (d.d」…..z，其為字典；a se R KX1，其是以字典D為基底得到的圖像塊的 稀疏表示；SsII EiH-Das |g C1 Il ?Λ為正則約束項；λ和V s均為正則化參數；
[0021] 步驟b):在字典更新過程中，圖像μ的初值取隨機矩陣，之后圖像μ取上一次圖 像更新過程產生的更新結果，并假設圖像μ在本次字典更新過程中保持恒定，則目標函數 二的表達式簡化為簡式一，具體為：
[0022] 娜Σ|Ε'μ-Da』::
[0023] 其中，字典D的更新采用K-SVD算法獲得，稀疏表示a s的更新采用OMP算法獲得；
[0024] 步驟c):在圖像更新過程中，字典D的初值取離散余弦變換矩陣，稀疏表示α 3的 初值以圖像μ和字典D的初值為基礎通過OMP算法獲得，且之后字典D和稀疏表示α 3取 上一次字典更新過程產生的更新結果，則目標函數二的表達式簡化為簡式二，具體為：
[0026] 以二次可分離代理函數來替代簡式二中的數據保真項Σ；^(?5/2)([Αμ],.-^ ,則圖 像μ的最終表達式為： CN 105118078 A 說明書 3/8 頁
[0028] j = l，2，...，N2
[0029] 其中，t為迭代次數，j表示圖像μ的第j個元素，二次可分離代理函數為：
[0032] 步驟d):在正則化參數λ為正無窮大的條件下，重復步驟b)和c)，直至滿足迭代 終止條件，獲得初步重建圖像
[0033] 優選的是，所述的欠采樣的CT圖像重建方法，其中，步驟3)中，相對誤差δ 的 表達式為：
[0035] 優選的是，所述的欠采樣的CT圖像重建方法，其中，步驟4)中，所述正則化參數λ 的自適應選擇模型的表達式為：
[0037] 其中，f ( δ s)的表達式為： ^ _ Jl.74485名 +0.58883? -6.88253 (f Sc >1.96
[0038] ' <<；!；) ' } -0.21545<>；； + 1.08602〇;, - 0.32634 //. 〇 I.%
[0039] 其中，Sg=IO6S^0
[0040] 優選的是，所述的欠采樣的CT圖像重建方法，其中，步驟5)中，所述加權懲罰因子 Ws的表達式為：
[0042] 其中，μ s表示從圖像μ中抽取出的第s個圖像塊，且μ s= reshape (Es μ，N。，N。）， E = K丨e ，其表示從圖像μ中提取Ν。X Ν。維圖像塊的算子；reshape函數表示將E s μ 重新排列為二維圖像塊；I I μ」U為圖像塊的全變分，計算公式為Μν = Σ (Q1.).. ，G為 梯度算子，（Gx)iij= (XiAfXiJXiJ1-Xiij), X代表任意矩陣；mean(| I μ」|τν)表示所有圖 像塊全變分的平均值；ε > 〇。
[0043] 優選的是，所述的欠采樣的CT圖像重建方法，其中，步驟6)中，所述目標函數一的 表達式為：
[0045] 優選的是，所述的欠采樣的CT圖像重建方法，其中，步驟7)中，采用與獲得初步重 建圖像相同的方法對目標函數一進行最小化求解，得到的圖像μ的表達式為：
[0047] 其中，加權懲罰因子Ws的初始值取1。
[0048] 優選的是，所述的欠采樣的CT圖像重建方法，其中，步驟8)中，所述加權懲罰因子 Ws的更新公式為：
[0051] 優選的是，所述的欠采樣的CT圖像重建方法，其中，在步驟9)中，迭代終止條件為 兩次迭代的掃描數據相對誤差增量小于給定閾值T1，迭代終止條件的表達式為：
[0053] 優選的是，所述的欠采樣的CT圖像重建方法，其中，在步驟d)中，迭代終止條件為 兩次迭代的掃描數據相對誤差增量小于給定閾值T2，迭代終止條件的表達式為：
[0055] 本發明的有益效果是：本案能夠很好的適應投影數據不充分的欠采樣條件以滿足 降低輻射劑量的需求，在數據采樣率下降到一定程度時依然能夠高質量重建CT圖像，對低 劑量引起的投影噪聲也有很好的魯棒性。
【附圖說明】
[0056] 圖1為欠采樣的CT圖像重建方法的流程示意圖。
[0057] 圖2為采用不同方法獲得的實物CT圖像的比對圖。
【具體實施方式】
[0058] 下面結合附圖對本發明做進一步的詳細說明，以令本領域技術人員參照說明書文 字能夠據以實施。
[0059] 本案提供一實施例的欠采樣的CT圖像重建方法，包括以下步驟：
[0060] 步驟1)加載待處理的原始CT掃描數據；
[0061] 步驟2)基于字典學習的迭代重建算法，假設正則化參數λ為正無窮大，得到此條 件下的初步重建圖像μL+ ;
[0062] 步驟3)計算出初步重建圖像的前向投影數據與原始CT掃描數據的相對誤差 3 λ -?〇〇 i
[0063] 步驟4)構建一個正則化參數λ的自適應選擇模型，該自適應選擇模型通過以相 對誤差S 為自變量來獲取與原始CT掃描數據相匹配的正則化參數λ的值；
[0064] 步驟5)構建一個基于圖像塊全變分的加權懲罰因子ws;
[0065] 步驟6)將該加權懲罰因子Ws與字典學習的約束項相結合，構建一個基于加權字 典學習迭代重建模型的目標函數一；
[0066] 步驟7)對該目標函數一最小化求解，得到圖像μ ;
[0067] 步驟8)根據圖像μ更新加權懲罰因子ws;
[0068] 步驟9)重復執行步驟7)和步驟8)，直至滿足迭代終止條件，輸出最終的CT重建 圖像μ。
[0069] 具體的圖像重建方法如下：
[0070] 步驟2)中的初步重建圖像通過以下方法獲得：
[0071] 步驟a):構建一個基于字典學習迭代重建模型的目標函數二，該目標函數二的表 達式為：
[0073] 其中，A = KWrf，其為CT成像的系統矩陣，aij表示系統矩陣A中第i行第j列 的元素，Rrf代表大小為I行N2列的矩陣集合，！為原始投影的總數目，N2表示圖像塊的 大小；? =H..，ξ)τ ，其為排成列向量的原始投影數據，Rixi代表長度為I的列向量； 咚= exp(-〖/||?||2)，其為統計權重；μ = 4)'，其為待求圖像排成的列向量，T代表轉 置；細/2)([Al4 -/>為數據保真項；A = ?<,) e ，.其為從圖像中提取NqXNq維小圖 像塊的算子，s代表圖像塊的編號，即第s個小圖像塊，nj代表抽取的小圖像塊以圖像中第η 行第j列的元素作為小圖像塊的左上角元素，nj遍歷整個圖像；Β = (?1:,?12,...、（）εΚΑ?，其為 字典，d代表原子，與圖像塊的大小相同，K代表原子的數目；a se R ΚΧ1，其是以字典D為基 底得到的圖像塊的稀疏表示；Σ』ε#-Dogi+Σλ ΙΚ ΙΙ。