Medical CT image metal artifact removing method based on X-ray multi-mode imaging

Abstract

The invention provides a medical CT (Computed Tomography) image metal artifact removing method based on X-ray multi-mode imaging, which comprises the following steps of: S1, obtaining projection data of an object through traditional X-ray energy loss measurement, and reconstructing a two-dimensional image of a detection plane through an ordered subset maximum expected value image reconstruction algorithm OSEM; s2, a position sensitive detector is used for collecting the incident position and the emergent position of the X-ray, so that scattering angle information is obtained, and a two-dimensional image is obtained through a closest point PoCA algorithm; and S3, fusing the two low-resolution images which are obtained in the step S1 and the step S2 and are staggered by half pixels in the linear array direction and the vertical linear array direction by using a sub-pixel super-resolution image algorithm to obtain a high-resolution image. According to the method, the problem of metal artifacts in the current CT imaging field can be solved to a certain extent, and the accuracy and reliability of diagnosis are improved.

Description

Technical Field

[0001] The present invention relates to the field of medical physics, and in particular to a method for removing metal artifacts from medical CT images based on X-ray multimodal imaging. Background Art

[0002] In clinical imaging diagnosis, metal objects such as dental fillings or silver clips in organs will produce obvious shadows and stripe artifacts during CT scanning. These artifacts not only seriously interfere with the imaging of tissues around the metal objects, but also have a great impact on the accurate judgment of the fault structure, resulting in a significant reduction in image quality and failure to meet the precise requirements of medical diagnosis. Therefore, in order to improve the accuracy and reliability of diagnosis, effective technical means must be used to suppress the negative impact of metal artifacts on CT diagnosis.

[0003] In view of this, the present invention proposes a method for removing metal artifacts from medical CT images based on X-ray multimodal imaging. Two images are obtained through X-ray scattering and transmission imaging methods, and then the two images are fused using a sub-pixel super-resolution reconstruction algorithm to obtain a high-resolution image, thereby reducing the impact of metal artifacts on CT diagnosis. Summary of the invention

[0004] The purpose of the present invention is to propose a method for removing metal artifacts from medical CT images based on X-ray multimodal imaging, which utilizes the scattering phenomenon between X-rays and irradiated objects and combines it with traditional CT imaging technology. It can solve the metal artifact problem currently encountered in the field of CT imaging to a certain extent and improve the accuracy and reliability of diagnosis.

[0005] To achieve the above object, the technical solution of the present invention is: a method for removing metal artifacts from medical CT images based on X-ray multimodal imaging, which specifically includes the following steps:

[0006] Step S1: obtaining projection data of the object by conventional X-ray energy loss measurement, and then reconstructing a two-dimensional image of the detection plane by using an ordered subset maximum expectation image reconstruction algorithm (OSEM);

[0007] Step S2: using a position sensitive detector to collect the positions of X-ray exit and incident, thereby obtaining scattering angle information, and obtaining a two-dimensional image through a closest point PoCA algorithm; wherein the two-dimensional image of the scattering imaging in step S2 and the two-dimensional image of the transmission imaging in step S1 are respectively offset by half a pixel in the linear array direction and the push-scan direction;

[0008] Step S3: The two low-resolution images obtained in step S1 and step S2, which are offset by half a pixel in the linear array direction and the vertical linear array direction, are fused using a sub-pixel super-resolution image algorithm to obtain a high-resolution image.

[0009] Preferably, the step S1 is specifically as follows: in each iteration process, the projection data is divided into L subsets, and each subset is used to correct all pixel values ​​of the reconstructed image to update the reconstructed image, and one iteration is completed when the reconstructed image update operation of all subsets is completed; the above iteration process is repeated until the change amount of the image update in each iteration is less than a set threshold, and the change amount is the root mean square error ∈ of the pixel value between the current iteration image and the previous iteration image; the iteration process is expressed as:

[0010]

[0011]

[0012]

