CEST image preprocessing method and storage medium

Abstract

The invention discloses a CEST image preprocessing method and a storage medium, and relates to the technical field of image processing. The CEST image preprocessing method comprises the following steps: acquiring a CEST image and a structure image corresponding to the CEST image; performing coarse segmentation on the structure image, marking a coarse segmentation result, and performing registration on the CEST image by using the structure image; mapping the marking result and the CEST image based on the registration result; and obtaining a fine segmentation result and a signal quantification result of the CEST image according to a mapping result. According to the method, the accuracy and reliability of CEST image segmentation are improved, processing is efficient, the automation degree is high, meanwhile, standardization of the processing flow is achieved, and the method has wide applicability and popularization value.

Description

Technical Field

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

[0002] In the field of Magnetic Resonance Image (MRI), Chemical Exchange Saturation Transfer (CEST) imaging is an emerging imaging technology. Since CEST can be used to evaluate the molecular properties in biological tissues, it is an imaging technology with great potential. Different from traditional MRI that only images free hydrogen atoms, CEST applies different saturation pulses to these hydrogen atoms according to the different resonance frequencies of hydrogen atoms on different groups, and indirectly obtains information about different molecules in the solute by detecting the change in the intensity of the water signal caused by saturation transfer, which can improve the sensitivity of traditional magnetic resonance imaging and achieve accurate judgment of diseases.

[0003] However, CEST imaging also faces many challenges in clinical applications. For example, the low resolution during acquisition makes the clinical application of CEST imaging difficult, and traditional CEST analysis methods can only manually delineate the Region of Interest (ROI), which not only takes a long time but may also lead to errors caused by human subjective factors. In particular, in the rat brain, due to many fine brain regions and the disadvantage of low resolution during CEST acquisition, large errors will inevitably occur when manually delineating the ROI. Therefore, a method for automatically segmenting CEST images is needed to solve the above problems. Currently, many researchers have proposed methods based on the segmentation of structural images. Since structural images have high resolution and short acquisition time, many researchers use structural images for segmentation. For example: (1) A method for segmenting the cerebral cortex and subcortical regions of the rat brain based on the T2-FLAIR sequence; (2) A digital stereotaxic rat brain atlas based on the fine anatomical contours of the Paxinos space and its automated application; (3) A method for individualizing the segmentation of the rat brain based on structural images. However, most of the above methods are based on high-resolution structural images for segmentation, and the segmentation results are relatively rough. For example, they can only roughly segment gray matter, white matter, and cerebrospinal fluid based on structural images and cannot register them on images of other modalities, which greatly limits the application of segmentation methods. Summary of the Invention

[0004] 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 preprocessing method and a storage medium for CEST images, so as to improve the accuracy and reliability of CEST image segmentation and realize the standardization of the processing flow.

[0005] In a first aspect, an embodiment of the present invention provides a preprocessing method for CEST images, including: acquiring a CEST image and a structural image corresponding to the CEST image; performing rough segmentation on the structural image, marking the rough segmentation result, and registering the CEST image by using the structural image; mapping the marking result to the CEST image based on the registration result; and obtaining a fine segmentation result and a signal quantification result of the CEST image according to the mapping result.

[0006] According to an embodiment of the present invention, the performing rough segmentation on the structural image and marking the rough segmentation result includes: performing rough segmentation on the structural image by using a preset template to obtain a rough segmentation result including multiple regions, and respectively assigning a region label of the region to which each voxel in each region belongs.

[0007] According to an embodiment of the present invention, the registering the CEST image by using the structural image includes: determining a first target corresponding layer of the CEST image in the structural image; downsampling the first target corresponding layer image to the resolution of the CEST image to obtain a first target image; initializing a deformation matrix; transforming the CEST image by using the deformation matrix; calculating the similarity between the first target image and the transformed CEST image; judging whether the iteration terminates according to the similarity; if so, using the transformed CEST image as the registered CEST image, otherwise updating the deformation matrix according to the similarity and returning to the step of transforming the CEST image by using the deformation matrix.

[0008] According to an embodiment of the present invention, the similarity between the first target image and the transformed CEST image is calculated by the following formula:

[0009]

[0010] where MI(I, J) represents the similarity between the first target image I and the transformed CEST image J, p(i) and p(j) respectively represent the marginal probability distributions of the images I and J, represents the joint probability distribution function, H(i, j) represents the number of voxels with value i in the image I and voxels with value j in the image J, and L represents the resolution size of the joint distribution histogram of the images I and J.

[0011] According to an embodiment of the present invention, determining whether to terminate the iteration based on the similarity includes:

[0012] Calculating a first mutual information gradient through the following formula:

[0013]

[0014] where θ represents the transformation parameter corresponding to the deformation matrix, represents the partial derivative of p(i,j) with respect to θ, represents the partial derivative of log p(i,j) with respect to θ;

[0015] If the first mutual information gradient is less than a first preset gradient threshold, it is determined that the first round of iteration terminates.

[0016] According to an embodiment of the present invention, determining whether to terminate the iteration based on the similarity includes: constructing an objective function based on the similarity after the termination of the first round of iteration; calculating the Jacobian matrix of the objective function, and calculating an approximate Hessian matrix based on the Jacobian matrix; calculating a second mutual information gradient based on the approximate Hessian matrix and the Jacobian matrix; if the second mutual information gradient is less than a second preset gradient threshold, it is determined that the second round of iteration terminates.

[0017] According to an embodiment of the present invention, updating the deformation matrix based on the similarity includes: in the first round of iteration, obtaining the updated deformation matrix through the following formula:

