阅读说明：本技术 一种基于曲面全变差的低剂量ct投影域去噪及重建方法 (Low-dose CT projection domain denoising and reconstructing method based on curved surface total variation ) 是由 江颖 刘婷 刘伟锋 吴锐帆 于 2020-12-17 设计创作，主要内容包括：本发明涉及一种基于曲面全变差的低剂量CT投影域去噪及重建方法,步骤一：根据投影域数据噪声的产生机制,建立优化泛函数据项Ψ；建立包括基于l-(1/2)范数的曲面全变差正则项的优化泛函正则项Υ,将优化泛函数据项和优化泛函正则项γ之和作为去噪优化泛函；步骤二：对步骤一中的去噪优化泛函进行求解,得到降噪后的CT投影域数据；步骤三：基于连续积分方程模型,根据步骤二中降噪后的CT投影域数据进行图像重建。本方法在进行投影域去噪的时候,建立基于l-(1/2)范数的曲面全变差正则项来提高去噪的精度。在完成投影域数据去噪后,采用连续积分方程模型来进行图像重建,加快了图像重建的速度和提高图像重建的精度。(The invention relates to a low-dose CT projection domain denoising and reconstructing method based on curved surface total variation, which comprises the following steps: establishing an optimized functional data item psi according to a generation mechanism of projection domain data noise; the establishing includes based on 1/2 The optimization functional regular term gamma of the norm surface total variation regular term is used as a denoising optimization functional, and the sum of the optimization functional data term and the optimization functional regular term gamma is used as a denoising optimization functional; step two: solving the denoising optimization functional in the first step to obtain denoising CT projection domain data; step three: and based on the continuous integral equation model, carrying out image reconstruction according to the CT projection domain data subjected to noise reduction in the second step. The method is based on l when the projection domain denoising is carried out 1/2 And the regular term of the total variation of the curved surface of the norm is used for improving the denoising precision. After the projection domain data denoising is finished, the continuous integral equation model is adopted to reconstruct the image, so that the image reconstruction speed is increased and the image reconstruction precision is improved.)

1. A low-dose CT projection domain denoising and reconstructing method based on curved surface total variation is characterized by comprising the following steps:

the method comprises the following steps: establishing an optimized functional data item psi according to a generation mechanism of projection domain data noise; the establishing includes based on_1/2The optimization functional regular term gamma of the norm surface total variation regular term is used as a denoising optimization functional, and the sum of the optimization functional data term gamma and the optimization functional regular term gamma is used as a denoising optimization functional;

step two: solving the denoising optimization functional in the first step to obtain denoising CT projection domain data;

step three: and based on the continuous integral equation model, carrying out image reconstruction according to the CT projection domain data subjected to noise reduction in the second step.

2. The method as claimed in claim 1, wherein in the step one, optimizing the functional regularization term γ further comprises l-based denoising and reconstructing₁And a norm curvelet transform coefficient regular term.

3. The method as claimed in claim 2, wherein the method of establishing the regularization term of the total variation of the curved surface is to consider the CT projection domain data as CT projection domain dataSurfaces in space, constraining the jump of the surface in the direction of principal curvature, grading the direction of principal curvatureThe norm is used as a surface total variation regular term; the concrete expression is as follows:

in the formula, Y is input projection domain data with noise; b is a staggered difference operator; s is a principal curvature direction matrix;is a transverse difference operator;a longitudinal difference operator; l is the lipschitz continuity coefficient; m is the number of rows of Y.

4. The method for denoising and reconstructing a low-dose CT projection domain based on total variation of a curved surface as claimed in claim 3, wherein the method for establishing the regularization term of the curvelet transform coefficient is to perform discrete curvelet transform on the CT projection domain data Y to obtain the coefficient CY, and the regularization term of the curvelet transform coefficient is expressed as: gamma ray₂(Y)＝||CY||₁。

5. The method for denoising and reconstructing a low-dose CT projection domain based on total variation of a curved surface as claimed in claim 4, wherein the optimization functional data item Ψ in the first step is specifically:

wherein Y is input projection domain data with noise; q is the number of photons reaching the detection plate; p is the number of photons received by the detection plate; i is₀Is the X-ray incident intensity.

6. The method for denoising and reconstructing a low-dose CT projection domain based on total variation of a curved surface as claimed in claim 5, wherein the denoising optimization functional in the first step is specifically:

in the formula I₀Intensity of X-ray incidence; lambda [ alpha ]₁And λ₂Taking the penalty coefficients of the two regularization terms as positive real numbers; c is a curvelet transform operator; s is a main curvature direction operator; and B is an alternating difference operator.

7. The method for denoising and reconstructing a low-dose CT projection domain based on total variation of a curved surface as claimed in claim 2, wherein the specific process in the third step is:

s3.1: the optimized functional for image reconstruction is established as follows:

wherein f is an image reconstruction function; g is a parallel beam projection domain function;is Radon transformation; λ is the regularization term coefficient; phi is a constant greater than zero and L-lipschitz is a continuous function; b represents a linear transformation;

s3.2: solving an optimized functional of image reconstruction;

s3.3: and (4) carrying out iterative solution on the optimized functional of image reconstruction to obtain an image reconstruction function f.

8. The method for denoising and reconstructing a low-dose CT projection domain based on total variation of curved surface as claimed in claim 7, wherein in said step S3.2, an optimization functional of image reconstruction is solved by using a neighboring gradient method.

9. The method for denoising and reconstructing a low-dose CT projection domain based on total variation of a curved surface according to claim 8, wherein in the step S3.3, an iterative formula is adopted to iteratively solve an optimization functional of image reconstruction, and the iterative formula specifically is as follows:

wherein the content of the first and second substances,representing a fourier transform;representing an inverse fourier transform; prox_ΦIs the neighboring gradient operator of Φ; l is the lipschitz continuity coefficient of phi; e represents the multiplication of vectors by elements; b^(k+1)Represents the intermediate result of the iteration of step k +1, z^(k+1)Representing the image function at step k + 1.

10. The method for denoising and reconstructing a low-dose CT projection domain based on total surface variation according to any one of claims 1-9, wherein in the second step, the optimization functional in the first step is solved by using an ADMM optimization algorithm with random step size.

Technical Field

The invention relates to the field of CT image processing methods, in particular to a low-dose CT projection domain denoising and reconstructing method based on curved surface total variation.

Background

Computed Tomography (CT) can obtain an anatomical image of the inside of a human body in a non-invasive manner, clearly display structural tissues and lesions of each part, and help a doctor to perform clinical diagnosis, and is an imaging technique widely applied to clinical examination. However, during the CT scan, the X-ray is required to penetrate the human body to acquire projection data of various angles. Many studies have shown that X-ray radiation generated during CT scanning induces metabolic abnormalities and even carcinogenic risks in human body, and its influence on human health is not negligible. To address the problem of iatrogenic radiation damage from standard CT scans, Naidich et al introduced Low Dose CT (LDCT) in the 90's of the 20 th century. LDCT reduces the radiation dose of X-rays by lowering the milliampere-second (mAs) setting in the CT scan setting, i.e. by lowering the current through the X-ray tube or by shortening the exposure time. However, lowering the mAs setting degrades the distribution obeyed by photons arriving at the detector from a gaussian distribution to a poisson distribution. The raw data generation process of LDCT is contaminated with a relatively higher level of noise than standard CT.

Total Variation (TV) regularization widely adopted in the field of image noise reduction at present assumes that a real image is a slicing constant, and for example, a low-dose energy spectrum CT image denoising method is disclosed in a patent publication with publication number "CN 104574416A" and publication date 2015, 4 and 29. However, the CT projection domain data usually does not satisfy this assumption, so that a direct application of the total variation regularization method to the projection domain noise reduction problem easily generates a step artifact, which affects the noise reduction accuracy. Meanwhile, the image needs to be reconstructed after the image is denoised, most of the current image reconstruction methods are based on discrete model algorithms, the solved image function solution belongs to the category of a fragmentation constant function, and the required reconstruction time is longer. Furthermore, the optimization model derived using this method naturally has model errors due to the error introduced by the discretization process.