[0013] In the formula, f (k) (i,j),f (k+1) (i,j),f l (i, j) represents the value of pixel (i, j) in the image of the kth iteration, the value of pixel (i, j) in the image of the k+1th iteration, and the value of pixel (i, j) in the image of the lth subset; (i, j) represents a pixel position in the image, i is the row index of the pixel, and j is the column index of the pixel; f(i, j) represents the pixel value at position (i, j) in the image; k represents the number of iterations, and the increase of k indicates that the algorithm has entered the next iteration. In each iteration, the algorithm will update the pixel value of image f and calculate the error of the image update to decide whether to continue the iteration; w ijmn is the reconstructed pixel f(i,j) at the mth detector unit P under projection angle n nm The weight factor, P l is the lth subset of the projection data partition, is the expected value of the reconstructed pixel at the mth detector unit under the projection angle n corresponding to each subset l; M is the total number of pixels in the image.

[0014] Preferably, the threshold is set to 10 -4 .

[0015] Preferably, the calculation of the scattering angle information is specifically as follows:

[0016] The target space straight line to be fitted Written in matrix equation form:

[0017]

[0018] Among them, x, y, and z are the three-dimensional position information of the particles arriving at each layer of detectors. 0 ,y 0 , a, b are the parameters to be fitted;

[0019] For the case of N measurement points, the corresponding least squares fitting parameters satisfy the following equation:

[0020]

[0021] In the equation, x 1 …x N ,y 2 …y N 、z 1 …z N It refers to the position coordinates of the same incident particle or the same emitted particle passing through N detectors in sequence. The coefficient matrix is ​​solved as follows:

[0022]

[0023] By obtaining various parameters through the above formula, the scattering projection angle of the cumulative scattering angle on the xz and yz planes can be obtained by the following formula:

[0024]

[0025] In the formula, θ x and θ y Respectively represent the scattering angle on the xz and yz planes. The subscripts out and in are used to distinguish the outgoing and incoming angles. out 、a in 、b out 、b in They represent the parameters of the straight line equations of the incident and exiting trajectories of the same particle respectively;

[0026] The scattering density λ is deduced based on the calculated scattering projection angle Lrad , the formula is as follows:

[0027]

[0028] Where L is the path length of the particle through the medium.

[0029] Preferably, the method of obtaining a two-dimensional image by the closest point PoCA algorithm specifically comprises: dividing the three-dimensional detection area into a plurality of cuboids, each cuboid being called a voxel; taking the intersection of the incident X-ray ray and the outgoing X-ray ray on the two-dimensional space xz or yz plane as the scattering point where scattering occurs; and obtaining the scattering density λ of the scattering point. Lrad When the scattering density λ LradThe voxels at the scattering points are assigned a value of zero, and the voxels that are not scattering points and are passed by the track are assigned a value of zero. The voxels that are not passed by the track are not assigned a value. When enough incident events are counted, almost all the voxels in the detection area are projected. By observing the xz or yz plane, a two-dimensional grayscale map of the scattering density distribution is obtained, that is, a scattering imaging picture is obtained.

[0030] Preferably, step S3 is specifically as follows: establishing a grid map of size 2m×2n for image fusion, assigning the average pixel value of two corresponding pixels in the two-dimensional image A and the two-dimensional image B to the corresponding pixel in the grid map, and obtaining a fused high-resolution image Q. The pixel value in the high-resolution image Q is calculated as follows:

[0031]

[0032] Among them, i and j are used to traverse the rows and columns of the two-dimensional image A and the two-dimensional image B, and i∈{1,2,…,m}, j∈{1,2,…,n}, q(,), a(,), b(,) represent the pixel values ​​of the pixel points corresponding to the corresponding rows and columns of the high-resolution image Q, the two-dimensional image A and the two-dimensional image B, respectively, and m and n are the number of rows and columns of the two low-resolution images A and B in sub-pixel imaging.

[0033] Compared with the prior art, the present invention has the following beneficial effects:

[0034] According to the method proposed in the present invention, two images of human tissue containing metal substances obtained after X-ray irradiation are reconstructed by sub-pixel full interpolation method to obtain a high-resolution image, so as to remove the metal artifact phenomenon generated in X-ray imaging. The present invention proposes to utilize the scattering phenomenon that occurs between X-rays and irradiated objects, and combine it with traditional CT imaging technology, which can solve the metal artifact problem encountered in the current CT imaging field to a certain extent and improve the accuracy and reliability of diagnosis. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] Figure 1 is a flow chart of the method of the present invention;

[0036] Figure 2 A grid diagram in which images A and B are alternately arranged in the present invention;

[0037] Figure 3 The figure is a schematic diagram of the placement of three layers of flat panel detectors in an embodiment of the present invention. DETAILED DESCRIPTION

[0038] The following is combined with Figure 1-3 , the technical solution of the present invention is specifically described.

[0039] like Figure 1As shown, the present invention proposes a method for removing metal artifacts from medical CT images based on X-ray multimodal imaging, which specifically includes the following steps:

[0040] Step S1: obtaining projection data of the object by conventional X-ray energy loss measurement, and then reconstructing a two-dimensional image of the detection plane by using an ordered subset maximum expectation image reconstruction algorithm (OSEM);

[0041] Step S2: utilizing the scattering physical process between X-rays and matter, using a position-sensitive detector to collect the exit and incident positions of the X-rays, thereby obtaining scattering angle information, and obtaining a two-dimensional image through a closest point PoCA algorithm; wherein the two-dimensional image of the scattering imaging in step S2 and the two-dimensional image of the transmission imaging in step S1 are respectively offset by half a pixel in the linear array direction and the push-scan direction;

[0042] Step S3: The two low-resolution images obtained in step S1 and step S2, which are offset by half a pixel in the linear array direction and the vertical linear array direction, are fused using a sub-pixel super-resolution image algorithm to obtain a high-resolution image, so as to reduce the influence of metal artifacts on the image.

[0043] The working principle of the method for removing metal artifacts from medical CT images of X-ray multimodal imaging according to the present invention is described in detail below.

[0044] First of all, CT imaging is to form an image by measuring the results of X-ray beams at different angles. In order to realize computed tomography reconstruction of human body imaging, it is necessary to determine the energy loss mechanism of X-rays:

[0045] I=I 0 e -μ·Δx

[0046] In the formula, I 0 is the incident X-ray intensity, I is the outgoing X-ray intensity, Δx is the thickness, and μ is the linear attenuation coefficient of the material. Generally speaking, μ is not a constant value, it changes with the energy of the X-ray, and different materials have different μ values.

[0047] After obtaining parallel projections of an object at different angles, traditional X-ray projection imaging uses a Fourier transform filtered back-projection algorithm to calculate the line integral of the density function for each projection ray along a certain projection direction, and then obtain the projection value on the ray. Calculate all projection values ​​in the projection direction, and then obtain the projection function in the projection direction. The Fourier slice theorem states that the one-dimensional Fourier transform function of the projection function of the density function in a certain direction is the value of the two-dimensional Fourier transform function of the original function on the plane along the same direction and on the straight line passing through the source point. The Fourier slice theorem points out the possibility of reconstructing images from projections. If enough projection function data is obtained at different angles and their Fourier transforms are performed, then the transformed data will fill the entire plane. Once all the values ​​of the frequency domain function are obtained, perform an inverse Fourier transform on it to obtain the original density function, which is the image we want to reconstruct.

[0048] The present invention will use an iterative reconstruction algorithm, namely the ordered subset maximum expectation reconstruction algorithm (OSEM). It is an improved maximum likelihood expectation method (MLEM). The OSEM method belongs to a type of "block" iterative method. In each iteration, the projection data is divided into L subsets. Generally, the similarity of the subsets is required to be as low as possible, so as to speed up the reconstruction of the image to approach the true value. In each subset, the values ​​of all pixels in the reconstructed image are corrected, and the reconstructed image is updated once. All subset operations are completed, which is called an iterative process. Repeat the above iterative process until the change in the image update in each iteration is less than the set threshold. The change is the root mean square error ∈ of the pixel value between the current iteration image and the previous iteration image. We set the threshold to 10 -4 In the MLEM method, the reconstructed image is updated only once after all the projection data are corrected in one iteration. In the OSEM method, the reconstructed image is updated L times in one iteration, which accelerates the convergence of the iteration process. This process can be described by the following formula:

[0049]

[0050]

[0051] Here, imaging is considered as an algebraic problem Wf=P, where f is the vector of the reconstructed image, P is the signal received by the detector, and W is the contribution factor of each pixel of the reconstructed image to each ray. (k) (i,j),f (k+1) (i,j),f l(i, j) represents the value of pixel (i, j) in the image of the kth iteration, the value of pixel (i, j) in the image of the k+1th iteration, and the value of pixel (i, j) in the image of the lth subset; (i, j) represents a pixel position in the image, i is the row index of the pixel, and j is the column index of the pixel; f(i, j) represents the pixel value at position (i, j) in the image; k represents the number of iterations, and the increase of k indicates that the algorithm has entered the next iteration. In each iteration, the algorithm will update the pixel value of image f and calculate the error of the image update to decide whether to continue the iteration; w ijmn is the reconstructed pixel f(i,j) at the mth detector unit P under projection angle n nm The weight factor is M, which is the total number of pixels in the image. l is the lth subset of the projection data partition, is the expected value of the detector m-th unit at the projection angle n corresponding to each subset l of the reconstructed pixel:

[0052]

[0053] The OSEM method has better reconstructed image quality and shorter calculation time, and conventional X-ray transmission imaging images are obtained after its algorithm.

[0054] Secondly, the present invention proposes to utilize the scattering process in the process of X-ray imaging. In the simulation process of Monte Carlo software geant4, we generally set three layers of position-sensitive detectors to receive the position coordinates of the X-ray outgoing and incident objects. Among them, the parameters of the three layers of flat-panel detectors are consistent, and they are placed in parallel above and below the detection object, with a spacing of 10 cm between each other. The external track data before and after the X-ray passes through the inspected object is the basis of X-ray scattering imaging. In the imaging system, the incident track and the outgoing track of each particle are respectively given by three three-dimensional coordinate points. In this embodiment, when the same particle passes through the three layers of position-sensitive detectors, the position (xi, yi, zi), i∈{1,6} can be obtained. These three points are often not collinear in space due to the errors of the detectors and background noise. It is necessary to use a certain method to fit and reconstruct the external track and calculate the corresponding scattering angle.

[0055] The present invention adopts the least square method to perform track reconstruction, and the method is described as follows:

[0056] The target space straight line to be fitted Written in matrix equation form:

[0057]

[0058] Among them, x, y, and z are the three-dimensional position information of the particles arriving at each layer of detectors.0 ,y 0 , a, and b are the parameters to be fitted.

[0059] For the case of N measurement points, the corresponding least squares fitting parameters satisfy the following equation:

[0060]

[0061] In the equation, x 1 …x N ,y 2 …y N 、z 1 …z N It refers to the position coordinates of the same incident particle or the same emitted particle passing through N detectors in sequence. In this embodiment, since the experimental simulation uses three layers of detectors on the upper and lower layers to collect particle (muon) position information, N is 3. Its coefficient matrix can be solved as follows:

[0062]

[0063] By obtaining various parameters through the above method, the scattering projection angle of the cumulative scattering angle on the xz and yz planes can be obtained by the following formula:

[0064]

[0065] In the formula, θ x and θ y Respectively represent the scattering angle on the xz and yz planes. The subscripts out and in are used to distinguish the outgoing and incoming angles. out 、a in 、b out 、b in They represent the parameters of the straight line equations of the incident and exiting trajectories of the same particle respectively;

[0066] The scattering density λ is deduced based on the calculated scattering projection angle Lrad , the formula is as follows:

[0067]

[0068] Where L is the path length of the particle (muon) through the medium.

[0069] The above method can be used to obtain the incident and outgoing track equations of each particle, and calculate the corresponding cumulative scattering angle.

[0070] Furthermore, the scattering angle of each X-ray is obtained by the above method, and then we need to obtain the scattering imaging image through the corresponding imaging algorithm. The present invention selects the closest point of approach (PoCA) algorithm commonly used in particle (muon) imaging. The algorithm uses an imaging framework similar to TGS to divide the three-dimensional detection area into small cuboids, which are called voxels, and each voxel is applicable to the above-mentioned imaging principle. The intersection of the X-ray incident ray and the outgoing ray on the xz or yz plane in two-dimensional space is used as the scattering point where scattering occurs. After obtaining the scattering density λ of the point Lrad When the incident events are counted, the voxels at the scattering point are assigned with the projection. Other voxels that the track passes through also have an impact on the scattering, but they are not scattering points, so they are assigned with a value of zero. Voxels that are not passed through by the track are not assigned with a value. When enough incident events are counted, almost all voxels in the detection area are projected. By observing the xz or yz plane, a two-dimensional grayscale image of the scattering density distribution can be obtained, that is, a scattering imaging picture can be obtained.