[0018]

[0019] where T rigid represents the updated deformation matrix in the first round of iteration, represents the T rigid corresponding transformation parameter, θ represents the transformation parameter corresponding to the deformation matrix before update, and α represents the update step size of gradient descent;

[0020] In the second round of iteration, obtaining the updated deformation matrix through the following formula:

[0021]

[0022] where represents the updated deformation matrix in the second round of iteration, represents corresponding transformation parameter, θ represents the transformation parameter corresponding to the deformation matrix before update, Δθ = -(H) -1 J T r, representing the second mutual information gradient, H = JT J + λI, which represents the approximate Hessian matrix, represents the Jacobian matrix between image I and image J, represents the partial derivative of E(θ) with respect to the j-th component in the transformation matrix, λ represents the optimization factor, I represents the identity matrix, r represents the residual vector, and E(θ) represents the objective function.

[0023] According to an embodiment of the present invention, determine the second target corresponding layer of the CEST image in the structural image, where the second target corresponding layer includes the first target corresponding layer; perform weighted averaging on the labeling results of the structural image of the second target corresponding layer to obtain the first target labeling result; downsample the region corresponding to the first target labeling result to the resolution of the registered CEST image and perform voxel-by-voxel mapping with the registered CEST image to obtain the mapping result.

[0024] According to an embodiment of the present invention, obtaining the segmentation result and signal quantification result of the CEST image from the mapping result includes: obtaining a fine segmentation result including multiple CEST segmentation regions from the mapping result; for different CEST segmentation regions, perform magnetization transfer asymmetry quantification and / or Lorentzian difference parameter quantification to obtain the signal quantification result.

[0025] In a second aspect, an embodiment of the present invention provides a computer-readable storage medium, on which a computer program is stored, characterized in that when the computer program is executed by a processor, it implements the preprocessing method of the CEST image described in the first aspect embodiment above.

[0026] The preprocessing method and storage medium of the CEST image according to the embodiments of the present invention first obtain the CEST image and the corresponding structural image of the CEST image, then perform rough segmentation on the structural image, label the rough segmentation result, and register the CEST image using the structural image; then map the labeling result with the CEST image based on the registration result, and obtain the fine segmentation result and signal quantification result of the CEST image according to the mapping result. Thus, the method improves the accuracy and reliability of CEST image segmentation, and is efficient and highly automated in processing, and at the same time realizes the standardization of the processing process, and has wide applicability and promotion value.

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

[0028] Figure 1 is a flowchart of the preprocessing method of the CEST image according to the embodiment of the present invention;

[0029] Figure 2 It is the preprocessing flowchart of the CEST image of a specific embodiment of the present invention;

[0030] Figure 3 It is the flowchart of rough segmentation and fine labeling of brain regions in an embodiment of the present invention;

[0031] Figure 4 It is the specific flowchart of the gradient descent algorithm based on high-density mutual information metric in an embodiment of the present invention;

[0032] Figure 5 It is the registration result and Dice coefficient statistical chart of different MI-SS parameters in an embodiment of the present invention;

[0033] Figure 6 It is the specific flowchart of optimizing the transformation matrix parameters in an embodiment of the present invention;

[0034] Figure 7 It is the flowchart of brain region segmentation of the CEST image of a rat brain in an embodiment of the present invention;

[0035] Figure 8 It is the segmentation result diagram of different brain regions of a rat brain in an embodiment of the present invention;

[0036] Figure 9 It is the quantitative result diagram of the hippocampal region of a rat brain at different frequency offsets of the CEST signal in an embodiment of the present invention;

[0037] Figure 10 It is the comparison result diagram between normal rats and AD rats in an embodiment of the present invention. Detailed implementation manners

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

[0039] The preprocessing method and storage medium of the CEST image according to the embodiments of the present invention will be described below with reference to the accompanying drawings.

[0040] Figure 1 It is the flowchart of the preprocessing method of the CEST image according to the embodiment of the present invention.

[0041] As Figure 1 shown, the preprocessing method of the CEST image includes:

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

[0043] Taking the acquisition of CEST images of the rat brain and their corresponding anatomical images as an example, when acquiring data, a 9.4T nuclear magnetic resonance scanner (Bruker) was used to perform magnetic resonance imaging on the rat brain. The scanning sequences were the water saturation imaging (WaterSaturation Shift Referencing, WASSR) sequence, the CEST sequence (RARE), and the T2 sequence (TurboRARE). The anatomical images were used for registering and segmenting the corresponding CEST images.

[0044] Among them, the CEST imaging parameters were set as follows: the saturation time was t sat = 2500 ms, the saturation field strength was 0.7 uT, 31 non-uniform samplings were performed within -10 ppm to 10 ppm, the repetition time (Repetition Time, TR) = 5500 ms, the echo time (Echo Time, TE) = 3.5 ms. Among them, the unsaturated image (S0) was obtained at a frequency offset of -200 ppm, the slice thickness was 1.5 mm, and the acquired matrix size was 72×96. The WASSR sequence was used for magnetic field correction, and 21 frequency offsets within the range of ±1 ppm were collected. The acquisition parameters and geometry were the same as those of other CEST experiments. The imaging parameter settings of the anatomical image sequence were as follows: the saturation time was t sat = 2500 ms, the repetition time TR = 4500 ms, the echo time TE = 36 ms, the flip angle = 180°, the slice thickness was 0.5 mm, and a total of 41 slices were acquired.

