CN104778667A

CN104778667A - Level-set-based correction method for cupping artifact in cone-beam CT

Info

Publication number: CN104778667A
Application number: CN201510176350.8A
Authority: CN
Inventors: 谢世朋; 丁铭晨; 李海波; 葛琦; 闫瑞菊
Original assignee: Nanjing Post and Telecommunication University
Current assignee: Nanjing Post and Telecommunication University
Priority date: 2015-04-14
Filing date: 2015-04-14
Publication date: 2015-07-15
Anticipated expiration: 2035-04-14
Also published as: CN104778667B

Abstract

The invention discloses a method for correcting cup artifacts in cone-beam CT based on a level set, and the method is applied to correction of cone-beam CT slice images. This method can self-adaptively correct the cupping artifact of cone beam CT, and the correction can be completed automatically without manual intervention. This method does not need to scan the measured object repeatedly; it does not increase the complexity of the cone beam CT system; for the reconstructed slice images, it can be directly oriented to the user without any modification to the original equipment of the original cone beam CT. After completing the correction work, this method can efficiently correct the cupping artifact of cone beam CT, and can also improve the contrast of the image.

Description

A Correction Method of Cupping Artifact in Cone Beam CT Based on Level Set

技术领域 technical field

本发明涉及一种使用水平集算法的锥束CT切片图像中杯状伪影的校正方法，属于图像处理技术领域。 The invention relates to a method for correcting cupping artifacts in cone-beam CT slice images using a level set algorithm, and belongs to the technical field of image processing.

背景技术 Background technique

锥束CT具有很好的扫描速度和辐射利用率，能有效的减少X射线管的负载输出，降低扫描成本，也能快速获得高分辨率三维断层图像数据。 Cone beam CT has good scanning speed and radiation utilization rate, can effectively reduce the load output of X-ray tubes, reduce scanning costs, and can quickly obtain high-resolution three-dimensional tomographic image data.

目前，影响锥束CT重建图像质量的因素有很多，如：X线散射、噪声、几何误差、能谱、探测单元响应不一致等。但由于锥束平板CT使用大范围的X射线平板探测器，这使得成像质量与传统CT相比较更易受到X射线散射及射束硬化的影响。因散射和射束硬化而形成的伪影、CT数的不准确等严重影响对重建图像的分析与判断。在针对人体的锥束CT重建图像中，这些伪影主要表现为CT值不均匀性伪影，该伪影的大多表现为杯状的伪影。这些伪影对于基于阈值的可视化显示方面和基于阈值的锥束CT图像分割方面影响非常严重。因此，针对锥束CT杯状伪影的校正显得尤为必要。 At present, there are many factors affecting the image quality of cone-beam CT reconstruction, such as: X-ray scattering, noise, geometric error, energy spectrum, inconsistent response of detection units, etc. However, because cone-beam flat-panel CT uses a wide range of X-ray flat-panel detectors, the imaging quality is more susceptible to X-ray scattering and beam hardening compared with traditional CT. Artifacts caused by scattering and beam hardening, inaccurate CT numbers, etc. seriously affect the analysis and judgment of reconstructed images. In cone-beam CT reconstructed images for the human body, these artifacts are mainly manifested as CT value inhomogeneity artifacts, most of which are cup-shaped artifacts. These artifacts have a serious impact on threshold-based visualization and threshold-based cone-beam CT image segmentation. Therefore, the correction of cupping artifacts in cone beam CT is particularly necessary.

为了减少杯状伪影(即：CT值不均匀性伪影)的影响，目前，大多数的校正方法主要是考虑因散射造成的杯状伪影，并且现有技术或文献资料的研究主要集中在对投影图像上的散射校正。可以将这些方法分成两类：一类是基于软件，另一类是基于硬件。大多数以软件为基础的方法是基于蒙特卡罗仿真，它是一种对锥束CT散射校正非常有效的方法。但它是特别耗时。近年来，人们已经提出一些改进的蒙特卡洛仿真算法，如：基于GPU的方法等。但是，即使采用快速算法，繁重的计算还是阻碍了其实际的应用。许多基于硬件的校正方法是也不错的，例如：部分射线遮挡、初级调制等。 In order to reduce the influence of cupping artifacts (i.e. CT value inhomogeneity artifacts), at present, most correction methods mainly consider the cupping artifacts caused by scattering, and the research of existing technologies or literature mainly focuses on Scattering correction on projected images. These methods can be divided into two categories: one is software-based, and the other is hardware-based. Most software-based methods are based on Monte Carlo simulation, which is a very effective method for scatter correction in cone-beam CT. But it is extremely time consuming. In recent years, people have proposed some improved Monte Carlo simulation algorithms, such as methods based on GPU. However, even with a fast algorithm, heavy computation hinders its practical application. Many hardware-based correction methods are also good, such as: partial ray occlusion, primary modulation, etc.

现有技术的缺点主要包括： The disadvantages of the prior art mainly include:

(1)现有技术主要是针对投影图像校正，没有直接针对重建后切片图像的校正。 (1) The existing technology is mainly aimed at the projection image correction, and there is no direct correction for the reconstructed slice image.

(2)现有技术大多数集中在因散射而造成的伪影校正的方法中，并且大多需要添加硬件设备，如：专利200710019084和201310039298，这两个专利都需要在昂贵的锥束CT设备上添加硬件设备，增加了操作的复杂性和对设备造成潜在的安全风险。特别是专利200710019084需要两次扫描被测物体，这样无疑增加的被测物的辐射量。 (2) Most of the existing technologies focus on the method of artifact correction caused by scattering, and most of them need to add hardware equipment, such as: patents 200710019084 and 201310039298, both of which require expensive cone-beam CT equipment Adding hardware devices increases the complexity of operations and poses potential security risks to devices. In particular, patent 200710019084 needs to scan the measured object twice, which undoubtedly increases the radiation dose of the measured object.

综上所述，在现有技术或文献资料的方法中，蒙特卡洛模拟方法非常耗费时间，初级射线调制方法中校正结果受限于调制板自身的结构，基于部分散射射线测量方法，有些需要增加照射剂量，有些方法对散射分布的估计准确度不高。而本发明能够很好地解决上面的问题。 To sum up, among the methods in the prior art or literature, the Monte Carlo simulation method is very time-consuming, and the correction result in the primary ray modulation method is limited by the structure of the modulation plate itself. Based on the partially scattered ray measurement method, some require Increasing the irradiation dose, some methods can not estimate the scattering distribution with high accuracy. And the present invention can well solve the above problems.

发明内容 Contents of the invention

本发明目的在于提出了一种基于水平集的锥束CT中杯状伪影的校正方法，该方法应用于锥束CT切片图像校正。该方法能够自适应的进行锥束CT的杯状伪影校正，无需人工干预就可以自动完成校正。该方法不需要重复扫描被测物体；不增加锥束CT系统的复杂度；针对重建后的切片图像，能够直接面向用户，不需要对原有锥束CT的原有设备进行任何改动，就可以完成校正工作，该方法能够高效地进行锥束CT的杯状伪影校正，同时还能够提高图像的对比度。 The object of the present invention is to propose a method for correcting cup artifacts in cone-beam CT based on level sets, and the method is applied to correction of cone-beam CT slice images. This method can self-adaptively correct the cupping artifact of cone beam CT, and the correction can be completed automatically without manual intervention. This method does not need to scan the measured object repeatedly; it does not increase the complexity of the cone beam CT system; for the reconstructed slice images, it can be directly oriented to the user without any modification to the original equipment of the original cone beam CT. After completing the correction work, this method can efficiently correct the cupping artifact of cone beam CT, and can also improve the contrast of the image.

本发明解决其技术问题所采取的技术方案是：一种基于水平集的锥束CT中杯状伪影的校正方法，该方法包括如下步骤： The technical scheme that the present invention solves its technical problem is: a kind of correction method of cupping artifact in the cone-beam CT based on level set, this method comprises the steps:

步骤1:获取带有伪影的锥束CT切片数据； Step 1: Obtain cone beam CT slice data with artifacts;

步骤2:根据公式(7)和公式(8)计算全局量函数ε(φ,f_s,p)和水平集约束方程F(φ,p,f_s)； Step 2: Calculate the global quantity function ε(φ,f _s ,p) and the level set constraint equation F(φ,p,f _s ) according to formula (7) and formula (8);

步骤3:对带有伪影的锥束CT切片数据，根据公式(10)，固定p和f_s，使用有限差分方法迭代演化出 Step 3: For the cone-beam CT slice data with artifacts, according to formula (10), fix p and f _s , and use the finite difference method to iteratively evolve

步骤4:根据公式(11)，固定φ和f_s，计算p变量值的估计值 Step 4: According to formula (11), fixing φ and f _s , calculate the estimated value of the p variable value

步骤5:根据公式(12)，固定φ和p，计算f_s变量值的估计值 Step 5: According to equation (12), fixing φ and p, calculate the estimated value of f _s variable value

步骤6:若不收敛或达未到迭代次数，则令且回到步骤3； Step 6: If does not converge or does not reach the number of iterations, then let And go back to step 3;