[0071] Finally, through the sub-pixel super-resolution reconstruction method, the two transmission and scattering imaging images mentioned above are fused with two low-resolution images that are offset by half a pixel in the linear array direction and the push-scan direction to obtain a high-resolution image.

[0072] Two low-resolution images with half a pixel offset in the linear array direction and the vertical linear array direction are represented by A and B. The size of images A and B is m×n, and a grid of size 2m×2n is established. The data images obtained by sampling A and B are as follows Figure 2 As shown, they are arranged in a grid in an interlaced manner. At this time, there are 2×(m×n) known data in the grid. That is, the remaining 2×m×n data of the A and B image data are unknown. It is necessary to estimate the value based on the correlation of the location point field data and reconstruct a high-resolution image.

[0073] The present invention adopts the full interpolation method in the interpolation method for value estimation. The interpolation process is to determine the numerical operation of the function between multiple adjacent sampling value points. Its essence is to express a certain continuous function in a discrete form of the linear operation of the basis function, which is an approximate theory in mathematics.

[0074] The traditional sub-pixel reconstruction method is to directly assign the pixel values ​​of two m×n low-resolution images A and B to the grid of the reconstructed high-resolution image, and the remaining pixel values ​​in the high-resolution image are calculated using different interpolation methods based on the neighborhood pixel values. If the resolution of the detector is a×a, then a(i,j) and b(i,j) are grayscale values ​​of an area of ​​a×a. If the pixel values ​​in A and B are directly assigned to the high-resolution image, the grayscale values ​​of an area of ​​a×a are used to directly replace the grayscale values ​​of an area of ​​a / 2×a / 2, and the result is not accurate enough.

[0075] The following are the specific steps of the full interpolation method: Here, the reconstructed high-resolution image is named Q. The area where the pixel q(2,2) of Q is located is the intersection of the areas where a(1,1) and b(1,1) are located in the low-resolution graphics A and B. Similarly, the area where each pixel in Q is located is the intersection of the areas where the two points in A and B are located. The average value of the corresponding two pixels in A and B is assigned to the corresponding pixel in Q. The calculation result will be closer to the true value than the traditional sub-pixel interpolation reconstruction method.

[0076] The calculation formula of the full interpolation method is:

[0077]

[0078] Among them, i∈{1,2,…,m}, j∈{1,2,…,n}, m and n are the number of rows and columns of the two low-resolution images A and B in sub-pixel imaging.

[0079] The above are preferred embodiments of the present invention. Any changes made according to the technical solution of the present invention, as long as the resulting functions do not exceed the scope of the technical solution of the present invention, belong to the protection scope of the present invention.

Claims

1. A method for removing metal artifacts from medical CT images based on X-ray multimodal imaging. It is characterized in that The specific steps include: Step S1: obtaining projection data of the object by conventional X-ray energy loss measurement, and then reconstructing a two-dimensional image of the detection plane by using an ordered subset maximum expectation image reconstruction algorithm (OSEM); Step S2: using a position sensitive detector to collect the positions of X-ray exit and incident, thereby obtaining scattering angle information, and obtaining a two-dimensional image through a closest point PoCA algorithm; wherein the two-dimensional image of the scattering imaging in step S2 and the two-dimensional image of the transmission imaging in step S1 are respectively offset by half a pixel in the linear array direction and the push-scan direction; Step S3: The two low-resolution images obtained in step S1 and step S2, which are offset by half a pixel in the linear array direction and the vertical linear array direction, are fused using a sub-pixel super-resolution image algorithm to obtain a high-resolution image.

