Image processing method and device, electronic equipment and computer readable storage medium

Abstract

The application provides an image processing method and device, electronic equipment and a computer readable storage medium. The image processing method comprises: acquiring voxel values of a plurality of voxels in a to-be-corrected attenuation image; acquiring a convolution kernel and normalizing the convolution kernel; performing deconvolution operation on the voxel values of the attenuation image by using the normalized convolution kernel to obtain voxel values of a corrected attenuation image. According to the embodiment of the application, the attenuation image is corrected, so that the attenuation image produces a certain "internal shrinkage" deformation while producing a "blurring" effect; the attenuation correction coefficient is calculated by using the corrected attenuation image, so that the attenuation correction coefficient is more accurate, and thus the image artifacts are eliminated.

Description

[0001] The present application relates to the field of data processing, in particular to an image processing method and device, electronic equipment and computer readable storage medium. BACKGROUND

[0002] Attenuation correction is a very important step in PET imaging. Attenuation correction needs to obtain an attenuation image first, in which each voxel value represents the linear attenuation coefficient at that voxel position. Then the integral value of linear attenuation coefficient along the line connecting each pair of detectors is calculated to obtain the attenuation correction coefficient, as shown in equation (1).

[0003]

[0004] where AC(i) represents the attenuation correction coefficient of the ith pair of detectors, A i and B i are the coordinates of the two detectors, considering that the detector usually has a certain volume, the coordinates are usually taken as the center point of the detector surface, and μ(x) represents the distribution of the linear attenuation coefficient in space with the coordinate x.

[0005] Equation (1) can also be expressed in a discrete form, as shown in equation (2).

[0006]

[0007] where μ j represents the jth voxel value of the attenuation image, c ij represents the length of the line connecting the ith pair of detectors A i B i through the jth voxel of the attenuation image.

[0008] The attenuation correction coefficient AC(i) is used to correct the scan data before image reconstruction, and the correction process can be represented by equation (3).

[0009]

[0010] where Y(i) is the number of coincidence events detected by the ith pair of detectors, i.e. the number of coincidence events without correction, r(i) is the number of random events of the ith pair of detectors, s(i) is the number of scatter events of the ith pair of detectors, N(i) is the normalization correction coefficient of the ith pair of detectors, Y AC (i) is the number of coincidence events after attenuation correction.

[0011] At present, the attenuation correction coefficient is mainly used in image reconstruction, for example, the classic reconstruction algorithm maximum likelihood expectation maximization (MLEM) algorithm:

[0012]

[0013] wherein, is the value of the jth voxel of the reconstructed image at the nth iteration, P ij is the probability that a gamma photon emitted by the jth voxel is detected by the ith pair of detectors, and is a constant known in advance.

[0014] The prior art usually obtains the attenuation image from a CT image or MR image registered with the PET image. There are many techniques that can obtain very accurate attenuation images, i.e. the obtained attenuation image can accurately reflect the linear attenuation coefficient at each position in space. However, the prior art has a calculation error when calculating the integral value of the linear attenuation coefficient with respect to the distance on the line connecting the pairs of detectors.

[0015] Figure 1a The line connecting the positions of a pair of detectors. Figure 1a The middle of the figure is the accurate spatial distribution of the linear attenuation coefficient, with dark color indicating high density and light color indicating low density. According to equation (2), when calculating the line integral, the straight line will pass through a longer distance in the high-density region. However, the flight path of the gamma photon detected by a pair of detectors can be different from Figure 1a . Figure 1b As shown in Figure 1c or as shown in Figure 1d . Figure 1b This is because the detector has a certain volume, Figure 1c because the gamma photon can penetrate the adjacent detector, Figure 1d because the gamma photon can first enter the adjacent detector and then be scattered back. Among them, Figure 1c the case shown in Figure 1c is the most common. Since the path of the gamma photon in does not pass through the high-density region, an error is caused in the calculation of the attenuation correction coefficient.

[0016] Figure 2a There are two main methods in the prior art to solve the problem of error in calculating the attenuation correction coefficient. As shown in Figure 2a , one method is to set an average depth when positioning the line connecting the detectors when calculating the attenuation correction coefficient, so that the path for calculating the attenuation correction coefficient will be moved from the dashed line to the solid line in the figure. Figure 1b The method shown in Figure 2b can take into account the effect of Figure 2b to some extent, but cannot comprehensively consider all effects. As shown in Figure 2bThe method shown can simulate paths of multiple gamma photons, but this method cannot set optimal numbers of straight lines, positions of straight lines and weights of each straight line. Currently, numbers and positions of straight lines are set according to experience, and weights of each straight line are usually the same.

[0017] Therefore, the attenuation correction coefficients calculated by the above two methods are not accurate enough, which further causes artifacts in the reconstructed image. SUMMARY

[0018] The present application provides an image processing method and device, electronic equipment and computer readable storage medium, to solve the problem of calculation error when calculating the attenuation correction coefficient.

[0019] According to an aspect of the present application, an image processing method is provided, which includes: obtaining voxel values of a plurality of voxels in an attenuation image to be corrected; obtaining a convolution kernel and normalizing the convolution kernel; performing deconvolution operation on the voxel values of the attenuation image by using the normalized convolution kernel to obtain voxel values of a corrected attenuation image.

[0020] According to some embodiments of the present application, the convolution kernel is calculated by formula (5):

[0021]

[0022] wherein (jx, jy, jz) is a coordinate corresponding to the jth voxel, (x, y, z) takes any value in the imaging field of view, κ j (x, y, z) is a value of the convolution kernel of the jth voxel at the spatial coordinate (x, y, z), is a 2D convolution kernel of the x-y plane, κ a (z) is a 1D convolution kernel in the axial direction.

[0023] According to some embodiments of the present application, the 2D convolution kernel of the x-y plane has rotational symmetry around the origin, and the 2D convolution kernel of the x-y plane satisfies formula (6):

[0024]

[0025] wherein, x' = x cos θ + y sin θ, y' = -x sin θ + y cos θ, wherein, is a 2D convolution kernel of the voxel in the x-y plane (r, 0) coordinate, (jx, jy, jz) is a coordinate corresponding to the voxel j, and (x, y, z) is any value in the imaging field of view.

[0026] According to some embodiments of the present application, the 2D kernel of the x-y plane in (r, 0) coordinates satisfies the PiG model, which satisfies equation (9):

[0027]

[0028] where A, B are the amplitudes of the Gaussian function, σ x0 ,σ y0 are the spreads of the Gaussian function in x and y directions on the side close to the origin, σ x1 ,σ y1 are the spreads of the Gaussian function in x and y directions on the side far from the origin, and x0 is the offset of the 2D kernel in the x-axis.

[0029] According to some embodiments of the present application, A = B, σ x0 = σ x1 = σ x ,σ y0 = σ y1 = σ y , the 2D kernel of the x-y plane in (r, 0) coordinates satisfies the SiG model, which satisfies equation (10):

[0030]

[0031] According to some embodiments of the present application, the A, B, σ x0 ,σ y0 ,σ x1 ,σ y1 , and x0 are obtained by scanning a line source at multiple different positions and reconstructing images; accumulating the reconstructed images to an x-y plane; and fitting the PiG or SiG model to the image of the x-y plane to obtain the A, B, σ x0 ,σ y0 ,σ x1 ,σ y1 , and x0.

[0032] According to some embodiments of the present application, the diameter of the line source is less than or equal to the edge length of the voxels of the attenuation image; the length of the line source is greater than or equal to 1 mm, preferably the length is greater than 10 mm; and the line source is parallel to the z-axis.

[0033] According to some embodiments of the present application, the A, B, σ x0 ,σ y0 ,σ x1 ,σ y1 , and x0 of the adjacent positions are calculated by using a 1D linear interpolation algorithm. x0 y0 x1 ,σ​​y1 x0.

[0034] According to some embodiments of the present application, the axial 1D kernel is obtained by scanning a plurality of face sources at different positions and reconstructing images; accumulating the reconstructed images to the z axis to obtain reconstructed image data; and averaging the reconstructed image data at different positions to obtain the axial 1D kernel.

[0035] According to some embodiments of the present application, the thickness of the face source is less than or equal to the edge length of the voxel of the attenuation image; the area of the face source is greater than or equal to 1 mm 2 , preferably greater than 10 mm 2 ; and the face source is perpendicular to the z axis.

[0036] According to some embodiments of the present application, the kernel of each voxel is normalized by formula (12):

[0037]

[0038] wherein κ′ j (x, y, z) is the normalized value of the kernel of the jth voxel at the position (x, y, z), κ j (x, y, z) is the kernel of the jth voxel.

[0039] According to some embodiments of the present application, the voxel value of the attenuation image is deconvoluted by using the normalized kernel by formula (13) to obtain the voxel value of the corrected attenuation image:

[0040] μ corr (x, y, z) = ∑ j κ′ j (x, y, z) μ j Formula (13)

[0041] wherein μ corr (x, y, z) is the voxel value of the corrected attenuation image at the position (x, y, z), and μ j is the voxel value of the uncorrected attenuation image at the jth voxel.

[0042] According to some embodiments of the present application, the image processing method further comprises: calculating an attenuation correction coefficient by using the voxel value of the corrected attenuation image; and performing image reconstruction by using the attenuation correction coefficient.

[0043] According to some embodiments of the present application, the attenuation image is obtained by a CT image or an MR image.

[0044] According to an aspect of the present application, an image processing apparatus is provided, comprising: a voxel value obtaining unit configured to obtain voxel values of a plurality of voxels in an attenuation image to be corrected; a convolution kernel calculating unit configured to obtain a convolution kernel and normalize the convolution kernel; and an attenuation image correcting unit configured to perform deconvolution operation on the voxel values of the attenuation image by using the normalized convolution kernel to obtain a corrected attenuation image.

[0045] According to an aspect of the present application, an electronic device is provided, comprising: one or more processors; a storage device configured to store a computer program; and when the computer program is executed by the one or more processors, the one or more processors are caused to implement the method according to any one of the above.

[0046] According to an aspect of the present application, a computer readable storage medium is provided, having stored thereon program instructions, which when executed implement the method according to any one of the above.

[0047] According to the example embodiments of the present application, by correcting the attenuation image, the attenuation image produces a "blurring" effect while also producing a certain "shrinkage" deformation; by using the corrected attenuation image to calculate the attenuation correction coefficient, the attenuation correction coefficient is more accurate, thereby eliminating image artifacts. Moreover, the entire correction process does not need to consider whether the detector has a certain volume, transmission and scattering of gamma photons between detectors, and the like. BRIEF DESCRIPTION OF DRAWINGS

[0048] In order to more clearly illustrate the technical solutions in the embodiments of the present application, the drawings needed in the embodiment description will be briefly introduced as follows.

[0049] Figure 1a A flight path of a gamma photon detected by a pair of detectors is shown.

[0050] Figure 1b Another flight path of a gamma photon detected by a pair of detectors is shown.

[0051] Figure 1c Another flight path of a gamma photon detected by a pair of detectors is shown.

[0052] Figure 1d Another flight path of a gamma photon detected by a pair of detectors is shown.

[0053] Figure 2a A method for calculating a linear attenuation coefficient of a pair of detectors by setting an average depth is shown.

[0054] Figure 2bA method of calculating a pair of linear attenuation coefficients of a detector is shown by integrating linear attenuation coefficients over a plurality of lines.

[0055] Figure 3 A flow chart of an image processing method according to an example embodiment of the application is shown.

[0056] Figure 4 A schematic diagram of an image processing according to an example embodiment of the application is shown.

[0057] Figure 5 A schematic diagram of a per voxel kernel space decomposition process according to an example embodiment of the application is shown.

[0058] Figure 6a A schematic diagram of a calculation of an x-y plane 2D kernel according to an example embodiment of the application is shown.

[0059] Figure 6b A schematic diagram of a calculation of an axial kernel according to an example embodiment of the application is shown.

[0060] Figure 7a A cross-sectional material profile of a phantom 1 with uniform activity distribution.

[0061] Figure 7b A coronal material profile of a phantom 2 with uniform activity distribution.

[0062] Figure 8a A reconstructed image of phantom 1 using a normal linear attenuation coefficient integration method is shown.

[0063] Figure 8b A reconstructed image of phantom 1 using a Figure 2a method is shown.

[0064] Figure 8c A reconstructed image of phantom 1 using a Figure 2b method is shown.

[0065] Figure 8d A reconstructed image of phantom 1 using a Figure 3 method is shown.

[0066] Figure 9a A reconstructed image of phantom 2 using a normal linear attenuation coefficient integration method is shown.

[0067] Figure 9b A reconstructed image of phantom 2 using a Figure 2a method is shown.

[0068] Figure 9c A reconstructed image of phantom 2 using a Figure 2b method is shown.

[0069] Figure 9d It is shown that the prosthesis 2 adopts Figure 3 the reconstructed image of the method shown.

[0070] Figure 10 It is shown that the image processing device block diagram according to an embodiment of the application.

[0071] Figure 11 It is shown that the image processing device block diagram according to another embodiment of the application. DETAILED DESCRIPTION

[0072] Example embodiments now will be described more fully hereinafter with reference to the accompanying drawings. Example embodiments, however, can be implemented in many different forms and should not be construed as limited to the embodiments set forth herein; rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of example embodiments to those skilled in the art. Like reference numerals refer to like elements throughout the several views.

[0073] The described features, structures, or characteristics can be combined in any suitable manner in one or more embodiments. In the following description, numerous specific details are provided to give a thorough understanding of embodiments of the disclosure. One skilled in the relevant art will recognize, however, that the

[0074] The flow charts shown in the figures are examples only and are not necessarily implemented in the order as shown. For example, one or more operations / steps can be eliminated, or one or more operations / steps can be combined or partially combined, and the order of operations / steps can be changed, depending on the actual implementation.

[0075] The terms "first", "second", and the like, in the description and in the claims of the present specification, as well as above-mentioned drawings, are used to distinguish between similar objects and are not necessarily used to describe a specific sequential or chronological order. Moreover, the terms "comprises", "comprising", "includes", "including", "has", "having" or the like are intended to encompass the inclusion of one or more steps or elements without limitation. For example, a process, method, system, product, or apparatus that comprises a list of steps or elements is not necessarily limited to those steps or elements which are recited.

[0076] Hereinafter, specific embodiments according to the present application will be described in detail with reference to the accompanying drawings.

[0077] Figure 3A flow chart of an image processing method according to an embodiment of the present application is shown. The following refers to Figure 3 An image processing method according to an embodiment of the present application is described in detail.

[0078] According to some example embodiments of the present application, Figure 3 The attenuation image in the image processing method shown is from a CT image or an MR image registered with the PET image.

[0079] Figure 4 A schematic diagram of correcting an attenuation image according to the image processing method shown is shown. Figure 3 A schematic diagram of correcting an attenuation image according to the image processing method shown is shown. Figure 3 The image processing method shown is to correct the accurate attenuation image, and then reconstruct the image using the corrected image according to the conventional integral algorithm for calculating the linear attenuation coefficient, such as formulas (1)-(4).