[0045] In some embodiments of the present invention, a magnetic field displacement map corresponding to the CEST image can also be obtained, and the magnetic field inhomogeneity correction of the CEST image can be performed using the magnetic field displacement map.

[0046] Specifically, the magnetic field displacement map is generated by calculating the frequency difference between the reference image and the water saturation image of the CEST sequence. Among them, the reference image refers to the CEST image acquired without applying the water saturation pulse (i.e., the unsaturated image S0), and the water saturation image is the CEST image acquired after applying the water saturation pulse to saturate water molecules. During the acquisition of the CEST image, water molecules have a specific resonance frequency, called the water signal frequency. If the main magnetic field is inhomogeneous, it will cause the water signal frequency to shift. Therefore, the WASSR technique can be used to estimate the shift of the B0 field by acquiring multiple images with different water signal frequencies and using the signals of these images. In actual operation, two types of images need to be acquired: one is the reference image, which is acquired without applying the water saturation pulse and serves as the benchmark for the water signal frequency; the other is the water saturation image, which is acquired after applying the water saturation pulse to saturate water molecules. By calculating the signal difference between these two types of images, the change in the water signal frequency can be estimated, thereby generating the magnetic field displacement map (B0 map). The B0 map can be used to characterize the shift amount of the main magnetic field and provide a basis for correcting the magnetic field inhomogeneity of the CEST image. Correcting the magnetic field inhomogeneity of the CEST image according to the B0 map can reduce the influence of the main magnetic field inhomogeneity on the quality and quantitative accuracy of the CEST image.

[0047] S12, coarsely segment the structural image, label the result of the coarse segmentation, and register the CEST image using the structural image.

[0048] In some embodiments of the present invention, coarsely segmenting the structural image and labeling the result of the coarse segmentation include: coarsely segmenting the structural image using a preset template to obtain a coarse segmentation result containing multiple regions, and assigning a region label of the region to which each voxel in each region belongs.

[0049] Taking the preprocessing of the CEST image of the rat brain region as an example, as Figure 2 、 Figure 3 shown, using an existing segmentation program (such as the SPM (Statistical Parametric Mapping) toolbox, which includes a developed standard brain template) to preliminarily segment the rat brain region in the structural image, the coarse segmentation results of the main tissues such as gray matter, white matter, and cerebrospinal fluid in the rat structural image can be obtained. At the same time, a transformation matrix (non-rigid transformation) can be generated during the segmentation process, and this transformation matrix can realize the spatial registration of the structural image and the standard brain template. Subsequently, taking the coarse segmentation result and the developed rat brain template as prior information, combining the anatomical features of different brain regions in the template, using an automated algorithm to traverse each voxel in the three-dimensional structural image, and assigning the brain region label to which each voxel belongs, finally 131 fine brain regions can be obtained, providing a basis for subsequent analysis.

[0050] Specifically, referring to Figure 3 , assume that the resolution of the structural image collected is [hx, hy, hs] (the size of each layer is hx * hy, and there are hs layers in total), and the resolution of the developed rat brain template is [140, 90, 98] (the size of each layer is 140 * 90, and there are 98 layers in total). After rough segmentation using the developed SPM toolbox, rough segmentation results of main tissues such as gray matter, white matter, and cerebrospinal fluid are generated (with the same resolution of [140, 90, 98]). Subsequently, the rough segmentation results are automatically labeled in combination with the developed brain template, and a total of 131 precise brain regions with a resolution of [140, 90, 98] are generated.

[0051] In some embodiments of the present invention, registering the CEST image using the structural image includes: determining a first target corresponding layer of the CEST image in the structure; downsampling the first target corresponding layer image to the resolution of the CEST image to obtain a first target image; initializing a deformation matrix; transforming the CEST image using the deformation matrix; calculating the similarity between the first target image and the transformed CEST image; determining whether to terminate the iteration according to the similarity; if so, taking the transformed CEST image as the registered CEST image, otherwise updating the deformation matrix according to the similarity and returning to the step of transforming the CEST image using the deformation matrix.

[0052] Specifically, an interlayer selection algorithm can be adopted to determine the first target corresponding layer according to the position parameters collected by the CEST image and the position parameters collected by the structural image. The first target corresponding layer image is the layer image closest to the CEST image. Then, the bilinear interpolation downsampling method or the bicubic interpolation downsampling method can be used to downsample each corresponding layer of the high-resolution structural image to the resolution of the CEST image, so as to achieve alignment and consistency in spatial resolution between the structural image and the CEST image.

[0053] In some examples, taking the bilinear interpolation downsampling method as an example, for a high-resolution structural image I and a structural image with the target resolution (the resolution of the CEST image) I', determine the corresponding position (x, y) of the voxel (x', y') in the target image in the high-resolution structural image. Due to the resolution difference, each voxel in the target image usually corresponds to a floating-point position in the high-resolution structural image. To estimate the value of this position, four adjacent pixel points are selected as references, and these four points form a rectangular area:

[0054]

[0055]

[0056]

[0057]

[0058] Calculate the corresponding weights according to the distances from the voxel (x’, y’) in the target image to these four adjacent points:

[0059] The weights in the horizontal direction are: , 1 - w x

[0060] The weights in the vertical direction are: , 1 - w y