步骤7:根据公式(13)，计算校正后的锥束CT切片图像。 Step 7: Calculate the corrected cone-beam CT slice image according to formula (13).

有效效果： Effective effect:

1、本发明能够直接针对重建后的切片图像的杯状伪影校正。 1. The present invention can directly correct the cupping artifact of the reconstructed slice image.

2、本发明计算量相对较小，能够高效地进行锥束CT切片图像杯状伪影校正的同时，也能够很好地提高图像的对比度。 2. The calculation amount of the present invention is relatively small, and it can efficiently correct the cupping artifact of the cone-beam CT slice image and at the same time can well improve the contrast of the image.

3、本发明能够直接面向CT切片需求用户，不需要对原有锥束CT原有设备进行任何改动就能完成校正工作。 3. The present invention can directly face users who need CT slices, and can complete the correction work without any modification to the original cone-beam CT equipment.

附图说明 Description of drawings

图1为本发明的方法流程图。 Fig. 1 is a flow chart of the method of the present invention.

图2为图像切片的示意图。 Figure 2 is a schematic diagram of image slices.

标识说明：(a)表示CTP486原始切片图；(b)表示使用本发明方法校正后的CTP486切片图；(c)表示CTP Top原始切片图；(d)表示使用本发明方法校正后的CTP Top切片图。 Identification description: (a) indicates the original slice image of CTP486; (b) indicates the corrected CTP486 slice image using the method of the present invention; (c) indicates the original slice image of CTP Top; (d) indicates the corrected CTP Top image using the method of the present invention slice graph.

图3为图像切片的水平剖面图。 Figure 3 is a horizontal cross-sectional view of an image slice.

标识说明：(a)、(b)表示中间行的水平剖面图。 Marking instructions: (a) and (b) represent the horizontal section of the middle row.

图4为CTP486杯状伪影指标τ_cup计算区域的示意图。 FIG. 4 is a schematic diagram of the calculation area of the CTP486 cupping artifact index τ _cup .

图5为CTP Top杯状伪影指标τ_cup计算区域。 Figure 5 shows the calculation area of the CTP Top cupping artifact index τ _cup .

图6为灰度不均匀校正前图像(左侧)和灰度不均匀校正后图像(右侧)的示意图。 FIG. 6 is a schematic diagram of an image before (left) correction of gray scale unevenness and an image after correction of gray scale unevenness (right side).

图7为图6头部锥束CT切片的1维水平剖面图。 Fig. 7 is a one-dimensional horizontal cross-sectional view of the cone-beam CT slice of the head in Fig. 6 .

图8为乳房锥束CT切片图像灰度不均匀校正实验，左列为校正前图像(a,c),右列为校正后图像(b,d)的示意图。 Figure 8 is a schematic diagram of the correction experiment of gray scale unevenness of breast cone-beam CT slice images. The left column is the image before correction (a, c), and the right column is the schematic diagram of the corrected image (b, d).

具体实施方式 Detailed ways

下面结合说明书附图对本发明创造作进一步的详细说明。 The invention will be described in further detail below in conjunction with the accompanying drawings.

如图1所示，一种基于水平集的锥束CT中杯状伪影的校正方法，该方法包括如下步骤： As shown in Figure 1, a method for correcting cupping artifacts in cone beam CT based on level sets, the method includes the following steps:

本发明的校正方法具体包括： The correction method of the present invention specifically comprises:

重建切片图像分解包括： The reconstructed slice image decomposition consists of:

锥束CT中重建算法是以FDK算法为基础的，重建图像集f可以写成 The reconstruction algorithm in cone beam CT is based on the FDK algorithm, and the reconstructed image set f can be written as

$\begin{matrix} f f = = \frac{11}{{44 π π}^{22}} {&Integral; &Integral;}_{00}^{22 π π} \frac{{d d}_{so so}^{22}}{{(({d d}_{so so} + + r r \cdot &Center Dot; s the s))}^{22}} {&Integral; &Integral;}_{- - \infty \infty}^{\infty \infty} \frac{{d d}_{so so}}{{(({d d}_{so so}^{22} + + {t t}^{22} + + {z z}^{22}))}^{11 / / 22}} \\ \cdot &Center Dot; {I I}_{33 D D.} ((t t,, z z ((r r)),, φ φ)) \cdot &Center Dot; {&Integral; &Integral;}_{- - \infty \infty}^{\infty \infty} | | ω ω | | {e e}^{j j 22 πω πω ((t t ((r r)) - - t t))} dωdtdφ dωdtdφ \end{matrix} - - - - - - ((11))$