[0080] In step S101, the voxel value of each voxel in the attenuation image to be corrected is obtained.

[0081] The voxel value of each voxel in the attenuation image represents the linear attenuation coefficient at the voxel position, and according to some embodiments, the voxel value is obtained by a CT image or an MR image.

[0082] In step S103, the convolution kernel is obtained, and the convolution kernel is normalized.

[0083] According to some example embodiments, the convolution kernel is obtained by setting an additional line source or surface source for scanning experiments when the attenuation image is obtained, and the convolution kernel of each voxel is a 3D convolution kernel, which has spatial decomposability, i.e., can be decoupled into the product of an x-y plane 2D convolution kernel and an axial convolution kernel, as shown in formula (5).

[0084]

[0085] where (jx, jy, jz) is the coordinate corresponding to the jth voxel, (x, y, z) can be any value within the imaging field of view (FOV), and each voxel corresponds to a 3D convolution kernel. For the jth voxel, its own coordinate is (jx, jy, jz), and κ j (x, y, z) is the value of the 3D convolution kernel of the jth voxel at the spatial coordinate (x, y, z), is the 2D convolution kernel in the x-y plane, and κ a (z) is the 1D convolution kernel in the axial direction, i.e., the z-axis direction; and κ a (z-jz) uses z-jz to reflect the axial 1D convolution kernel invariance with the axial translation of the voxel j, and "jz" offsets the offset amount of the voxel j in the z-axis.

[0086] Figure 5 FIG. 1 shows a schematic diagram of a per-voxel kernel space decomposition process according to an example embodiment of the present application.

[0087] According to some embodiments, when (x, y, z) is far away from (jx, jy, jz), such as 8mm-15mm, the per-voxel kernel K j (x, y, z) = 0, so that the computation can be saved, instead of computing for all (x, y, z).

[0088] According to some embodiments, for the 2D kernel in x-y plane, which depends on the x-axis and y-axis coordinates jxand jyof the jth voxel.

[0089] According to some embodiments, K a (z) is the 1D kernel in z-axis, which is not affected by the z-axis coordinate of the jth voxel.

[0090] According to some example embodiments of the present application, the 2D kernel in x-y plane K has rotational symmetry around the origin, as shown in equation (6).

[0091]

[0092] wherein, is the 2D kernel of the voxel in the (jx, jy) coordinates of the x-y plane, is the 2D kernel of the voxel in the (r, 0) coordinates of the x-y plane, (jx, jy, jz) is the coordinate corresponding to the jth voxel, and (x, y, z) can take any value within the FOV.

[0093] According to some embodiments, depends on the x-axis coordinate of the jth voxel and the radius r of the jth voxel from the origin, wherein, x' and y' can be calculated by equations (7) and (8).

[0094] x' = x cos θ + y sin θ, y' = -x sin θ + y cos θ equation (7)

[0095]

[0096] wherein, θ is an intermediate variable, which can be determined by equation (8).

[0097] According to some example embodiments of the present application, satisfy an asymmetric 2D piece-wise Gaussian function, i.e., a PiG (piece-wise Gaussian) model, as shown in equation (9).

[0098]

[0099] where A, B are the amplitudes of the Gaussian functions, σ x0 ,σ y0 are the spreads of the Gaussian functions in the x and y directions on the side close to the origin, σ x1 ,σ y1 are the spreads of the Gaussian functions in the x and y directions on the side far from the origin, and x0 is the shift of the 2D kernel in the x axis.

[0100] According to some embodiments, σ x0 ,σ y0 ,σ x1 ,σ y1 represents the "blurring" effect on the attenuation image in the image processing procedure. x0 is the shift of the kernel in the x axis, which represents the "shrinkage" effect on the attenuation image in the image processing procedure.

[0101] When the kernel is narrow, for example, the kernel is non-zero only within 1-2 mm adjacent to voxel j, the PiG model is greatly affected by noise. According to some embodiments, when the kernel is narrow, a degenerated model SiG (single Gaussian) model of the PiG model is used. At this time, the parameters in equation (9) satisfy A = B, σ x0 = σ x1 = σ x ,σ y0 = σ y1 = σ y , and the SiG model as shown in equation (10) is obtained.

[0102]

[0103] Figure 6a A schematic diagram of calculating the x-y plane kernel according to an embodiment of the present application is shown in FIG. 6a. The parameters of the PiG model and / or the SiG model are obtained by scanning a plurality of line sources at different positions and reconstructing images, and then accumulating the reconstructed images to the x-y plane, and finally fitting A, B, σ x0 ,σ y0 ,σ x1 ,σ y1 , and x0.

[0104] According to some embodiments of the present application, the fitted parameters A, B, σ x0 ,σ y0,σ x1 ,σ y1 , the relationship between x0 and r, and observe whether the parameters change smoothly with r. For locations with large jitter, the data can be fitted again using the SiG model and the fitting result can be used to replace the fitting result of the PiG model.

[0105] According to some embodiments, the diameter of the line source is less than or equal to the side length of the voxel of the attenuation image, the length of the line source is greater than or equal to 1 mm, the line source passes through the x-axis and is parallel to the z-axis, and the coordinate of the line source passing through the x-axis corresponds to The r in.