Disclosure of Invention

The invention provides a low-dose CT projection domain denoising and reconstructing method based on curved surface total variation, aiming at solving the problems of low denoising precision and long reconstruction time caused by step artifacts in low-dose CT denoising in the prior art, and realizing high-precision denoising and accelerating image reconstruction speed.

In order to solve the technical problems, the invention adopts the technical scheme that: a low-dose CT projection domain denoising and reconstructing method based on curved surface total variation comprises the following steps:

the method comprises the following steps: establishing an optimized functional data item psi according to a generation mechanism of projection domain data noise; the establishment includes based onThe optimization functional regular term gamma of the norm surface total variation regular term takes the sum of the optimization functional data term and the optimization functional regular term gamma as an optimization functional;

step two: solving the optimized functional in the first step to obtain the CT projection domain data after noise reduction;

step three: and based on the continuous integral equation model, carrying out image reconstruction according to the CT projection domain data subjected to noise reduction in the second step.

In the above technical solution, the noise in the input projection domain data is removed by establishing an optimization functional including a regularization term and a data term and performing optimization solution on the optimization functional, so as to implement denoising of the CT projection domain. The main curvature direction contains first-order gradient information and second-order gradient information, the direction with large change intensity of a certain point of an image and the smooth direction are more accurate than the first-order gradient direction, artifacts are avoided in the image, and the image denoising is more accurate. A continuous integral equation model is adopted to design a CT image reconstruction algorithm during image reconstruction, and the continuous integral equation is based on an actual physical model and has higher precision and reconstruction speed.

Preferably, in said first step, optimizing functional regularization term further comprises optimizing functional regularization term based on yAnd a norm curvelet transform coefficient regular term. According to the sparsity of curvelet coefficients of CT projection domain data, based onThe regular term of the curvelet transform coefficient of the norm can improve the denoising effect.

Preferably, the method for establishing the curved surface total variation regularization term comprises the following steps: treating CT projection domain data asSurfaces in space, constraining the jump of the surface in the direction of principal curvature, grading the direction of principal curvatureThe norm is used as a surface total variation regular term, and is specifically expressed as:

in the formula, Y is input projection domain data with noise; b is a staggered difference operator; s is a principal curvature direction matrix;is a transverse difference operator;a longitudinal difference operator; l is the lipschitz continuity coefficient; m is the number of rows of Y.

In particular, a principal curvature direction matrix is definedExpressed as:

in the formula, v₁And v₂Are respectively provided withA maximum principal curvature and a minimum principal curvature direction of Y; d_YρIs the principal curvature absolute difference.

Preferably, the method for establishing the regularization term of the curvelet transform coefficient is to perform discrete curvelet transform on the CT projection domain data Y to obtain the coefficient CY, which is specifically expressed as: gamma ray₂(Y)＝||CY||₁。

Preferably, the optimization functional data item Ψ in the step one is specifically:

wherein Y is input projection domain data with noise; q is the number of photons reaching the detection plate; p is the number of photons received by the detection plate; i is₀Is the X-ray incident intensity.

Preferably, the optimization functional in the step one is specifically:

in the formula I₀Intensity of X-ray incidence; lambda [ alpha ]₁And λ₂Taking the penalty coefficients of the two regularization terms as positive real numbers; c is a curvelet transform operator; s is a main curvature direction operator; and B is an alternating difference operator.

Preferably, the specific process in the third step is as follows:

s3.1: the optimized functional for image reconstruction is established as follows:

wherein f is an image reconstruction function; g is a parallel beam projection domain function;is Radon transformation; λ is the regularization term coefficient; phi is a constant greater than zero and L-lipschitz is a continuous function; b represents a linear transformation;

s3.2: solving an optimized functional of image reconstruction;

s3.3: and (4) carrying out iterative solution on the optimized functional of image reconstruction to obtain an image reconstruction function f.

On the basis of a continuous integral equation model, Radon transformation is introduced into the establishment of an optimization functional of image reconstruction, and the subsequent iterative algorithm is accelerated through the convolution property of the Radon transformation, so that the image reconstruction speed is accelerated, and the reconstruction time is shorter.