其中d_so表示源到旋转轴的距离，I_3D(t,z(r),φ)表示投影图像的序列。这里，投影图像I_3D分解为如下： where d _so denotes the distance from the source to the axis of rotation, and I _3D (t,z(r),φ) denotes the sequence of projected images. Here, the projected image I _3D is decomposed as follows:

I_3D＝P_3D+S_3D+n (2) I _3D ＝P _3D +S _3D +n (2)

其中P_3D是从没有伪影的真实图像，S_3D是由散射和射束硬化造成的伪影部分，n是噪声。P_3D是物体切片对象的固有物理属性，如果有N中物质，那么可以假设它被分为N个灰度(即：CT值)不变的区域。根据观察，本发明能发现锥束CT切片伪影S_3D是缓慢变化的杯状伪影。重写重建图像公式如下， where P _3D is the real image without artifacts, S _3D is the artifact part caused by scattering and beam hardening, and n is the noise. P _3D is the inherent physical property of the object slice object. If there are N kinds of substances, it can be assumed that it is divided into N regions with constant gray levels (ie: CT values). According to the observation, the present invention can find that the cone-beam CT slice artifact S _3D is a slowly changing cupping artifact. Rewrite the reconstructed image formula as follows,

$\begin{matrix} f f = = \frac{11}{{44 π π}^{22}} {&Integral; &Integral;}_{22}^{22 π π} \frac{{d d}_{so so}^{22}}{{(({d d}_{so so} + + r r \cdot &Center Dot; s the s))}^{22}} {&Integral; &Integral;}_{- - \infty \infty}^{\infty \infty} \frac{{d d}_{so so}}{{(({d d}_{so so}^{22} + + {t t}^{22} + + {z z}^{22}))}^{11 / / 22}} \cdot \cdot \\ (({P P}_{33 D D.} ((t t,, z z ((r r)),, φ φ)) + + {S S}_{33 D D.} ((t t,, z z ((r r)),, φ φ)) + + n no)) \cdot \cdot {&Integral; &Integral;}_{- - \infty \infty}^{\infty \infty} | | ω ω | | {e e}^{j j 22 πω πω ((t t ((t t)) - - t t))} dωdtdφ dωdtdφ \\ = = {f f}_{p p} + + {f f}_{s the s} + + {f f}_{n no} \end{matrix} - - - - - - ((33))$

其中 in

$\begin{matrix} {f f}_{s the s} = = \frac{11}{{44 π π}^{22}} {&Integral; &Integral;}_{00}^{22 π π} \frac{{d d}_{so so}^{22}}{{(({d d}_{so so} + + r r \cdot &Center Dot; s the s))}^{22}} {&Integral; &Integral;}_{- - \infty \infty}^{\infty \infty} \frac{{d d}_{so so}}{{(({d d}_{so so}^{22} + + {t t}^{22} + + {z z}^{22}))}^{11 / / 22}} \\ \cdot &Center Dot; {S S}_{33 D D.} ((t t,, z z ((r r)),, φ φ)) \cdot &Center Dot; {&Integral; &Integral;}_{- - \infty \infty}^{\infty \infty} | | ω ω | | {e e}^{j j 22 πω πω ((t t ((r r)) - - t t))} dωdtdφ dωdtdφ \end{matrix},,$

也就是说重建图像可以表示为三个独立成分相加。 That is to say, the reconstructed image can be expressed as the addition of three independent components.

定义邻域内能量函数包括： Defining the energy function in the neighborhood includes:

本发明考虑一个半径ρ的圆形的区域Θ_y∈Ω，并定义Θ_y＝{x:|x-y|≤ρ}，y∈Ω。整个域切片图像记为Ω，Ω可分为N个子区域，分别为由于f_s同样是一个缓慢变化的图像，所以在圆形域Θ_y中，所有X的值f_s(x)接近于值f_s(y)。因此，在每一子区域Θ_y∩Ω_i，强度f_p(x)+f_s(x)接近于系数p_i+f_s(y)。 The present invention considers a circular region Θ _y ∈ Ω with a radius ρ, and defines Θ _y ={x:|xy|≤ρ}, y ∈ Ω. The entire domain slice image is recorded as Ω, and Ω can be divided into N sub-regions, respectively Since f _s is also a slowly changing image, all X values f _s (x) are close to the value f _s (y) in the circular field Θ _y . Therefore, in each subregion Θ _y ∩Ω _i , the intensity f _p (x)+f _s (x) is close to the coefficient p _i +f _s (y).

因此，对于区域Θ_y∩Ω_i本发明有： Therefore, for the region Θ _y ∩ Ω _i the present invention has:

f(x)≈f_s(y)+p_i， (4) f(x)≈f _s (y)+p _i , (4)

对邻域Θ_y内，本发明定义能量函数为 In the neighborhood Θ _y , the present invention defines the energy function as

${ϵ ϵ}_{{Θ Θ}_{y the y}} = = \underset{{Ω Ω}_{i i}}{&Integral; &Integral;} K K ((y the y - - x x)) {| | ff ff ((x x)) - - (({f f}_{s the s} ((y the y)) + + {p p}_{i i})) | |}^{22} dx dx,, - - - - - - ((55))$

其中α是一个归一化常数，使得∫K(u)＝1，δ是高斯函数的标准方差(或尺度参数)，ρ是Θ_y圆形域的半径。 in α is a normalization constant such that ∫K(u)=1, δ is the standard deviation (or scale parameter) of the Gaussian function, and ρ is the radius of the circular domain of _Θy .

定义全局能量函数及水平集函数包括： Defining the global energy function and level set function includes:

对于整个图像区域，能量函数可以表示为： For the entire image area, the energy function can be expressed as:

$\begin{matrix} ϵ ϵ = = &Integral; &Integral; {ϵ ϵ}_{{Θ Θ}_{y the y}} dy dy \\ = = &Integral; &Integral; (({Σ Σ}_{i i = = 11}^{N N} \underset{{Ω Ω}_{i i}}{&Integral; &Integral;} K K ((y the y - - x x)) {| | f f ((x x)) - - (({f f}_{s the s} ((y the y)) + + {p p}_{i i})) | |}^{22} dx dx)) dy dy \end{matrix} - - - - - - ((66))$

本发明使用水平集来表示能量函数，以N＝2的水平集为例，其他多水平集情况可类似推导。为方便起见，本发明用向量p＝(p₁,…p_N)代表常数p_1,…p_N。因此，水平集函数φ，向量p，和杯状f_s是能量ε的变量，故ε可以写成ε(φ,f_s,p)的如下形式： The present invention uses a level set to represent the energy function, taking the level set of N=2 as an example, other multi-level set situations can be derived similarly. For the sake of convenience, the present invention uses the vector p=(p ₁ ,...p _N ) to represent the constants p _1, ...p _N . Therefore, the level set function φ, vector p, and cup f _s are variables of energy ε, so ε can be written as ε(φ,f _s ,p) in the following form:

$\begin{matrix} ϵ ϵ ((φ φ,, {f f}_{s the s},, p p)) = = &Integral; &Integral; (({Σ Σ}_{i i = = 11}^{22} &Integral; &Integral; K K ((y the y - - x x)) {| | f f ((x x)) - - (({f f}_{s the s} ((y the y)) + + {p p}_{i i})) | |}^{22} {M m}_{i i} ((φ φ ((x x)))) dx dx)) dy dy \\ = = &Integral; &Integral; (({Σ Σ}_{i i = = 11}^{22} &Integral; &Integral; K K ((y the y - - x x)) {| | f f ((x x)) - - (({f f}_{s the s} ((y the y)) + + {p p}_{i i})) | |}^{22} dy dy)) {M m}_{i i} ((φ φ ((x x)))) dx dx \end{matrix},, - - - - - - ((77))$

其中，M₁(φ)＝H(φ)和M₂(φ)＝1-H(φ)，H是Heaviside函数，而φ是一个水平集函数，其零水平集为C₀＝{x:φ(x)＝0}，它分图像域为两个不相交的Ω₁＝{x:φ(x)＞0}和Ω₂＝{x:φ(x)＜0}区域。对于2个水平集，C₀＝{x:φ(x)＝0}和C_k＝{x:φ_k(x)＝0}则把图像分为三个区域，即N＝3。 Among them, M ₁ (φ)=H(φ) and M ₂ (φ)=1-H(φ), H is a Heaviside function, and φ is a level set function, and its zero level set is C ₀ ={x: φ(x)=0}, which divides the image domain into two disjoint Ω ₁ ={x:φ(x)>0} and Ω ₂ ={x:φ(x)<0} regions. For 2 level sets, C ₀ ={x:φ(x)=0} and C _k ={x: _φk (x)=0} divide the image into three regions, ie N=3.

为了约束零水平集函数，本发明对整个水平集函数添加约束项。新水平集方程定义为： In order to constrain the zero level set function, the present invention adds constraint items to the whole level set function. The new level set equation is defined as:

F(φ,p,f_s)＝ε(φ,f_s,p)+vL(φ)+μR_q(φ) (8) F(φ,p,f _s )＝ε(φ,f _s ,p)+vL(φ)+μR _q (φ) (8)