[0106] According to some embodiments, the line source scan data at each position is reconstructed using the MLEM algorithm.

[0107] According to some embodiments, the line source is parallel to the z-axis and placed at multiple positions on the x-axis for scanning. The parameters of the PiG model and / or SiG model that are not at the aforementioned position are calculated based on the known A, B, σ of the adjacent positions. x0 ,σ y0 ,σ x1 ,σ y1 ,x0 is calculated using 1D linear interpolation algorithm.

[0108] Formula (11) shows a method for calculating parameter A using a 1D linear interpolation algorithm.

[0109]

[0110] Among them, A(r) is the value of parameter A when the radius is r, r k The x-axis coordinate where the kth line source is placed during scanning, k = 1, 2, ..., N (N is the total number of line source scanning positions). The calculation of the remaining parameters is similar to that of parameter A and will not be repeated here.

[0111] According to some embodiments, the rotational symmetry of the transaxial 2D convolution kernel around the origin in the xy plane is used to calculate other voxels that are not on the x-axis.

[0112] Figure 6b FIG. 4 shows a schematic diagram of calculating an axial convolution kernel according to an embodiment of the present application, as shown in FIG. Figure 6b As shown in the figure, for the axial 1D convolution kernel, by scanning the surface sources at multiple different positions and reconstructing the image, the reconstructed image is accumulated to the z-axis to obtain the reconstructed image data, and the reconstructed image data at different positions are averaged to obtain the axial 1D convolution kernel.

[0113] According to some embodiments, the thickness of the surface source is less than or equal to the side length of a voxel in the attenuation image, and / or the area of ​​the surface source is greater than or equal to 1 mm. 2 , the surface source is perpendicular to the z-axis.

[0114] According to some embodiments, the MLEM algorithm is used to reconstruct the surface source scanning data at multiple different positions.

[0115] According to some embodiments, the reconstructed image is accumulated to the z-axis according to the position of the surface source to obtain reconstructed image data, the reconstructed image data is aligned and averaged, and the averaged value is the 1D convolution kernel κ a (z).

[0116] According to some embodiments, the convolution kernel of each voxel is normalized by formula (12).

[0117]

[0118] In step S105 , a deconvolution operation is performed on the voxel values ​​of the attenuation image using the normalized convolution kernel to obtain a corrected attenuation image.

[0119] According to some example embodiments of the present application, the voxel values ​​of the attenuation image are corrected using the normalized convolution kernel through formula (13).

[0120] μ corr (x,y,z)=∑ j κ′ j (x,y,z)μ j Formula (13)

[0121] Among them, μ corr (x, y, z) is the voxel value of the corrected attenuation image at the position (x, y, z), μ j is the voxel value of the jth voxel in the attenuation image before correction, κ′ j (x, y, z) is the normalized value of the convolution kernel of the j-th voxel at the (x, y, z) position, κ j (x,y,z) is the convolution kernel of the j-th voxel.

[0122] According to some embodiments, after obtaining the voxel values ​​of the corrected attenuation image, an attenuation correction coefficient is calculated using formula (1) or formula (2). The calculated attenuation correction coefficient is used to perform image reconstruction using formulas (3) and (4) to obtain a reconstructed image with higher resolution.

[0123] pass Figure 3 The image processing method shown in the figure corrects the attenuation image without modifying the integral formula for calculating the linear attenuation coefficient. The attenuation correction coefficient is calculated using the corrected attenuation image, so that the attenuation correction coefficient will not be calculated over a long distance in the high-density area. At the same time, there is a certain "blurring" effect on the attenuation image, which causes the attenuation image to produce an "inward-shrinking" deformation.

[0124] Figure 3 The entire image processing process shown does not need to consider whether the detector has a certain volume, transmission and scattering of gamma photons between the detectors, and the like, so that Figure 1b 、 Figure 1c and Figure 1d The situation shown is comprehensively considered, Figure 1b 、 Figure 1c and Figure 1d The problems in the above methods are solved.

[0125] Figure 7a It is a cross-sectional material distribution map of a prosthesis 1 with uniform activity distribution. Figure 7b It is a coronal plane material distribution map of a prosthesis 2 with uniform activity distribution. As shown in Figure 7a and Figure 7b , the water and polytetrafluoroethylene (PTFE) distribution are marked with numerical values in the figure. The linear attenuation coefficient of water is about 0.096 cm -1 , and the linear attenuation coefficient of polytetrafluoroethylene is about 0.18 cm -1 , which is close to the linear attenuation coefficient of bone.

[0126] Figure 8a The reconstructed image of the prosthesis 1 using the normal linear attenuation coefficient integration method is shown. Figure 8b The reconstructed image of the prosthesis 1 using the average depth setting method is shown. Figure 8c The reconstructed image of the prosthesis 1 using the multiple linear attenuation coefficient integration method is shown. Figure 8d The reconstructed image of the prosthesis 1 using the correction method shown in Figure 3 is shown.

[0127] Figure 9a The reconstructed image of the prosthesis 2 using the normal linear attenuation coefficient integration method is shown. Figure 9b The reconstructed image of the prosthesis 2 using the average depth setting method is shown. Figure 9c The reconstructed image of the prosthesis 2 using the multiple linear attenuation coefficient integration method is shown. Figure 9d The reconstructed image of the prosthesis 2 using the correction method shown in Figure 3 is shown.

[0128] From the reconstructed images shown in Figures 8a to 8d and Figures 9a to 9d , it can be seen that the image obtained by using the image processing method described in the embodiments of the present application is the most uniform, and there is basically no visible artifact. The reconstructed images of the other three methods still have visible artifacts. Therefore, the effect of the reconstructed image method using the example embodiments of the present application is better than that of the other three methods.

[0129] Figure 10 A block diagram of an image processing apparatus is shown according to an embodiment of the present application. As shown in Figure 10 The image processing apparatus shown according to an embodiment of the present application comprises a voxel value obtaining unit 1001, a convolution kernel calculating unit 1003 and an attenuation image correcting unit 1005.

[0130] The voxel value obtaining unit 1001 is configured to obtain voxel values of a plurality of voxels in an attenuation image to be corrected, the convolution kernel calculating unit 1003 is configured to obtain a convolution kernel and normalize the convolution kernel, and the attenuation image correcting unit 1005 is configured to perform a deconvolution operation on the voxel values of the attenuation image by using the normalized convolution kernel to obtain a corrected attenuation image.

[0131] Figure 11 A block diagram of another image processing apparatus is shown according to an embodiment of the present application. Figure 11 The image processing apparatus shown is merely an example and should not bring any limitation to the functions and usage range of the embodiments of the present application.

[0132] As shown in Figure 11 The image processing apparatus is shown in the form of a general computing device. The components of the image processing apparatus can include, but are not limited to, at least one processor 210, at least one memory 220, a bus 230 connecting different system components including the memory 220 and the processor 210, a display unit 240, etc. The memory 220 stores program codes which can be executed by the processor 210 to make the processor 210 perform the methods according to various exemplary embodiments of the present application described in the present specification. For example, the processor 210 can perform the method shown in Figure 3 .

[0133] The memory 220 can include a readable medium in the form of a volatile storage unit, such as a random access memory (RAM) 2201 and / or a cache memory unit 2202, and can further include a read-only memory (ROM) 2203.

[0134] The memory 220 can further include program / utilities 2204 having a set of (at least one) program modules 2205, such as an operating system, one or more application programs, other program modules, and program data, each of which can include implementation of a network environment or some combination of these examples.

[0135] The bus 230 can represent one or more of several types of bus structures, including a storage unit bus or bus controller, a peripheral bus, a graphics acceleration port, a processing unit bus, or a local bus using any of a variety of bus architectures.

[0136] The image processing device can also communicate with one or more external devices 300 such as a keyboard, a pointing device, a Bluetooth device, etc., and can also communicate with one or more devices that enable a user to interact with the image processing device, and / or any devices (e.g., routers, modems, etc.) that enable the image processing device to communicate with one or more other computing devices. Such communication can occur via an input / output (I / O) interface 250. Still yet, the image processing device can communicate with one or more networks (such as a local area network (LAN), a wide area network (WAN), and / or the public network, such as the Internet) via a network adapter 260. The network adapter 260 can communicate with the other components of the image processing device via the bus 230. It should be appreciated that although not shown, other hardware and / or software modules could be used in connection with the image processing device. Such modules include, but are not limited to, microcode, device drivers, redundant processing units, external disk drive arrays, RAID systems, tape drives, and data archival storage systems, etc.

[0137] Those skilled in the art will readily understand that the example embodiments described herein can be implemented by software and / or by software in combination with the necessary hardware. The technical solutions according to the embodiments of the present application can be embodied in the form of a software product, which can be stored in a computer readable storage medium (which can be a CD-ROM, a U disk, a mobile hard disk, etc.) or a network, and includes a plurality of computer program instructions to make a computer device (which can be a personal computer, a server, or a network device, etc.) execute the above-mentioned methods according to the embodiments of the present application.

[0138] The software product can adopt any combination of one or more readable media. The readable medium can be a readable signal medium or a readable storage medium. The readable storage medium may, for example, be but is not limited to an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, device or apparatus, or any combination of the above. More specific examples (non-exhaustive list) of the readable storage medium include an electrical connection having one or more wires, a portable disk, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or flash memory), an optical fiber, a portable compact disk read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the above.

[0139] The computer readable storage medium can include a computer-readable medium in baseband or propagated as a carrier wave in a propagated signal, wherein the latter encompasses a baseband propagation. The propagated signal can take a wide variety of forms including, but not limited to, electro-magnetic, optical, or any suitable combination thereof. A computer readable storage medium can be any medium (medium or media) suitable for storing or transmitting program code in the form of instructions or data structures. Program Code embodied on a computer-readable medium can be transmitted using any appropriate medium, including but not limited to wireless, wired, optical fiber cable, RF, etc., or any suitable combination thereof.

[0140] The program code can be implemented in any of various ways, including procedure-based, object-based, or component-based technologies, and the program code can be implemented across several devices or spread across network coupled devices, which can use various languages including C, C++, Java, and others.

[0141] The above computer readable medium stores one or more program instructions, when the one or more program instructions are executed by a device, the computer readable medium enables the device to implement the above functions.

[0142] Those skilled in the art can understand that the above modules can be distributed in the device according to the description of the embodiments, and can also be changed to be in one or more devices different from the embodiments. The modules of the above embodiments can be combined into one module, and one module can be further split into multiple sub-modules.

[0143] From the above description of the embodiments, those skilled in the art can easily understand that the example embodiments described herein can be implemented by software, or by software in combination with necessary hardware. The technical solutions according to the embodiments of the present application can be embodied in the form of a software product, which can be stored in a computer readable storage medium (which can be a CD-ROM, a U disk, a mobile hard disk, etc.) or a network, and includes a plurality of computer program instructions to enable a computing device (which can be a personal computer, a server, or a network device, etc.) to execute the above-mentioned methods according to the embodiments of the present application.

[0144] The software product can employ any combination of one or more readable media. The readable media can be a readable signal medium or a readable storage medium. The readable storage medium, for example, can be but is not limited to an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any suitable combination of the foregoing. More specific examples (a non-exhaustive list) of the readable storage medium include an electrical connection having one or more wires, a portable disk, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or flash memory), an optical fiber, a portable compact disc read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the foregoing.

[0145] The computer readable storage medium can include a data signal transported, for example, over a baseband or a carrier wave by a contract bridge, where the contract bridge can be any data gate or physical layer that conveys the program code. The data signal can be transmitted by using any appropriate medium, including but not limited to the use of electronic, electromagnetic, optical, and radio frequencies, or any suitable combination of the foregoing.

[0146] The program code for carrying out operations of the present application can be written in any combination of one or more programming languages, including an object oriented programming language such as Java, C++, or the like, and conventional procedural programming languages, such as the "C" programming language or similar programming languages. The program code can execute entirely on the user's computing device, partly on the user's computing device, as a stand-alone software package, partly on the user's computing device and partly on a remote computing device or entirely on the remote computing device or server. In the latter scenario, the remote computing device can be connected to the user's computing device through any type of network, including a local area network (LAN) or a wide area network (WAN), or the connection can be made to an external computing device, such as through the Internet using an Internet Service Provider.

[0147] The above computer readable medium carries one or more program instructions, when the one or more program instructions are executed by a device, the computer readable medium enables the foregoing functions.

[0148] Those skilled in the art can understand that the above modules can be distributed in the device according to the description of the embodiments, and can also be changed in one or more devices different from the embodiments. The modules of the above embodiments can be combined into one module, and one module can be further split into multiple sub-modules.

[0149] According to some example embodiments of the present application, by correcting the attenuation image, the attenuation correction coefficient is calculated using the corrected attenuation image, so that when calculating the attenuation correction coefficient, it will not pass through the high-density area for too long a distance, and at the same time, the attenuation image also has a certain "blurring" effect, so that the attenuation image produces a kind of "internal shrinkage" deformation. Moreover, the whole correction process does not need to consider whether the detector exists a certain volume, the transmission and scattering of gamma photons between the detectors and other problems.

[0150] Although the present application provides the method operation steps as described in the above embodiments or flowcharts, more or less operation steps can be included in the method based on conventional or non-creative labor. In steps that do not have necessary causal relationship in logic, the execution order of these steps is not limited to the execution order provided in the embodiments of the present application.

[0151] Each of the embodiments in the specification is described in a progressive manner, and the same or similar parts between each embodiment can be referred to each other. Each embodiment focuses on the difference from other embodiments.

[0152] The above describes the embodiments of the present application in detail, and the specific examples are applied to explain the principles and implementation modes of the present application. The above embodiment description is only used to help understand the method of the present application and its core idea. Meanwhile, the changes or deformations made by the skilled in the art according to the idea of the present application, based on the specific implementation mode and application range of the present application, all belong to the protection scope of the present application. In summary, the content of the specification should not be understood as a limitation of the present application.

Claims

1. An image processing method, characterized by, The image processing method comprises: acquiring voxel values of a plurality of voxels in an attenuation image to be corrected; obtaining a convolution kernel by performing a scanning experiment in which an additional line source or area source is arranged when the attenuation image is acquired, and normalizing the convolution kernel; performing deconvolution operation on the voxel values of the attenuation image by using the normalized convolution kernel to obtain voxel values of a corrected attenuation image; wherein the convolution kernel can be decoupled into multiplication of an x-y plane 2D convolution kernel and an axial convolution kernel, the x-y plane 2D convolution kernel has rotational symmetry around an origin, the 2D convolution kernel of the x-y plane in (r, 0) coordinates satisfies a PiG model, and the PiG model satisfies formula (9): where A, B are the amplitude of the Gaussian function, σ x0 ,σ y0 are the spread of the Gaussian function in the x and y direction on the side close to the origin, σ x1 ,σ y1 are the spread of the Gaussian function in the x and y direction on the side far from the origin, x0 is the offset of the 2D kernel on the x axis, x' = xcosθ + ysinθ, y' = -xsinθ + ycosθ, is the 2D kernel of the voxel on the x-y plane (r, 0) coordinate, 2. The image processing method of claim 1, wherein, the convolution kernel is calculated by formula (5): where (jx, jy, jz) is the coordinate corresponding to the jth voxel, (x, y, z) takes any value within the imaging field of view, and Kj(x, y, z) is the value of the convolution kernel of the jth voxel at the spatial coordinate (x, y, z), is a 2D convolution kernel in the x-y plane, and K a(z) is a 1D convolution kernel in the axial direction.

3. The image processing method of claim 2, wherein, the 2D convolution kernel of the x-y plane satisfies formula (6):

4. The image processing method of claim 1, wherein, A = B, σx0 = σx1 = σx, σy0 = σy1 = σy, the 2D convolution kernel of the x-y plane in (r, 0) coordinates satisfies a SiG model, and the SiG model satisfies formula (10):

5. The image processing method of claim 4, wherein, the A, B, σx0, σy0, σx1, σy1, x0 are obtained in the following manner: scanning a plurality of line sources at different positions and reconstructing images; accumulating the reconstructed images to an x-y plane; fitting the PiG or the SiG model by using the images of the x-y plane to obtain the A, B, σx0, σy0, σx1, σy1, x0.

6. The image processing method according to claim 5, wherein: a diameter of the line source is less than or equal to a voxel edge length of the attenuation image; a length of the line source is greater than or equal to 1 mm; the line source is parallel to a z axis.

7. The image processing method according to claim 4, wherein: the A, B, σx0, σy0, σx1, σy1, x0 are calculated by using a 1D linear interpolation algorithm according to known A, B, σx0, σy0, σx1, σy1, x0 of adjacent positions.

8. The image processing method of claim 2, wherein, the axial 1D convolution kernel is obtained in the following manner: scanning a plurality of area sources at different positions and reconstructing images; accumulating the reconstructed images to a z axis to obtain reconstructed image data; averaging the reconstructed image data at different positions to obtain the axial 1D convolution kernel.

9. The image processing method according to claim 8, wherein: a thickness of the area source is less than or equal to a voxel edge length of the attenuation image; The area of the area source is greater than or equal to 1 mm 2 ; the area source is perpendicular to the z axis.

10. The image processing method of claim 2, wherein, the convolution kernel of each voxel is normalized by formula (12): where κj(x, y, z) is the normalized value of the convolution kernel of the jth voxel at the position (x, y, z), and κj(x, y, z) is the convolution kernel of the jth voxel. j (x,y,z) is the normalized value of the convolution kernel of the jth voxel at the position (x, y, z), and κj(x, y, z) is the convolution kernel of the jth voxel.

11. The image processing method of claim 10, wherein, deconvolution operation is performed on the voxel values of the attenuation image by using the normalized convolution kernel by formula (13) to obtain voxel values of a corrected attenuation image: μcorr(x, y, z) = ∑ jκj (x, y, z) μj formula (13) wherein μcorr(x, y, z) is a voxel value of the corrected attenuation image at a position (x, y, z), and μj is a voxel value of the attenuation image before correction at a jth voxel.

12. The image processing method of claim 10, wherein, The image processing method further comprises: calculating an attenuation correction coefficient by using the voxel values of the corrected attenuation image; performing image reconstruction by using the attenuation correction coefficient.

13. The image processing method according to claim 1, wherein: An attenuation image is acquired by a CT image or an MR image.

14. An image processing apparatus characterized by comprising: The image processing apparatus comprises: a voxel value acquisition unit configured to acquire voxel values of a plurality of voxels in an attenuation image to be corrected; a convolution kernel calculation unit configured to acquire a convolution kernel by performing a scanning experiment in which an additional line source or surface source is arranged in acquiring the attenuation image, and normalize the convolution kernel; an attenuation image correction unit configured to perform a deconvolution operation on the voxel values of the attenuation image using the normalized convolution kernel to obtain a corrected attenuation image; wherein the convolution kernel can be decoupled into multiplication of an x-y plane 2D convolution kernel and an axial convolution kernel, the x-y plane 2D convolution kernel has rotational symmetry around an origin, and a 2D convolution kernel on a (r, 0) coordinate of the x-y plane satisfies a PiG model, the PiG model satisfies formula (9): where A, B are the amplitude of the Gaussian function, σχ0, σγ0are the spread of the Gaussian function in x and y direction on the side close to the origin, σχ1, σγ1are the spread of the Gaussian function in x and y direction on the side far from the origin, x0is the offset of the 2D kernel on the x-axis, x' = xcosθ + ysinθ, y' = -xsinθ + ycosθ, ' is the 2D kernel of the voxel on the x-y plane (r, 0) coordinate, ​ 15. An electronic device, comprising: comprising: one or more processors; a storage device configured to store a computer program; when the computer program is executed by the one or more processors, the one or more processors are caused to implement the method according to any one of claims 1-13.

16. A computer-readable storage medium, characterized in that, a program instruction is stored thereon, and the program instruction is executed to implement the method according to any one of claims 1-13.

Images