Preferably, in step S3.2, an optimization functional of image reconstruction is solved by using a neighboring gradient method, and the specific process is as follows:

first of all, calculateOf the gradient of (c). Is provided withThenWhereinIs the conjugate operator of the Radon transform. By usingThe convolution property of (a) and the fourier transform accomplish an accelerated computation.

Computing Is f at L²A certain set of basis functions in space, z_iAre the coefficients of the basis functions and,is the detection range of the detection plate,thereby can be solvedWherein the acceleration can be achieved by using the convolution propertyCalculating, setting a matrixElements thereofSince W is a circulant matrix, thereforeW is the first row of the matrix W, z ═ z₁,z₂,…,z_L]^T。

Computing Wherein

Representing the fourier transform of g (θ, u) with respect to the u component.

Preferably, in step S3.3, an iterative formula is used to iteratively solve the optimization functional of image reconstruction, where the iterative formula specifically includes:

wherein the content of the first and second substances,representing a fourier transform;representing an inverse fourier transform; prox_ΦIs the neighboring gradient operator of Φ; l is the lipschitz continuity coefficient of phi; e represents the multiplication of vectors by elements; b^(k+1)Represents the intermediate result of the iteration of step k +1, z^{(k +1)}Representing the image function at step k + 1.

Preferably, in the second step, the optimization functional in the first step is solved by using an ADMM optimization algorithm with random step sizes. The convergence speed of the optimization algorithm can be accelerated by adopting the random step length.

Compared with the prior art, the invention has the beneficial effects that: the method is based on the establishment of the projection and the denoisingThe method comprises the steps of introducing a curved surface total variation regular term of norm, replacing a gradient direction with a principal curvature direction, and avoiding the occurrence of lithography in an image by combining the properties of a projection domain, thereby improving the denoising precision. After the projection domain data denoising is finished, the continuous integral equation model is adopted to reconstruct the image, so that the image reconstruction speed is increased and the image reconstruction precision is improved.

Drawings

FIG. 1 is a schematic overall flow chart of a low-dose CT projection domain denoising and reconstructing method based on curved surface total variation according to the present invention;

FIG. 2 is a noisy CT image;

FIG. 3 is a high dose CT image;

fig. 4 is a CT image of fig. 2 processed by the present method.

Detailed Description

The drawings are for illustrative purposes only and are not to be construed as limiting the patent; for the purpose of better illustrating the embodiments, certain features of the drawings may be omitted, enlarged or reduced, and do not represent the size of an actual product; it will be understood by those skilled in the art that certain well-known structures in the drawings and descriptions thereof may be omitted. The positional relationships depicted in the drawings are for illustrative purposes only and are not to be construed as limiting the present patent.

The technical scheme of the invention is further described in detail by the following specific embodiments in combination with the attached drawings:

example 1

Fig. 1 shows an embodiment of a low-dose CT projection domain denoising and reconstructing method based on total variation of curved surface, which includes the following steps:

the method comprises the following steps: establishing an optimized functional data item psi according to a generation mechanism of projection domain data noise; the establishment includes based onThe total variation regularization term of the norm on the curved surfaceThe optimization functional regular term gamma of the standard curve-wave transformation coefficient regular term takes the sum of the optimization functional data term and the optimization functional regular term gamma as an optimization functional;

the method for establishing the curved surface total variation regular term comprises the following steps: treating CT projection domain data asSurfaces in space, constraining the jump of the surface in the direction of principal curvature, grading the direction of principal curvatureThe norm is used as a surface total variation regular term, and is specifically expressed as:

in the formula, Y is input projection domain data with noise; b is a staggered difference operator; s is a principal curvature direction matrix;is a transverse difference operator;a longitudinal difference operator; l is the lipschitz continuity coefficient; m is the number of rows of Y.

The method for establishing the regularization term of the curvelet transform coefficient comprises the following steps of carrying out discrete curvelet transform on CT projection domain data Y to obtain a coefficient CY, wherein the coefficient CY is specifically represented as: gamma ray₂(Y)＝||CY||₁。