[0061] These weights reflect the spatial relationship between the voxel in the target image and its adjacent points. The closer the points are, the higher their weights, thus more significantly affecting the estimated value of the voxel in the target image. Subsequently, use the horizontal direction weight w x and the vertical direction weight w y to perform weighted summation on the four adjacent points within the rectangular region. Specifically, first perform weighted interpolation on two adjacent points along the horizontal direction, and then perform secondary weighted interpolation on the result of the horizontal interpolation along the vertical direction to finally obtain the estimated value (x’, y’) of the voxel in the target image in the structural image:

[0062] I′(x′,y′)=(1 - w x )(1 - w y )Q 11 +w x (1 - w y )Q 21 +(1 - w x )w y Q 12 +w x w y Q 22

[0063] In the high - resolution structural image, the fine segmentation results of specific brain regions have been obtained through rough segmentation and labeling. Since the resolution of CEST images is relatively low, direct point - to - point mapping may lead to loss of segmentation information or error accumulation. Therefore, a method that can maintain the integrity of spatial information during the resolution conversion process is required. The bilinear interpolation downsampling method can effectively smooth the detailed information in the high - resolution image while ensuring that each pixel in the low - resolution image can reflect the local features of the original image. Therefore, using bilinear interpolation can not only quickly complete the downsampling process but also maximize the retention of the spatial accuracy of the brain region segmentation results when the resolution is reduced.

[0064] In some other examples, bicubic interpolation downsampling is taken as an example. First, for the voxels in the target image, the voxel data of two rows above and below it are selected, and cubic polynomial interpolation calculations are respectively performed along the horizontal direction. The expression form of the cubic polynomial is:

[0065] P(x) = ax 3 + bx 2 + cx + d

[0066] where a, b, c, and d are the coefficients of the polynomial, determined by the values of the surrounding known data points and their derivatives (or differences), and x represents the relative position of the target interpolation point between adjacent voxels. After completing the interpolation in the horizontal direction, based on the obtained interpolation result, cubic polynomial interpolation is applied again along the vertical direction to calculate the final value of the target voxel. This process can ensure that the downsampled image after interpolation reaches a high precision in terms of spatial resolution and gray continuity.

[0067] In some embodiments of the present invention, the similarity between the first target image and the transformed CEST image is calculated by the following formula:

[0068]

[0069] where MI(I, J) represents the similarity between the first target image I and the transformed CEST image J, p(i) and p(j) respectively represent the marginal probability distributions of images I and J, represents the joint probability distribution function, H(i, j) represents the number of voxels with value i in image I and voxels with value j in image J, and L represents the resolution size of the joint distribution histogram of images I and J.

[0070] Specifically, to achieve the precise registration of the structural image and the CEST image after resolution alignment, a two-dimensional array H of size L * L is created to construct the joint distribution histogram of the two images. Each element H(i, j) in the array represents the number of pairs of voxels with value i in image I and voxels with value j in image J. By traversing the two images pixel by pixel, each pair of pixel pairs that meet the above conditions (i.e., J(x, y) from the CEST image and I(x, y) from the first target image) is recorded in the array H. Subsequently, H is normalized and converted into the above joint probability distribution p(i, j). Based on the normalized joint probability distribution p(i, j), the mutual information between the two images is further calculated by the above formula MI(I, J) to obtain the similarity. Based on this similarity, the registration process can be further optimized.

[0071] In some embodiments of the present invention, it is determined whether to terminate the iteration according to the similarity, including: calculating the first mutual information gradient by the following formula:

[0072]

[0073] Among them, θ represents the transformation parameter corresponding to the deformation matrix, including the rotation angle θ r , the translation vector t x and t y , represents the partial derivative of p(i, j) with respect to θ, represents the partial derivative of log p(i, j) with respect to θ; if the first mutual information gradient is less than the first preset gradient threshold, it is determined that the first round of iteration terminates.

[0074] Furthermore, during the first round of iteration, the updated deformation matrix is obtained through the following formula:

[0075]

[0076] Among them, T rigid represents the updated deformation matrix during the first round of iteration, represents the T rigid corresponding transformation parameter, θ represents the transformation parameter corresponding to the deformation matrix before update, and α represents the update step size of gradient descent, which is used to control the amplitude of parameter update.

[0077] In some embodiments of the present invention, after the termination of the first round of iteration, the transformation parameters can be further optimized to improve the registration accuracy, and at this time, the second round of iteration is started. During the second round of iteration, it is judged whether to terminate the iteration according to the similarity, including: constructing an objective function based on the similarity after the termination of the first round of iteration; calculating the Jacobian matrix of the objective function, and calculating the approximate Hessian matrix according to the Jacobian matrix; calculating the second mutual information gradient according to the approximate Hessian matrix and the Jacobian matrix; if the second mutual information gradient is less than the preset gradient threshold, it is determined that the second round of iteration terminates.

[0078] Furthermore, in the second round of iteration, the updated deformation matrix is obtained through the following formula:

[0079]

[0080] Among them, represents the updated deformation matrix during the second round of iteration, represents corresponding transformation parameter, θ represents the transformation parameter corresponding to the deformation matrix before update, Δθ = -(H) -1 J T r, represents the second mutual information gradient, H = J T J + λI, represents the approximate Hessian matrix, represents the Jacobian matrix between image I and image J, denotes the partial derivative of E(θ) with respect to the j-th component in the transformation matrix, λ denotes the optimization factor (by setting different λ values, a dynamic balance can be found between the gradient descent algorithm and the Gauss-Newton method), I denotes the identity matrix, r denotes the residual vector, and E(θ) denotes the objective function.

[0081] Specifically, in the second-round iteration process, for the transformation parameter θ between the obtained CEST image and the first target image, an objective function E(θ) based on mutual information is defined as follows:

[0082] E(θ) = -MI(I, J′)

[0083] For mutual information gradient descent registration, the goal is to maximize the mutual information MI, while the goal of the LM (Levenberg-Marquardt) algorithm is to minimize the error function. Therefore, in the above formula, θ is the transformation parameter obtained by the high-density mutual information gradient descent algorithm and serves as the starting value of the optimization algorithm, J’ is the transformed CEST image obtained after preliminary registration (i.e., after the end of the first-round iteration), and MI(I, J’) represents the mutual information between the two images.

[0084] By calculating the contribution of geometric transformation information such as the rotation angle and translation distance of each voxel in the CEST image to the rate of change of the objective function, the Jacobian matrix J of the objective function is established (as above) to quantify the sensitivity of the objective function when the parameters change. The construction of the Jacobian matrix provides a mathematical basis for obtaining gradient information in the optimization process, thus supporting accurate parameter updates.

[0085] During the process of optimizing the parameters, an approximate calculation of the Hessian matrix is introduced to further improve the optimization accuracy through second-order information. Specifically, the approximate Hessian matrix is used to describe the local curvature characteristics of the objective function near the current parameter value, and can more accurately guide the iteration direction and step size at each voxel position. By combining the Jacobian matrix and the approximate Hessian matrix, the optimization algorithm can more efficiently approach the minimum point of the objective function, while significantly accelerating the convergence speed and avoiding being trapped in the dilemma of local extrema.

[0086] After that, the obtained not only contains information on the rotation angle and translation distance, but also comprehensively combines the optimization direction and step size adjustment of the current iteration to ensure the accuracy and stability of parameter updates. Based on the latest transformation parameters, a new transformation matrix is calculated

[0087] Finally, use the updated transformation matrix Perform a rigid registration operation on the CEST image to achieve its precise alignment with the high-resolution T2 structural image. Through this method, the spatial deviation between the CEST image and the structural image can be effectively corrected, ensuring the spatial consistency of the two types of images during subsequent analysis, thereby providing a reliable alignment basis for brain segmentation operations and quantitative analysis. The flowchart for optimizing the transformation matrix using the LM algorithm is as Figure 6 shown.

[0088] Thus, by iteratively updating the parameter vector, the sum of the squares of the errors between the predicted value and the observed value is minimized, thereby optimizing the objective function. When the current solution is far from the optimal solution, the algorithm behaves like the gradient descent method and has strong convergence; while when the current solution is close to the optimal solution, the algorithm behaves like the Gauss-Newton method and can quickly approach the local optimal solution. During the registration process of the CEST image and the structural image, on the one hand, there is a significant correlation between the structural information and texture features of the two types of images (for example, the increase in signal intensity of the thalamus relative to other regions in the CEST image and the structural image, and the decrease in signal intensity of the corpus callosum relative to the cortex and hippocampal regions in the CEST image and the structural image), which provides good initial matching conditions for the optimization based on the objective function; on the other hand, since the initial transformation matrix from the CEST image to the T2 structural image has been obtained through the previous (i.e., the first round of iteration) calculation, this matrix can better describe the preliminary spatial relationship between the two types of images and avoid the problem of local error accumulation caused by the low resolution and relatively poor signal-to-noise ratio of the CEST image. Therefore, it provides reliable initial conditions for the efficient convergence of the LM algorithm. In addition, by adjusting the damping factor, it is possible to flexibly switch between the gradient descent method and the Gauss-Newton method, thereby dynamically balancing the convergence speed and the stability of the solution during the optimization process. Especially during the registration process of the CEST image and the T2 structural image, this flexibility can effectively cope with the optimization difficulties brought by local noise or texture discontinuities in the image.

[0089] To further improve the registration accuracy, a sample size parameter of mutual information measurement (Mutual Information-Sample Size, MI-SS) can be introduced into the registration algorithm. This parameter determines the number of voxel samples used to calculate the similarity between the CEST image and the first target image during the registration process. Especially in regions with rich texture changes in the image, by increasing the number of sampling points, the quality of the mutual information estimation can be significantly improved, thereby improving the accuracy and robustness of the registration. During the parameter calculation process, the numerical gradient of each parameter is calculated through the update formula, and when the gradient of the objective function is small enough (for example, the preset gradient threshold is set to Epsilon = 10 -5 ), it means that it has approached the local optimal point, and the iterative optimization is stopped. The specific registration process of the CEST image and the structural image is as Figure 4 shown.

[0090] Due to the significant structural information and texture feature similarities between the anatomical image and the CEST image, this similarity is mainly reflected in the boundaries of anatomical structures (for example, the differences between the corpus callosum and the hippocampal region, and between the corpus callosum and the sensory cortex are very obvious in both the CEST image and the anatomical image), signal intensities (for example, the signals of the hippocampal region and the thalamus are higher than those of other regions in both the anatomical image and the CEST image), and the local consistency of tissue distribution. Since these similarities can be reflected as a higher mutual information in the joint distribution of the images, mutual information can be used to evaluate the spatial alignment degree of the two images. In the present invention, to achieve the precise registration of the CEST image and the anatomical image, a gradient descent algorithm based on mutual information is adopted. By maximizing the mutual information value between the images, the registration process is optimized, which can make full use of the common features of the two types of images and achieve a registration effect with high robustness and high precision.