其中用来约束零水平轮廓线的光滑性，H是Heaviside函数, in Used to constrain the smoothness of the zero-level contour line, H is the Heaviside function,

$R_{p} (φ) = &Integral; q (&dtri; φ) dx$ 用来维持符号距离属性，其中 $q = \frac{{(s - 1)}^{2}}{2} .$ $R_{p} (φ) = &Integral; q (&dtri; φ) dx$ Used to maintain the signed distance property, where $q = \frac{{(the s - 1)}^{2}}{2} .$

7.2.4估计变量φ,p,f_s 7.2.4 Estimated variables φ,p,f _s

固定p和f_s，最小化F(φ,p,f_s)使用标准梯度下降方法，解梯度流方程得到： Fix p and f _s , minimize F(φ,p,f _s ) using standard gradient descent methods, and solve the gradient flow equation get:

$\begin{matrix} \frac{&PartialD; &PartialD; φ φ}{&PartialD; &PartialD; t t} = = - - {Σ Σ}_{i i = = 11}^{N N} \frac{{&PartialD; &PartialD; M m}_{i i}}{&PartialD; &PartialD; φ φ} {e e}_{i i} + + vδ vδ ((φ φ)) div div ((\frac{&dtri; &dtri; φ φ}{| | &dtri; &dtri; φ φ | |})) \\ + + μ μ div div ((dq dq ((&dtri; &dtri; φ φ)) &dtri; &dtri; φ φ)) \end{matrix},, - - - - - - ((99))$

这里e_i＝∫K(y-x)|f(x)-(f_s(y)+p_i)|²dy，是梯算子，div(·)是散度算子，δ为狄拉克δ函数。 Here e _i =∫K(yx)|f(x)-(f _s (y)+p _i )| ² dy, is the ladder operator, div(·) is the divergence operator, and δ is the Dirac δ function.

我们使用有限差分方法迭代演化出φ，结果用表示。方法描述如下：设为φ离散化的值，为公式(9)右式近似值。则有进而得到迭代公式： We iteratively evolve φ using the finite difference method, resulting in express. The method is described as follows: is the discretized value of φ, is the approximate value of the right-hand formula of formula (9). then there is Then get the iterative formula:

${φ φ}_{i i,, j j}^{k k + + 11} = = {φ φ}_{i i,, j j}^{k k} + + ΔtL ΔtL (({φ φ}_{i i,, j j}^{k k})),, - - - - - - ((1010))$

其中令c₀为大于0的常数。 Which order c ₀ is a constant greater than 0.

固定φ和f_s，用表示p的优化后值。最小化ε(φ,f_s,p)，则 $\frac{&PartialD; ϵ (φ, f_{s}, p)}{&PartialD; p} = 0,$ 得到 Fixing φ and f _s , use Indicates the optimized value of p. Minimize ε(φ,f _s ,p), then $\frac{&PartialD; ϵ (φ, f_{the s}, p)}{&PartialD; p} = 0,$ get

${\overset{^^}{p p}}_{i i} = = \frac{&Integral; &Integral; ((&Integral; &Integral; ((K K ((y the y - - x x)) ((f f ((x x)) - - {f f}_{s the s} ((y the y)))))) dy dy)) {M m}_{i i} ((φ φ ((x x)))) dx dx}{&Integral; &Integral; ((&Integral; &Integral; K K ((y the y - - x x)) dy dy)) {M m}_{i i} ((φ φ ((x x)))) dx dx},, - - - - - - ((1111))$

固定φ和p。最小化ε(φ,f_s,p)，则得到 Fix φ and p. Minimize ε(φ,f _s ,p), then get

${\overset{^^}{f f}}_{s the s} = = {Σ Σ}_{i i = = 11}^{N N} \frac{&Integral; &Integral; ((&Integral; &Integral; ((K K ((y the y - - x x)) ((f f ((x x)) - - {p p}_{i i})))) dy dy)) {M m}_{i i} ((φ φ ((x x)))) dx dx}{&Integral; &Integral; ((&Integral; &Integral; K K ((y the y - - x x)) dy dy)) {M m}_{i i} ((φ φ ((x x)))) dx dx} - - - - - - ((1212))$

通过迭代法，每一次均最小化F(φ,p,f_s)。本次迭代中使用上一次的估计结果进行迭代，进行直至数据收敛或达到迭代次数为止。最后，本发明得到校正后图像为： By iterative method, F(φ,p,f _s ) is minimized each time. In this iteration, the last estimated result is used to iterate until the data converges or the number of iterations is reached. Finally, the corrected image obtained by the present invention is:

${f f}_{p p} = = f f - - {\overset{^^}{f f}}_{s the s} - - - - - - ((1313))$

本发明的实验过程和结果具体包括： Experimental process of the present invention and result specifically comprise:

定量分析指标定义包括： Quantitative analysis index definitions include:

本发明定义杯状伪影τ_cup＝100(u_M,edge-u_M,center)/u_M,edge，其中u_M,center和u_M,edge是模体中心和边缘的CT值(即：灰度值)。 The present invention defines cupping artifact τ _cup =100(u _{M, edge} -u _{M, center} )/u _{M, edge} , wherein u _{M, center} and u _{M, edge} are the CT values of the phantom center and edge (ie: grayscale value).

均方根对比度表示为其中I_ij是二维图像的(i,j)位置像素值。 The root mean square contrast ratio is expressed as where I _ij is the pixel value at position (i, j) of the two-dimensional image.

对比度信噪比度CNR计算公式为CNR＝|u_M,1-u_M,2|/σ_M，其中u_M,1，u_M,2为两个对比区域的均值，σ_M,1,σ_M,2为对比区域的标准差，像素噪声σ_M定义为σ_M＝(σ_M,1+σ_M,2)/2。 The calculation formula of contrast signal-to-noise ratio CNR is CNR=|u _M,1 -u _M,2 |/σ _M , where u _M,1 , u _M,2 are the mean values of two contrast areas, σ _M,1 ,σ _M,2 is the standard deviation of the comparison area, and the pixel noise σ _M is defined as σ _M =(σ _M,1 +σ _M,2 )/2.

Catphan 500模体锥束CT切片杯状伪影校正实验。 Catphan 500 phantom cone beam CT slice cupping artifact correction experiment.

本发明的实验使用Catphan500模体进行实验，使用模体中CTP486和模体顶部(即：CTP Top)。从图2中能看出，本发明能够消除杯状伪影到人眼难以觉察水平；在图3中表示的是出对应图2中的1维水平剖面，从中能看出本发明能明显的消除切片图像的杯状伪影。对散射校正效果定量分析见表1和表2。从表中能得到，本发明提出的方法使得杯状伪影τ_cup平均下降约91.8％。 In the experiment of the present invention, the Catphan500 phantom was used for the experiment, and the CTP486 in the phantom and the top of the phantom (ie: CTP Top) were used. As can be seen from Fig. 2, the present invention can eliminate cupping artifacts to a level that is difficult for the human eye to detect; in Fig. 3, it is shown that the 1-dimensional horizontal section corresponding to Fig. 2, from which it can be seen that the present invention can obviously Remove cupping artifacts from sliced images. See Table 1 and Table 2 for the quantitative analysis of the scattering correction effect. It can be obtained from the table that the method proposed by the present invention reduces the cupping artifact τ _cup by about 91.8% on average.

表1：CTP486切片图像校正前后定量分析，校正前的切图表示为PI_NONE,校正后的切片图表示为PI_SC Table 1: Quantitative analysis of CTP486 slice image before and after correction, the slice image before correction is represented as PI_NONE, and the slice image after correction is represented as PI_SC

表2：CTP Top切片图像校正前后定量分析，校正前的切图表示为PI_NONE,校正后的切片图表示为PI_SC。 Table 2: Quantitative analysis of CTP Top slice images before and after correction. The slice before correction is denoted as PI_NONE, and the slice after correction is denoted as PI_SC.

人头骨的锥束CT切片杯状伪影校正实验 Correction of Cupping Artifacts in Cone Beam CT Slices of Human Skull

对于人头骨的锥束CT杯状伪影校正前后切片图如图6所示。图7为图6图像的中间行剖面图。未校正前的大脑组织出现灰度不均(图7左侧)，表现为边缘区域亮度的升高。校正后灰度均匀，达到了理想效果。 Figure 6 shows the cone-beam CT slices of the human skull before and after cupping artifact correction. FIG. 7 is a cross-sectional view of the middle row of the image of FIG. 6 . Uncorrected brain tissue showed gray scale unevenness (left side of Figure 7), manifested as an increase in the brightness of the edge area. After correction, the gray level is even, and the ideal effect is achieved.

乳房锥束CT杯状伪影校正实验 Correction of Cupping Artifacts in Breast Cone Beam CT

如图8所示，乳房锥束CT实验显示本发明的方法能把灰度不均匀性伪影降低到人眼觉察不到的状态，且本发明提升了图像的对比度。校正后的RMSC是校正前的约1.3倍。 As shown in FIG. 8 , breast cone-beam CT experiments show that the method of the present invention can reduce gray level non-uniformity artifacts to a state that is imperceptible to human eyes, and the present invention improves the contrast of the image. The RMSC after correction is about 1.3 times that before correction.