Optimizing the functional into optimizing functional data item psi and curved surface total variation regular item gamma₁(Y) and the Curve-wave transform coefficient regular term γ₂(Y) the sum of the three is as follows:

wherein Y is input projection domain data with noise; q is the number of photons reaching the detection plate; p is the number of photons received by the detection plate; i is₀Is the X-ray incident intensity.

Step two: solving the optimization functional in the first step by using an ADMM optimization algorithm with random step length to obtain the CT projection domain data after noise reduction;

step three: and based on a continuous integral equation model, reconstructing an image according to the CT projection domain data subjected to noise reduction in the second step, wherein the specific process is as follows:

s3.1: the optimized functional for image reconstruction is established as follows:

wherein f is a drawingAn image reconstruction function; g is a parallel beam projection domain function;is Radon transformation; λ is the regularization term coefficient; phi is a constant greater than zero and L-lipschitz is a continuous function; b represents a linear transformation;

s3.2: solving an optimized functional of image reconstruction by adopting a neighboring gradient method;

s3.3: and (3) carrying out iterative solution on the optimized functional of image reconstruction to obtain an image reconstruction function f, wherein a specific iterative formula is as follows:

wherein the content of the first and second substances,representing a fourier transform;representing an inverse fourier transform; prox_ΦIs the neighboring gradient operator of Φ; l is the lipschitz continuity coefficient of phi; e represents the multiplication of vectors by elements; b^(k+1)Representing the intermediate result of the iteration of the (k + 1) th step; z is a radical of^{(k +1)}Representing the image function at step k + 1.

Specifically, in step S3.2, an optimization functional of image reconstruction is solved by using a neighboring gradient method, and a specific process is as follows:

first of all, calculateOf the gradient of (c). Is provided withThenWhereinIs the conjugate operator of the Radon transform. By usingThe convolution property of (a) and the fourier transform accomplish an accelerated computation.

Computing Is f at L²A certain set of basis functions in space, z_iAre the coefficients of the basis functions and,is the detection range of the detection plate,thereby can be solvedWherein the acceleration can be achieved by using the convolution propertyCalculating, setting a matrixElements thereofSince W is a circulant matrix, thereforeW is the first row of the matrix W, z ═ z₁,z₂,…,z_L]^T。

Computing Wherein Representing the fourier transform of g (θ, u) with respect to the u component.

The working principle of the embodiment is as follows: and removing noise in the input projection domain data by establishing an optimization functional comprising a regular term and a data item and performing optimization solution on the optimization functional, so as to realize the denoising of the CT projection domain. The main curvature direction contains first-order gradient information and second-order gradient information, the direction with large change intensity of a certain point of an image and the smooth direction are more accurate than the first-order gradient direction, artifacts are avoided in the image, and the image denoising is more accurate. A continuous integral equation model is adopted to design a CT image reconstruction algorithm during image reconstruction, and the continuous integral equation is based on an actual physical model and has higher precision and reconstruction speed.

The method has the beneficial effects that the method is established based on the projection and denoising in the process of projection and denoisingThe method comprises the steps of introducing a curved surface total variation regular term of norm, replacing a gradient direction with a principal curvature direction, and avoiding the occurrence of lithography in an image by combining the properties of a projection domain, thereby improving the denoising precision. After the projection domain data denoising is finished, the continuous integral equation model is adopted to reconstruct the image, so that the image reconstruction speed is increased and the image reconstruction precision is improved. As shown in FIG. 4, the low-dose noisy CT image of FIG. 2 is subjected to de-noising reconstruction by the method to obtain a CT image close to FIG. 3. The CT image of FIG. 3 is a high dose CT image, which is a low dose CT imageDose (noisy) images are denoised to reconstruct the object.

It should be understood that the above-described embodiments of the present invention are merely examples for clearly illustrating the present invention, and are not intended to limit the embodiments of the present invention. Other variations and modifications will be apparent to persons skilled in the art in light of the above description. And are neither required nor exhaustive of all embodiments. Any modification, equivalent replacement, and improvement made within the spirit and principle of the present invention should be included in the protection scope of the claims of the present invention.