[0091] In addition, in the gradient descent algorithm based on mutual information measurement, a mutual information measurement sample size parameter MI-SS is introduced to optimize the quality of mutual information calculation in the registration process. Specifically, the accuracy of mutual information calculation depends to a large extent on the number of samples used to statistically calculate the joint distribution. When MI-SS increases, the coverage range of the sample points participating in the calculation is wider. Especially in regions with rich texture changes, the increased sampling points can more comprehensively capture the correlation between the images, reducing the influence of noise and uneven distribution on the estimation of mutual information; and a larger MI-SS can also smooth the joint probability distribution and weaken the influence of local extrema. In the present invention, based on the traditional gradient descent algorithm based on mutual information measurement, by introducing and increasing the value of MI-SS, the accuracy of mutual information estimation is significantly improved, thus overcoming the possible local maximum problem and the risk of falling into the local extremum trap in the traditional mutual information algorithm.

[0092] The unregistered CEST image, and the registered CEST images obtained by using the gradient descent algorithm based on different MI-SS values are as Figure 5 shown. See Figure 5, the red channel is the registered CEST image, and the green channel is the T2 image of the target (i.e., the target images of all corresponding layers mentioned above). It can be seen that as the MI-SS increases, the registration effect also improves. At the same time, the evaluation index Dice coefficient is also used. By constructing the binary regions of the pre-registration, post-registration, and target images, the size of the intersection of the two regions is judged to evaluate the accuracy of the algorithm. From the registration results of different MI-SS parameters, as the MI-SS increases, the Dice coefficient also increases, which proves the role of the MI-SS parameter in this algorithm. Moreover, when comparing the Dice coefficients of the registration results of different experimental MI-SS = 500 and MI-SS = 150 on four rats (HP004, HP005, HP012, HP013), significant differences were found (p < 0.5, marked with * in Figure 5 ).

[0093] S13, map the labeling result to the CEST image based on the registration result.

[0094] In some embodiments of the present invention, mapping the labeling result to the CEST image based on the registration result includes: determining the second target corresponding layer of the CEST image in the structural image, where the second target corresponding layer includes the first target corresponding layer; performing weighted averaging on the labeling results of the structural images of the second target corresponding layer to obtain the first target labeling result; downsampling the region corresponding to the first target labeling result to the resolution of the registered CEST image and performing voxel-by-voxel mapping with the registered CEST image to obtain the mapping result.

[0095] Specifically, similar to determining the first target corresponding layer, an interlayer selection algorithm can be used to determine the second target corresponding layer according to the position parameters of the CEST image acquisition and the position parameters of the structural image acquisition. The second target corresponding layer has several layers, such as three layers, including the first target corresponding layer. Then, perform weighted averaging on the labeling results of the structural images of the second target corresponding layer to obtain the first target labeling result; downsample (the above bilinear interpolation downsampling method, bicubic interpolation downsampling method, etc. can be used) the corresponding regions of each layer of the second target corresponding layer of the first target labeling result to the resolution of the registered CEST image and perform voxel-by-voxel mapping with the registered CEST image to obtain the mapping result. Taking the above-mentioned rat brain region segmentation as an example, using 131 fine brain regions segmented in the structural image, after performing weighted averaging on the labeling results and downsampling to the resolution of the registered CEST image, and performing one-to-one voxel-by-voxel mapping with the registered CEST image, fine multi-brain region segmentation of the CEST image can be achieved. The flowchart of the entire segmentation process is as Figure 7 shown.

[0096] S14. Obtain the segmentation result and signal quantification result of the CEST image according to the mapping result.

[0097] In some embodiments of the present invention, obtaining the segmentation result and signal quantification result of the CEST image according to the mapping result includes: obtaining a fine segmentation result including multiple CEST segmentation regions according to the mapping result;

[0098] For different CEST segmentation regions, perform magnetization transfer asymmetry quantification and / or Lorentzian difference parameter quantification to obtain the signal quantification result.

[0099] Specifically, for the segmentation result of the CEST-MRI signal, automatically perform multi-parameter CEST-MRI signal quantification of different brain regions. The following formula can be used to automatically perform magnetization transfer asymmetry quantification:

[0100]

[0101] where S sat is the saturated CEST image, S 0 is the intensity of the water signal without the saturation pulse. Negative offset and positive offset respectively represent the negative frequency offset and positive frequency offset with the same distance from the reference water frequency under the action of the saturation pulse;

[0102] Use the following formula to automatically perform Lorentzian difference parameter quantification:

[0103]

[0104] where A is the amplitude of the peak in the Z-spectrum, Γ is the full width at half maximum of the peak, and Δω 0 is the resonance frequency offset. Calculate the Lorentzian difference value of each brain region through this formula and give the statistical results of different parameter quantification signals.

[0105] In some embodiments of the present invention, according to different requirements, an automatic data analysis method for single-group data and batch data is proposed to improve the efficiency of data processing and the reliability of the results. The following takes the preprocessing of the CEST image of the rat brain region as an example for illustration:

[0106] For a single set of data, after image registration is completed, Z-spectrum files for each brain region of each experimental rat can be automatically generated. The Z-spectrum is obtained by averaging the signal values of all pixel points in a certain brain region, and can effectively reflect the overall signal characteristics of that brain region. However, during the CEST image acquisition process, some image layers may not completely cover all brain regions, resulting in the appearance of NaN (Not a Number) values in the Z-spectrum data of some brain regions. To solve this problem, a data cleaning mechanism is introduced during the Z-spectrum generation process, which automatically filters out the brain regions containing NaN values, ensuring that the finally generated Z-spectrum files only include the brain regions existing in the actually acquired layers. At the same time, different quantitative methods are used to quantify the CEST signals of each brain region.

