一種改進的全變分降噪算法在低劑量工業(yè)CT重建圖像中的應(yīng)用

摘""" 要：針對低劑量CT重建中的噪聲抑制問題，文章將全變分降噪模型應(yīng)用到CT重建的圖像投影域，并介紹了ROF模型及能量泛函的建立，以及在低劑量CT圖像處理中的應(yīng)用。提出一種改進的ROF模型的實現(xiàn)方法。首先，利用梯度下降算法對投影圖進行迭低降噪處理。其次，使用濾波反投影算法對工業(yè)CT圖像進行重建。最后，通過CUDA并行運算實現(xiàn)了整個改進算法的過程，提高了算法的運行時間。通過模擬噪聲和鋰電池快速在線工業(yè)CT設(shè)備圖像處理結(jié)果證明：所提出的方法能有效降低重建圖像的噪聲，提高圖像峰值信噪比。

關(guān)鍵詞：全變分；低劑量；CT重建

中圖分類號：TP338.6""""""" 文獻標志碼：A"""""""" 文章編號：1009-5128（2024）05-0083-05

收稿日期：2024-01-15

基金項目：陜西省科技廳工業(yè)攻關(guān)項目：SOC及模數(shù)混合大規(guī)模集成電路直流參數(shù)測試系統(tǒng)（2016GY-117）

作者簡介：葛春平，男，甘肅慶陽人，渭南師范學(xué)院物理與電氣工程學(xué)院副教授，理學(xué)碩士，主要從事數(shù)字圖像算法及CT系統(tǒng)重建算法研究；何冰，男，陜西合陽人，渭南師范學(xué)院物理與電氣工程學(xué)院教授，工學(xué)博士，主要從事圖像特征提取、計算機視覺研究。

CT技術(shù)自20世紀 70年代問世以來，在醫(yī)學(xué)影像診斷領(lǐng)域取得了巨大成就，但由于較高的X射線劑量對人體有較大的輻射危害，因而低劑量CT重建方法在醫(yī)療CT中率先展開研究。近年來在工業(yè)檢測中，CT也從以往的實驗室抽檢發(fā)展為在線全檢測。工業(yè)CT中一般都采用較快的旋轉(zhuǎn)速度和較高的取像幀率來提高掃描速度和增加投影數(shù)［1］，增加旋轉(zhuǎn)速度和取像幀率雖然都能降低單幀投影圖像的輻射劑量，但會導(dǎo)致投影圖像噪聲水平升高。相對解析重建算法、迭代重建算法在噪聲抑制方面具有較好的表現(xiàn)，但由于迭代重建算法具有計算復(fù)雜性高、算法并行設(shè)計困難等缺點，所以濾波反投影等解析重建算法仍然是重建算法中最常用和有效的方法之一［2］。

Rudin等［3］提出的經(jīng)典全變分（Total varivation，TV）模型在圖像處理中已獲得極大成功，被廣泛應(yīng)用在圖像復(fù)原［4］、圖像降噪［5］等方面。ROF全變分模型建立在以下四個假設(shè)基礎(chǔ)之上：一是稀疏性假設(shè)， ROF模型基于圖像的稀疏性假設(shè)［6］，即圖像中的信號在某種變換域中是稀疏的。通常這種假設(shè)可以通過對圖像梯度進行稀疏表示來實現(xiàn)。ROF模型利用了圖像的梯度信息，因為許多自然圖像的邊緣是稀疏的。二是全變分假設(shè)，ROF模型引入全變分作為正則項，以促使圖像在梯度方向上保持稀疏性。全變分是圖像梯度的L1范數(shù)，用于測量圖像在空間中的變化。通過最小化全變分，ROF模型能夠去除圖像中的噪聲，同時保持圖像的邊緣結(jié)構(gòu)。三是噪聲假設(shè)，ROF模型假設(shè)圖像中的噪聲是平穩(wěn)且均勻分布的，并且噪聲水平較低。這種假設(shè)有助于ROF模型在去噪過程中更好地保持圖像的結(jié)構(gòu)信息。四是平滑性假設(shè)，ROF模型在圖像去噪的同時也具有平滑圖像的效果。這是因為全變分正則項促使圖像在梯度方向上保持稀疏性，從而限制了圖像的變化。

本文通過分析CT投影數(shù)據(jù)的噪聲統(tǒng)計特性，使用ROF全變分模型建立能量泛函，提出了一種新的ROF模型離散化方法，根據(jù)歐拉-拉格朗日方程（E-L方程），使用該離散化方法，通過梯度下降算法對能量泛函進行極小值迭代求解，實現(xiàn)了對投影圖像的降噪處理。通過對原始重建圖像和降噪后重建圖像的比對分析了算法的有效性。

1"" 全變分降噪算法

1.1"" ROF模型及能量泛函的建立

假設(shè)成像系統(tǒng)是線性移不變系統(tǒng)，可用模型[uo=][k]*[u+n]來描述退化的圖像。其中：[u]表示理想圖像，[k]表示成像系統(tǒng)的點擴散函數(shù)（PSF），[n]為加性噪聲。通過[uo]求解[u]是一個不適定的數(shù)學(xué)問題，可使用Tikhonv正則化模型，但是該模型在去噪的同時，容易造成圖像邊緣模糊，為克服這種弊端，Rudin等提出了經(jīng)典的全變分模型。

針對高斯噪聲，ROF建立了式（1）所示的全變分模型。其中：[u]是待恢復(fù)的圖像，[Ω]表示圖像的空間域，[f]是觀測圖像（含噪聲），λ是正則化參數(shù)，用于控制全變分項的影響，?[u] 是圖像[u]的梯度。式（1）分為前后兩項的和，分別為正則項和保真項。

""""""""""""""""""""""""""""""""""""""""""""""""""""" [minuE]（[u]）[=Ω?u][ + λ2][Ω（u-f）][2]。""""""""""""""""""""""""" """"""""""""""""""""""""""""""（1）

ROF模型的提出是建立在降噪后圖像的總梯度之和變小的假設(shè)之上，所以式（1）中的正則項采用了計算圖像全域范圍內(nèi)的梯度大小。但是，毫無保留地追求小的梯度會引來兩個問題：一是圖像整體變暗，因為像素整體灰度的降低也就意味著像素的梯度變小；二是當梯度減小時，高頻細節(jié)信息的損失。因此在ROF模型中，在以總梯度為能量表示的同時，也增加了第2項保真項，用于擬合觀測數(shù)據(jù)。這一項對比觀測圖像[f]和恢復(fù)的圖像[u]，使得它們在某種程度上相似。

最小化能量泛函的過程實際上就是在找到一幅圖像[u]，使得全變分正則項和數(shù)據(jù)擬合項的加權(quán)和最小。這樣的最優(yōu)圖像在去噪的同時，由于全變分正則項的存在，會在梯度方向上保持稀疏性，從而有助于保持圖像的邊緣結(jié)構(gòu)。要推導(dǎo) ROF 模型的E-L方程，需要對 E （[u]）進行變分［7］，即找到能量泛函E （[u]）的變分導(dǎo)數(shù)等于0的條件。定義式（2）所示函數(shù)F，求泛函F最小值需要滿足式（3）所示E-L方程，分別求式（2）中F 對u、[ux]、[uy]的偏導(dǎo)后，式（3）可表示為式（4），將式（4）使用散度表示可簡寫為式（5）。這個方程描述了ROF 模型的平衡條件。在實際求解中，通常使用迭代算法來求解方程，例如梯度下降算法等。ROF 模型的解對應(yīng)于能量泛函取極小值的圖像，具有去噪效果。""""""""""""""""""""""""""""""""""""""""""""""" [F=] [λ2] （[u-f]）[2+（?u?x）2+（?u?y）2]，""""""""""""""""""""""""""""""""""""""""nbsp;"""""""""""""" （2）""""" """""""""""""""""""""""""""""""""""""""""""""""[?F?u=??x]（[?F?ux]）[+??y] （[?F?uy]），"""""""""""""""nbsp;"""""""""""""""""""""""""""""""""""""""""""""" （3）""""""""""""""""""""""""""""""""""""""""""""""""""" [?F?u=??x]（[?u?x?u]）[+][??y]（[?u?y?u]），"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" （4）""""""""""""""""""""""""""""""""""""""""""""""""""""" [diν]（[?u?u]）[-λ]（[u-f]）[=0]。"""""" """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""" （5）

1.2"" 改進的ROF模型

為了在離散圖像上實現(xiàn) ROF 模型，需要對E-L方程進行離散化。這通常涉及對空間和變量的離散化，本文主要采用差分方法來逼近導(dǎo)數(shù)和積分過程。

1.2.1"" 空間離散化

圖像為離散的像素網(wǎng)格，記為[ui，j]，其中[i]，[j]分別表示像素的坐標。圖像梯度[?u]可通過差分來表示為""""""""""""""""""""""""""""""" [?u=]（[?u?x]，[?u?y]）[=]（[ui+1，j-ui-1，j2x]，[ui，j+1-ui-1，j-12y]），

其中：[x]和[y]是像素之間的水平和垂直距離。

1.2.2"" ROF模型的離散化及求解

只有將連續(xù)的E-L方程離散化，才能采用數(shù)值優(yōu)化算法進行求解，通常采用差分近似。一個常見的方法是使用有限差分來近似梯度和散度運算。將ROF模型中的梯度和散度運算使用中心差分表示，可得到離散化后的式（6）"""""""""""""""""""""""""""""""""""" [minE （u）=λi，j（][?ui+1，j?ui+1，j-?ui，j?ui，j）+12i，j（ui，j-fi，j）2]。"" """""""""""""""""""""""""""（6）

在求解ROF模型約束優(yōu)化問題時，梯度下降法是一種簡單易行的迭代方式［5］，根據(jù)式（6）可以給出式（7）所示的公式，通過式（7）對圖像進行迭代更新，給出了使用第[k]次圖像迭代結(jié)果計算第[k]+1次圖像的過程，其中[η]表示學(xué)習(xí)率，選擇適當?shù)膶W(xué)習(xí)率對收斂至關(guān)重要。通常需要進行調(diào)試和測試以找到合適的值。迭代次數(shù)和收斂準則的選擇也是重要的參數(shù)。可以通過觀察能量泛函的變化或者圖像的變化來確定停止的條件。""""""""" """""""""""""""""""""[ui，j（k+1）=] [ui，j（k）-η] ［ [λi，j]（[?u（k）i+1，j?u（k）i+1，j-?u（k）i，j?u（k）i，j]）[+（u（k）i，j-fi，j）] ］。"""""""""""""""""""""""""""" （7）

1.3"" ROF模型在低劑量CT圖像處理中的應(yīng)用

ROF全變分模型的離散化及迭代算法應(yīng)用到低劑量CT重建時有兩種選擇：一是對原始噪聲投影圖像重建后，在重建圖像上進行降噪處理；二是先對投影圖像進行降噪后再進行重建操作。圖1（a）和圖1（b）分別為同一疊片鋰電池圖像先重建后降噪和先降噪后重建的圖像。通過對比可以發(fā)現(xiàn)圖1中對重建圖像降噪，如果迭代次數(shù)小，降噪效果不明顯；迭代次數(shù)多，圖像易出現(xiàn)蠟質(zhì)化而損失細節(jié)［8］。所以在低劑量重建時，一般都采用對投影圖像先進行降噪處理后再重建。但是對投影圖像降噪意味著所有投影角度圖像需要全部進行一遍迭代處理，計算量巨大，在算法設(shè)計中可以采用CUDA并行實現(xiàn)來加速算法。

1.4"" 使用CUDA加速ROF模型的求解算法

CUDA（Compute Unified Device Architecture，CUDA）是NVIDIA提供的并行計算平臺和編程模型，用于在NVIDIA GPU上加速計算。使用CUDA可顯著提高ROF模型求解算法的計算速度，特別適用于CT重建中對投影圖像進行降噪的這種大規(guī)模圖像數(shù)據(jù)的處理［9］。

使用CUDA提高算法的并行性一般采用以下幾個步驟：

步驟1：待處理的圖像數(shù)據(jù)傳輸?shù)紾PU內(nèi)存中。

步驟2：將ROF模型的核心算法編寫為CUDA核函數(shù)。

步驟3：根據(jù)ROF模型的梯度下降迭代步驟，設(shè)計GPU上的并行計算。

步驟4：在CUDA核函數(shù)執(zhí)行完成后，使用cudaMemcpy函數(shù)將處理后的圖像數(shù)據(jù)傳輸回主機。

對于迭代算法可采用將迭代過程運行在CPU上，將單次迭代的圖像處理算法運行在GPU上，這樣可以兼顧整體算法設(shè)計的邏輯性和時效性問題。

2"" 實驗結(jié)果與分析

2.1"" 模擬噪聲實驗結(jié)果分析

我們對標準測試Lena圖像引入均值為0、方差為0.02的高斯噪聲，分別通過本文方法和11*11的中值濾波算法進行降噪處理，并計算兩種算法的峰值信噪比PSNR，實驗結(jié)果如圖2所示。本文方法的圖像保邊性和PSNR值都要優(yōu)于中值濾波算法。這種保邊性正好是低劑量CT重建時最重要的特性，因為在濾波反投影等解析重建算法中，高通濾波可以通過減小低頻分量來去除偽影，提高圖像的邊緣清晰度。重建偽影通常是由圖像中低頻成分的疊加所致，而這種保邊性可以很大程度保持重建圖像的對比度。

2.2"" 鋰離子疊片電池CT圖像實驗分析

為了分析本文降噪算法在低劑量CT重建中的實際應(yīng)用效果，我們在日聯(lián)科技鋰離子電池快速CT專用設(shè)備LX9800上取了疊片鋰電池的實際CT投影圖像進行分析對照。設(shè)備取圖電壓115 kV，電流300 μA，通過5 s旋轉(zhuǎn)時間取了360張投影圖像進行CT重建，重建方法為濾波反投影算法。

圖3分別表示投影圖像、無降噪重建、降噪迭代5次重建、降噪迭代15次重建的切層圖像。對照圖3（a）和圖3（b）可以發(fā)現(xiàn)，濾波反投影（FBP）算法重建結(jié)果相對投影圖像會放大圖像噪聲，這是因為FBP算法中普遍采用的S-L等濾波器為加窗斜坡濾波器，放大了高頻噪聲信號。圖3（c）和圖3（d）分別表示經(jīng)過5次迭代和15次迭代后的重建切片圖像，正負極電極對比度均比無降噪重建圖像圖3（a）有明顯提升。

3"" 結(jié)語

本文推導(dǎo)了求解ROF模型的E-L方程，并對其進行了離散化方法分析，在此基礎(chǔ)上介紹了將其應(yīng)用到快速CT重建算法的投影圖像降噪上［10］。仿真和設(shè)備實驗數(shù)據(jù)結(jié)果表明：本文的算法針對快速CT低劑量重建可以很好地抑制噪聲和偽影，實現(xiàn)圖像的優(yōu)質(zhì)重建，有實際應(yīng)用意義。同時也對比分析了低劑量CT重建時先降噪后重建和先重建后降噪的圖像差異。

參考文獻：

［1］" LIU Y P，ZHANG F S.Cone_beam Computed Tomography System［M］.Beijing：Intellectual Property Publishing House，2008.

［2］" ZHU K F，JIANG W，WANG D F，et al.New kind of clarity evaluation function of image［J］.Infrared Laser and Engineering，2005，34（4）：464-468.

［3］" RUDIN L，OSHER S，F(xiàn)ATEMI E，et al.Nonlinear total variation based noise removal algorithms［J］.Physica D Nonlinear Phenomena，1992，60（1）：259-268.

［4］" 陳劍鋒，陳錦釗.水下圖像增強和復(fù)原方法研究［J］.閩江學(xué)院學(xué)報，2023（5）：73-80.

［5］" 雷戈，周冬明，楊浩，等.結(jié)合殘差塊和MLP卷積的真實圖像去噪網(wǎng)絡(luò)［J］.電子器件，2023（5）：1325-1332.

［6］" 王麗艷.斷層圖像稀疏性重建模型與算法研究［D］.南京：南京理工大學(xué)，2024.

［7］" 郭從洲，秦志遠，時文俊.基于能量泛函和視覺特性的全變分圖像降噪模型［J］.中國圖象圖形學(xué)報，2014（9）：1282-1287.

［8］" 李智，富文軍.基于小波變換和全變分的圖像去噪［J］.電子設(shè)計工程，2020（2）：126-129.

［9］" JALALZAI K.Some remarks on the staircasing phenomenon in total-variation based image denoising［J］.Journal of Mathematical Imaging amp; Vision，2016，54（2）：256-268.

［10］" 牛善洲.基于變分正則化的低劑量CT成像方法研究［D］.廣州：南方醫(yī)科大學(xué)，2015.

【責(zé)任編輯""nbsp; 牛懷崗】

Application of an Improved Total Variation Denoising Algorithm in Low-Dose Industrial Computed Tomography Reconstruction Images

GE Chunping1，2，HE Bing1，2，YUAN Wei1，2，LIN Guancheng1，2

（1. School of mathematics and Physics， Weinan Normal University， Weinan 714099，China；

2. Engineering Research Center of X-ray Imaging and Detection Shaanxi University， Weinan 714099，China）

Abstract： Aiming at the noise suppression problem in low-dose CT reconstruction， the total variation noise reduction model is applied to the image projection domain of CT reconstruction， and the establishment of ROF model and energy functional is introduced， as well as its application in low-dose CT image processing. An improved implementation method of ROF model is proposed. First， the gradient descent method is used to perform iteration and noise reduction processing on the projection image. Then， the filtered back-projection algorithm is used to reconstruct the industrial CT image. Finally， the entire process of improving the algorithm is implemented through CUDA parallel operation， and the operation of the algorithm is improved time. The results of fast online industrial CT equipment image processing using simulated noise and lithium batteries prove that the proposed method can effectively reduce the noise level of the reconstructed image and improve the peak signal-to-noise ratio of the image.

Key words：total variation； low dose； CT reconstruction