CN107292846B - A recovery method of incomplete CT projection data under circular orbit - Google Patents

A recovery method of incomplete CT projection data under circular orbit Download PDF

Info

Publication number
CN107292846B
CN107292846B CN201710501987.9A CN201710501987A CN107292846B CN 107292846 B CN107292846 B CN 107292846B CN 201710501987 A CN201710501987 A CN 201710501987A CN 107292846 B CN107292846 B CN 107292846B
Authority
CN
China
Prior art keywords
projection
formula
frequency domain
data
iteration
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201710501987.9A
Other languages
Chinese (zh)
Other versions
CN107292846A (en
Inventor
徐圆
马健晖
周凌宏
丁振强
黄杰星
罗宁
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Southern Medical University
Original Assignee
Southern Medical University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Southern Medical University filed Critical Southern Medical University
Priority to CN201710501987.9A priority Critical patent/CN107292846B/en
Publication of CN107292846A publication Critical patent/CN107292846A/en
Application granted granted Critical
Publication of CN107292846B publication Critical patent/CN107292846B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00Image enhancement or restoration
    • G06T5/77Retouching; Inpainting; Scratch removal
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10072Tomographic images
    • G06T2207/10081Computed x-ray tomography [CT]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20048Transform domain processing
    • G06T2207/20056Discrete and fast Fourier transform, [DFT, FFT]

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Image Analysis (AREA)

Abstract

The invention relates to a method for recovering incomplete CT projection data under a circular orbit, which firstly deduces an accurate consistency condition under circular orbit scanning data according to a John equation and realizes a corresponding projection data recovery method according to the condition. The method comprises the following steps: firstly, obtaining the position information of a bad pixel point; then, initializing the bad pixel point projection by using a spatial interpolation method; respectively carrying out Fourier transform on the initialized projection and the projections at two adjacent angles in front and at the back of the initialized projection to obtain transformed projection frequency data; carrying out Fourier inversion on the weighted projection frequency domain data to generate corrected projection data, and obtaining a first iteration final value; and taking the final result of the last iteration as an initial value, and repeating the third step to the fifth step until the requirement of the iteration times is met. The invention has the advantages of image quality, simple realization and applicability.

Description

一种圆轨道下不完全CT投影数据的恢复方法A recovery method of incomplete CT projection data under circular orbit

技术领域technical field

本发明涉及一种医学图像处理,具体涉及图像的增强或复原,特别是涉及一种圆轨道下不完全CT投影数据的恢复方法。The invention relates to medical image processing, in particular to image enhancement or restoration, and in particular to a method for restoring incomplete CT projection data under circular orbits.

背景技术Background technique

锥形束计算机断层扫描是一种可以快速且低辐射剂量获得人体的三维解剖图像的医学成像方法。数字平面探测器的出现使得真正的锥形束计算机断层扫描(Cone-beamComputed Tomography,CBCT)在临床应用成为可能。Cone beam computed tomography is a medical imaging method that can obtain three-dimensional anatomical images of the human body quickly and with low radiation dose. The advent of digital flat detectors makes it possible for true cone beam computed tomography (Cone-beamComputed Tomography, CBCT) to be applied in clinic.

由于CBCT采用锥形束扫描结构,所采集的投影数据具有高度冗余性,想要利用数据的冗余性来达到高效率的目的就必须满足一致性条件(Consistency Conditions)。约翰于1938年从数学的角度推导出了一个超双曲线偏微分方程(Ultrahyperbolic PartialDifferential Equation,PDE)作为线积分的一致性条件,故此方程称为约翰方程。在70年代,随着CT的出现和兴起,约翰方程(John's Equation)被广泛地应用在CT领域并进行了更深入地推导和延伸。Since CBCT adopts a cone beam scanning structure, the acquired projection data is highly redundant. Consistency Conditions must be satisfied if the data redundancy is to be used to achieve high efficiency. In 1938, John deduced a hyperbolic partial differential equation (Ultrahyperbolic PartialDifferential Equation, PDE) from a mathematical point of view as the consistency condition of the line integral, so the equation is called John's equation. In the 1970s, with the emergence and rise of CT, John's Equation was widely used in the field of CT and was further derived and extended.

时至今日,已经有许多学者提出了多种基于约翰方程的一致性条件,如Patch提出一种基于第三代螺旋CT的一致性条件(Sarah K.Patch,Computation of UnmeasuredThird-Generation VCT Views From Measured Views,2002,IEEE Trans.Med.Imaging),在此文献中,作者通过使用第三代螺旋CT参数来替换基于笛卡尔坐标系的约翰方程中的参数,并在此基础上进行了一系列地深入推导,将二阶偏微分方程转换成了一组一阶常微分方程(first-order ordinary differential equation),对此方程进行傅里叶变换,然后应用其频率域的形式可以从获得的投影中计算出未测量的投影数据。根据此理论,文献中使用“旋转+螺旋+旋转”的射线源轨道,螺距为0.128m,每转一圈采集984个投影,然后针对获得的投影数据使用推导方程对应的傅立叶变换形式进行外插,从而获得未测量的投影数据。但是Patch推导出的一致性条件是基于螺旋CT的几何结构推导而成,螺旋结构要求变量z沿着纵轴在不停地变化,这会在公式中引入一个变化量,使得公式计算复杂,解决问题困难。文献中由于采用“旋转+螺旋+旋转”扫描方式,需要在“螺旋”这个扫描阶段前后分别进行工作台移动速度的增加和减小,以保证在这三个扫描阶段数据采集的快速转换和平滑过渡,但是在实际应用中,一般螺旋CT的工作台移动速度都是恒定的,因此,文献中的一致性条件很难具体应用。此外,Hao Yan等人在使用移动的波束停止阵列(beam stop array,BSA)来进行散射校正时采用了Patch基于三代螺旋CT的一致性条件(Hao Yan,XuanqinMou,Shaojie Tang,Qiong Xu and Maria Zankl,Projection movingbeam stop array,2010,Phys.Med.Biol.)。波束停止阵列是一个遮挡板,上面有规律地排列着遮挡点。X射线经过遮挡板照射到物体上,其中经过遮挡板上遮挡点的X射线会被屏蔽掉。如果在探测器遮挡点对应位置出现数值,则这个数值是散射值。由于使用了BSA,得到的投影数据是不完整的,因此需要使用一致性条件对其进行恢复。此文献采用Patch基于三代螺旋CT的推导公式[Patch(2002)],并结合空间插值方法,提出一种PC_SI方法用于散射校正领域中。但如前面提到的那样,由于螺旋CT的z轴在不停地变化,导致此公式很难具体应用,所以此文献参照Defrise等人的做法,用探测器的纵坐标来近似代替z坐标,但是这种近似做法,会导致精度的下降。从几何结构上来看,只需要令z=0,螺旋CT就可以退化成圆轨道CBCT,但是文献中的公式是对z求导,当z为常数时导数没有意义,所以此文献中所用的方法并不适用于圆轨道CBCT。Up to now, many scholars have proposed a variety of consistency conditions based on John's equation, such as Patch proposed a consistency condition based on the third generation spiral CT (Sarah K. Views, 2002, IEEE Trans.Med.Imaging), in this paper, the author replaces the parameters in the John's equation based on the Cartesian coordinate system by using the third-generation helical CT parameters, and based on this, a series of In-depth derivation, converting a second-order partial differential equation into a set of first-order ordinary differential equations, Fourier transforming this equation, and then applying its frequency-domain form can be obtained from the projections obtained Unmeasured projection data is calculated. According to this theory, the "rotation + helix + rotation" ray source orbit is used in the literature, the pitch is 0.128m, 984 projections are collected per revolution, and then the obtained projection data is extrapolated using the Fourier transform corresponding to the derivation equation. , to obtain unmeasured projection data. However, the consistency condition derived by Patch is derived based on the geometric structure of the helical CT. The helical structure requires the variable z to change continuously along the vertical axis, which will introduce a change in the formula, which makes the calculation of the formula complicated and solves the problem. Problem is difficult. In the literature, due to the "rotation + helix + rotation" scanning method, it is necessary to increase and decrease the moving speed of the table before and after the "spiral" scanning stage to ensure the fast conversion and smoothness of data acquisition in these three scanning stages. However, in practical applications, the moving speed of the working table of general helical CT is constant, so the consistency conditions in the literature are difficult to be applied concretely. In addition, Hao Yan et al. adopted the consistency condition of Patch based on three generations of helical CT when using a moving beam stop array (BSA) for scatter correction (Hao Yan, Xuanqin Mou, Shaojie Tang, Qiong Xu and Maria Zankl). , Projection movingbeam stop array, 2010, Phys.Med.Biol.). A beam stop array is a shielding plate with regularly arranged shielding points. The X-rays are irradiated onto the object through the shielding plate, and the X-rays passing through the shielding points on the shielding plate will be shielded. If a value appears at the corresponding position of the detector occlusion point, this value is the scattering value. Due to the use of BSA, the resulting projection data is incomplete and needs to be recovered using a consistency condition. This document adopts Patch's derivation formula based on three-generation helical CT [Patch (2002)], combined with the spatial interpolation method, and proposes a PC_SI method for the field of scatter correction. However, as mentioned above, since the z-axis of helical CT is constantly changing, this formula is difficult to apply specifically, so this document refers to the practice of Defrise et al., and uses the ordinate of the detector to approximately replace the z-coordinate. However, this approximation will lead to a decrease in accuracy. From the geometrical point of view, it is only necessary to set z=0, and the spiral CT can be degenerated into a circular orbit CBCT, but the formula in the literature is to derive the z derivative. When z is a constant, the derivative is meaningless, so the method used in this document Not suitable for circular orbit CBCT.

因此,针对现有技术不足,提供一种圆轨道下不完全CT投影数据的恢复方法以克服现有技术不足甚为必要。Therefore, in view of the deficiencies of the prior art, it is necessary to provide a method for restoring incomplete CT projection data under a circular orbit to overcome the deficiencies of the prior art.

发明内容SUMMARY OF THE INVENTION

本发明所要解决的技术问题是提供一种圆轨道下不完全CT投影数据的恢复方法,和传统的空间插值方法相比,本发明的方法能够很好地恢复丢失信息中的高频成分,进而大幅度减少图像中的条形伪影,显著提高图像质量。在圆轨道下使用Patch的方法需要先使用探测器的纵坐标来近似代替z坐标,而本方法采用的一致性条件是基于圆轨道CBCT几何结构推导得到的,能够避免因为近似z坐标而导致的精度下降,并且Patch的方法中z轴在不停变化,使得此方法具体应用困难。与之相反,本发明方法只是一个简单的一阶求导公式,实现简单,具有可应用性。The technical problem to be solved by the present invention is to provide a method for restoring incomplete CT projection data under a circular orbit. Compared with the traditional spatial interpolation method, the method of the present invention can restore the high-frequency components in the lost information well, and further Dramatically reduces striping artifacts in images, significantly improving image quality. The method of using Patch in the circular orbit needs to use the ordinate of the detector to approximate the z coordinate first, and the consistency condition adopted by this method is derived based on the CBCT geometry of the circular orbit, which can avoid the approximation caused by the z coordinate. The accuracy decreases, and the z-axis in the Patch method is constantly changing, making the specific application of this method difficult. On the contrary, the method of the present invention is only a simple first-order derivation formula, which is simple to implement and has applicability.

本发明采用如下技术方案:The present invention adopts following technical scheme:

一种圆轨道下不完全CT投影数据的恢复方法,通过如下步骤进行:A method for recovering incomplete CT projection data under a circular orbit is carried out through the following steps:

第一步,获得坏像素点位置信息;The first step is to obtain the location information of bad pixels;

第二步,对每个投影中的坏像素点使用水平一维三次样条插值方法进行线性插值,以恢复坏像素点的低频信息,经过此步骤得到每个角度下的初始化投影;In the second step, the horizontal one-dimensional cubic spline interpolation method is used to linearly interpolate the bad pixels in each projection to restore the low-frequency information of the bad pixels. After this step, the initial projection at each angle is obtained;

第三步,令N=1,进入第四步,The third step, let N=1, enter the fourth step,

第四步,令M=N,M为当前迭代次数;The fourth step, let M=N, M is the current iteration number;

将每个角度下的初始化投影进行二维傅里叶变换,获得每个角度下的初始化投影在频率域中的频域数据;Perform a two-dimensional Fourier transform on the initialized projection at each angle to obtain the frequency domain data of the initialized projection at each angle in the frequency domain;

第五步,分别将前后两个相邻投影的频率域形式带入到一致性条件对应的公式中计算投影结果;In the fifth step, the frequency domain forms of the two adjacent projections before and after are respectively brought into the formula corresponding to the consistency condition to calculate the projection result;

第六步,对加权后的投影进行傅里叶反变换,得到第M次迭代终值;The sixth step is to perform inverse Fourier transform on the weighted projection to obtain the final value of the Mth iteration;

第七步,判断当前迭代次数M是否大于等于迭代终止次数S,如果是,则进入第九步;否则进入第八步;The seventh step is to judge whether the current iteration number M is greater than or equal to the iteration termination number S, and if so, enter the ninth step; otherwise, enter the eighth step;

第八步,将第M次迭代终值作为初始化投影,令N=N+1,进入第三步;The eighth step, take the final value of the Mth iteration as the initialization projection, let N=N+1, and enter the third step;

第九步;以第M次迭代终值为最终恢复的数据结果。The ninth step: take the final value of the M-th iteration as the final recovered data result.

优选的,上述的圆轨道下不完全CT投影数据的恢复方法,Preferably, the recovery method of incomplete CT projection data under the above-mentioned circular orbit,

第五步中使用的一致性条件,具体形式如下:The consistency condition used in the fifth step is as follows:

Figure BDA0001333901850000031
Figure BDA0001333901850000031

其中,其中ρ是射线源到旋转中心的距离,d是旋转中心到探测器的距离,(α12)是平板探测器坐标,θ是方位角。where ρ is the distance from the ray source to the rotation center, d is the distance from the rotation center to the detector, (α 1 , α 2 ) is the flat panel detector coordinate, and θ is the azimuth angle.

优选的,上述第五步具体是使用相同的权重因子ω对前后两个投影结果进行加权,得到加权更新后的投影结果。Preferably, the fifth step above specifically uses the same weighting factor ω to weight the two projection results before and after, to obtain a weighted and updated projection result.

优选的,上述的圆轨道下不完全CT投影数据的恢复方法,Preferably, the recovery method of incomplete CT projection data under the above-mentioned circular orbit,

第五步,具体使用一致性条件即公式(28)对每一个初始化投影的频域数据进行恢复,公式(28)的具体形式如下:The fifth step is to restore the frequency domain data of each initialized projection using the consistency condition, that is, formula (28). The specific form of formula (28) is as follows:

u*(θ+dθ)=u*(θ)+dθ·U,k2≠0……公式(28);u * (θ+dθ)=u * (θ)+dθ·U,k 2 ≠0…Formula (28);

Figure BDA0001333901850000032
Figure BDA0001333901850000032

其中,u*(θ+dθ)和u*(θ)分别表示θ+dθ和θ角度下初始化投影的频域数据,u*(θ)为已知量,u*(θ+dθ)是要求的未知量,dθ代表相邻两个角度的间隔,k1、k2是平板探测器的频率域坐标,i代表复数的虚部;Among them, u*(θ+dθ) and u*(θ) represent the frequency domain data of the initialized projection at θ+dθ and θ angles, respectively, u*(θ) is a known quantity, and u*(θ+dθ) is the requirement The unknown quantity of , dθ represents the interval between two adjacent angles, k 1 and k 2 are the frequency domain coordinates of the flat panel detector, and i represents the imaginary part of the complex number;

公式(28)的具体使用方法如下:The specific usage of formula (28) is as follows:

首先将u*(θ+dθ)代入到公式(28)中计算u*(θ):First, substitute u*(θ+dθ) into formula (28) to calculate u*(θ):

u*(θ)=u*(θ+dθ)-dθ·U;u * (θ)=u * (θ+dθ)-dθ·U;

然后用u*(θ-dθ)代替公式(28)中的u*(θ),u*(θ)代替公式(28)中的u*(θ+dθ)计算得到:Then use u*(θ-dθ) to replace u*(θ) in formula (28), and u*(θ) to replace u*(θ+dθ) in formula (28) to calculate:

u*(θ)=u*(θ-dθ)+dθ·U;u * (θ)=u * (θ-dθ)+dθ·U;

由于选取前后两个相邻投影频域数据的任一个都可以实现投影数据恢复,为了平衡恢复效果,使用权重因子ω,ω∈[0,1],对上述两个公式的计算结果进行加权平均,得到加权后的结果uR*(θ),即:Since the projection data recovery can be achieved by selecting any one of the two adjacent projection frequency domain data before and after, in order to balance the recovery effect, the weighting factor ω, ω∈[0,1] is used to perform a weighted average of the calculation results of the above two formulas , the weighted result u R* (θ) is obtained, namely:

uR*(θ)=ω(u*(θ-dθ)+dθ·U)+(1-ω)(u*(θ+dθ)-dθ·U),k2≠0;u R* (θ)=ω(u * (θ-dθ)+dθ·U)+(1-ω)(u * (θ+dθ)-dθ·U),k 2 ≠0;

由于公式(28)只有在k2≠0时才有效,所以对频域数据u*(θ)中k2=0的部分,直接使用θ角度下初始化投影频域数据中的相应部分来补偿;Since formula (28) is only valid when k 2 ≠0, for the part where k 2 =0 in the frequency domain data u*(θ), directly use the corresponding part in the initialized projected frequency domain data under the angle θ to compensate;

第六步,具体是对uR*(θ)进行傅里叶反变换,得到第M次迭代终值uR(θ);The sixth step is to perform inverse Fourier transform on u R* (θ) to obtain the final value u R (θ) of the Mth iteration;

uR(θ)=F-1(uR*(θ))……(1);u R (θ)=F -1 (u R* (θ))...(1);

其中公式(1)中uR(θ)表示θ角度下加权后的频域数据uR*(θ)的反傅里叶变换结果,F-1代表傅里叶反变换;Wherein u R (θ) in formula (1) represents the inverse Fourier transform result of the weighted frequency domain data u R* (θ) under the angle θ, and F -1 represents the inverse Fourier transform;

第七步,具体是判断当前迭代次数M是否大于等于迭代终止次数S,如果是,则进入第九步;否则进入第八步;The seventh step is to determine whether the current iteration number M is greater than or equal to the iteration termination number S, and if so, enter the ninth step; otherwise, enter the eighth step;

第八步,具体是将第M次迭代终值uR(θ)作为初始化投影,令N=N+1,进入第三步;The eighth step is to use the final value u R (θ) of the Mth iteration as the initialization projection, let N=N+1, and enter the third step;

第九步,具体是以第M次迭代终值uR(θ)为最终恢复的数据结果。The ninth step, specifically, takes the final value u R (θ) of the Mth iteration as the final recovered data result.

优选的,上述的圆轨道下不完全CT投影数据的恢复方法,第一步,获得坏像素点位置信息具体是将实际采集的投影数据转换成对应的弦图,即:弦图的x轴为投影的像素点位置,y轴为该像素位置对应的投影角度位置;对某个探测器单元存在损坏的条件下,则在弦图中对应的像素位置上会出现一条竖线,根据此理论在弦图中获得坏像素点的位置信息,并记录坏像素点位置。Preferably, in the above-mentioned method for recovering incomplete CT projection data under a circular orbit, the first step is to obtain the position information of bad pixels, specifically converting the actually collected projection data into a corresponding chord diagram, that is, the x-axis of the chord diagram is The pixel position of the projection, the y-axis is the projection angle position corresponding to the pixel position; under the condition that a certain detector unit is damaged, a vertical line will appear on the corresponding pixel position in the chord diagram. According to this theory, in The position information of bad pixels is obtained from the chord diagram, and the positions of bad pixels are recorded.

其中本方法中应用的公式(28)推导过程如下:The derivation process of formula (28) applied in this method is as follows:

由于本发明提出的一致性条件是基于约翰方程,所以先对约翰方程进行简单介绍和分析。FJohn于1938年在文献《THE ULTRAHYPERBOLIC DIFFERENTIAL EQUATION WITHFOURINDEPENDENT VARIABLES》中提出一个超双曲线二阶偏微分方程,作为线积分的一致性条件,具体形式如下:Since the consistency condition proposed by the present invention is based on the John equation, the John equation is briefly introduced and analyzed first. In 1938, FJohn proposed a hyperbolic second-order partial differential equation in the document "THE ULTRAHYPERBOLIC DIFFERENTIAL EQUATION WITHFOURINDEPENDENT VARIABLES" as the consistency condition of the line integral. The specific form is as follows:

Figure BDA0001333901850000041
Figure BDA0001333901850000041

公式中i,j为下标;In the formula, i and j are subscripts;

u是函数f穿过ε和η两点的线积分:u is the line integral of the function f through ε and η:

u(ε;η)=∫f(ε+t(η-ε))dt (3);u(ε; η)=∫f(ε+t(η-ε))dt (3);

以下为圆轨道下基于约翰方程的一致性条件推导步骤:The following are the derivation steps of the consistency condition based on the John equation under the circular orbit:

首先采用圆轨道CBCT坐标参数对ε和η进行参数化:First, ε and η are parameterized using the circular orbit CBCT coordinate parameters:

Figure BDA0001333901850000042
Figure BDA0001333901850000042

Figure BDA0001333901850000043
Figure BDA0001333901850000043

其中ρ是射线源到旋转中心的距离,d是旋转中心到探测器的距离,(α12)是平板探测器坐标,θ是方位角。where ρ is the distance from the ray source to the center of rotation, d is the distance from the center of rotation to the detector, (α 1 , α 2 ) are the coordinates of the flat panel detector, and θ is the azimuth angle.

分别求ε和η的每个分量求偏导:Find the partial derivatives for each component of ε and η separately:

Figure BDA0001333901850000051
Figure BDA0001333901850000051

Figure BDA0001333901850000052
Figure BDA0001333901850000052

Figure BDA0001333901850000053
Figure BDA0001333901850000053

Figure BDA0001333901850000054
Figure BDA0001333901850000054

Figure BDA0001333901850000055
Figure BDA0001333901850000055

Figure BDA0001333901850000056
Figure BDA0001333901850000056

为书写方便,定义如下算子:For the convenience of writing, the following operators are defined:

Figure BDA0001333901850000057
Figure BDA0001333901850000057

Figure BDA0001333901850000058
Figure BDA0001333901850000058

Figure BDA0001333901850000059
Figure BDA0001333901850000059

Figure BDA00013339018500000510
Figure BDA00013339018500000510

Figure BDA00013339018500000511
Figure BDA00013339018500000511

把公式(6)-(11)代入到公式(2)中得到:Substitute formulas (6)-(11) into formula (2) to get:

Figure BDA00013339018500000512
Figure BDA00013339018500000512

Figure BDA00013339018500000513
Figure BDA00013339018500000513

Figure BDA00013339018500000514
Figure BDA00013339018500000514

由公式(13)-(16)对H1,H2,

Figure BDA00013339018500000515
的定义,代入到公式(17)-(19)中,得到:By formulas (13)-(16) for H 1 , H 2 ,
Figure BDA00013339018500000515
The definition of , and substituting it into formulas (17)-(19), we get:

Figure BDA0001333901850000061
Figure BDA0001333901850000061

从上述公式中,得出:From the above formula, we get:

Figure BDA0001333901850000062
Figure BDA0001333901850000062

对上面公式(21)移项合并同类项:For the above formula (21), shift items and merge similar items:

Figure BDA0001333901850000063
Figure BDA0001333901850000063

Figure BDA0001333901850000064
Figure BDA0001333901850000064

求解公式(23),最终推导结果如下:Solving formula (23), the final derivation result is as follows:

Figure BDA0001333901850000065
Figure BDA0001333901850000065

公式(24)即为我们推导的圆轨道下基于约翰方程的一致性条件。对公式(24)进行傅里叶变换,得到对应的频率域的表达形式如下:Equation (24) is the consistency condition based on the John equation under the circular orbit we deduce. Fourier transform is performed on formula (24), and the corresponding expression in the frequency domain is obtained as follows:

Figure BDA0001333901850000066
Figure BDA0001333901850000066

在上述公式中,uθ *表示u对θ求导后的傅里叶变换,(k1,k2)是(α12)的对应频率域形式。In the above formula, u θ * represents the Fourier transform of u after derivation of θ, and (k 1 , k 2 ) is the corresponding frequency domain form of (α 1 , α 2 ).

为了方便表达和应用,将公式(25)右边式子用U代替:For the convenience of expression and application, the right-hand side of formula (25) is replaced by U:

Figure BDA0001333901850000071
Figure BDA0001333901850000071

考虑到求导公式的一般定义,可以将公式(25)改下成如下形式:Considering the general definition of the derivation formula, formula (25) can be changed into the following form:

Figure BDA0001333901850000072
Figure BDA0001333901850000072

即:which is:

u*(θ+dθ)=u*(θ)+dθ·U,k2≠0 (28);u * (θ+dθ)=u * (θ)+dθ·U,k 2 ≠0 (28);

公式(28)是本发明改进后的一致性条件,本发明使用公式(28)进行投影数据的恢复。The formula (28) is the improved consistency condition of the present invention, and the present invention uses the formula (28) to recover the projection data.

传统的空间插值技术使用同一个角度的投影的其他完整像素点数据来恢复丢失的数据,所使用重建信息只局限于丢失数据所在角度的投影,所以插值方法只能恢复丢失数据中的低频信息,对丢失的高频信息基本没有恢复。本发明方法应用的一致性条件,即不仅基于同一个投影中的数据来恢复丢失数据,还利用相邻角度投影中的对应位置数据信息。因此本发采用一致性条件可以很好地恢复丢失数据中的高频信息,能基本消除使用插值技术获得的重建图像中的条形伪影,显著提高图像质量。The traditional spatial interpolation technology uses other complete pixel data projected at the same angle to restore the lost data, and the reconstruction information used is limited to the projection at the angle of the lost data, so the interpolation method can only restore the low-frequency information in the lost data, There is basically no recovery of the lost high-frequency information. The consistency condition applied by the method of the present invention is not only to recover the lost data based on the data in the same projection, but also to use the corresponding position data information in the adjacent angle projections. Therefore, by adopting the consistency condition, the present invention can restore the high-frequency information in the lost data well, basically eliminate the stripe artifacts in the reconstructed image obtained by using the interpolation technology, and significantly improve the image quality.

附图说明Description of drawings

图1为实例中圆轨道CBCT扫描几何示意图及参数示意图,图中,ρ是射线源到旋转中心的距离,d是旋转中心到探测器的距离,(α12)是平板探测器坐标,θ是方位角。Figure 1 is a schematic diagram of the geometry and parameters of the circular orbit CBCT scanning in the example. In the figure, ρ is the distance from the ray source to the rotation center, d is the distance from the rotation center to the detector, and (α 1 , α 2 ) are the coordinates of the flat panel detector. , θ is the azimuth angle.

图2是投影数据对应的弦图,横坐标x轴为像素点位置,纵坐标y轴为投影角度。Figure 2 is the chord diagram corresponding to the projection data, the x-axis of the abscissa is the pixel position, and the y-axis of the ordinate is the projection angle.

图3为本发明的数据恢复方法流程图。FIG. 3 is a flow chart of the data recovery method of the present invention.

图4为三种方法的仿真结果示意图,图4(a)为理想图像,图4(b)为未经任何处理的原始图像,图4(c)为线性差值获得的结果,图4(d)、(e)、(f)为本发明方法一次迭代、二次迭代、三次迭代后的结果图,图(g)、(h)、(i)为Patch方法一次迭代、二次迭代、三次迭代后的结果图。Figure 4 is a schematic diagram of the simulation results of the three methods, Figure 4(a) is the ideal image, Figure 4(b) is the original image without any processing, Figure 4(c) is the result obtained by linear difference, Figure 4( d), (e), (f) are the results of the first iteration, the second iteration, and the third iteration of the method of the present invention. Figures (g), (h), and (i) are the first iteration, the second iteration, and the third iteration of the Patch method. Result plot after three iterations.

具体实施方式Detailed ways

实施例1。Example 1.

一种圆轨道下不完全CT投影数据的恢复方法,通过如下步骤进行:A method for recovering incomplete CT projection data under a circular orbit is carried out through the following steps:

第一步,获得坏像素点位置信息,将实际采集的投影数据转换成对应的弦图,即:弦图的x轴为投影的像素点位置,y轴为该像素位置对应的投影角度位置;对某个探测器单元存在损坏的条件下,则在弦图中对应的像素位置上会出现一条竖线,根据此理论在弦图中获得坏像素点的位置信息,并记录坏像素点位置。The first step is to obtain the bad pixel position information, and convert the actually collected projection data into the corresponding chord diagram, that is: the x-axis of the chord diagram is the projected pixel position, and the y-axis is the projection angle position corresponding to the pixel position; Under the condition of damage to a certain detector unit, a vertical line will appear on the corresponding pixel position in the chord diagram. According to this theory, the position information of the bad pixel point is obtained in the chord diagram, and the position of the bad pixel point is recorded.

第二步,对每个投影中的坏像素点使用水平一维三次样条插值方法进行线性插值,以恢复坏像素点的低频信息,经过此步骤得到每个角度下的初始化投影。In the second step, the horizontal one-dimensional cubic spline interpolation method is used to linearly interpolate the bad pixels in each projection to restore the low-frequency information of the bad pixels. After this step, the initial projection at each angle is obtained.

第三步,令N=1,进入第四步,The third step, let N=1, enter the fourth step,

第四步,令M=N,M为当前迭代次数;The fourth step, let M=N, M is the current iteration number;

将每个角度下的初始化投影进行二维傅里叶变换,获得每个角度下的初始化投影在频率域中的频域数据;Perform a two-dimensional Fourier transform on the initialized projection at each angle to obtain the frequency domain data of the initialized projection at each angle in the frequency domain;

第五步,分别将前后两个相邻投影的频率域形式带入到一致性条件对应的公式中计算投影结果;具体是使用相同的权重因子ω对前后两个投影结果进行加权,得到加权更新后的投影结果。In the fifth step, the frequency domain forms of the two adjacent projections before and after are respectively brought into the formula corresponding to the consistency condition to calculate the projection result; specifically, the same weight factor ω is used to weight the two projection results before and after, and the weighted update is obtained. The projection result after.

第五步中使用的一致性条件,形式如下:The consistency condition used in the fifth step is of the form:

Figure BDA0001333901850000081
Figure BDA0001333901850000081

其中,其中ρ是射线源到旋转中心的距离,d是旋转中心到探测器的距离,(α12)是平板探测器坐标,θ是方位角,如图1所示。where ρ is the distance from the ray source to the center of rotation, d is the distance from the center of rotation to the detector, (α 1 , α 2 ) is the coordinates of the flat panel detector, and θ is the azimuth angle, as shown in Figure 1.

第五步,可具体使用一致性条件即公式(28)对每一个初始化投影的频域数据进行恢复,公式(28)的具体形式如下:In the fifth step, the consistency condition, that is, formula (28), can be used to restore the frequency domain data of each initialized projection. The specific form of formula (28) is as follows:

u*(θ+dθ)=u*(θ)+dθ·U,k2≠0……公式(28);u * (θ+dθ)=u * (θ)+dθ·U,k 2 ≠0…Formula (28);

Figure BDA0001333901850000082
Figure BDA0001333901850000082

其中,u*(θ+dθ)和u*(θ)分别表示θ+dθ和θ角度下初始化投影的频域数据,u*(θ)为已知量,u*(θ+dθ)是要求的未知量,dθ代表相邻两个角度的间隔,k1、k2是平板探测器的频率域坐标,i代表复数的虚部;Among them, u*(θ+dθ) and u*(θ) represent the frequency domain data of the initialized projection at θ+dθ and θ angles, respectively, u*(θ) is a known quantity, and u*(θ+dθ) is the requirement The unknown quantity of , dθ represents the interval between two adjacent angles, k 1 and k 2 are the frequency domain coordinates of the flat panel detector, and i represents the imaginary part of the complex number;

公式(28)的具体使用方法如下:The specific usage of formula (28) is as follows:

首先将u*(θ+dθ)代入到公式(28)中计算u*(θ):First, substitute u*(θ+dθ) into formula (28) to calculate u*(θ):

u*(θ)=u*(θ+dθ)-dθ·U;u * (θ)=u * (θ+dθ)-dθ·U;

然后用u*(θ-dθ)代替公式(28)中的u*(θ),u*(θ)代替公式(28)中的u*(θ+dθ)计算得到:Then use u*(θ-dθ) to replace u*(θ) in formula (28), and u*(θ) to replace u*(θ+dθ) in formula (28) to calculate:

u*(θ)=u*(θ-dθ)+dθ·U;u * (θ)=u * (θ-dθ)+dθ·U;

由于选取前后两个相邻投影频域数据的任一个都可以实现投影数据恢复,为了平衡恢复效果,使用权重因子ω,ω∈[0,1],对上述两个公式的计算结果进行加权平均,得到加权后的结果uR*(θ),即:Since the projection data recovery can be achieved by selecting any one of the two adjacent projection frequency domain data before and after, in order to balance the recovery effect, the weighting factor ω, ω∈[0,1] is used to perform a weighted average of the calculation results of the above two formulas , the weighted result u R* (θ) is obtained, namely:

uR*(θ)=ω(u*(θ-dθ)+dθ·U)+(1-ω)(u*(θ+dθ)-dθ·U),k2≠0;u R* (θ)=ω(u * (θ-dθ)+dθ·U)+(1-ω)(u * (θ+dθ)-dθ·U),k 2 ≠0;

由于公式(28)只有在k2≠0时才有效,所以对频域数据u*(θ)中k2=0的部分,直接使用θ角度下初始化投影频域数据中的相应部分来补偿;Since formula (28) is only valid when k 2 ≠0, for the part where k 2 =0 in the frequency domain data u*(θ), directly use the corresponding part in the initialized projected frequency domain data under the angle θ to compensate;

第六步,具体是对uR*(θ)进行傅里叶反变换,得到第M次迭代终值uR(θ);The sixth step is to perform inverse Fourier transform on u R* (θ) to obtain the final value u R (θ) of the Mth iteration;

uR(θ)=F-1(uR*(θ))……(1);u R (θ)=F -1 (u R* (θ))...(1);

其中公式(1)中uR(θ)表示θ角度下加权后的频域数据uR*(θ)的反傅里叶变换结果,F-1代表傅里叶反变换。In the formula (1), u R (θ) represents the inverse Fourier transform result of the weighted frequency domain data u R* (θ) at the angle θ, and F −1 represents the inverse Fourier transform.

第七步,具体是判断当前迭代次数M是否大于等于迭代终止次数S,如果是,则进入第九步;否则进入第八步。The seventh step is to determine whether the current iteration number M is greater than or equal to the iteration termination number S, and if so, enter the ninth step; otherwise, enter the eighth step.

第八步,具体是将第M次迭代终值uR(θ)作为初始化投影,令N=N+1,进入第三步。The eighth step is to take the final value u R (θ) of the Mth iteration as the initialization projection, and let N=N+1, and enter the third step.

第九步,具体是以第M次迭代终值uR(θ)为最终恢复的数据结果。The ninth step, specifically, takes the final value u R (θ) of the Mth iteration as the final recovered data result.

其中本方法中应用的公式(28)推导过程如下:The derivation process of formula (28) applied in this method is as follows:

由于本发明提出的一致性条件是基于约翰方程,所以先对约翰方程进行简单介绍和分析。F John于1938年在文献《THE ULTRAHYPERBOLIC DIFFERENTIAL EQUATION WITHFOUR INDEPENDENT VARIABLES》中提出一个超双曲线二阶偏微分方程,作为线积分的一致性条件,具体形式如下:Since the consistency condition proposed by the present invention is based on the John equation, the John equation is briefly introduced and analyzed first. In 1938, F John proposed a hyperbolic second-order partial differential equation in the document "THE ULTRAHYPERBOLIC DIFFERENTIAL EQUATION WITHFOUR INDEPENDENT VARIABLES" as the consistency condition of the line integral. The specific form is as follows:

Figure BDA0001333901850000091
Figure BDA0001333901850000091

u是函数f穿过ε和η两点的线积分:u is the line integral of the function f through ε and η:

u(ε;η)=∫f(ε+t(η-ε))dt (3);u(ε; η)=∫f(ε+t(η-ε))dt (3);

以下为圆轨道下基于约翰方程的一致性条件推导步骤:The following are the derivation steps of the consistency condition based on the John equation under the circular orbit:

首先采用圆轨道CBCT坐标参数对ε和η进行参数化:First, ε and η are parameterized using the circular orbit CBCT coordinate parameters:

Figure BDA0001333901850000092
Figure BDA0001333901850000092

Figure BDA0001333901850000093
Figure BDA0001333901850000093

其中ρ是射线源到旋转中心的距离,d是旋转中心到探测器的距离,(α12)是平板探测器坐标,θ是方位角。where ρ is the distance from the ray source to the center of rotation, d is the distance from the center of rotation to the detector, (α 1 , α 2 ) are the coordinates of the flat panel detector, and θ is the azimuth angle.

分别求ε和η的每个分量求偏导:Find the partial derivatives for each component of ε and η separately:

Figure BDA0001333901850000101
Figure BDA0001333901850000101

Figure BDA0001333901850000102
Figure BDA0001333901850000102

Figure BDA0001333901850000103
Figure BDA0001333901850000103

Figure BDA0001333901850000104
Figure BDA0001333901850000104

Figure BDA0001333901850000105
Figure BDA0001333901850000105

Figure BDA0001333901850000106
Figure BDA0001333901850000106

为书写方便,定义如下算子:For the convenience of writing, the following operators are defined:

Figure BDA0001333901850000107
Figure BDA0001333901850000107

Figure BDA0001333901850000108
Figure BDA0001333901850000108

Figure BDA0001333901850000109
Figure BDA0001333901850000109

Figure BDA00013339018500001010
Figure BDA00013339018500001010

Figure BDA00013339018500001011
Figure BDA00013339018500001011

把公式(6)-(11)代入到公式(2)中得到:Substitute formulas (6)-(11) into formula (2) to get:

Figure BDA00013339018500001012
Figure BDA00013339018500001012

Figure BDA00013339018500001013
Figure BDA00013339018500001013

Figure BDA00013339018500001014
Figure BDA00013339018500001014

由公式(13)-(16)对H1,H2,

Figure BDA00013339018500001015
的定义,代入到公式(17)-(19)中,得到:By formulas (13)-(16) for H 1 , H 2 ,
Figure BDA00013339018500001015
The definition of , and substituting it into formulas (17)-(19), we get:

Figure BDA0001333901850000111
Figure BDA0001333901850000111

从上述公式中,得出:From the above formula, we get:

Figure BDA0001333901850000112
Figure BDA0001333901850000112

对上面公式(21)移项合并同类项:For the above formula (21), shift items and merge similar items:

Figure BDA0001333901850000113
Figure BDA0001333901850000113

Figure BDA0001333901850000114
Figure BDA0001333901850000114

求解公式(23),最终推导结果如下:Solving formula (23), the final derivation result is as follows:

Figure BDA0001333901850000115
Figure BDA0001333901850000115

公式(24)即为我们推导的圆轨道下基于约翰方程的一致性条件。对公式(24)进行傅里叶变换,得到对应的频率域的表达形式如下:Equation (24) is the consistency condition based on the John equation under the circular orbit we deduce. Fourier transform is performed on formula (24), and the corresponding expression in the frequency domain is obtained as follows:

Figure BDA0001333901850000116
Figure BDA0001333901850000116

在上述公式中,uθ *表示u对θ求导后的傅里叶变换,(k1,k2)是(α12)的对应频率域形式。In the above formula, u θ * represents the Fourier transform of u after derivation of θ, and (k 1 , k 2 ) is the corresponding frequency domain form of (α 1 , α 2 ).

为了方便表达和应用,将公式(25)右边式子用U代替:For the convenience of expression and application, the right-hand side of formula (25) is replaced by U:

Figure BDA0001333901850000121
Figure BDA0001333901850000121

考虑到求导公式的一般定义,可以将公式(25)改下成如下形式:Considering the general definition of the derivation formula, formula (25) can be changed into the following form:

Figure BDA0001333901850000122
Figure BDA0001333901850000122

即:which is:

u*(θ+dθ)=u*(θ)+dθ·U,k2≠0 (28);u * (θ+dθ)=u * (θ)+dθ·U,k 2 ≠0 (28);

公式(28)是本发明改进后的一致性条件,本发明使用公式(28)进行投影数据的恢复。The formula (28) is the improved consistency condition of the present invention, and the present invention uses the formula (28) to recover the projection data.

传统的空间插值技术使用同一个角度的投影的其他完整像素点数据来恢复丢失的数据,所使用重建信息只局限于丢失数据所在角度的投影,所以插值方法只能恢复丢失数据中的低频信息,对丢失的高频信息基本没有恢复。本发明方法应用的一致性条件,即不仅基于同一个投影中的数据来恢复丢失数据,还利用相邻角度投影中的对应位置数据信息。因此本发采用一致性条件可以很好地恢复丢失数据中的高频信息,能基本消除使用插值技术获得的重建图像中的条形伪影,显著提高图像质量。The traditional spatial interpolation technology uses other complete pixel data projected at the same angle to restore the lost data, and the reconstruction information used is limited to the projection at the angle of the lost data, so the interpolation method can only restore the low-frequency information in the lost data, There is basically no recovery of the lost high-frequency information. The consistency condition applied by the method of the present invention is not only to recover the lost data based on the data in the same projection, but also to use the corresponding position data information in the adjacent angle projections. Therefore, by adopting the consistency condition, the present invention can restore the high-frequency information in the lost data well, basically eliminate the stripe artifacts in the reconstructed image obtained by using the interpolation technology, and significantly improve the image quality.

和传统的空间插值方法相比,本发明的方法能够很好地恢复丢失信息中的高频成分,进而大幅度减少图像中的条形伪影,显著提高图像质量。在圆轨道下使用Patch的方法需要先使用探测器的纵坐标来近似代替z坐标,而本方法采用的一致性条件是基于圆轨道CBCT几何结构推导得到的,能够避免因为近似z坐标而导致的精度下降,并且Patch的方法中z轴在不停变化,使得此方法具体应用困难。与之相反,本发明方法只是一个简单的一阶求导公式,实现简单,具有可应用性。Compared with the traditional spatial interpolation method, the method of the present invention can restore the high-frequency components in the lost information well, thereby greatly reducing the stripe artifacts in the image and significantly improving the image quality. The method of using Patch in the circular orbit needs to use the ordinate of the detector to approximate the z coordinate first, and the consistency condition adopted by this method is derived based on the CBCT geometry of the circular orbit, which can avoid the approximation caused by the z coordinate. The accuracy decreases, and the z-axis in the Patch method is constantly changing, making the specific application of this method difficult. On the contrary, the method of the present invention is only a simple first-order derivation formula, which is simple to implement and has applicability.

实施例2。Example 2.

本实例分别使用传统空间插值方法、Patch的方法(背景技术中有详细说明)、以及本方法对一组不完整的投影数据进行恢复,不完整的投影数据由一台平板探测器具有坏像素点的CBCT采集得到,其CT扫描参数具体如下:射线源到旋转中心的距离ρ为500mm,旋转中心到探测器的距离d为500mm,平板探测器大小为850×200mm,像素点尺寸为1×1mm,围绕物体旋转360度采集1080个投影。This example uses the traditional spatial interpolation method, the Patch method (detailed in the background art), and the present method to recover a set of incomplete projection data. The incomplete projection data is generated by a flat panel detector with bad pixels. The CT scanning parameters are as follows: the distance ρ from the ray source to the rotation center is 500mm, the distance d from the rotation center to the detector is 500mm, the size of the flat panel detector is 850×200mm, and the pixel size is 1×1mm , rotate 360 degrees around the object to collect 1080 projections.

本发明技术的方法具体实施方式如下:The specific embodiment of the method of the present technology is as follows:

第一步,获得坏像素点位置信息。将实际采集的投影数据转换成对应的弦图,如图2所示。图中黑色竖线所对应的横坐标即为坏像素点位置,通过此方式可以找到所有坏像素点的位置信息,并记录。The first step is to obtain the location information of bad pixels. Convert the actual collected projection data into the corresponding chord diagram, as shown in Figure 2. The abscissa corresponding to the black vertical line in the figure is the position of the bad pixel. In this way, the position information of all the bad pixels can be found and recorded.

第二步,对每一个投影中的坏像素点使用水平一维三次样条插值方法线性插值,并将插值后的结果作为该坏像素点的初始值。In the second step, use the horizontal one-dimensional cubic spline interpolation method to linearly interpolate the bad pixels in each projection, and use the interpolation result as the initial value of the bad pixel.

第三步,将插值后的投影以及前后两个相邻投影进行傅里叶变换,得到对应的频率域形式,以方便应用本发明推导出的一致性条件。In the third step, Fourier transform is performed on the interpolated projection and the two adjacent projections before and after to obtain the corresponding frequency domain form, so as to facilitate the application of the consistency condition derived by the present invention.

第四步,将上一步得到的前后两个相邻投影的频率域形式代入到公式(28)中,恢复中间投影的丢失信息(尤其是高频成分信息),并使用同一个权重因子ω对两个结果进行加权,由于没有更多的约束条件表明前后两个相邻投影哪个对中间投影的贡献更大,所以本实例令ω=0.5。然后用计算加权后的结果修正该投影初始值。The fourth step is to substitute the frequency domain form of the two adjacent projections obtained in the previous step into formula (28) to restore the lost information (especially the high-frequency component information) of the intermediate projection, and use the same weight factor ω to The two results are weighted, and since there are no more constraints to indicate which of the two adjacent projections before and after has a greater contribution to the middle projection, this example sets ω = 0.5. The projected initial value is then corrected with the calculated weighted result.

第五步,对修正后的值进行傅里叶反变换,得到第一次迭代最终值。The fifth step is to perform inverse Fourier transform on the corrected value to obtain the final value of the first iteration.

第六步,将上一次迭代的最终结果作为初值,重复上面第三步到第五步,本实例共迭代三次。In the sixth step, the final result of the previous iteration is used as the initial value, and the third to fifth steps above are repeated. This example is iterated three times in total.

具体步骤流程图如图3所示。The flow chart of the specific steps is shown in Figure 3.

仿真结果如图4所示,图4(a)为理想图像,图4(b)为未经任何处理的原始图像,图4(c)为线性差值获得的结果,图4(d)、(e)、(f)为本发明方法一次迭代、二次迭代、三次迭代后的结果图,图(g)、(h)、(i)为Patch方法一次迭代、二次迭代、三次迭代后的结果图。The simulation results are shown in Fig. 4. Fig. 4(a) is the ideal image, Fig. 4(b) is the original image without any processing, Fig. 4(c) is the result obtained by linear difference, Fig. 4(d), (e) and (f) are the results of the first iteration, the second iteration and the third iteration of the method of the present invention. Figures (g), (h) and (i) are the first iteration, the second iteration and the third iteration of the Patch method. result graph.

使用绝对误差(Absolute Error)对重建图像进行分析:Analyze the reconstructed image using Absolute Error:

Figure BDA0001333901850000131
Figure BDA0001333901850000131

其中N为像素点个数,ri是重建图像第i个点对应的灰度值,ti是理想图像第i个点对应的灰度值。结果如下:where N is the number of pixels, ri is the gray value corresponding to the ith point of the reconstructed image, and t i is the gray value corresponding to the ith point of the ideal image. The result is as follows:

方法method 第101层的MAEMAE at Tier 101 线性空间插值Linear Space Interpolation 1.3368×10<sup>+6</sup>1.3368×10<sup>+6</sup> Patch方法一次迭代Patch method one iteration 3.0641×10<sup>+5</sup>3.0641×10<sup>+5</sup> Patch方法二次迭代Patch method second iteration 3.1477×10<sup>+5</sup>3.1477×10<sup>+5</sup> Patch方法三次迭代Patch method three iterations 2.6800×10<sup>+5</sup>2.6800×10<sup>+5</sup> 本方法一次迭代One iteration of this method 2.9724×10<sup>+5</sup>2.9724×10<sup>+5</sup> 本方法二次迭代Second iteration of this method 2.9658×10<sup>+5</sup>2.9658×10<sup>+5</sup> 本方法三次迭代This method iterates three times 2.6187×10<sup>+5</sup>2.6187×10<sup>+5</sup>

从重建图像效果来看,使用空间插值方法获得的重建图像条形伪影严重,严重影响了图像诊断信息的提取,而Patch方法和本方法则可以显著地改善图像质量。From the perspective of the reconstructed image effect, the reconstructed image obtained by using the spatial interpolation method has serious stripe artifacts, which seriously affects the extraction of image diagnostic information, while the Patch method and this method can significantly improve the image quality.

使用绝对误差对重建图像进行定量分析可以发现:和Patch方法相比,本发明方法获得的图像质量更接近理想图像质量,每一次的迭代效果都优于Patch方法,而且后面的迭代效果优于前一次的迭代效果。Using the absolute error to quantitatively analyze the reconstructed image, it can be found that: compared with the Patch method, the image quality obtained by the method of the present invention is closer to the ideal image quality, and the effect of each iteration is better than that of the Patch method, and the latter iteration effect is better than the previous one. One iteration effect.

最后应当说明的是,以上实施例仅用以说明本发明的技术方案而非对本发明保护范围的限制,尽管参照较佳实施例对本发明作了详细说明,本领域的普通技术人员应当理解,可以对本发明的技术方案进行修改或者等同替换,而不脱离本发明技术方案的实质和范围。Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit the protection scope of the present invention. Although the present invention has been described in detail with reference to the preferred embodiments, those of ordinary skill in the art should The technical solutions of the present invention may be modified or equivalently replaced without departing from the spirit and scope of the technical solutions of the present invention.

Claims (3)

1.一种圆轨道下不完全CT投影数据的恢复方法,其特征在于,通过如下步骤进行:1. the recovery method of incomplete CT projection data under a circular track, is characterized in that, is carried out by the following steps: 第一步,获得坏像素点位置信息;The first step is to obtain the location information of bad pixels; 第二步,对每个投影中的坏像素点使用水平一维三次样条插值方法进行线性插值,以恢复坏像素点的低频信息,经过此步骤得到每个角度下的初始化投影;In the second step, the horizontal one-dimensional cubic spline interpolation method is used to linearly interpolate the bad pixels in each projection to restore the low-frequency information of the bad pixels. After this step, the initial projection at each angle is obtained; 第三步,令N=1,进入第四步;The third step, let N=1, enter the fourth step; 第四步,令M=N,M为当前迭代次数;The fourth step, let M=N, M is the current number of iterations; 将每个角度下的初始化投影进行二维傅里叶变换,获得每个角度下的初始化投影在频率域中的频域数据;Perform a two-dimensional Fourier transform on the initialized projection at each angle to obtain the frequency domain data of the initialized projection at each angle in the frequency domain; 第五步,分别将前后两个相邻投影的频率域形式带入到一致性条件对应的公式中计算投影结果;使用相同的权重因子ω对前后两个投影结果进行加权,得到加权更新后的投影结果;In the fifth step, the frequency domain forms of the two adjacent projections before and after are respectively brought into the formula corresponding to the consistency condition to calculate the projection result; the same weighting factor ω is used to weight the two projection results before and after, and the weighted update is obtained. projection result; 第六步,对加权后的投影结果进行傅里叶反变换,得到第M次迭代终值;The sixth step is to perform inverse Fourier transform on the weighted projection result to obtain the final value of the Mth iteration; 第七步,判断当前迭代次数M是否大于等于迭代终止次数S,如果是,则进入第九步;否则进入第八步;The seventh step is to judge whether the current iteration number M is greater than or equal to the iteration termination number S, and if so, enter the ninth step; otherwise, enter the eighth step; 第八步,将第M次迭代终值作为初始化投影,令N=N+1,进入第三步;The eighth step, take the final value of the Mth iteration as the initialization projection, let N=N+1, and enter the third step; 第九步,以第M次迭代终值为最终恢复的数据结果;The ninth step, with the final value of the M-th iteration as the final recovered data result; 第五步中使用的一致性条件,具体形式如下:The consistency condition used in the fifth step is as follows:
Figure FDA0002673576870000011
Figure FDA0002673576870000011
其中,ρ是射线源到旋转中心的距离,d是旋转中心到探测器的距离,(α12)是平板探测器坐标,θ是方位角;Among them, ρ is the distance from the ray source to the rotation center, d is the distance from the rotation center to the detector, (α 1 , α 2 ) are the coordinates of the flat panel detector, and θ is the azimuth angle; u是目标函数f穿过ε和η两点的线积分:u is the line integral of the objective function f through two points ε and η: u(ε;η)=∫f(ε+t(η-ε))dt;t为射线f穿过ε和η两点的积分变量符号。u(ε; η)=∫f(ε+t(η-ε))dt; t is the integral variable sign of the ray f passing through the two points ε and η.
2.根据权利要求1所述的圆轨道下不完全CT投影数据的恢复方法,其特征在于:对第五步中的一致性条件进行傅里叶变换,得到对应的频率域的表达形式如下:2. the recovery method of incomplete CT projection data under the circular orbit according to claim 1, is characterized in that: carry out Fourier transform to the consistency condition in the 5th step, obtain the expression form of corresponding frequency domain as follows:
Figure FDA0002673576870000012
Figure FDA0002673576870000012
其中,uθ *表示u对θ求导后的傅里叶变换,(k1,k2)是(α12)的对应频率域形式;Among them, u θ * represents the Fourier transform of u after derivation of θ, (k 1 , k 2 ) is the corresponding frequency domain form of (α 1 , α 2 ); 将公式(25)右边式子用U代替,得到公式(26):Substitute the right-hand side of formula (25) with U to obtain formula (26):
Figure FDA0002673576870000021
Figure FDA0002673576870000021
其中,u*(θ+dθ)和u*(θ)分别表示θ+dθ和θ角度下初始化投影的频域数据,u*(θ)为已知量,u*(θ+dθ)是要求的未知量,dθ代表相邻两个角度的间隔,(k1,k2)是平板探测器的频率域坐标,i代表复数的虚部;Among them, u*(θ+dθ) and u*(θ) represent the frequency domain data of the initialized projection at θ+dθ and θ angles, respectively, u*(θ) is a known quantity, and u*(θ+dθ) is the requirement The unknown quantity of , dθ represents the interval between two adjacent angles, (k 1 , k 2 ) is the frequency domain coordinate of the flat panel detector, i represents the imaginary part of the complex number; 将公式(25)改写成如下形式:Rewrite formula (25) into the following form:
Figure FDA0002673576870000022
Figure FDA0002673576870000022
对公式(27)进行等式变形,得到如下形式:Equationally transforming formula (27), the following form is obtained: u*(θ+dθ)=u*(θ)+dθ·U,k2≠0 (28);u * (θ+dθ)=u * (θ)+dθ·U,k 2 ≠0 (28); 第五步,具体使用公式(28)对每一个初始化投影的频域数据进行恢复;The 5th step, specifically uses formula (28) to restore the frequency domain data of each initialization projection; 公式(28)的具体使用方法如下:The specific usage of formula (28) is as follows: 首先将u*(θ+dθ)代入到公式(28)中计算u*(θ):First, substitute u*(θ+dθ) into formula (28) to calculate u*(θ): u*(θ)=u*(θ+dθ)-dθ·U;u * (θ)=u * (θ+dθ)-dθ·U; 然后用u*(θ-dθ)代替公式(28)中的u*(θ),u*(θ)代替公式(28)中的u*(θ+dθ)计算得到:Then use u*(θ-dθ) to replace u*(θ) in formula (28), and u*(θ) to replace u*(θ+dθ) in formula (28) to calculate: u*(θ)=u*(θ-dθ)+dθ·U;u * (θ)=u * (θ-dθ)+dθ·U; 由于选取前后两个相邻投影频域数据的任一个都可以实现投影数据恢复,为了平衡恢复效果,使用权重因子ω,ω∈[0,1],对上述两个公式的计算结果进行加权平均,得到加权后的结果uR*(θ),即:Since the projection data recovery can be achieved by selecting any one of the two adjacent projection frequency domain data before and after, in order to balance the recovery effect, the weighting factor ω, ω∈[0,1] is used to perform a weighted average of the calculation results of the above two formulas , the weighted result u R* (θ) is obtained, namely: uR*(θ)=ω(u*(θ-dθ)+dθ·U)+(1-ω)(u*(θ+dθ)-dθ·U),k2≠0;u R* (θ)=ω(u * (θ-dθ)+dθ·U)+(1-ω)(u * (θ+dθ)-dθ·U),k 2 ≠0; 由于公式(28)只有在k2≠0时才有效,所以对频域数据u*(θ)中k2=0的部分,直接使用θ角度下初始化投影频域数据中的相应部分来补偿;Since formula (28) is only valid when k 2 ≠0, for the part where k 2 =0 in the frequency domain data u*(θ), directly use the corresponding part in the initialized projected frequency domain data under the angle θ to compensate; 第六步,具体是对uR*(θ)进行傅里叶反变换,得到第M次迭代终值uR(θ);The sixth step is to perform inverse Fourier transform on u R* (θ) to obtain the final value u R (θ) of the Mth iteration; uR(θ)=F-1(uR*(θ)) (1);u R (θ)=F -1 (u R* (θ)) (1); 其中公式(1)中uR(θ)表示θ角度下加权后的频域数据uR*(θ)的反傅里叶变换结果,F-1代表傅里叶反变换;Wherein u R (θ) in formula (1) represents the inverse Fourier transform result of the weighted frequency domain data u R* (θ) under the angle θ, and F -1 represents the inverse Fourier transform; 第七步,具体是判断当前迭代次数M是否大于等于迭代终止次数S,如果是,则进入第九步;否则进入第八步;The seventh step is to determine whether the current iteration number M is greater than or equal to the iteration termination number S, and if so, enter the ninth step; otherwise, enter the eighth step; 第八步,具体是将第M次迭代终值uR(θ)作为初始化投影,令N=N+1,进入第三步;The eighth step is to use the final value u R (θ) of the Mth iteration as the initialization projection, let N=N+1, and enter the third step; 第九步,具体是以第M次迭代终值uR(θ)为最终恢复的数据结果。The ninth step, specifically, takes the final value u R (θ) of the Mth iteration as the final recovered data result.
3.根据权利要求2所述的圆轨道下不完全CT投影数据的恢复方法,其特征在于,3. the recovery method of incomplete CT projection data under circular orbit according to claim 2, is characterized in that, 第一步,获得坏像素点位置信息具体是将实际采集的投影数据转换成对应的弦图,即:弦图的x轴为投影的像素点位置,y轴为该像素位置对应的投影角度位置;对某个探测器单元存在损坏的条件下,则在弦图中对应的像素位置上会出现一条竖线,根据此理论在弦图中获得坏像素点的位置信息,并记录坏像素点位置。The first step, obtaining the bad pixel position information is to convert the actually collected projection data into the corresponding chord diagram, that is: the x-axis of the chord diagram is the projected pixel position, and the y-axis is the projection angle position corresponding to the pixel position. ; Under the condition that a certain detector unit is damaged, a vertical line will appear on the corresponding pixel position in the chord diagram. According to this theory, the position information of the bad pixel point is obtained in the chord diagram, and the position of the bad pixel point is recorded. .
CN201710501987.9A 2017-06-27 2017-06-27 A recovery method of incomplete CT projection data under circular orbit Active CN107292846B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710501987.9A CN107292846B (en) 2017-06-27 2017-06-27 A recovery method of incomplete CT projection data under circular orbit

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710501987.9A CN107292846B (en) 2017-06-27 2017-06-27 A recovery method of incomplete CT projection data under circular orbit

Publications (2)

Publication Number Publication Date
CN107292846A CN107292846A (en) 2017-10-24
CN107292846B true CN107292846B (en) 2020-11-10

Family

ID=60098750

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710501987.9A Active CN107292846B (en) 2017-06-27 2017-06-27 A recovery method of incomplete CT projection data under circular orbit

Country Status (1)

Country Link
CN (1) CN107292846B (en)

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6173030B1 (en) * 1998-11-25 2001-01-09 General Electric Company Almost-everywhere extrapolation using 2D transforms from cone-beam data
EP1769460B1 (en) * 2004-07-13 2008-01-16 Philips Intellectual Property & Standards GmbH A computed tomography method for the reconstruction of object images from real and fictitious measured values
CN101627919A (en) * 2009-08-20 2010-01-20 浙江大学 PET concentration reestablishing method based on Kalman filtration in limited sampling angle
CN102680948A (en) * 2012-05-15 2012-09-19 东南大学 Method for estimating modulation frequency and starting frequency of linear frequency-modulated signal
CN104021582A (en) * 2014-05-28 2014-09-03 山东大学 CT (Computed Tomography) iterative image reconstruction method
CN104050631A (en) * 2013-11-25 2014-09-17 中国科学院上海应用物理研究所 Low-dose CT image reconstruction method
CN105590332A (en) * 2015-12-24 2016-05-18 电子科技大学 Rapid algebraic reconstruction technique applied to computed tomography imaging
CN106612427A (en) * 2016-12-29 2017-05-03 浙江工商大学 Method for generating spatial-temporal consistency depth map sequence based on convolution neural network

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6173030B1 (en) * 1998-11-25 2001-01-09 General Electric Company Almost-everywhere extrapolation using 2D transforms from cone-beam data
EP1769460B1 (en) * 2004-07-13 2008-01-16 Philips Intellectual Property & Standards GmbH A computed tomography method for the reconstruction of object images from real and fictitious measured values
CN101627919A (en) * 2009-08-20 2010-01-20 浙江大学 PET concentration reestablishing method based on Kalman filtration in limited sampling angle
CN102680948A (en) * 2012-05-15 2012-09-19 东南大学 Method for estimating modulation frequency and starting frequency of linear frequency-modulated signal
CN104050631A (en) * 2013-11-25 2014-09-17 中国科学院上海应用物理研究所 Low-dose CT image reconstruction method
CN104021582A (en) * 2014-05-28 2014-09-03 山东大学 CT (Computed Tomography) iterative image reconstruction method
CN105590332A (en) * 2015-12-24 2016-05-18 电子科技大学 Rapid algebraic reconstruction technique applied to computed tomography imaging
CN106612427A (en) * 2016-12-29 2017-05-03 浙江工商大学 Method for generating spatial-temporal consistency depth map sequence based on convolution neural network

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Consistency conditions upon 3D CT data and the wave equation;S K Patch 等;《Physics in Medicine and Biology》;20021231(第47期);第1-14页 *
Improved 2D rebinning of helical cone-beam CT data using John"s equation;M. Defrise 等;《2002 IEEE Nuclear Science Symposium Conference Record》;20031027;第1465-1469页 *
WE-AB-207A-02:John’s Equation Based Consistency Condition and Incomplete Projection Restoration Upon Circular Orbit CBCT;J. Ma 等;《Medical Physics》;20160630;第43卷(第6期);参见第3797-3798页WE-AB-207A-02部分 *
基于投影域数据恢复的低剂量CT稀疏角度重建;刘文磊 等;《CT理论与应用研究》;20130731;第22卷(第3期);第421-428页 *

Also Published As

Publication number Publication date
CN107292846A (en) 2017-10-24

Similar Documents

Publication Publication Date Title
US8731266B2 (en) Method and system for correcting artifacts in image reconstruction
CN109146994B (en) Metal artifact correction method for multi-energy spectrum X-ray CT imaging
Zbijewski et al. Characterization and suppression of edge and aliasing artefacts in iterative x-ray CT reconstruction
Park et al. Characterization of metal artifacts in X‐ray computed tomography
CN103479379B (en) A kind of image rebuilding method of tilting screw scanning and device
CN103413280A (en) Low-dose X-ray CT image reconstruction method
CN101980302A (en) Non-local Average Low Dose CT Reconstruction Method Guided by Projection Data Restoration
CN103106676B (en) A kind of X ray CT image rebuilding method based on the filtering of low dosage data for projection
Yan et al. Projection correlation based view interpolation for cone beam CT: primary fluence restoration in scatter measurement with a moving beam stop array
Zhang et al. Dynamic cone-beam CT reconstruction using spatial and temporal implicit neural representation learning (STINR)
US20070253529A1 (en) Systems and methods for digital volumetric laminar tomography
Yang et al. Shading correction assisted iterative cone-beam CT reconstruction
Li et al. Sparse CT reconstruction based on multi-direction anisotropic total variation (MDATV)
CN110211193B (en) 3D CT inter-slice image interpolation restoration and super-resolution processing method and device
Friot et al. Iterative tomographic reconstruction with TV prior for low-dose CBCT dental imaging
CN103793890A (en) Method for recovering and processing energy spectrum CT images
CN105844678A (en) Low dose X-ray CT image reconstruction method based on completely generalized variational regularization
CN107292846B (en) A recovery method of incomplete CT projection data under circular orbit
CN105354868B (en) Limited angle conical beam CT image rebuilding method based on several picture square
CN112656438A (en) Low-dose CT projection domain denoising and reconstructing method based on curved surface total variation
CN112085811A (en) Method and device for CT local reconstruction
CN110264536A (en) A method of high-low resolution projection relation is calculated in the reconstruction of parallel beam oversubscription
Pengpan et al. Cone Beam CT using motion-compensated algebraic reconstruction methods with limited data
Je et al. Dental cone-beam CT reconstruction from limited-angle view data based on compressed-sensing (CS) theory for fast, low-dose X-ray imaging
Je et al. Simulation and experimental studies of three-dimensional (3D) image reconstruction from insufficient sampling data based on compressed-sensing theory for potential applications to dental cone-beam CT

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant