WO2022213321A1

WO2022213321A1 - Wavelet fusion-based pet image reconstruction method

Info

Publication number: WO2022213321A1
Application number: PCT/CN2021/085939
Authority: WO
Inventors: 万丽雯; 李彦明; 郑海荣; 刘新; 任溥弦
Original assignee: 深圳高性能医疗器械国家研究院有限公司
Priority date: 2021-04-08
Filing date: 2021-04-08
Publication date: 2022-10-13

Abstract

Disclosed in the present invention is a wavelet fusion-based PET image reconstruction method. The method comprises: performing wavelet decomposition on a first PET image to obtain a corresponding first high-frequency image and a corresponding first low-frequency image, and performing wavelet decomposition on a second PET image to obtain a corresponding second high-frequency image and a corresponding second low-frequency image, wherein the first PET image and the second PET image are PET reconstructed images obtained by means of different modes; within different frequency channels, processing a wavelet coefficient by using a set fusion rule according to the wavelet coefficient, and performing inverse wavelet transformation on a processed new wavelet coefficient to obtain a fused image; and evaluating the fused image and iterating to obtain an optimal solution, and using the obtained fused image that meets a set evaluation criterion as a reconstructed image. By means of the reconstructed image obtained in the present invention, the contour and details of an image can be well preserved.

Description

A PET Image Reconstruction Method Based on Wavelet Fusion

technical field

The invention relates to the technical field of medical image processing, and more particularly, to a PET image reconstruction method based on wavelet fusion.

Background technique

Positron emission tomography (PET) is a highly sensitive and versatile molecular imaging method for metabolic imaging and is widely used in clinical practice. Using PET images, it is possible to find the early stage of tumors, investigate and deal with cardiovascular diseases and neurological diseases. However, since PET scans need to inject contrast agents into patients, long-term dynamic scans and high doses of drugs can make PET images clearer, but they increase the burden on patients. Therefore, it is essential to use the low-dose (As Low As Reasonably Achievable, ALARA) principle reasonably under the premise of satisfying the clinical diagnosis, and to minimize the radiation dose to patients. However, since the PET image reconstruction decreases with the dose, the image quality will be poor, which will also affect the accuracy of the Patlak map and some other data. Therefore, research and development of new low-dose PET image reconstruction methods can reduce the amount of drug injected to patients under the effect of PET image reconstruction and Patlak map, which is easier for patients to accept and reduce the risk of high doses for clinical medical methods. application prospects.

In the prior art, there are many kinds of PET image reconstruction algorithms, each with its own advantages and disadvantages. For the reconstruction of parametric images, although EM (expectation maximization) can reconstruct the image very well, the noise is large and the count is too high, while MR Guided Kernel EM can reduce noise and reconstruct the image at lower counts (the number of photons in opposite directions generated by positive and negative e-cigarette extinguishing), but the details are not good enough, and for the part of interest Ki in the parametric image, the net inflow rate is lower than the actual value.

SUMMARY OF THE INVENTION

The purpose of the present invention is to overcome the above-mentioned defects of the prior art, and to provide a PET image reconstruction method based on wavelet fusion. Based on the wavelet transform, the images reconstructed by Kernel EM and EM are decomposed respectively, and then features are extracted to perform frame-by-frame fusion, so that the PET image is decomposed. Reconstruction results in relatively clear images at reduced doses.

The technical solution of the present invention is to provide a PET image reconstruction method based on wavelet fusion. The method includes the following steps:

Perform wavelet decomposition on the first PET image to obtain a corresponding first high-frequency image and a first low-frequency image, and perform wavelet decomposition on the second PET image to obtain a corresponding second high-frequency image and second low-frequency image, wherein the first PET image The image and the second PET are PET reconstructed images obtained in different ways;

In different frequency channels, the wavelet coefficients are processed according to the set fusion criteria according to the wavelet coefficients, and the processed new wavelet coefficients are subjected to wavelet inverse transformation to obtain the fused image;

The fused image is evaluated and the optimal solution is obtained iteratively, and the obtained fused image that meets the set evaluation criteria is used as the reconstructed image.

Compared with the prior art, the advantage of the present invention is that, in order to solve the technical problem of improving the quality of the low-dose image, the present invention uses its respective advantages on the basis of the reconstruction algorithm to convert it from the time domain to the time domain through wavelet decomposition. The frequency domain extracts the advantages of their respective reconstruction effects in the form of weights and fuses them to improve the quality of the reconstructed image, and uses wavelets to extract high and low frequency information and feature fusion. The decision of feature fusion, extracting the advantages of the results of different image reconstruction methods, improves the image quality.

Other features and advantages of the present invention will become apparent from the following detailed description of exemplary embodiments of the present invention with reference to the accompanying drawings.

Description of drawings

The accompanying drawings, which are incorporated in and constitute a part of this specification, illustrate embodiments of the invention and together with the description serve to explain the principles of the invention.

1 is a flowchart of a method for reconstructing a PET image based on wavelet fusion according to an embodiment of the present invention;

2 is a schematic process diagram of a method for reconstructing a PET image based on wavelet fusion according to an embodiment of the present invention;

FIG. 3 is a frequency distribution diagram after performing three-layer wavelet transform on an image according to an embodiment of the present invention;

FIG. 4 is a schematic process diagram of a method for reconstructing a PET image based on wavelet fusion according to another embodiment of the present invention.

Detailed ways

Various exemplary embodiments of the present invention will now be described in detail with reference to the accompanying drawings. It should be noted that the relative arrangement of components and steps, the numerical expressions and numerical values set forth in these embodiments do not limit the scope of the invention unless specifically stated otherwise.

The following description of at least one exemplary embodiment is merely illustrative in nature and is in no way intended to limit the invention, its application, or uses.

Techniques, methods, and apparatus known to those of ordinary skill in the relevant art may not be discussed in detail, but where appropriate, such techniques, methods, and apparatus should be considered part of the specification.

In all examples shown and discussed herein, any specific values should be construed as illustrative only and not limiting. Accordingly, other instances of the exemplary embodiment may have different values.

It should be noted that like numerals and letters refer to like items in the following figures, so once an item is defined in one figure, it does not require further discussion in subsequent figures.

1 and 2, the PET image reconstruction method based on wavelet fusion provided by the present invention includes the following steps.

Step S1, reconstruct the PET image to obtain a reconstructed image.

Specifically, the likelihood function composed of dataset A is expressed as:

L(A|X)=∏ _i ∏ _j p(A _ij |β=H _ij *X _j ) (1)

where β represents the expectation of

E(β)=H _ij *X _j (2)

The probability of the Poisson distribution when the variable takes A _ij , H _ij is a constant, which represents the system response matrix, which is expressed as the probability that the gamma photon emitted by the jth pixel is detected by the ith LOR, X is a dimensional array to represent a two-dimensional image.

At this time, the probability distribution function of A _ij is expressed as:

After taking logarithmic simplification, we can get

Among them, const represents some other quantity unrelated to X.

In order to eliminate the random variable A _ij in the above formula, the original random variable A _ij is replaced by A _ij under the condition that the expected number X of the number of photons emitted by all pixels facing all directions obtained by the maximum likelihood is known, so that achieve the purpose of iterative update.

Because an iterative update is required to reconstruct the image, in an iteration, let the expectation of the current A _ij be the A _ij in the next iteration. X ^current is the current image, and c is used to represent current.

which is

Find the expectation formula from the binomial distribution: if r～B(r,p), then E(r)=np, it can be known that _n is Pi and p is

and β _ij =H _ij *X _j .

thus obtainable

For the obtained logarithmically simplified formula L(A|X)=∑ _i ∑ _j (A _ij ln(H _ij X _j )-H _ij X _j ), find the partial derivative to get its extreme point :

Substitute the obtained A _ij into, let the partial derivative be 0 to find its extreme point:

An updated X _j is obtained each time X _j is iterated as X ^current . That is to get the iterative formula:

In step S2, the dynamic frame in the dynamic PET is fused into the forward projection model as prior information to obtain a reconstructed image.

Specifically, the original data forward projection model

y=PKα+r (12)

Among them, α represents the coefficient vector of K, P represents the system matrix, r represents the scattering of photons and random numbers, K is the key, and K represents the relationship between the neighbors f _j of each pixel f _i using the KNN method.

K _ij =K(fi , _{f j} ₎ (13)

KNN (k-nearest neighbor algorithm) is sought based on its Euclidean distance, and dynamic frames in dynamic PET are fused into the forward projection model as a priori information, thereby changing the PET image reconstruction of low-count scans. Step S1 is used as the basis for reconstructing the image, and the time frame is divided into three time periods to be accumulated to obtain features, which are used as prior information to reconstruct the PET image.

Similar to step S1, combine the projection model

y=PKα+r (14)

Use the Poisson-likelihood function to estimate the kernel's coefficient image.

In the above formula, bP(a) is the penalty function, and b is its regularization parameter. After obtaining the kernel-based coefficient image, the reconstructed image can be obtained according to the projection model:

The reconstructed image x has been regularized by K, so when estimating the maximum likelihood function, the regularization parameter b can be set to 0, and the simplified formula is:

Similar to step S1, after performing the desired calculation and logarithmic simplification, the iterative formula based on K can be obtained:

To sum up, the present invention reconstructs dynamic PET images, obtains each slice image corresponding to the first group of time frames, reconstructs the PET image based on the kernel method, and obtains each slice image corresponding to the second group of time frames, wherein the first group of time frames corresponds to each slice image. And the second group of image reconstruction results are PET reconstructed images obtained according to different ways.

In step S3, the images obtained in step S1 and step S2 are extracted frame by frame using wavelet transform for fusion.

Wavelet transform is to decompose the image into sub-images of each frequency segment in the frequency domain, which are used to represent each feature component of the original image. That is, the image is decomposed into a low-frequency image and a high-frequency image after being decomposed by wavelet transform, and the coefficients of low-frequency and high-frequency are obtained at the same time. The low-frequency image contains low-frequency components, and the high-frequency image contains high-frequency components. , the high frequency components are the details of the image. Then, in different frequency channels, different fusion criteria are used to process the wavelet coefficients according to the difference of the wavelet coefficients. The processed new wavelet coefficients preserve more frequency band features. Finally, inverse wavelet transform is performed on the new wavelet coefficients to obtain the fused image.

Specifically, after each wavelet transform, the image is decomposed into four frequency bands: LL, HL, LH, and HH.

LL band, this band maintains the original image content information, the energy of the image is concentrated in this band:

HL band, which preserves high-frequency information in the horizontal direction of the image:

LH band, which preserves high-frequency information in the vertical direction of the image:

HH band, this band maintains the high frequency information of the image in the diagonal direction:

Among them, f(k,l) represents the details of the image projected at each scale. (m,n) represents the pixels of the image, and j represents the number of layers of low-pass and high-pass decomposition coefficients.

For example, when the image is subjected to three-layer wavelet transform, the frequency distribution diagram is shown in Figure 3. It should be understood that more layers of wavelet transform may also be employed.

In one embodiment, taking the three-layer wavelet transform as an example, for the high-frequency H1(j,k) and H2(j,k) of the two images, take their maximum value as the high-frequency information H( j, k), j=1, 2, 3, j is the number of decomposition layers. For the low-frequency information L1(j,k) and L2(j,k) of the two images, the low-frequency information L(j,k) of the fused image is obtained by averaging.

In step S4, the fusion result is evaluated and an optimal solution is obtained iteratively.

For example, as shown in FIG. 4 , for the obtained fused image, the mean value and the standard deviation are used to judge the fusion result. Among them, the mean value can reflect the average brightness of the image, and the standard deviation reflects the degree of dispersion of the image gray level relative to the mean value. A larger standard deviation means that most values are more different from their mean; a smaller standard deviation means that these values are closer to the mean. That is, the smaller the standard deviation, the more accurate the data.

For example, the fusion results are evaluated using the information entropy of the image, which defines the average amount of information obtained when observing the output of a single source symbol. The maximum entropy occurs when the probability of occurrence of each symbol of the source is equal, and the source provides the maximum possible average information amount of the source symbol at this time. The greater the information entropy, the greater the amount of information contained in the image, that is, the better the fusion effect.

For example, the average gradient can be used to represent an improvement in the sharpness of the image. Image gradients can also reflect small detail contrast and texture transformation characteristics in the image:

Among them, A and B represent the size of the image, Δ _x f(i, j) and Δ _y f(i, j) are the first-order differences of the pixel (i, j) in the x and y directions, respectively

Δ _x f(i,j)=f(i+1,j)-f(i,j) (24)

Δy _f (i,j)=f(i,j+1)-f(i,j) (25)

Each frame of dynamic PET imaging is iteratively fused to obtain a wavelet-fused dynamic PET image.

It should be noted that, in addition to using the images obtained in steps S1 and S2 for fusion, the present invention can also use other reconstruction algorithms to perform fusion according to the reconstruction effect after proper deformation, for example, filtered back-projection reconstruction method, based on least squares PET image reconstruction, etc.

To sum up, the present invention performs frame-by-frame fusion for details of high-frequency components and contours of low-frequency components by considering the differences in the results of different PET image reconstruction algorithms, thereby reducing image noise after image reconstruction and improving image details. In addition, on the basis of reducing the count of reconstructed PET images, the reconstructed images based on low-dose images can preserve the contours and details well, and improve the fitting degree between the net inflow of the tumor area and the true value.

The present invention may be a system, method and/or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for causing a processor to implement various aspects of the present invention.

A computer-readable storage medium may be a tangible device that can hold and store instructions for use by the instruction execution device. The computer-readable storage medium may be, for example, but not limited to, an electrical storage device, a magnetic storage device, an optical storage device, an electromagnetic storage device, a semiconductor storage device, or any suitable combination of the foregoing. More specific examples (non-exhaustive list) of computer readable storage media include: portable computer disks, hard disks, random access memory (RAM), read only memory (ROM), erasable programmable read only memory (EPROM) or flash memory), static random access memory (SRAM), portable compact disk read only memory (CD-ROM), digital versatile disk (DVD), memory sticks, floppy disks, mechanically coded devices, such as printers with instructions stored thereon Hole cards or raised structures in grooves, and any suitable combination of the above. Computer-readable storage media, as used herein, are not to be construed as transient signals per se, such as radio waves or other freely propagating electromagnetic waves, electromagnetic waves propagating through waveguides or other transmission media (eg, light pulses through fiber optic cables), or through electrical wires transmitted electrical signals.

The computer readable program instructions described herein may be downloaded to various computing/processing devices from a computer readable storage medium, or to an external computer or external storage device over a network such as the Internet, a local area network, a wide area network, and/or a wireless network. The network may include copper transmission cables, fiber optic transmission, wireless transmission, routers, firewalls, switches, gateway computers, and/or edge servers. A network adapter card or network interface in each computing/processing device receives computer-readable program instructions from a network and forwards the computer-readable program instructions for storage in a computer-readable storage medium in each computing/processing device .

The computer program instructions for carrying out the operations of the present invention may be assembly instructions, instruction set architecture (ISA) instructions, machine instructions, machine dependent instructions, microcode, firmware instructions, state setting data, or instructions in one or more programming languages. Source or object code written in any combination, including object-oriented programming languages, such as Smalltalk, C++, Python, etc., and conventional procedural programming languages, such as the "C" language or similar programming languages. The computer readable program instructions may execute entirely on the user's computer, partly on the user's computer, as a stand-alone software package, partly on the user's computer and partly on a remote computer, or entirely on the remote computer or server implement. In the case of a remote computer, the remote computer may be connected to the user's computer through any kind of network, including a local area network (LAN) or a wide area network (WAN), or may be connected to an external computer (eg, using an Internet service provider through the Internet connect). In some embodiments, custom electronic circuits, such as programmable logic circuits, field programmable gate arrays (FPGAs), or programmable logic arrays (PLAs), can be personalized by utilizing state information of computer readable program instructions. Computer readable program instructions are executed to implement various aspects of the present invention.

Aspects of the present invention are described herein with reference to flowchart illustrations and/or block diagrams of methods, apparatus (systems) and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and/or block diagrams, and combinations of blocks in the flowchart illustrations and/or block diagrams, can be implemented by computer readable program instructions.

These computer readable program instructions may be provided to a processor of a general purpose computer, special purpose computer or other programmable data processing apparatus to produce a machine that causes the instructions when executed by the processor of the computer or other programmable data processing apparatus , resulting in means for implementing the functions/acts specified in one or more blocks of the flowchart and/or block diagrams. These computer readable program instructions can also be stored in a computer readable storage medium, these instructions cause a computer, programmable data processing apparatus and/or other equipment to operate in a specific manner, so that the computer readable medium on which the instructions are stored includes An article of manufacture comprising instructions for implementing various aspects of the functions/acts specified in one or more blocks of the flowchart and/or block diagrams.

Computer readable program instructions can also be loaded onto a computer, other programmable data processing apparatus, or other equipment to cause a series of operational steps to be performed on the computer, other programmable data processing apparatus, or other equipment to produce a computer-implemented process , thereby causing instructions executing on a computer, other programmable data processing apparatus, or other device to implement the functions/acts specified in one or more blocks of the flowcharts and/or block diagrams.

The flowchart and block diagrams in the Figures illustrate the architecture, functionality, and operation of possible implementations of systems, methods and computer program products according to various embodiments of the present invention. In this regard, each block in the flowchart or block diagrams may represent a module, segment, or portion of instructions, which comprises one or more functions for implementing the specified logical function(s) executable instructions. In some alternative implementations, the functions noted in the blocks may occur out of the order noted in the figures. For example, two blocks in succession may, in fact, be executed substantially concurrently, or the blocks may sometimes be executed in the reverse order, depending upon the functionality involved. It is also noted that each block of the block diagrams and/or flowchart illustrations, and combinations of blocks in the block diagrams and/or flowchart illustrations, can be implemented in dedicated hardware-based systems that perform the specified functions or actions , or can be implemented in a combination of dedicated hardware and computer instructions. It is well known to those skilled in the art that implementation in hardware, implementation in software, and implementation in a combination of software and hardware are all equivalent.

Various embodiments of the present invention have been described above, and the foregoing descriptions are exemplary, not exhaustive, and not limiting of the disclosed embodiments. Numerous modifications and variations will be apparent to those of ordinary skill in the art without departing from the scope and spirit of the described embodiments. The terminology used herein was chosen to best explain the principles of the embodiments, the practical application or technical improvement in the marketplace, or to enable others of ordinary skill in the art to understand the embodiments disclosed herein. The scope of the invention is defined by the appended claims.

Claims

A PET image reconstruction method based on wavelet fusion, comprising the following steps:

Step S11, performing wavelet decomposition on the first PET image to obtain a corresponding first high-frequency image and a first low-frequency image, and performing wavelet decomposition on the second PET image to obtain a corresponding second high-frequency image and a second low-frequency image, wherein The first PET image and the second PET are PET reconstructed images obtained according to different methods;

Step S12, in different frequency channels, process the wavelet coefficients according to the set fusion criteria according to the wavelet coefficients, and perform wavelet inverse transformation on the processed new wavelet coefficients to obtain a fused image;

In step S13, the fused image is evaluated and an optimal solution is obtained iteratively, and the obtained fused image satisfying the set evaluation standard is used as a reconstructed image.
The method according to claim 1, wherein step S11 comprises:

Multi-layer wavelet decomposition is performed on the first PET image and the second PET image respectively, and the obtained first high-frequency image is expressed as H1(j,k), the first low-frequency image is expressed as L1(j,k), and the second high-frequency image is expressed as L1(j,k). The low-frequency image is represented as H2(j,k), and the second low-frequency image is represented as L1(j,k), where j represents the layer index of wavelet decomposition.
The method according to claim 2, wherein step S12 comprises:

For the first high-frequency image H1(j,k) and the second high-frequency image H2(j,k), take the maximum value as the high-frequency image information H(j,k) of the fused image;

For the first low-frequency image L1(j,k) and the second low-frequency image L2(j,k), the low-frequency image information L(j,k) of the fused image is obtained by taking an average.
The method according to claim 1, wherein the fused image is evaluated by one or more of mean, standard deviation and mean gradient, wherein the mean is used to represent the average brightness of the fused image, and the standard deviation It reflects the degree of dispersion of the gray level of the fusion image relative to the mean value, and the average gradient reflects the small detail contrast and texture transformation characteristics in the fusion image.
The method according to claim 4, wherein the average gradient is expressed as:

Among them, A and B represent the size of the image, and Δ x f(i, j) and Δ y f(i, j) are the first-order differences of the pixel point (i, j) in the x and y directions, respectively.
The method of claim 1, wherein the first PET image is a PET reconstructed image obtained according to a Kernel EM algorithm.
The method of claim 1, wherein the second PET image is obtained according to the steps of:

Combine the projection model y=PKα+r with the Poisson likelihood function to estimate the coefficient image of the kernel:

After the kernel-based coefficient image is obtained, the reconstructed image is obtained according to the projection model, which is expressed as:

Let the regularization parameter b be 0, when estimating the maximum likelihood function, it is expressed as:

After performing the desired calculations and logarithmic simplification, the iterative formula based on K is obtained:

Among them, α represents the coefficient vector of K, P represents the system matrix, r represents the scattering of photons and random numbers, K represents the relationship between the neighbors f j of each pixel f i , bP(a) is the penalty function, and b is its regularization parameter, X current represents the current image.
The method of claim 2, wherein three-layer wavelet decomposition is performed on the first PET image and the second PET image respectively.
A computer-readable storage medium having stored thereon a computer program, wherein the program, when executed by a processor, implements the steps of the method according to any one of claims 1 to 8.
A computer device, comprising a memory and a processor, a computer program that can be run on the processor is stored in the memory, and characterized in that, when the processor executes the program, any one of claims 1 to 8 is implemented The steps of the method described in item.