[0107] For experimental data that needs to be processed in batches, an efficient batch processing mechanism is provided. Through this mechanism, the data of multiple experimental rats can be processed at one time, and the Z-spectrum and LD spectral lines of each brain region of each rat can be automatically generated after the processing is completed. These results are systematically stored in the root directory of each rat, facilitating subsequent rapid comparison and comprehensive analysis between groups. The batch processing function significantly reduces the workload of manual operations, while ensuring the consistency of the data analysis process between different rats, providing strong support for large-scale experimental research.

[0108] The beneficial effects of the preprocessing method of the CEST images according to the embodiments of the present invention are illustrated by the following experiments:

[0109] Six rats were selected as the research objects in the experiment, three of which were normal control group rats (numbered 002, 004, 005), and the other three were rats with Alzheimer's disease (AD) (numbered 010, 012, 013). A complete processing flow of multi-modal image registration, multi-brain region segmentation, and multi-signal quantitative analysis was performed on all rats. The segmentation results are as Figure 8 shown, where the grayscale image is the whole rat brain, and the colored part is the segmentation result of each brain region. It should be noted that only the segmentation results of six rats in three typical brain regions (for each rat, from left to right are the hippocampus, corpus callosum, and sensory cortex) are shown here.

[0110] In addition, based on the characteristics of CEST images at different frequency offsets, the Lorentzian difference method was used to quantitatively analyze the signals. Specifically, at five key frequency offsets (+3.1 ppm for Amine, +3.5 ppm for Amide Proton Transfer (APT), +2.1 ppm for GAmine, -3.5 ppm for Nuclear Overhauser Effect (NOE), +2.6 ppm for Phosphocreatine (pCr)), the signals in each brain region were processed respectively to generate the quantitative maps and quantitative values of each brain region at the five frequency offsets (as Figure 9 shown). These results clearly demonstrated the differences in metabolite contents among different brain regions. For example, in the LD map of Figure 9 , the signal value at the -3.5 ppm frequency offset was significantly higher than the quantitative maps at the other four frequency offsets, indicating a large deposition of macromolecular lipids at -3.5 ppm. These quantitative maps provided important basis for subsequent biological analysis.

[0111] Furthermore, the Z-spectrum curves of the hippocampal regions of rats in the control group and the AD group were compared experimentally. The comparison results are as Figure 10 shown. Figure 10 The left figure of Figure 10 is the Z-spectrum curves of rats in the AD group and the control group, where the dashed line is the Z-spectrum of the AD group and the solid line is the Z-spectrum of the control group. It can be observed from the figure that the signal of rats in the AD group increased significantly at +3.5 ppm, and the Amide Proton Transfer (APT) calculated by the MTRasym quantitative method also decreased significantly in the hippocampal region of AD rats (p < 0.05, marked with * in the right figure of

[0112] ). This phenomenon is highly consistent with the research conclusions in the existing literature, further verifying the practicability of the present invention. In summary, the preprocessing method of CEST images proposed by the present invention, which covers key steps such as image registration, brain region segmentation, and signal quantitative analysis, can achieve the following beneficial effects:

[0113] 1) The accuracy and reliability of the segmentation results are significantly improved: By combining the registration algorithm with high-density mutual information metric, the segmentation strategy guided by the standard template, and the bilinear interpolation downsampling method, the texture features and structural information between the two types of images are fully utilized to achieve high-precision geometric alignment, laying a solid foundation for subsequent analysis; at the same time, the subjectivity problem existing in the process of traditional manual delineation of the region of interest is overcome, significantly improving the accuracy and consistency of the brain region segmentation results and reducing the influence of human errors;

[0114] 2) High efficiency and high automation in processing: Throughout the implementation process, the execution efficiency has always been an important consideration. By means of automated image registration and segmentation algorithms, the traditional manual operation process has been greatly simplified. Compared with manually delineating the region of interest, this method significantly shortens the data processing time, can quickly complete the analysis of a large number of experimental data, thereby improving the research efficiency and being applicable to high-throughput experimental requirements;

[0115] 3) Standardization of the processing flow: Based on existing standard templates and segmentation tools, a set of standardized processing flows for multi-region segmentation and multi-parameter signal quantification of CEST images have been developed. Through experimental verification on data collected multiple times, this flow performs excellently in terms of accuracy and robustness, ensuring the stability and repeatability of experimental results. At the same time, the standardized processing method effectively reduces the errors caused by differences in experimental conditions, providing guarantee for the consistency of data analysis;

[0116] 4) Wide applicability and promotion value: The method of the present invention can not only efficiently solve the key problems in CEST image processing, but also provide technical support for the application of CEST technology in a wider range of fields. The standardized processing flow reduces the technical threshold of experimental operations, helping to promote the further application and popularization of CEST technology in neuroscience research and disease diagnosis.

[0117] Based on the preprocessing method of CEST images in the above embodiments, the present invention also proposes a computer-readable storage medium.

[0118] In the embodiments of the present invention, a computer program is stored on the computer-readable storage medium. When the computer program is executed by a processor, the preprocessing method of CEST images in the above embodiments is implemented.

[0119] It should be noted that the logic and / or steps represented in the flowchart or otherwise described 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, apparatuses, 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 combination with an instruction execution system, apparatus, or device. More specific examples (a non-exhaustive list) of the computer-readable medium 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 medium on which the program can be printed, because the program can be obtained electronically, for example, by optically scanning the paper or other medium, followed by editing, interpretation, or otherwise processing as appropriate, and then storing it in a computer memory.

