CEST image motion correction method and storage medium

Abstract

The invention discloses a CEST image motion correction method and a storage medium, and relates to the technical field of image processing. The method comprises the following steps: acquiring a CEST image and a structure image corresponding to the CEST image; obtaining a reference image according to the CEST image and the structure image; a transformation matrix between the reference image and a floating image is determined, and the floating image is a target frequency offset image of the CEST image; and performing registration on each frequency offset image of the CEST image according to the reference image and the transformation matrix to obtain the CEST image after motion correction. The method can improve the accuracy and reliability of a processing result, has the characteristics of high processing efficiency, wide application scene and the like, and is beneficial to popularization.

Description

Technical Field

[0001] The present invention relates to the technical field of image processing, and particularly relates to a motion correction method and a storage medium for CEST images. Background Art

[0002] CEST (Chemical Exchange Saturation Transfer) imaging is an important medical imaging technique for evaluating molecular characteristics in biological tissues. This technique obtains a Z-spectrum by detecting the change in water signal intensity caused by saturation transfer at different resonance frequencies, and the multi-solute molecular information quantified therefrom can provide assistance for clinical applications such as disease diagnosis and glycogen detection.

[0003] The analysis of CEST images generally requires accurate Z-spectrum data as a prerequisite. However, in the actual clinical image acquisition process, due to irresistible factors such as easy signal interference and long scanning time, the directly obtained Z-spectrum data is not ideal. Therefore, it is very necessary to perform a series of reasonable preprocessing on the acquired image data to improve the reliability of the quantitative results before data analysis.

[0004] Currently, the main challenges faced by CEST technology in the clinical acquisition registration process are as follows: (1) image deformation caused by magnetic field inhomogeneity or tissue movement; (2) in order to comprehensively understand the lesion situation, it is usually necessary to combine CEST images with other multi-modal imaging. However, the contrast mechanisms of CEST images and other imaging modalities are different, and there are significant differences in the gray-scale distribution and feature performance of the images; (3) in clinical practice, physiological movements such as the patient's breathing and heartbeat are inevitable, which will cause changes in the position and shape of CEST images at different time points.

[0005] Therefore, the following solutions have been proposed in related technologies: (1) Main magnetic field inhomogeneity correction: methods for correcting the water saturation frequency shift of the B0 field through water saturation shift referencing (WASSR) and correcting the B0 field by combining a two-pool model with Lorentz fitting; (2) Multi-modal registration: using image information (such as gray scale, edges, contours, etc.) or based on a spatial transformation model, and matching CEST images with other modal images through an optimization algorithm; (3) Motion correction: a robust principal component analysis (PRCA) motion correction method for image motion correction by decomposing the image data affected by motion into a low-rank part and a sparse part to separate the background and motion components.

[0006] However, most of these methods are only applicable to experimental subjects with simple structures such as the human brain. For body regions with relatively complex structures such as the liver, they have not been widely tested and applied. Therefore, developing a motion correction method applicable to CEST body magnetic resonance images is of great significance for promoting the clinical application of CEST technology. Summary of the Invention

[0007] The present invention aims to solve at least one of the technical problems in the related art to some extent. To this end, the object of the present invention is to provide a motion correction method and a storage medium for CEST images, so as to improve the accuracy and reliability of CEST image motion correction and facilitate popularization.

[0008] In a first aspect, an embodiment of the present invention provides a motion correction method for CEST images, including: acquiring a CEST image and a structural image corresponding to the CEST image; obtaining a reference image according to the CEST image and the structural image; determining a transformation matrix between the reference image and a floating image, where the floating image is a target frequency offset image of the CEST image; registering each frequency offset image of the CEST image according to the reference image and the transformation matrix to obtain a motion-corrected CEST image.

[0009] According to an embodiment of the present invention, the obtaining a reference image according to the CEST image and the structural image includes: determining a mapping relationship between the levels of the CEST image and the levels of the structural image; performing alignment processing on the structural image and the CEST image based on the mapping relationship using a gradient descent optimization algorithm; downsampling the aligned structural image to a size corresponding to the CEST image using a bicubic interpolation algorithm; and fusing the target frequency offset image with the downsampled structural image to obtain the reference image.

[0010] According to an embodiment of the present invention, the fusing the target frequency offset image with the downsampled structural image includes: determining a target layer image corresponding to the target frequency offset image in the downsampled structural image according to the mapping relationship; and performing pixel-by-pixel linear weighting of the pixel values of the target frequency offset image and the target layer image.

[0011] According to an embodiment of the present invention, the determining a transformation matrix between the reference image and the floating image includes: in the (n + 1)-th iteration, updating the coefficients of the B-spline control knots through the following formula:

[0012]