2. The method for removing metal artifacts from medical CT images based on X-ray multimodal imaging according to claim 1, It is characterized in that The step S1 is specifically as follows: in each iteration process, the projection data is divided into L subsets, and each subset is used to correct all pixel values ​​of the reconstructed image to update the reconstructed image. The completion of the reconstructed image update operation of all subsets is considered as the completion of one iteration; the above iteration process is repeated until the change amount of the image update in each iteration is less than the set threshold value, and the change amount is the root mean square error ∈ of the pixel value between the current iteration image and the previous iteration image; the iteration process is expressed as: In the formula, f (k) (i, j), f (k+1) (i, j), f l (i, j) represents the value of the pixel (i, j) in the image of the kth iteration, the value of the pixel (i, j) in the image of the k+1th iteration, and the value of the pixel (i, j) in the image of the lth subset; (i, j) represents a pixel position in the image, i is the row index of the pixel, and j is the column index of the pixel; f(i, j) represents the pixel value at position (i, j) in the image; k represents the number of iterations, and the increase of k indicates that the algorithm has entered the next iteration. In each iteration, the algorithm will update the pixel value of image f and calculate the error of the image update to decide whether to continue the iteration; w ijmn is the reconstructed pixel f(i, j) at the mth detector unit P under projection angle n nm The weight factor, P l is the lth subset of the projection data partition, is the expected value of the reconstructed pixel at the mth detector unit under the projection angle n corresponding to each subset l; M is the total number of pixels in the image.

3. The method for removing metal artifacts from medical CT images based on X-ray multimodal imaging according to claim 2, It is characterized in that The threshold is set to 10 -4 .

4. The method for removing metal artifacts from medical CT images based on X-ray multimodal imaging according to claim 1, It is characterized in that The calculation of the scattering angle information is specifically as follows: The target space straight line to be fitted Written in matrix equation form: Among them, x, y, z are the three-dimensional position information of the particles arriving at each layer of detectors, and x0, y0, a, and b are the parameters to be fitted; For the case of N measurement points, the corresponding least squares fitting parameters satisfy the following equation: In the equation, x 1 …x N ,y 2 …y N 、z 1 …z N It refers to the position coordinates of the same incident particle or the same emitted particle passing through N detectors in sequence. The coefficient matrix is ​​solved as follows: By obtaining various parameters through the above formula, the scattering projection angle of the cumulative scattering angle on the xz and yz planes can be obtained by the following formula: In the formula, θ x and θ y Respectively represent the scattering angle on the xz and yz planes. The subscripts out and in are used to distinguish the outgoing and incoming angles. out 、a in 、b out 、b in They represent the parameters of the straight line equations of the incident and exiting trajectories of the same particle respectively; The scattering density λ is deduced based on the calculated scattering projection angle Lrad , the formula is as follows: Where L is the path length of the particle through the medium.

5. The method for removing metal artifacts from medical CT images based on X-ray multimodal imaging according to claim 4, It is characterized in that The method of obtaining a two-dimensional image by the closest point PoCA algorithm specifically comprises the following steps: dividing the three-dimensional detection area into a plurality of cuboids, each of which is called a voxel; taking the intersection of the incident X-ray ray and the outgoing X-ray ray on the two-dimensional space xz or yz plane as the scattering point where scattering occurs; and obtaining the scattering density λ of the scattering point. Lrad When the scattering density λ Lrad The voxels at the scattering points are assigned a value of zero, and the voxels that are not scattering points and are passed by the track are assigned a value of zero. The voxels that are not passed by the track are not assigned a value. When enough incident events are counted, almost all the voxels in the detection area are projected. By observing the xz or yz plane, a two-dimensional grayscale map of the scattering density distribution is obtained, that is, a scattering imaging picture is obtained.

6. The method for removing metal artifacts from medical CT images based on X-ray multimodal imaging according to claim 1, It is characterized in that The step S3 is specifically as follows: a grid map of size 2m×2n is established for image fusion, and the average pixel value of two corresponding pixels in the two-dimensional image A and the two-dimensional image B is assigned to the corresponding pixel in the grid map to obtain a fused high-resolution image Q. The pixel value in the high-resolution image Q is calculated as follows: Among them, i and j are used to traverse the rows and columns of the two-dimensional image A and the two-dimensional image B, and i∈{1, 2, …, m}, j∈{1, 2, …, n}, q(,), a(,), b(,) represent the pixel values ​​of the pixel points corresponding to the corresponding rows and columns of the high-resolution image Q, the two-dimensional image A and the two-dimensional image B, respectively, and m and n are the number of rows and columns of the two low-resolution images A and B in sub-pixel imaging.

Images