[0120] 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-described embodiments, 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 by 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), and the like.

[0121] In the description of this specification, the description referring to terms such as "one embodiment", "some embodiments", "example", "specific example", or "some examples", etc. 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 can be combined in any one or more embodiments or examples in a suitable manner.

[0122] In the description of the present invention, it should be understood that the orientation or positional relationship indicated by the terms "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 relationship 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 therefore should not be construed as a limitation on the present invention.

[0123] 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 the features. In the description of the present invention, "a plurality of" means at least two, such as two, three, etc., unless otherwise specifically defined.

[0124] In the present invention, unless otherwise clearly defined and limited, the terms "mounted", "connected", "coupled", "fixed", etc. 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.

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

[0126] 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 a limitation on 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 CEST image preprocessing method, characterized in that: include: Acquire a CEST image and a structural image corresponding to the CEST image; Performing rough segmentation on the structural image, marking the rough segmentation result, and registering the CEST image using the structural image; Mapping the labeling result with the CEST image based on the registration result; The fine segmentation result and signal quantification result of the CEST image are obtained according to the mapping result.

2. The CEST image preprocessing method according to claim 1, characterized in that: The roughly segmenting the structure image and marking the rough segmentation result includes: The structural image is roughly segmented using a preset template to obtain a rough segmentation result including multiple regions, and a region label of the region to which the voxel belongs is respectively assigned to each voxel in each region.

3. The CEST image preprocessing method according to claim 1, characterized in that: The registering the CEST image by using the structural image includes: Determine a first target corresponding layer of the CEST image in the structural image; Downsampling the first target corresponding layer image to the resolution of the CEST image to obtain a first target image; Initialize the deformation matrix; transforming the CEST image using the deformation matrix; Calculating the similarity between the first target image and the transformed CEST image; Determine whether to terminate the iteration according to the similarity; If so, the transformed CEST image is used as the registered CEST image; otherwise, the deformation matrix is ​​updated according to the similarity, and the process returns to the step of transforming the CEST image using the deformation matrix.

4. The CEST image preprocessing method according to claim 3, characterized in that: The similarity between the first target image and the transformed CEST image is calculated by the following formula: Among them, MI(I,J) represents the similarity between the first target image I and the transformed CEST image J, p(i) and p(j) represent the marginal probability distribution of images I and J respectively. represents the joint probability distribution function, H(i,j) represents the number of voxels with value i in image I and the number of voxels with value j in image J, and L represents the resolution of the joint distribution histogram of images I and J.

5. The CEST image preprocessing method according to claim 4, characterized in that: The step of judging whether to terminate the iteration according to the similarity comprises: The first mutual information gradient is calculated by the following formula: Among them, θ represents the transformation parameter corresponding to the deformation matrix, represents the partial derivative of p(i,j) with respect to θ, represents the partial derivative of the logarithm p(i,j) with respect to θ; If the first mutual information gradient is less than a first preset gradient threshold, it is determined that the first round of iteration is terminated.

6. The CEST image preprocessing method according to claim 5, characterized in that: The step of judging whether to terminate the iteration according to the similarity comprises: Construct the objective function based on the similarity after the first round of iterations ends; Calculating the Jacobian matrix of the objective function, and calculating the approximate Hessian matrix according to the Jacobian matrix; Calculate a second mutual information gradient according to the approximate Hessian matrix and the Jacobian matrix; If the second mutual information gradient is less than a second preset gradient threshold, it is determined that the second round of iteration is terminated.

7. The CEST image preprocessing method according to claim 6, characterized in that: The updating of the deformation matrix according to the similarity comprises: In the first round of iteration, the updated deformation matrix is ​​obtained by the following formula: Among them, T rigid represents the updated deformation matrix during the first round of iterations, Indicates the T rigid The corresponding transformation parameters, θ represents the transformation parameters corresponding to the deformation matrix before updating, and α represents the update step size of gradient descent; In the second round of iteration, the updated deformation matrix is ​​obtained by the following formula: in, represents the updated deformation matrix during the second round of iteration, express The corresponding transformation parameters, θ represents the transformation parameters corresponding to the deformation matrix before updating, Δθ = -(H) -1 J T r, represents the second mutual information gradient, H = J T J+λI, represents the approximate Hessian matrix, represents the Jacobian matrix between image I and image J, represents the partial derivative of E(θ) with respect to the j-th component in the transformation matrix, λ represents the optimization factor, I represents the identity matrix, r represents the residual vector, and E(θ) represents the objective function.

8. The CEST image preprocessing method according to claim 3, characterized in that: Mapping the marking result with the CEST image based on the registration result includes: Determine a second target corresponding layer of the CEST image in the structural image, wherein the second target corresponding layer includes the first target corresponding layer; Taking a weighted average of the marking results of the structural image of the layer corresponding to the second target, to obtain a first target marking result; The area corresponding to the first target labeling result is downsampled to the resolution of the registered CEST image, and is voxel-by-voxel mapped with the registered CEST image to obtain the mapping result.

9. The CEST image preprocessing method according to claim 1, characterized in that: The step of obtaining the segmentation result and signal quantification result of the CEST image according to the mapping result includes: Obtaining a fine segmentation result including a plurality of CEST segmentation regions according to the mapping result; For different CEST segmentation regions, magnetization transfer asymmetry quantification and / or Lorentz difference parameter quantification are performed to obtain the signal quantification result.

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