[0013] where denotes the coefficient of the B-spline control node at the index position (i, j, k) in the n-th iteration, where i, j, and k represent the index values in three dimensions. denotes the gradient information of the similarity with respect to the coefficient a n denotes the step size in the n-th iteration, H n denotes the approximate Hessian matrix calculated based on the historical gradient information and its corresponding step size; using perform a B-spline transformation on the floating image, and calculate the similarity between the reference image and the transformed floating image; if at least one of the following conditions holds: the similarity is less than a preset similarity threshold, n reaches a preset iteration number threshold, and the norm of g n is less than a preset convergence threshold, then take H n as the transformation matrix, otherwise continue to the next iteration.

[0014] According to an embodiment of the present invention, the similarity between the reference image and the transformed floating image is calculated by the following formula:

[0015]

[0016] where MI(X, Y) represents the similarity between the reference image X and the transformed floating image Y, p(x, y) represents the joint probability distribution function of X and Y, characterizing the probability of pixel pairs with gray value x in X and gray value y in Y, p(x) and p(y) respectively represent the marginal probability distribution functions of X and Y, and p(x, y) characterizes.

[0017] According to an embodiment of the present invention, p(x, y), p(x), and p(y) are calculated based on the downsampling points obtained by randomly downsampling the reference image X and the transformed floating image Y.

[0018] According to an embodiment of the present invention, the approximate Hessian matrix H n is calculated based on the gradient information and its corresponding step size in the most recent m iterations, where the value range of m is 1 - 7.

[0019] According to an embodiment of the present invention, the target frequency offset image is the image at the highest frequency offset of the signal value. Registering each frequency offset image of the CEST image according to the reference image and the transformation matrix includes: registering the target frequency offset image using the reference image and the transformation matrix; updating the reference image to the registered target frequency offset image, and updating the target frequency offset image to the previous adjacent frequency offset image of the target frequency offset image; returning to the step of registering the target frequency offset image using the reference image and the transformation matrix until all frequency offset images are registered.

[0020] According to an embodiment of the present invention, before obtaining the reference image based on the CEST image and the structural image, the method further includes: acquiring a magnetic field displacement map corresponding to the CEST image; and performing magnetic field inhomogeneity correction on the CEST image by using the magnetic field displacement map.

[0021] In a first aspect, an embodiment of the present invention provides a computer-readable storage medium, on which a computer program is stored. When the computer program is executed by a processor, the method described in the first aspect is implemented.

[0022] The motion correction method and storage medium for the CEST image according to the embodiment of the present invention first acquire the CEST image and the corresponding structural image of the CEST image, then obtain the reference image based on the CEST image and the structural image, then determine the transformation matrix between the reference image and the target frequency offset image of the CEST image, and finally register each frequency offset image of the CEST image according to the reference image and the transformation matrix to obtain the motion-corrected CEST image. Thereby, the accuracy and reliability of the processing result can be improved, and it has the characteristics of high processing efficiency and wide application scenarios, which is conducive to popularization.

[0023] Additional aspects and advantages of the present invention will be given in part in the following description, become apparent in part from the following description, or be learned through the practice of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0024] Figure 1 is a flowchart of the motion correction method for the CEST image according to the embodiment of the present invention;

[0025] Figure 2 is a flowchart of the preprocessing of the multimodal fusion image according to an embodiment of the present invention;

[0026] Figure 3 is a flowchart of obtaining the transformation matrix according to an embodiment of the present invention;

[0027] Figure 4 is a flowchart of the adjacent frequency offset registration according to an embodiment of the present invention;

[0028] Figure 5 is an intuitive display diagram of the registration coincidence degree according to an embodiment of the present invention;

[0029] Figure 6 is a comparison diagram of the pixel spectrograms before and after registration according to an embodiment of the present invention;

[0030] Figure 7 is a comparison diagram of the quantification of the magnetization transfer rate asymmetry before and after registration according to an embodiment of the present invention;

[0031] Figure 8 is a statistical comparison diagram of the quantification results according to an embodiment of the present invention;

[0032] Figure 9 It is a similarity comparison diagram of repeated acquisition experiments in an embodiment of the present invention;

[0033] Figure 10 It is a statistical chart of parameter selection experiments in an embodiment of the present invention. Detailed implementation manners

[0034] The embodiments of the present invention will be described in detail below. Examples of the embodiments are shown in the drawings, where the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described by referring to the drawings below are exemplary and are intended to explain the present invention, and should not be construed as limiting the present invention.

[0035] The motion correction method and storage medium of the CEST image according to the embodiment of the present invention will be described below with reference to the drawings.

[0036] Figure 1 It is a flowchart of the motion correction method of the CEST image according to the embodiment of the present invention.

[0037] As Figure 1 shown, the motion correction method of the CEST image includes:

[0038] S11, obtaining a CEST image and a structural image corresponding to the CEST image.

[0039] Taking the acquisition of 3D liver magnetic resonance image data as an example, the subject needs to fast for 6-8 hours before magnetic resonance imaging (MRI) examination. During the scan, the subject lies in the supine position, and a 3.0T MR system is used in combination with a 16-channel phased array body coil for scanning. Cross-sectional imaging is used for routine liver MRI examination, including CEST images (3D-APT) and T1-weighted imaging (T1WI). 3D-APT imaging is acquired by the fast spin echo Dixon technique (Turbo Spin Echo Dixon technique, Dixon-TSE), and the specific parameters are as follows: repetition time (TR) = 2.8 milliseconds; echo time (TE) = 7.3 milliseconds; field of vision (FOV) = 300×225 mm²; matrix size = 80×80; slice thickness = 7 mm; number of slices = 41; saturation frequency offset = ±6, ±4.8, ±3.6, ±2.4, ±1.2, 0, and -1560 ppm; the imaging acquisition time is about 20 minutes. The subject received respiratory training before the scan and was required to hold their breath for several seconds when the machine made a sound. The T1w structural image is acquired by the fast field echo Dixon technique (Fast Field Echo Dixon technique, FFE-Dixon), and the specific parameters are as follows: TR = 3.6 milliseconds; TE = 1.1 milliseconds; FOV = 400×350 mm²; matrix size = 672×672; slice thickness = 5 mm; number of slices = 91.

[0040] In some embodiments of the present invention, the magnetic field displacement map corresponding to the CEST image can also be obtained, and the CEST image is corrected for magnetic field inhomogeneity by using the magnetic field displacement map.

[0041] Specifically, the magnetic field displacement map can be automatically generated by the Dixon method. The Dixon method mainly realizes water-fat separation based on the precession frequency difference between water and fat signals in the magnetic field. This method uses multiple echo times (TE) to collect data, and separates the signal components of water and fat by processing the signals collected at different TEs. In this process, the B0 field inhomogeneity will affect the signal phases of water and fat, so to a certain extent, the information of the B0 field can be indirectly reflected by the result of water-fat separation. This method is very useful for imaging tissues with a high fat content or when it is necessary to accurately distinguish water and fat signals. For example, in imaging of the abdomen, pelvis and other parts, the organ structure can be clearly displayed, and at the same time, the information after water-fat separation can be used to infer the situation of the B0 field. This B0 map (i.e., the magnetic field displacement map) can characterize the main offset and is used to guide the magnetic field inhomogeneity correction of CEST images to reduce the influence of the main magnetic field inhomogeneity on the quality and accuracy of CEST images.

[0042] S12. Obtain a reference image according to the CEST image and the structural image.

[0043] In some embodiments of the present invention, obtaining a reference image according to the CEST image and the structural image includes: determining the mapping relationship between the levels of the CEST image and the levels of the structural image; performing alignment processing on the structural image and the CEST image based on the mapping relationship using the gradient descent optimization algorithm; using the bicubic interpolation algorithm to downsample the aligned structural image to the size corresponding to the CEST image; fusing the target frequency offset image with the downsampled structural image to obtain a reference image.

[0044] Exemplarily, fusing the target frequency offset image with the downsampled structural image includes: determining the target layer image corresponding to the target frequency offset image in the downsampled structural image according to the mapping relationship; performing pixel-by-pixel linear weighting of the pixel values of the target frequency offset image and the target layer image.

[0045] Specifically, as Figure 2 shown, assuming the structural image ( Figure 2The T1w in it has M layers, and the CEST image has N layers. By combining the slice thickness and spatial starting positioning parameters of the two images, the mapping relationship between the M-layer structural image and the N-layer CEST image can be obtained. This mapping relationship can be used to determine the structural image layers corresponding to each layer of the CEST image. Then, according to the mapping relationship, the gradient descent optimization algorithm is used to register the structural image to the CEST image to ensure the spatial alignment of the corresponding layers of the structural image and the CEST image. Then, the aligned structural image is downsampled to the matrix size (i.e., the dimension) corresponding to the CEST image using the bicubic interpolation algorithm, and based on the corresponding layer image with ω% weight (i.e., the above-mentioned target layer image) combined with the target frequency offset image with (1 - ω%) weight (which can be the CEST image at the highest signal frequency offset, the frequency offset farthest from the water frequency offset, usually the frequency offset of the first acquisition), a pixel-by-pixel numerical linear combination is performed to obtain the fused reference image, where ω ranges from 0 to 50 and can be dynamically adjusted according to the image characteristics. Specifically, the reference image can be obtained through the following formula:

[0046] O(x,y) = ω% * I 1 (x,y) + (1 - ω%) * I 2 (x,y)

[0047] where O(x,y) represents the pixel value at the pixel position (x,y) in the reference image O, and I 1 (x,y) represents the pixel value at the pixel position (x,y) in the target layer image I 1 and I 2 (x,y) represents the pixel value at the pixel position (x,y) in the target frequency offset image.

[0048] Since the high-signal-value frequency offset CEST image has the best image quality at a specific frequency offset and can well reflect the chemical exchange characteristics of tissues, while the structural image clearly shows the anatomical structure of tissues, the reference image obtained by fusing the two can provide richer and more comprehensive image features for subsequent registration.

[0049] In some examples, the above-mentioned bicubic interpolation algorithm is as follows:

[0050] A cubic polynomial function is used to perform a weighted average on the surrounding 16 pixels. Let I be the original image, (x,y) be a position in the target image, and the value of I(x,y) is calculated by the following formula:

[0051]

[0052] where W represents the cubic interpolation kernel function, usually in the following form:

[0053]

[0054] Among them, a is usually taken as -0.5 or -0.75.

[0055] In some examples, the above gradient descent optimization algorithm is as follows:

[0056] Define the transformation matrix:

[0057]

[0058] Among them, θ represents the rotation angle.

[0059] Translation vector represents the translation amounts in the x and y directions. For a point on the floating image (target frequency offset image) (where x 1 is the coordinate point along the horizontal direction, and x 2 is the coordinate point along the vertical direction), the point x' after the rigid transformation can be expressed as:

[0060] x' = Rx + t

[0061] Define the loss function:

[0062] Let I R (x) represent the pixel value of the reference image at point x, and I F (x') represent the pixel value of the floating image at point x' after transformation. The loss function J(R, t) can be defined as:

[0063]

[0064] Among them, Ω represents the domain of the image, and N represents the number of pixels in the image.

[0065] Taking the two-dimensional case as an example, the gradient calculation for the translation vector t is as follows:

[0066]

[0067] Update the transformation parameters:

[0068]

[0069] Among them, k represents the number of iterations, and α is the learning rate, which is used to control the update step size.

[0070] S13. Determine the transformation matrix between the reference image and the floating image, where the floating image is the target frequency offset image of the CEST image.

[0071] In some embodiments of the present invention, determining the transformation matrix between the reference image and the floating image includes: in the (n + 1)-th iteration, updating the coefficients of the B-spline control knots through the following formula:

[0072]

[0073] Among them, represents the coefficient of the B-spline control node at the index position (i, j, k) in the nth iteration, where i, j, and k represent the index values in three dimensions. represents the gradient information of the similarity with respect to the coefficient a n represents the step size in the nth iteration, and H n represents an approximate Hessian matrix calculated based on historical gradient information (such as the gradient information of the most recent m iterations) and their corresponding step sizes. Optionally, m can take values in the range of 1 - 7 and can be dynamically adjusted according to the image characteristics. By approximating the Hessian matrix through the most recent m gradient information, it is not necessary to store the complete Hessian matrix, which can accelerate the optimization process while effectively utilizing historical gradient information.

[0074] After that, use to perform a B-spline transformation on the floating image and calculate the similarity between the reference image and the transformed floating image; if at least one of the following conditions is met: the similarity is less than a preset similarity threshold, n reaches a preset iteration number threshold, or the norm of g n is less than a preset convergence threshold, then take H n as the transformation matrix, otherwise continue the next iteration.

[0075] Among them, the similarity between the reference image and the transformed floating image (which can measure the degree of dependence between two random variables) can be calculated by the following formula:

[0076]

[0077] Among them, MI(X, Y) represents the similarity between the reference image X and the transformed floating image Y, p(x, y) represents the joint probability distribution function of X and Y, which characterizes the probability of pixel pairs with gray value x in X and gray value y in Y, and p(x) and p(y) respectively represent the marginal probability distribution functions of X and Y. p(x, y) is characterized. Exemplarily, p(x, y), p(x), and p(y) are calculated based on the downsampling points obtained by randomly downsampling the reference image X and the transformed floating image Y.

[0078] Specifically, as Figure 3 shown, when obtaining the transformation matrix, use the L - BFGS - 2 optimization algorithm of historical gradients to optimize the coefficient c ijk of the B-spline control node, and the optimization objective is to minimize the difference between the reference image I R (x, y, z) and the floating image I FThe similarity between (x′, y′, z′). When calculating the similarity, the downsampling point method obtained by ω′% (the value range of ω′ can be 0 - 100 and can be dynamically adjusted according to the image characteristics) random sampling can be adopted. For example, when the value of ω′ is 50, that is, half of the pixel points are randomly selected from the reference image and the floating image respectively. Thus, while maintaining the image feature information and ensuring the accuracy of the similarity, the calculation amount can be reduced, thereby improving the calculation efficiency. The specific process of obtaining the transformation matrix is as follows:

[0079] In each iteration, first, according to the coefficient c of the current B-spline control knot ijk , perform B-spline transformation on the floating image I F , and obtain the transformed image I F (x′, y′, z′). Through the linear combination of B-spline basis functions, the local transformation of the image is accurately simulated, providing the transformed image data for subsequent similarity calculation. Randomly sample ω′% from the reference image I R and the transformed floating image I F respectively. Use a random number generator to generate random indices within the image pixel index range and select the pixel points corresponding to the indices. Based on these sampled points, calculate the joint probability distribution function p(x, y) and the marginal probability distribution functions p(x), p(y), and then calculate the similarity. After that, the numerical differentiation or analytical differentiation method can be used to calculate the gradient information g ijk of the similarity (S) with respect to the coefficient c n . Store the calculated gradient information g n in the historical gradient queue. When the number of stored gradient information exceeds m, remove the earliest gradient information to ensure that the queue always stores the latest m historical gradient information, providing data support for the calculation of the approximate Hessian matrix. Then, according to the update formula of the L - BFGS - 2 algorithm, that is, use the current gradient g n , step size a n and the approximate Hessian matrix H n to update the coefficient c ijk . The step size a n can adopt a fixed step size or an adaptive step size strategy. The adaptive step size strategy can dynamically adjust the step size according to the change of gradient information and the objective function to accelerate convergence and avoid the algorithm falling into local optimum. The approximate Hessian matrix H n can be calculated according to the stored historical gradient information and the corresponding step size through a specific quasi - Newton update formula.

[0080] After each iteration, check whether the iteration termination condition is met. If the similarity is less than the preset similarity threshold, or the current iteration count reaches the preset iteration count threshold, or the gradient norm is less than the preset convergence threshold, it is considered that the algorithm converges, that is, the iteration termination condition is met. At this time, end the training and use the current approximate Hessian matrix as the transformation matrix to be obtained. If the iteration termination condition is not met, continue the next iteration until the iteration termination condition is met.

[0081] S14. Register each frequency-offset image of the CEST image according to the reference image and the transformation matrix to obtain the motion-corrected CEST image.

[0082] In some embodiments of the present invention, the target frequency-offset image is the image at the frequency offset with the highest signal value. Registering each frequency-offset image of the CEST image according to the reference image and the transformation matrix includes: registering the target frequency-offset image using the reference image and the transformation matrix; updating the reference image to the registered target frequency-offset image, and updating the target frequency-offset image to the adjacent frequency-offset image before the target frequency-offset image; returning to the step of registering the target frequency-offset image using the reference image and the transformation matrix until the registration of all frequency-offset images is completed.

[0083] Specifically, as Figure 4 shown, after obtaining the transformation matrix, using the target frequency-offset image (such as the image at 6 ppm in Figure 4 ) as the starting point, perform the registration operation of adjacent frequency offsets in sequence according to the frequency-offset order (that is, each frequency-offset image is first registered as a floating image to the previous adjacent frequency offset, and the registered image is used as the reference image for registering the next frequency offset until the image at the farthest frequency offset on the opposite side of the water frequency offset is registered). For example, referring to Figure 4 , first perform the registration at 6 ppm frequency offset, and then perform the registration at 4.8 ppm frequency offset based on the registration result at 6 ppm frequency offset, and so on until the registration at -6 ppm frequency offset is completed, ending the registration process. In this way, a high image alignment degree and the continuity of the CEST Z-spectrum can be achieved simultaneously.

[0084] The following uses experimental data to illustrate the advantages of the present invention compared with the comparative method. The experiments include intuitive qualitative comparison experiments, objective quantitative comparison experiments, repeatability acquisition similarity comparison experiments, and optimal parameter selection experiments. In the experiments, for the asymmetric quantization map, the difference before and after processing the asymmetric quantization map can be compared. The asymmetric quantization map at a specific frequency signal can be solved by the following formula to obtain the information of the saturated frequency solute:

[0085]

[0086] where S 0 represents the intensity of the water signal without the saturation pulse, and respectively represent the signal intensities at the distances from the reference water frequency under the action of the saturation pulse and .

[0087] Intuitive qualitative comparison experiment: Blur and rotation (20°) interference are artificially added to a CEST image to obtain the CEST image to be corrected, and the method of the present invention is used to correct it. The images before and after correction are compared as Figure 5 shown. The first row is a three-dimensional display, and the second row is a two-dimensional plane display. In each row, the first column is the coincidence display of the reference image and the floating image before motion correction, and the second column is the coincidence display after motion correction. It can be seen that for the artificially added blur and rotation interference, the present invention can accurately register to the reference image, which proves the stability of the present invention in real human data.

[0088] Objective quantitative comparison experiment: The gradient descent rigid registration method and the RPCA method are selected as the comparison methods. As Figure 6 shown, (a) is the CEST water map at 6 ppm, and (b) and (c) are the B0 map and the T1w structural image corresponding to the CEST water map at 6 ppm respectively. The method of the present invention and the two comparison methods are respectively used to correct the image in (a). The pixel spectrograms of all frequency offsets before and after correction are as Figure 6 shown in (d) in. In (d), the red frame represents the pixel spectrogram of the corresponding row at all frequency offsets after the data of all frequency offsets in the row where the yellow dotted line is located in (a) are gathered together. Ideally, the pixel spectra of the same row at each frequency offset should be aligned vertically along a straight line. In this experiment, motion is artificially added to several frequency offsets of the original data, and the effect is as shown without correction. It can be seen that compared with the comparison methods, the present invention can more accurately align the corresponding organs at each frequency offset, thereby ensuring the authenticity of the subsequent quantitative results.

[0089] Figure 6 The quantitative results corresponding to (d) in Figure 7 are as Figure 7 shown. The Magnetization Transfer Ratio Asymmetry (MTRasym) is selected as the standard for evaluating the quantitative results. If the alignment between each frequency offset is not achieved, it will lead to the subtraction result of the symmetric frequency offset Z-spectrum values with respect to the water frequency offset (~0 ppm) being either high or low, and there will be higher or lower abnormal signal regions in the quantitative map, as

[0090] Figure 6The comparison of the objective quantification results corresponding to (d) in China is as follows Figure 8 as shown, which shows two evaluation criteria as follows:

[0091] (1) Neighboring Structural Similarity Index (SSIM):

[0092] Using this index can measure the improvement of the Z-spectrum before and after correction, and reflect the final motion correction effect. Different from general medical image processing, the gray levels and contrasts of CEST frequency offset maps are different. Therefore, calculating the floating image and the reference image separately to measure the registration result cannot fully reflect the registration situation. Since the motionless Z-spectrum signal is theoretically continuous and smooth, this prior knowledge can be used to place the floating image in the context of the sequence map to more accurately evaluate its motion correction situation. The present invention uses the structural similarity index between the corrected saturation-weighted image and the average image at two adjacent offsets to evaluate the registration quality. For one of the frequency offset images S i , the registration effect of the current layer is measured by evaluating the similarity between each layer and its adjacent layer. For the evaluation of the registration quality of a single frequency offset, the adjacent SSIM is used to quantify the correction effect, and the formula is as follows:

[0093]

[0094] where S i represents the current frequency offset image, and S i+1 and S i-1 respectively represent a frequency offset image before and after S i in the acquisition sequence. When the motion correction effect is better, the difference between the current frequency offset image and the adjacent frequency offset image is smaller, and the adjacent SSIM value is larger, and vice versa. Using this index to calculate all sequence maps and taking the average value can reflect the overall final motion correction effect of the CEST sequence images. The formula is as follows:

[0095]

[0096] When the overall continuity and similarity of the sequence maps are greater, the value of meanSSIM (average SSIM) is greater, and the motion correction effect is considered better, and vice versa.

[0097] As Figure 8 shown in (a) in China, the average SSIM value of the uncorrected image data is about 0.8. After correction, the average SSIM value of the method of the present invention is significantly improved compared with the comparative method, reaching 0.9, which proves the effectiveness of the method of the present invention.

[0098] (2) CEST Quantification Map Uniformity Index (UI):

[0099] This metric can be used to measure the improvement of the CEST quantification maps before and after correction. The CEST quantification maps within the same tissue region (such as the liver region) of normal subjects can be regarded as relatively uniform, and the specific uniformity calculation formula is as follows:

[0100]

[0101] where N represents the total number of pixels in the quantification map, S i represents the value of the i-th quantification map point, represents the average value of all values, R represents the scaling factor (which can be adjusted according to the research scenario and can be set to 5), represents the sum of the absolute values of the differences between all data points and the average value (measuring the degree of difference between data points). When the UI value is larger, it indicates higher consistency or uniformity, and the motion correction effect is considered better, and vice versa.

[0102] As Figure 8 shown in (b) below, the UI value is calculated by selecting the MTRasym quantification map at 3.5 ppm and the average MTRasym quantification map from 0.5 - 1.5 ppm of the subject. It can be seen that the average UI values of the two quantification maps are around 65 - 70 before correction. After correction, the UI value of the method of the present invention is significantly improved compared with the comparative method, reaching around 75 - 80, which once again proves the effectiveness of the method of the present invention.

[0103] Repeatability acquisition similarity comparison experiment:

[0104] In this experiment, image data of the same liver part of the same subject were collected twice (half an hour apart), as Figure 9 shown in acquisition 1 in (a) and acquisition 2 in (b) below. Five identical regions of interest (ROIs) were selected respectively. After subtracting the Z-spectrum data of these regions before and after acquisition and taking the absolute value, the distribution of the differences was statistically analyzed and the average value was marked. Theoretically, the average Z-spectra of all ROIs for the two acquisitions are approximately the same and there is almost no difference. Since the gradient descent rigid registration and RPCA have little effect (p < 0.05), Figure 9 only the results before correction and after correction using the method of the present invention are shown in the figure below. As can be seen from the figure below, compared with before correction, the overall difference of the method of the present invention is reduced and the distribution is more concentrated. Figure 9 As can be seen from the figure below, compared with before correction, the overall difference of the method of the present invention is reduced and the distribution is more concentrated.

[0105] Parameter selection experiment:

[0106] Figure 10Figure (a) shows the different sampling intervals and the number of historical gradients of different L-BFGS-2 optimization backtracking used to calculate the mutual information (i.e., similarity), the corresponding mutual information values, and the total iteration duration. It can be seen that after the sampling interval reaches 1 / 2, the obtained mutual information values tend to be stable. Generally, the mutual information values obtained by backtracking 6 or 7 historical gradients are not much different, and are much larger than the model design with backtracking 5 historical gradients. However, the time consumed by backtracking 7 historical gradients is much longer than that of backtracking 6 historical gradients. At the same time, considering the computational efficiency, setting the sampling interval to the 1 / 2 region and backtracking 6 historical gradients is a better choice.

[0107] Figure 10 Figure (b) shows the mutual information values corresponding to different reference images of 5 subjects. It can be seen that due to the low-resolution drawback of CEST images, the mutual information values obtained using CEST images as reference images are smaller than those obtained using the fused images of CEST images and T1w anatomical images as reference images.

[0108] Thus, the rationality and scientificity of the parameter selection and design of the method of the present invention are proved.

[0109] In summary, for the motion correction method of CEST images in the embodiments of the present invention, using a multi-modal fused image as the reference image, a high-signal frequency-offset CEST image as the floating image, introducing the L-BFGS-2 optimization algorithm and B-spline transformation to iteratively obtain the transformation matrix, and using the transformation matrix to perform registration for each frequency offset in the order of adjacent frequency offsets, and using the image after each registration as the reference image for the next adjacent frequency offset, the following beneficial effects can be achieved:

[0110] 1) The accuracy and reliability of the processing results are improved;

[0111] 2) High efficiency in processing: In the entire motion correction process, each step comprehensively considers the computing time as an important metric. On the premise of ensuring the image processing effect, the algorithm with the highest time efficiency is preferentially selected to ensure the high efficiency of the entire process. In addition, the standardized processing process also reduces the errors caused by experimental differences, facilitating the effect evaluation and comparison of future new methods;

[0112] 3) Wide application scenarios: The present invention is not only applicable to the liver region, but can also be migrated to similar multi-organ tissue environments, such as the chest cavity, pelvic cavity, etc., to improve the data quality. At the same time, this method is also applicable to target bodies with relatively simple structures, such as human brains, rats, etc. The tests on different data sets show that the present invention performs excellently in terms of accuracy and robustness, which is conducive to its wide application and popularization in the CEST field.

[0113] Based on the motion correction method of CEST images in the above embodiments, the present invention proposes a computer-readable storage medium.

[0114] In this embodiment, a computer program is stored on a computer-readable storage medium. When the computer program is executed by a processor, the motion correction method for CEST images in the above embodiment is implemented.

[0115] It should be noted that the logic and / or steps represented in the flowchart or described in other ways herein, for example, can be considered as a definite sequence list of executable instructions for implementing logical functions, and can be specifically implemented in any computer-readable medium for use by an instruction execution system, apparatus, or device (such as a computer-based system, a system including a processor, or other systems that can fetch and execute instructions from the instruction execution system, apparatus, or device), or in combination with these instruction execution systems, apparatus, or devices. For the purposes of this specification, a "computer-readable medium" can be any device that can contain, store, communicate, propagate, or transport a program for use by or in connection with an instruction execution system, apparatus, or device. More specific examples (non-exhaustive list) of computer-readable media include the following: an electrical connection portion with one or more wirings (electronic device), a portable computer diskette (magnetic device), a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or flash memory), an optical fiber device, and a portable compact disc read-only memory (CDROM). Additionally, the computer-readable medium can even be paper or other suitable media on which the program can be printed, because the program can be obtained electronically, for example, by optically scanning the paper or other media, followed by editing, interpretation, or otherwise processing as appropriate, and then storing it in a computer memory.

[0116] It should be understood that various parts of the present invention can be implemented by hardware, software, firmware, or a combination thereof. In the above embodiment, multiple steps or methods can be implemented by software or firmware stored in a memory and executed by a suitable instruction execution system. For example, if implemented in hardware, as in another embodiment, any one or a combination of the following techniques well known in the art can be used: discrete logic circuits having logic gate circuits for implementing logical functions on data signals, application specific integrated circuits having appropriate combinational logic gate circuits, programmable gate arrays (PGAs), field programmable gate arrays (FPGAs), etc.

[0117] In the description of this specification, the description with reference to terms such as "one embodiment", "some embodiments", "example", "specific example", or "some examples" means that the specific features, structures, materials, or characteristics described in connection with the embodiment or example are included in at least one embodiment or example of the present invention. In this specification, the schematic representations of the above terms do not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials, or characteristics described may be combined in any one or more embodiments or examples in a suitable manner.

[0118] In the description of the present invention, it should be understood that the orientation or positional relationships indicated by terms such as "center", "longitudinal", "transverse", "length", "width", "thickness", "upper", "lower", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer", "clockwise", "counterclockwise", "axial", "radial", "circumferential", etc. are based on the orientation or positional relationships shown in the drawings, and are only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and thus should not be construed as a limitation of the present invention.

[0119] In addition, the terms "first" and "second" are only used for descriptive purposes and should not be construed as indicating or implying relative importance or implicitly specifying the quantity of the indicated technical features. Thus, features defined with "first" and "second" may explicitly or implicitly include at least one of such features. In the description of the present invention, the meaning of "a plurality" is at least two, such as two, three, etc., unless otherwise specifically and clearly defined.

[0120] In the present invention, unless otherwise clearly specified and limited, terms such as "mounted", "connected", "connected to", "fixed" and the like should be understood in a broad sense. For example, it may be a fixed connection, a detachable connection, or integrated; it may be a mechanical connection or an electrical connection; it may be directly connected or indirectly connected through an intermediate medium, and it may be the internal communication of two elements or the interaction relationship between two elements, unless otherwise clearly limited. For those of ordinary skill in the art, the specific meanings of the above terms in the present invention can be understood according to specific circumstances.

[0121] In the present invention, unless otherwise clearly specified or limited, the first feature being "on" or "under" the second feature may mean that the first and second features are in direct contact, or the first and second features are indirectly in contact through an intermediate medium. Moreover, the first feature being "above", "over" and "on top of" the second feature may mean that the first feature is directly above or obliquely above the second feature, or merely indicates that the horizontal height of the first feature is higher than that of the second feature. The first feature being "under", "beneath" and "underneath" the second feature may mean that the first feature is directly below or obliquely below the second feature, or merely indicates that the horizontal height of the first feature is less than that of the second feature.