Claims

1. a method for correcting goblet artifacts in cone beam CT based on level set, is characterized in that, described method comprises the steps:

Step 1: obtain cone beam CT slice data with artifacts;

Step 2: Calculate the global quantity function ε(φ,f _s ,p) and the level set constraint equation F(φ,p,f _s ) according to formula (7) and formula (8);

Step 3: For the cone-beam CT slice data with artifacts, according to formula (10), fix p and f _s , and use the finite difference method to iteratively evolve

Step 4: According to formula (11), fixing φ and f _s , calculate the estimated value of the p variable value

Step 5: According to equation (12), fixing φ and p, calculate the estimated value of f _s variable value

Step 6: If does not converge or does not reach the number of iterations, then let And return to step 3 above;

Step 7: Calculate the corrected cone-beam CT slice image according to formula (13).

2. a kind of correction method of cupping artifact in the cone beam CT based on level set according to claim 1, is characterized in that, the formula (7) of described step 2 is:

\begin{matrix} ϵ ϵ ((φ φ,, {f f}_{s the s},, p p)) = = &Integral; &Integral; (({Σ Σ}_{i i = = 11}^{22} &Integral; &Integral; K K ((y the y - - x x)) | | f f ((x x)) - - (({f f}_{s the s} ((y the y)) + + {p p}_{i i})) {| |}^{22} {M m}_{i i} ((φ φ ((x x)))) dx dx)) dy dy \\ = = &Integral; &Integral; (({Σ Σ}_{i i = = 11}^{22} &Integral; &Integral; K K ((y the y - - x x)) | | f f ((x x)) - - (({f f}_{s the s} ((y the y)) + + {p p}_{i i})) {| |}^{22} dy dy)) {M m}_{i i} ((φ φ ((x x)))) dx dx \end{matrix};;

The formula (8) of described step 2 is:

F(φ,p,f _s )=ε(φ,f _s ,p)+vL(φ)+μR _q (φ).

3. the correction method of cupping artifact in a kind of cone beam CT based on level set according to claim 1, is characterized in that, the formula (10) of described step 3 is:

{φ φ}_{i i,, j j}^{k k + + 11} = = {φ φ}_{i i,, j j}^{k k} + + ΔtL ΔtL (({φ φ}_{i i,, j j}^{k k})) . .

4. the correction method of cupping artifact in a kind of level set based cone beam CT according to claim 1, is characterized in that, the formula (11) of described step 4 is:

{\overset{^^}{p p}}_{i i} = = \frac{&Integral; &Integral; ((&Integral; &Integral; ((K K ((y the y - - x x)) ((f f ((x x)) - - {f f}_{s the s} ((y the y)))))) dy dy)) {M m}_{i i} ((φ φ ((x x)))) dx dx}{&Integral; &Integral; ((&Integral; &Integral; K K ((y the y - - x x)) dy dy)) {M m}_{i i} ((φ φ ((x x)))) dx dx} . .

5. a kind of correction method of cupping artifact in the cone beam CT based on level set according to claim 1, is characterized in that, the formula (12) of described step 5 is:

{\overset{^^}{f f}}_{s the s} = = {Σ Σ}_{i i = = 11}^{N N} \frac{&Integral; &Integral; ((&Integral; &Integral; ((K K ((y the y - - x x)) ((f f ((x x)) - - {p p}_{i i})))) dy dy)) {M m}_{i i} ((φ φ ((x x)))) dx dx}{&Integral; &Integral; ((&Integral; &Integral; K K ((y the y - - x x)) dy dy)) {M m}_{i i} ((φ φ ((x x)))) dx dx} . .

6. the correction method of cupping artifact in a kind of cone beam CT based on level set according to claim 1, is characterized in that, the formula (13) of described step 7 is:

{f f}_{p p} = = f f - - {\overset{^^}{f f}}_{s the s} . .

7. The method for correcting cupping artifacts in a level set-based cone beam CT according to claim 1, wherein the method adaptively corrects the cupping artifacts of cone beam CT without artificial Intervention can automatically complete the correction; the method does not need to repeatedly scan the measured object; does not increase the complexity of the cone beam CT system; for the reconstructed slice image, it can directly face the user, without the need for the original cone beam CT original There is no equipment to make any changes, and the calibration work can be completed.

8 . The method for correcting cupping artifacts in cone beam CT based on level set according to claim 1 , wherein the method is applied to correction of cone beam CT slice images. 9 .