CN111652950B - Compressed X-ray tomosynthesis method of multi-light-source snapshot - Google Patents

Abstract

A multi-light source snapshot compression X-ray tomosynthesis method, based on the PRISM model to establish a nonlinear CXT reconstruction framework, and then propose an improved split Bregman algorithm to solve the nonlinear inverse reconstruction problem, so as to obtain the detected Three-dimensional object image; since the present invention allows multiple light sources to illuminate the object at the same time, compared with the traditional sequential illumination method, the scanning time can be effectively shortened, the problem of image distortion caused by the deformation of the object can be overcome, and the reconstruction algorithm can be greatly improved. Robustness to object motion or perturbation.

Description

A Compressed X-ray Tomosynthesis Method for Multiple Light Source Snapshots

technical field

The invention belongs to the technical field of computational imaging, and in particular relates to a compressed X-ray tomosynthesis method for snapshots of multiple light sources.

Background technique

At present, X-ray computed tomography (computed tomography, CT for short) has been widely used in many fields, such as medical diagnosis, safety inspection, industrial non-destructive testing and other fields. Traditional CT needs to collect complete measurement data around the detected object, and then use the filtered back-projection algorithm (filered back-projection gorithm, FBP for short) to reconstruct the attenuation coefficient inside the object. However, in the medical field, how to obtain high-resolution reconstructed images while reducing the radiation dose is a very challenging requirement. X-ray tomosynthesis (XT) is an alternative to CT that does not require complete projection measurements, and uses simple perceptual geometry to reconstruct the object under test from a set of incomplete projections. Since projections only need to be obtained from a limited angle, the XT technique has many advantages, such as greatly reducing the radiation dose, shortening the detection time, and simplifying the sensing geometry. The incomplete measurement of XT can reduce a lot of radiation, but it will inevitably cause the problem of underdetermined inverse reconstruction. As an emerging technology, compressive X-ray tomosynthesis (CXT) enables us to solve the underdetermined inverse reconstruction problem by using incomplete measurements and ensure better image reconstruction quality. Traditional CXT technology can reconstruct three-dimensional objects from two-dimensional projection data generated by a set of X-ray light sources. The coded aperture is used to modulate the illumination structure to reduce radiation dose, and to adjust the structure of the perception matrix to improve image reconstruction. quality. Specifically, the traditional CXT system uses a movable cone-beam X-ray light source or a set of cone-beam X-ray light sources fixed at different positions to illuminate the object from different angles in turn, and collect multiple two-dimensional projection information. In these methods, each detector unit only receives X-rays from one light source at a certain moment, so as to obtain non-overlapping measurement information, so the imaging model has the form of a linear equation. Then, based on the sparsity assumption and a linear compressive sensing (CS) framework, the attenuation coefficients of the detected objects can be numerically reconstructed from the 2D measurements. However, these measurement methods suffer from the following drawbacks. First of all, taking pictures in sequence from different angles will prolong the scanning time of the test. Some seriously ill patients may not be able to keep still for a long time, or they may feel uncomfortable during long-term testing. In addition, some living biological tissues may not be able to remain absolutely still within the detection time. However, using the traditional CXT method, once the object is deformed, the reconstruction result of the object will be distorted, even severe image distortion.

To sum up, traditional CXT cannot meet the needs of specific situations. In order to avoid the influence of object motion or disturbance and shorten the measurement time, we propose a nonlinear CXT technique using snapshots of multiple light sources.

SUMMARY OF THE INVENTION

In order to solve the above problems, the present invention provides a compressed X-ray tomosynthesis method for snapshots of multiple light sources, which can effectively shorten the scanning time and overcome the problem of image distortion caused by object deformation.

A compressed X-ray tomosynthesis method for snapshots of multiple light sources, comprising the following steps:

S1: The revised framework for constructing a refactoring framework based on the PRISM model is as follows:

表示X射线光源对三维物体进行扫描后得到的衰减系数向量的优化值，x表示X射线光源对三维物体进行扫描后得到的衰减系数向量的理论值，y_k表示第k次扫描中各X射线光源经编码孔径调制后在探测器上的实际成像模型，也即第k次扫描时探测器上采集的实际强度值，且k＝1,2,…,K，K表示扫描的总次数，f_k(x)表示第k次扫描中各X射线光源经编码孔径调制后在探测器上的模拟成像模型，也即第k次扫描时探测器上通过模拟计算得到的强度值，e_k表示y_k与f_k(x)的残差，λ_*表示低秩正则项的权重系数，||·||_*表示对矩阵秩操作的正则化核函数，T_T表示将辅助向量d转化成与三维物体的张量维度相同的张量，T_M(·)表示将张量转化成衰减系数矩阵，μ*和μ₁表示设定权重系数，λ₁表示稀疏正则项的权重系数，θ表示稀疏系数向量，v表示衡量辅助向量d和衰减系数向量x之间距离的辅助向量，w表示衡量辅助向量c和稀疏系数向量θ之间距离的辅助向量；in,

represents the optimal value of the attenuation coefficient vector obtained after the X-ray light source scans the three-dimensional object, x represents the theoretical value of the attenuation coefficient vector obtained after the X-ray light source scans the three-dimensional object, and y _k represents each X-ray in the kth scan The actual imaging model on the detector after the light source is modulated by the coded aperture, that is, the actual intensity value collected on the detector during the kth scan, and k=1,2,…,K, K represents the total number of scans, f _k (x) represents the simulated imaging model on the detector after each X-ray light source is modulated by the coded aperture in the kth scan, that is, the intensity value obtained by the simulation calculation on the detector in the kth scan, e _k represents y The residual of _k and f _k (x), λ _* represents the weight coefficient of the low-rank regular term, || · || _* represents the regularization kernel function of the matrix rank operation, T _T represents the transformation of the auxiliary vector d into a three-dimensional A tensor with the same tensor dimension of the object, T _M ( ) represents converting the tensor into a matrix of attenuation coefficients, μ* and μ ₁ represent the set weight coefficient, λ ₁ represents the weight coefficient of the sparse regular term, θ represents the sparse coefficient vector, v represents the auxiliary vector measuring the distance between the auxiliary vector d and the attenuation coefficient vector x, w represents the auxiliary vector measuring the distance between the auxiliary vector c and the sparse coefficient vector θ;

S2: The iterative formula for obtaining the attenuation coefficient vector x, the residual _ek , the auxiliary vector d, the auxiliary vector v, the auxiliary vector c and the auxiliary vector w in the revised framework is as follows:

v ⁿ⁺¹ =v ⁿ +x ⁿ⁺¹ -d ⁿ⁺¹

w ⁿ⁺¹ =w ⁿ +θ ⁿ⁺¹ -c ⁿ⁺¹

θ ⁿ⁺¹ = ^ψT x ⁿ⁺¹

Among them, n represents the number of iterations, ψ represents the sparse coefficient, and T represents the transpose;

S3: Set the initial value of each iteration formula in step S2 to 0, and repeat the iteration formula in step S2 continuously, wherein, after each iteration is completed, it is judged whether the attenuation coefficient vector x ⁿ⁺¹ obtained by this iteration meets the setting requirements , the setting requirements include that the difference between the attenuation coefficient vectors obtained from two adjacent iterations is less than the set precision or the number of iterations reaches the set upper limit. If satisfied, the attenuation coefficient vector x ⁿ⁺¹ obtained in this iteration is The optimal solution is taken as the image of the compressed X-ray tomosynthesis; if not, the iteration is continued until the attenuation coefficient vector x ⁿ⁺¹ meets the set requirements.

Further, the construction method of the described reconstruction framework based on the PRISM model is as follows:

S11: Construct a non-limiting optimization function for reconstructing the internal attenuation coefficient of the three-dimensional object according to K scans as follows:

S12: Expand the tensor of the three-dimensional object into a two-dimensional matrix, and then perform a low-rank restriction operation on the two-dimensional matrix to obtain the regular term R(x):

R(x)=λ _* ||T _M (χ)|| _* +λ ₁ ||Τ _S (x)|| ₁ =λ _* ||X|| _* +λ ₁ ||θ ₁

Wherein, T _M (x) represents converting tensor x into attenuation coefficient matrix X, and T _S (x) represents converting attenuation coefficient vector x into corresponding sparse coefficients;

S13: Build a reconstruction framework based on the PRISM model according to the unrestricted optimization function and regular term:

Further, the calculation method of the imaging model f _k (x) is:

表示第k次扫描时打开的X射线光源索引所组成的集合，b_p表示第p个光源对其下方的编码孔径进行扫描后得到的向量，⊙表示点乘运算，H_p表示第p个X射线光源对应的系统矩阵。in,

represents the set of indices of the X-ray light sources turned on in the kth scan, b _p represents the vector obtained after the pth light source scans the coded aperture below it, ⊙ represents the point multiplication operation, and H _p represents the pth X The system matrix corresponding to the ray source.

Further, the acquisition method of the vector b _p is specifically:

S14: Construct the cost function for optimizing the coded aperture as follows:

M为探测器元素的个数，N为物体体素的个数，

为经过编码孔径调制后第i个探测器元素的平均测量强度实际值，m₁为所有探测器元素的平均测量强度理论值，

为经过编码孔径调制后的第j个体素的平均感知信息实际值，m₂为每个体素平均感知信息量理论值，s_t为第i个探测器元素对应的复用编码，m₃为每个探测器元素对应的透光编码孔径元素数量的理论值；in,

M is the number of detector elements, N is the number of object voxels,

is the actual value of the average measured intensity of the i-th detector element after coded aperture modulation, m ₁ is the theoretical value of the average measured intensity of all detector elements,

is the actual value of the average perceptual information of the j-th voxel after coded aperture modulation, m ₂ is the theoretical value of the average perceptual information of each voxel, s _t is the multiplexing code corresponding to the i-th detector element, m ₃ is each The theoretical value of the number of light-transmitting coded aperture elements corresponding to each detector element;

S15: Use the DSB algorithm to solve the cost function to obtain the optimized b ₁ , . . . , b _p .

Further, a compressed X-ray tomosynthesis method for snapshots of multiple light sources, further comprising the following steps:

S4: Construct the objective function F according to the optimal solution x ⁿ + ¹ :

S5: The objective function F is optimized by the conjugate gradient method, and the obtained optimized value x* is used as an image of compressed X-ray tomography.

Beneficial effects:

A multi-light source snapshot compression X-ray tomosynthesis method, based on the PRISM model to establish a nonlinear CXT reconstruction framework, and then propose an improved split Bregman algorithm to solve the nonlinear inverse reconstruction problem, so as to obtain the detected Three-dimensional object image; since the present invention allows multiple light sources to illuminate the object at the same time, compared with the traditional sequential illumination method, the scanning time can be effectively shortened, the problem of image distortion caused by the deformation of the object can be overcome, and the reconstruction algorithm can be greatly improved. Robustness to object motion or perturbation.

Description of drawings

Fig. 1 is a flow chart of a compressed X-ray tomosynthesis method of a multi-light source snapshot;

2 is a schematic diagram of a traditional CXT system and a multi-light source snapshot CXT system proposed by the present invention;

Figure 3 is a schematic diagram of the initial object and the deformation of the object during irradiation;

Figure 4 shows the reconstruction results obtained by using the traditional sequential irradiation method when the object is deformed;

5 is an object reconstruction result obtained by using the nonlinear CXT method proposed by the present invention when a random coded aperture modulated light source is used;

6 is a schematic diagram of an optimization result of a coded aperture;

FIG. 7 is a schematic diagram of an object reconstruction result obtained by using the nonlinear CXT method proposed by the present invention when an optimized coded aperture is used.

Detailed ways

In order to make those skilled in the art better understand the solutions of the present application, the following will clearly and completely describe the technical solutions in the embodiments of the present application with reference to the accompanying drawings in the embodiments of the present application.

As shown in Figure 1, a method for compressed X-ray tomosynthesis of multi-light source snapshots includes the following steps:

S1: The revised framework for constructing the reconstruction framework based on the PRISM model (prior rank, intensity and sparsity model, prior low rank, intensity and sparsity model) is as follows:

表示X射线光源对三维物体进行扫描后得到的衰减系数向量的优化值，x表示X射线光源对三维物体进行扫描后得到的衰减系数向量的理论值，y_k表示第k次扫描中各X射线光源经编码孔径调制后在探测器上的实际成像模型，也即第k次扫描时探测器上采集的实际强度值，且k＝1,2,…,K，K表示扫描的总次数，f_k(x)表示第k次扫描中各X射线光源经编码孔径调制后在探测器上的模拟成像模型，也即第k次扫描时探测器上通过模拟计算得到的强度值，e_k表示对应于第k次扫描的辅助向量，实际为y_k与f_k(x)的残差，λ_*表示低秩正则项的权重系数，||·||_*表示对矩阵秩操作的正则化核函数，T_T表示将辅助向量d转化成与三维物体的张量维度相同的张量，T_M(·)表示将张量转化成衰减系数矩阵，μ_*和μ₁表示设定权重系数，λ₁表示稀疏正则项的权重系数，θ表示稀疏系数向量，v表示衡量辅助向量d和衰减系数向量x之间距离的辅助向量，w表示衡量辅助向量c和稀疏系数向量θ之间距离的辅助向量，且d≈x，c≈θ；in,

represents the optimal value of the attenuation coefficient vector obtained after the X-ray light source scans the three-dimensional object, x represents the theoretical value of the attenuation coefficient vector obtained after the X-ray light source scans the three-dimensional object, and y _k represents each X-ray in the kth scan The actual imaging model on the detector after the light source is modulated by the coded aperture, that is, the actual intensity value collected on the detector during the kth scan, and k=1,2,…,K, K represents the total number of scans, f _k (x) represents the simulated imaging model on the detector after each X-ray light source is modulated by the coded aperture in the kth scan, that is, the intensity value obtained by the simulation calculation on the detector in the kth scan, and e _k represents the corresponding The auxiliary vector in the kth scan is actually the residual of y _k and f _k (x), λ _* represents the weight coefficient of the low-rank regular term, || · || _* represents the regularization kernel function of the matrix rank operation , T _T represents the transformation of the auxiliary vector d into a tensor with the same dimension as the tensor of the three-dimensional object, T _M ( ) represents the transformation of the tensor into the attenuation coefficient matrix, μ _* and μ ₁ represent the set weight coefficient, λ ₁ represents the weight coefficient of the sparse regular term, θ represents the sparse coefficient vector, v represents the auxiliary vector that measures the distance between the auxiliary vector d and the attenuation coefficient vector x, w represents the auxiliary vector that measures the distance between the auxiliary vector c and the sparse coefficient vector θ, And d≈x, c≈θ;

It should be noted that the correction framework is to impose low-rank constraints on the attenuation coefficient vector x, and the subsequent MSB algorithm needs to separate the regular term, and it is difficult to directly derive the attenuation coefficient vector x with a low-rank term. For this problem, therefore, it is necessary to introduce an auxiliary vector d, and indirectly constrain the attenuation coefficient vector x by constraining the auxiliary vector d; for the same reason, it is also necessary to introduce an auxiliary vector c; in addition, it can be seen from the subsequent iterative formula that the auxiliary vector v is not The exact distance between the auxiliary vector d and the attenuation coefficient vector x, the auxiliary vector w is not the exact distance between the auxiliary vector c and the sparse coefficient vector θ, and the auxiliary vector v and the auxiliary vector w are also related to the previous iterative steps, are A cumulative variable, so v and w are just an auxiliary vector to measure distance.

S2: adopt the modified split Bregman algorithm (modified split Bregman algorithm, MSB for short) to obtain the attenuation coefficient vector x, the residual _ek , the auxiliary vector d, the auxiliary vector v, the auxiliary vector c and the The iterative formula of the auxiliary vector w is as follows:

v ⁿ⁺¹ =v ⁿ +x ⁿ⁺¹ -d ⁿ⁺¹

w ⁿ⁺¹ =w ⁿ +θ ⁿ⁺¹ -c ⁿ⁺¹

θ ⁿ⁺¹ = ^ψT x ⁿ⁺¹

Among them, n represents the number of iterations, ψ represents the sparse coefficient, and T represents the transpose;

S3: Set the initial value of each iteration formula in step S2 to 0, and repeat the iteration formula in step S2 continuously, wherein, after each iteration is completed, it is judged whether the attenuation coefficient vector x ⁿ⁺¹ obtained by this iteration meets the setting requirements , the setting requirements include that the difference between the attenuation coefficient vectors obtained from two adjacent iterations is less than the set precision or the number of iterations reaches the set upper limit. If satisfied, the attenuation coefficient vector x ⁿ⁺¹ obtained in this iteration is The optimal solution is taken as the image of the compressed X-ray tomosynthesis; if not, the iteration is continued until the attenuation coefficient vector x ⁿ⁺¹ meets the set requirements.

It should be noted that, in order to further improve the accuracy of the reconstructed image, the present invention proposes an MSB method. After an approximate solution is obtained by using the iterative steps in step S3, a subsequent processing method is used to make the attenuation coefficient vector closer to Optimal solution. The subsequent processing method is specifically:

S4: Construct the objective function F according to the optimal solution x ⁿ⁺¹ :

S5: The objective function F is optimized by the conjugate gradient method, and the obtained optimized value x ^* is used as the image of the compressed X-ray tomography.

It can be seen that the present invention can make the calculated measurement value closer to the actual measurement value, thereby making the attenuation coefficient vector closer to the attenuation coefficient of the real three-dimensional object.

The following describes the construction method of the reconstruction framework based on the PRISM model in detail, including the following steps:

表示三维物体的张量表示；令

表示对三维物体衰减系数进行扫描后得到的向量；将探测器栅格化为M_x×M_y个元素，其中M_x和M_y分别为沿着x和y轴的维度；设CXT系统中共存在P个X射线光源，CXT成像过程中总共对物体照射K次(K可以等于1)，每一次均可采用P个光源中的某一个照射物体，或者采用P个光源中的多个光源同时照射；

表示对第k(k＝1,2……K)次快照中探测器的测量强度值进行扫描后得到的向量。S101: As shown in Figure 2, rasterize the three-dimensional object into a data cube of N _x ×N _y ×N _z , where N _x , N _y and N _z are the dimensions of the object along the x, y and z axes, respectively; make

represents a tensor representation of a three-dimensional object; let

Represents the vector obtained after scanning the attenuation coefficient of the 3D object; rasterizes the detector into M _x ×M _y elements, where M _x and M _y are the dimensions along the x and y axes, respectively; let the CXT system coexist There are P X-ray light sources. During the CXT imaging process, the object is irradiated K times in total (K can be equal to 1). Each time, one of the P light sources can be used to illuminate the object, or multiple light sources of the P light sources can be used to illuminate the object at the same time. ;

Represents a vector obtained by scanning the measured intensity values of the detector in the kth (k=1, 2...K) snapshot.

That is to say, the present invention allows multiple X-ray light sources to scan the object at the same time, so compared with the traditional sequential irradiation method, the detection time can be effectively shortened; at the same time, in the traditional sequential irradiation, if the detected object is deformed, It will lead to distortion of the object reconstruction result, even serious image distortion; the present invention can alleviate or avoid this problem, and greatly improve the robustness of the reconstruction algorithm to object motion or disturbance.

表示。则第k次拍照的成像模型的离散形式可以表述为：S102: Set the set formed by the index of the light source turned on when taking the photo for the kth (k=1, 2...K) times as

express. Then the discrete form of the imaging model of the k-th photo can be expressed as:

where H _p is the system matrix corresponding to the pth X-ray light source.

考虑编码孔径的调制作用后，第k次拍照的成像模型可以表述为：S103: Consider placing an coded aperture under each light source, and the vector obtained after scanning the coded aperture under the pth light source is denoted as

After considering the modulation effect of the coded aperture, the imaging model of the k-th photo can be expressed as:

where ⊙ is the dot multiplication operation.

S104: Considering all K scans, the non-restrictive optimization problem for reconstructing the attenuation coefficient inside the object can be expressed as:

Among them, y _k represents the actual intensity value collected on the detector during the k-th photograph, and f _k (x) is the intensity value on the detector during the k-th photograph calculated by the imaging model, namely

S105: In order to further improve the reconstruction quality, a PRISM model is adopted, and a low-rank regular term and a sparse regular term are added to the objective function. The low-rank regularization term can effectively reduce image noise and artifacts; at the same time, the two regularization terms can also limit the solution space, so that the object image can be more accurately reconstructed from the underdetermined measurement model. To take advantage of the low-rank properties of 3D objects, we expand the 3D object tensor into a 2D matrix form, and then perform a low-rank restriction operation on this 2D matrix. Therefore, the regular term in the PRISM model can be expressed as:

R(x)=λ _* ||T _M (χ)|| _* +λ ₁ ||Τ _S (x)|| ₁ =λ _* ||X|| _* +λ ₁ ||θ ₁

where λ _* and λ ₁ are the weight coefficients of the low-rank regular term and the sparse regular term, respectively; T _M (χ) represents the transformation of the tensor χ into the attenuation coefficient matrix X; || · || _* represents the regularity of the matrix rank operation The kernel function is defined as the sum of all singular values of a certain matrix, namely ||X _i || _* =∑ _k σ _k ; Τ _S (x) represents the conversion of the attenuation coefficient vector x into the corresponding sparse coefficient, assuming that the test is to be The attenuation coefficient vector x of the object can be sparsely represented under a sparse basis Ψ, then the sparse coefficient vector obtained by the T _S operation is θ=Ψ ^T x.

S106: Combined with the unconstrained problem in step S104 and the regular term in step S105, the reconstruction framework based on the PRISM model can be expressed as:

Further, in the reconstruction problem, only a sub-optimal solution can be obtained using random coded apertures. Therefore, the optimization of the coded aperture is particularly important to improve the reconstruction quality of the image. In order to further improve the image reconstruction effect of nonlinear CXT and improve the reconstruction accuracy, the present invention designs a coded aperture optimization method based on uniform perception criteria. The three uniform perception criteria are as follows:

Criterion 1: Detector measures uniform perception of intensity. To fully utilize the information measured on the detector, each detector element should receive approximately the same measured intensity value, indicating approximately the same amount of information received. Since the information on each detector element is a superposition of information from multiple X-ray sources, this is a nonlinear model. Note that each row of the system matrix H corresponds to a particular detector element. It is assumed that the transmittance of the coded aperture is μ _T , so theoretically the detector can only receive μ _T × 100% of the amount of information on the detector without the coded aperture. When no coded aperture is added, and the object space is represented by an all-ones vector μ _N = [1, ^.

表示经过编码孔径调制后的实际平均测量强度值，则The theoretically uniform measured intensity value for all detector elements is therefore m ₁ =E(μ _T q), where E(·) is the operation to find the mathematical expectation. However, the actual measured intensity of the detector is the coded aperture modulated measured value. definition

represents the actual average measured intensity value after coded aperture modulation, then

The purpose of Criterion 1 is to try to homogenize the actual measured intensity on the detector so that it is close to the theoretical value m ₁ .

Criterion 2: Achieve uniform measurements in each data cube voxel. An X-ray should pass through a certain number of voxels. Thus, the definition vector r represents the perceived entropy of a data cube voxel when no coded aperture is added. Note that each column of the system matrix H represents the contribution coefficient of a particular voxel to the measured values of all detector units, then

是经过编码孔径调制后的平均体素感知信息，则where μ _M is an M-length all-ones vector, ie μ _M = [1, . . . , 1] ^T . Each element of the vector r represents the mean value of the system matrix columns corresponding to all P X-ray sources. Now considering the transmittance of the coded aperture, theoretically, the average perceived amount of information per voxel should be m ₂ =E(μ _T r). In order to achieve uniform measurements of cubic voxels, all elements of the vector r should have a minimum variance. definition

is the average voxel perception information after coded aperture modulation, then

The purpose of this criterion is to uniformly perceive the amount of information in the voxel, making it close to the theoretical value m ₂ .

是复用编码。因此，复用编码s的每一个元素值理论上应该接近于m₃。Criterion 3: Achieving uniform coding in multiplexed illumination of multiple light sources. Since multiple X-ray light sources illuminate the object to be measured at the same time, the attenuated X-rays hit the same detector, and the coded apertures under multiple light sources also affect the measurement value of the same detector. Therefore, uniform encoding should be achieved in the simultaneous illumination of multiple light sources. Considering P X-ray sources, over all coded apertures, the number of all light-transmitting coded aperture elements corresponding to a particular detector element should theoretically be m ₃ = μ _T P . definition vector

is a multiplexed code. Therefore, the value of each element of the multiplexed code s should theoretically be close to m ₃ .

Therefore, according to the previous three constraints, a cost function can be defined for optimizing the coded aperture, which includes the following steps:

S107: Construct a cost function for optimizing the coded aperture as follows:

由于第二项对应N个元素的和，因此

由于第三项对应M个元素的和，因此

M为探测器元素的个数，N为物体体素的个数，下标i,j,t表示对应向量的元素索引，具体的，

为经过编码孔径调制后第i个探测器元素的平均测量强度实际值，m₁为所有探测器元素的平均测量强度理论值，

为经过编码孔径调制后的第j个体素的平均感知信息实际值，m₂为每个体素平均感知信息量理论值，s_t为第i个探测器元素对应的复用编码，m₃为每个探测器元素对应的透光编码孔径元素数量的理论值，同时，m₁，m₂和m₃的值仅仅与X射线光源的数量和编码孔径的透过率有关；Among them, since the first term corresponds to the sum of M elements, so

Since the second term corresponds to the sum of N elements, so

Since the third term corresponds to the sum of M elements, so

M is the number of detector elements, N is the number of object voxels, and the subscripts i, j, t represent the element index of the corresponding vector. Specifically,

is the actual value of the average measured intensity of the i-th detector element after coded aperture modulation, m ₁ is the theoretical value of the average measured intensity of all detector elements,

is the actual value of the average perceptual information of the j-th voxel after coded aperture modulation, m ₂ is the theoretical value of the average perceptual information of each voxel, s _t is the multiplexing code corresponding to the i-th detector element, m ₃ is each The theoretical value of the number of light-transmitting coded aperture elements corresponding to each detector element, at the same time, the values of m ₁ , m ₂ and m ₃ are only related to the number of X-ray light sources and the transmittance of the coded aperture;

S108: Use a direct binary search method (direct binary search, DBS for short) to solve the cost function to obtain optimized b ₁ , . . . , b _p .

It should be noted that the DBS algorithm is a search method that iteratively evaluates the influence of each element change in a binary image on the cost function, and an optimized coded aperture pattern can be obtained. Using the optimized coded aperture to modulate the light intensity distribution of each X-ray light source, the measurement data can be obtained on the detector.

It can be seen that the existing CXT-based image reconstruction methods are all sequential illumination methods. The light sources at different positions are turned on in sequence to illuminate the object, and multiple times of information are collected to reconstruct the image. However, during the irradiation process, if the object is deformed due to the movement of the object or external disturbance, it will cause the distortion of the reconstruction result of the object, and even serious image distortion. According to the CXT model, the present invention allows multiple light sources to illuminate at the same time, thereby constructing a nonlinear snapshot system. Compared with the traditional reconstruction method of sequential illumination, the nonlinear CXT method of the multi-light source snapshots involved in the present invention greatly improves the robustness of the imaging of deformed objects, shortens the detection time, and improves the performance by optimizing the coded aperture. imaging accuracy.

The imaging effects of the present invention and the traditional method are described below in conjunction with the accompanying drawings:

Figure 3(a)-(d) are the four-layer slices of the initial object, (e)-(h) represent the tiny expansion inside the object, which is to simulate the tiny expansion of the lungs during inhalation, (i)-(l) It is the difference between the object with a slight expansion inside and the initial object, (m)-(p) represents a slight contraction inside the object, which is to simulate a slight contraction of the lungs when exhaling, and (q)-(t) is a slight contraction. The difference between the object and the original object. Figure 4(a)-(d) shows the reconstructed results when the object is slightly expanded and slightly contracted during the second and fourth illuminations, respectively, using the traditional sequential illumination method when using random coded apertures. The reconstructed average PSNR value is 10.21, and (e)-(h) represent the reconstructed results when the object is slightly expanded when there is only one irradiation, and the reconstructed average PSNR value is 12.84. Figures 5(a)-(d) show the results of object reconstruction using the method proposed in the present invention when a deformed object is irradiated with a random coded aperture, and the reconstructed average PSNR value is 39.95. In Fig. 6, (a) is the initial coding aperture, (b) is the coding aperture optimized by the method of the present invention, and (c) is the difference between the two. FIG. 7 shows the result of object reconstruction using the method proposed in the present invention with the optimized coded aperture, and the reconstructed average PSNR value is 41.58. Comparing FIG. 7 , FIG. 5 and FIG. 4 , it can be seen that the multi-light source snapshot method in the present invention can effectively improve the robustness of object reconstruction, and the coded aperture optimization can improve the reconstruction accuracy of the object.

Of course, the present invention can also have other various embodiments. Without departing from the spirit and essence of the present invention, those skilled in the art can of course make various corresponding changes and deformations according to the present invention, but these corresponding Changes and deformations should belong to the protection scope of the appended claims of the present invention.

Claims (5)

1. A method for compressive X-ray tomosynthesis of a multi-light-source snapshot is characterized by comprising the following steps:

s1: a correction framework for constructing a reconstruction framework based on the PRISM model is as follows:

wherein,

represents the optimized value of attenuation coefficient vector obtained after the X-ray light source scans the three-dimensional object, and the X tableShows the theoretical value, y, of the attenuation coefficient vector obtained after the X-ray source scans the three-dimensional object _k Represents the actual imaging model of each X-ray source on the detector after the coded aperture modulation in the K-th scanning, and K is 1,2, …, K represents the total number of scanning, f _k (x) Representing the simulated imaging model of each X-ray source on the detector after modulation by the coded aperture in the k-th scan, e _k Denotes y _k And f _k (x) Residual of (a), λ _* Weight coefficient representing low-rank regular term, | | · | | non-woven phosphor _* Regularizing kernel function, T, representing operation on matrix rank _T Representing the conversion of an auxiliary vector d into a tensor, T, of the same dimension as the tensor of a three-dimensional object _M (-) denotes the conversion of a tensor into a matrix of attenuation coefficients, μ _* And mu ₁ Indicates a set weight coefficient, λ ₁ Representing a weight coefficient of a sparse regular term, theta represents a sparse coefficient vector, v represents an auxiliary vector for measuring the distance between an auxiliary vector d and an attenuation coefficient vector x, and w represents an auxiliary vector for measuring the distance between an auxiliary vector c and a sparse coefficient vector theta;

s2: obtaining attenuation coefficient vector x and residual error e in the correction frame by adopting the corrected split Blackermann algorithm _k The iterative formulas of the auxiliary vector d, the auxiliary vector v, the auxiliary vector c and the auxiliary vector w are as follows:

v ⁿ⁺¹ ＝v ⁿ +x ⁿ⁺¹ -d ⁿ⁺¹

w ⁿ⁺¹ ＝w ⁿ +θ ⁿ⁺¹ -c ⁿ⁺¹

θ ⁿ⁺¹ ＝ψ ^T x ⁿ⁺¹

wherein n represents the number of iterations, ψ represents a sparse coefficient, and T represents transposition;

s3: setting the initial value of each iteration formula in the step S2 to be 0, and continuously repeating the iteration formula in the step S2, wherein the attenuation coefficient vector x obtained by the current iteration is judged after each iteration is finished ⁿ⁺¹ Whether a set requirement is met or not, wherein the set requirement comprises that the difference value between attenuation coefficient vectors obtained by two adjacent iterations is smaller than a set precision or the iteration frequency reaches a set upper limit, if so, the attenuation coefficient vector x obtained by the iteration is obtained ⁿ⁺¹ The optimal solution is used as an image of compressed X-ray tomosynthesis; if not, continuing the iteration until the attenuation coefficient vector x ⁿ⁺¹ The setting requirements are met.

2. The method of claim 1, wherein the construction method of the PRISM model-based reconstruction framework is as follows:

s11: constructing a non-limiting optimization function for reconstructing the attenuation coefficient inside the three-dimensional object from the K scans is as follows:

s12: unfolding tensor of the three-dimensional object into a two-dimensional matrix, and then carrying out low-rank limiting operation on the two-dimensional matrix to obtain a regular term R (x):

R(x)＝λ _* ||T _M (χ)|| _* +λ ₁ ||Τ _S (x)|| ₁ ＝λ _* ||X|| _* +λ ₁ ||θ ₁

wherein, T _M (X) watchTransforming tensor χ into attenuation coefficient matrix X, T _S (x) Representing the conversion of the attenuation coefficient vector x into corresponding sparse coefficients;

s13: constructing a reconstruction frame based on a PRISM model according to a non-limiting optimization function and a regular term:

3. the method of claim 1, wherein the simulation model f is a model of a multi-source snapshot compression X-ray tomosynthesis _k (x) The calculating method comprises the following steps:

wherein,

representing the set of indices of the X-ray source turned on at the k-th scan, b _p Indicates the vector obtained after the p-th light source scans the coded aperture below the p-th light source, indicates a dot product operation, H _p The system matrix corresponding to the p-th X-ray source is shown.

4. The method of claim 3, wherein the vector b is a compressed X-ray tomosynthesis of the multi-source snapshot _p The obtaining method specifically comprises the following steps:

s14: the cost function for optimizing the coded aperture is constructed as follows:

wherein,

m is the number of detector elements, N is the number of object voxels,

is the average measured intensity actual value of the ith detector element after coded aperture modulation, m ₁ The theoretical value of the measured intensity is averaged over all detector elements,

is the average actual value of the perception information of the jth voxel after the coded aperture modulation, m ₂ Average perceptual information content theoretical value, s, for each voxel _t For multiplexing of the i-th detector element, m ₃ Encoding the theoretical value of the number of aperture elements for the light transmission corresponding to each detector element;

s15: solving the cost function by adopting a DSB algorithm to obtain optimized b ₁ ,...,b _p 。

5. The method for compressive X-ray tomosynthesis of a multi-light-source snapshot of claim 1, further comprising the steps of:

s4: according to the optimal solution x ⁿ⁺¹ Constructing an objective function F:

s5: optimizing the objective function F by adopting a conjugate gradient method, and obtaining an optimized value x ^* As a compressed X-ray tomosynthesis image.

Images