[0122] Although the embodiments of the present invention have been shown and described above, it can be understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Those of ordinary skill in the art can make changes, modifications, substitutions and variations to the above embodiments within the scope of the present invention.

Claims

1. A motion correction method for a CEST image, characterized in that: include: Acquire a CEST image and a structural image corresponding to the CEST image; Obtaining a reference image according to the CEST image and the structural image; Determine a transformation matrix between the reference image and a floating image, wherein the floating image is a target frequency offset image of the CEST image; The frequency-offset images of the CEST image are registered according to the reference image and the transformation matrix to obtain a motion-corrected CEST image.

2. The method according to claim 1, characterized in that: The step of obtaining a reference image according to the CEST image and the structural image comprises: Determining a mapping relationship between the levels of the CEST image and the levels of the structural image; Based on the mapping relationship, the structural image and the CEST image are aligned using a gradient descent optimization algorithm; Using a bicubic interpolation algorithm to downsample the aligned structural image to a size corresponding to the CEST image; The target frequency offset image is fused with the downsampled structural image to obtain the reference image.

3. The method according to claim 2, characterized in that The step of fusing the target frequency offset image with the downsampled structure image includes: Determine a target layer image corresponding to the target frequency offset image in the downsampled structure image according to the mapping relationship; The target frequency deviation image and the target layer image are linearly weighted pixel by pixel.

4. The method according to claim 1, characterized in that: The determining of the transformation matrix between the reference image and the floating image comprises: In the n+1th iteration, the coefficients of the B-spline control nodes are updated by the following formula: in, Represents the coefficient of the B-spline control node at index position (i, j, k) in the nth iteration, where i, j, k represent the index values ​​in three dimensions. Represents the similarity with respect to the coefficient The gradient information of n represents the step size in the nth iteration, H n Represents the approximate Hessian matrix calculated based on historical gradient information and its corresponding step size; use Performing a B-spline transformation on the floating image, and calculating the similarity between the reference image and the transformed floating image; If the similarity is less than the preset similarity threshold, n reaches the preset iteration number threshold, g n If the norm of is less than at least one of the preset convergence thresholds, then H n As the transformation matrix, otherwise continue to the next iteration.

5. The method according to claim 4, characterized in that The similarity between the reference image and the transformed floating image is calculated by the following formula: Among them, MI(X,Y) represents the similarity between the reference image X and the transformed floating image Y, p(x,y) represents the joint probability distribution function of X and Y, representing the probability of a pixel pair with grayscale value x in X and grayscale value y in Y, p(x) and p(y) represent the edge probability distribution functions of X and Y respectively, and p(x,y) represents.

6. The method according to claim 5, characterized in that p(x,y), p(x) and p(y) are calculated based on the down-sampling points obtained by randomly down-sampling the reference image X and the transformed floating image Y.

7. The method according to claim 4, characterized in that The approximate Hessian matrix H n It is calculated based on the gradient information of the most recent m iterations and their corresponding step sizes, where the value of m ranges from 1 to 7.

8. The method according to claim 1, characterized in that: The target frequency deviation image is an image at the highest frequency deviation of the signal value, and the registering of the frequency deviation images of the CEST image according to the reference image and the transformation matrix includes: Registering the target frequency-shifted image using the reference image and the transformation matrix; Updating the reference image to the registered target frequency offset image, and updating the target frequency offset image to an adjacent frequency offset image before the target frequency offset image; Return to the step of registering the target frequency-offset image using the reference image and the transformation matrix until the registration of all frequency-offset images is completed.

9. The method according to any one of claims 1 to 8, characterized in that Before obtaining a reference image according to the CEST image and the structural image, the method further includes: Obtaining a magnetic field displacement map corresponding to the CEST image; The CEST image is corrected for magnetic field inhomogeneity using the magnetic field displacement map.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the method according to any one of claims 1 to 9 is implemented.

Images