Preprocessing method of cest images

Abstract

The application discloses a CEST image preprocessing method, comprising the following steps: acquiring a CEST sequence, acquiring structure item images corresponding to each CEST image in the CEST sequence, and determining a magnetic field displacement map corresponding to the CEST sequence; correcting the magnetic field non-uniformity of the CEST images in the CEST sequence based on the magnetic field displacement map by using a cubic spline interpolation method, and registering the corresponding CEST images based on each structure item image. The CEST image preprocessing method can improve the accuracy and reliability of the processing result, has the characteristics of high processing efficiency and wide application scenarios, and is beneficial to promotion.

Description

Technical Field

[0001] This invention relates to the field of image processing technology, and in particular to a preprocessing method for CEST images. Background Technology

[0002] CEST (Chemical Exchange Saturation Transfer) imaging is an important medical imaging technique for assessing molecular properties in biological tissues. This technique obtains the Z-spectrum by detecting changes in water signal intensity caused by saturation transfer at different resonance frequencies. Researchers can then quantify information about multiple solute molecules from this spectrum, ultimately providing support for clinical applications such as disease diagnosis and glycogen detection.

[0003] However, CEST imaging faces a number of challenges in clinical applications. Accurate Z-spectrum is crucial for obtaining valuable information; however, during data acquisition, human and environmental factors can cause data to deviate from ideal conditions. Therefore, performing a series of preprocessing steps on the acquired images is essential to improve the practical usability of the data before final analysis.

[0004] Currently, the uncertainties in CEST images mainly include the following aspects: (1) Motion artifacts: During the acquisition process, unpredictable motion of the subject can cause misalignment of tissues in the sequence images; (2) Z-spectrum shift: The inhomogeneity of the spatial magnetic field is ultimately reflected in the shift of the Z-spectrum, in which the main magnetic field (B0) affects the displacement of the Z-spectrum in the frequency direction, ultimately leading to frequency point distortion; (3) Misalignment of reference structure items: CEST images have the disadvantage of unclear tissue structures. Usually, the anatomical structure items of the registration are used for positioning. However, some body parts, such as the pelvic cavity, have complex tissue structures and large differences in grayscale and contrast with other modal images, which makes registration difficult.

[0005] To this end, the following solutions have been proposed in related technologies: (1) Motion correction scheme: The low-rank approximation of z-spectrum (LRAZ) and the two-stage registration scheme using robust principal component analysis plus principal component analysis (RPCA+PCA) both separate motion and contrast, obtain a reference image, and realize the pair registration of the sequence image and the reference image; (2) Principal magnetic field inhomogeneity correction scheme: Collect images near the water resonance frequency to display the offset of each voxel's Z spectrum, such as the WASSR (water saturation displacement) image, and then guide the Z spectrum to be corrected based on the image; (3) Image multimodal registration: such as feature point matching based on MI (mutual information), SSIM (structural similarity) and SIFT (scale invariant feature transformation).

[0006] However, while these methods have achieved good results on experimental subjects with simple structures such as the human brain, they have not yet been widely tested in more complex regions such as the pelvis. Different research teams may use different methods, leading to inconsistencies in processing results. Therefore, developing a specialized preprocessing workflow for body CEST images is of positive significance for the development of CEST imaging in clinical applications. Summary of the Invention

[0007] This invention aims to at least partially solve one of the technical problems in related technologies. Therefore, the purpose of this invention is to propose a preprocessing method for CEST images to improve the accuracy and reliability of CEST image processing results, achieve efficient processing, and facilitate widespread application and promotion.

[0008] To achieve the above objectives, this invention proposes a preprocessing method for CEST images, which involves obtaining a CEST sequence, obtaining the structural term image corresponding to each CEST image in the CEST sequence, and determining the magnetic field displacement map corresponding to the CEST sequence; using cubic spline interpolation based on the magnetic field displacement map to perform magnetic field non-uniformity correction on the CEST images in the CEST sequence, and performing structural term registration on the corresponding CEST images based on each structural term image.

[0009] In addition, the CEST image preprocessing method according to embodiments of the present invention may also have the following additional technical features:

[0010] According to one embodiment of the present invention, before performing main magnetic field inhomogeneity correction and structural term registration on the CEST images, the method further includes: performing motion correction on the CEST images in the CEST sequence based on Z-spectrum similarity.

[0011] According to an embodiment of the present invention, the motion correction of each CEST image in the CEST sequence based on Z-spectrum similarity includes: expanding each CEST image in the current CEST sequence into an m*n row and one column one-dimensional column vector, where m*n is the size of the CEST image; merging the one-dimensional column vectors in the column direction to obtain a Z-spectrum set matrix of all voxels in the CEST sequence, where the column direction of the Z-spectrum set matrix represents different acquisition frequencies and the row direction represents different voxels; using the entropy weighted topsis algorithm to obtain a Z-spectrum similarity map based on the Z-spectrum set matrix, and calculating the current sharpness of the Z-spectrum similarity map; if the current sharpness or the current iteration number satisfies a first preset condition, The current CEST sequence is then used as the motion-corrected CEST sequence. The first preset condition includes at least one of the following: the current sharpness is greater than or equal to a preset sharpness threshold, the difference between the current sharpness and the sharpness of the previous iteration is less than or equal to a first difference threshold, and the current iteration number is greater than or equal to a preset iteration number threshold. If the current sharpness and the current iteration number do not meet the preset conditions, the parameters of the current first deformation matrix are updated using the gradient descent method, and the updated first deformation matrix is ​​used to transform each CEST image in the current CEST sequence, proceeding to the step of expanding each CEST image in the current CEST sequence into an m*n row and column one-dimensional column vector.

[0012] According to one embodiment of the present invention, calculating the current sharpness of the Z-spectral similarity map includes: calculating the Z-spectral similarity using the Sobel operator. Figure X gradient G in direction x The gradients G in the (x,y) and Y directions y (x,y); The gradient value G(x,y) of each voxel point (x,y) in the Z-spectrum similarity map is obtained by the following formula:

[0013]

[0014] The current sharpness of the Z-spectrum similarity map is obtained by averaging the gradient values ​​of all voxel points.

[0015] According to an embodiment of the present invention, the first deformation matrix is ​​represented by the following formula:

[0016]

[0017] Where T1 represents the first deformation matrix, θ represents the rotation angle, and d x d represents displacement in the X-axis direction. y This indicates displacement along the Y-axis.

[0018] According to one embodiment of the present invention, the magnetic field displacement map is obtained by calculating the water signal frequency difference between the water saturation image and the reference image corresponding to the CEST sequence. The reference image is a CEST image acquired without a water saturation pulse in the CEST sequence, and the water saturation image is a CEST image acquired after water molecules are saturated with a water saturation pulse in the CEST sequence. The magnetic field non-uniformity correction of the CEST images in the CEST sequence based on the magnetic field displacement map using cubic spline interpolation includes: for each voxel X in the CEST sequence, for each pair of adjacent data points (ω) of voxel X... i y i ) and (ω i+1 y i+1 Construct a cubic polynomial S between ) i (w):

[0019] S i (ω)=a i (ω-ω i ) 3 +b i (ω-ω i ) 2 +c i (ω-ω i )+d i

[0020] Where i = 1, 2, ..., n-1, n represents the number of CEST images in the CEST sequence, a i b i c i d i The coefficients are represented by equation S. i (ω i )=y i S i (ω i+1 )=y i+1 S i (ω i+1 ) = S i+1 (ω i+1 ), S i (ω i+1 ) = S i+1 (ω i+1 S″1(ω1)=0, S n-1 (ω n ) = 0 is determined, ω i y represents the acquisition frequency of the i-th CEST image in the CEST sequence. i This indicates that the i-th CEST image in the CEST sequence is at frequency ω. i The grayscale value of voxel X; using the aforementioned cubic polynomial S i(ω) Obtain the interpolation sequence, and obtain the B0 field offset sequence based on the interpolation sequence and the magnetic field displacement map; perform Z-spectrum correction on the corresponding CEST image in the CEST sequence based on each B0 field offset in the B0 field offset sequence.

[0021] According to an embodiment of the present invention, the structural term registration of the corresponding CEST image based on each of the structural term images includes: performing coarse structural term registration of the CEST image and its corresponding structural term image based on normalized mutual information and the Powell optimization algorithm to obtain an aligned structural term map corresponding to the CEST image; obtaining a floating deformation field and a reference deformation field based on the CEST image and the aligned structural term map using a pre-trained ResNet network, obtaining a target deformation field based on the floating deformation field and the reference deformation field, and performing fine structural term registration of the CEST image using the target deformation field.

[0022] According to an embodiment of the present invention, the coarse registration of the CEST image and its corresponding structural item image based on normalized mutual information and the Powell optimization algorithm includes: transforming the structural item image corresponding to the CEST image using a second deformation matrix; calculating the current normalized mutual information between the CEST image and the transformed structural item image; if the current normalized mutual information satisfies a second preset condition, then using the transformed structural item image as the aligned structural item image of the CEST image, wherein the second preset condition includes: the current normalized mutual information is greater than or equal to a preset normalized mutual information threshold, and / or, the difference between the current normalized mutual information and the normalized mutual information of the previous iteration is less than a second difference threshold; if the current normalized mutual information does not satisfy the second preset condition, then optimizing the second deformation matrix using the Powell algorithm, and proceeding to the step of transforming the structural item image corresponding to the CEST image using the second deformation matrix.

[0023] According to one embodiment of the present invention, the second deformation matrix is ​​represented by the following formula:

[0024]

[0025] Where T2 represents the second deformation matrix, θ represents the rotation angle, and t x t represents displacement along the X-axis. y represents the displacement along the Y-axis, and s represents the scaling factor.

[0026] According to an embodiment of the present invention, the step of obtaining a floating deformation field and a reference deformation field based on the CEST image and the aligned structure term map using a pre-trained ResNet network, and obtaining a target deformation field based on the floating deformation field and the reference deformation field, includes: extracting features from the CEST image using a pre-trained ResNet50 network to obtain a first feature map; dividing the CEST image into a 2*2 grid to obtain (m / 2)*(n / 2) grids; and inputting the information of each voxel point on the first feature map into the corresponding grid. The floating deformation field is constructed by taking the midpoint of each cell as a control point; the alignment structure item map corresponding to the CEST image is extracted using a pre-trained ResNet50 network to obtain a second feature map, and the alignment structure item map is divided into 2*2 cells to obtain (m / 2)*(n / 2) cells; the information of each voxel point on the second feature map is input into the corresponding cell, and the midpoint of each cell is taken as a control point to construct the reference deformation field; for each control point D in the deformation field, multiple control points D around the control point D are selected. i A set of points to be matched is constructed, and the feature vector T(D) of the corresponding control point D is obtained on the CEST image. The feature vector T′(D) of each control point in the set of points to be matched is then obtained on the alignment structure item map. i ), calculate T(D) and T′(D) i The Euclidean distance between the control points D in the floating deformation field is determined, and the optimal matching point of the control point D in the reference deformation field is determined based on the Euclidean distance; the target deformation field is obtained based on all the optimal matching points.

[0027] The CEST image preprocessing method of this invention introduces cubic spline interpolation in the main magnetic field correction process; and applies ResNet neural network to extract feature descriptions in the structural term registration stage, thereby improving the accuracy and reliability of the processing results. It also features high processing efficiency and wide application scenarios, making it easy to promote.

[0028] Additional aspects and advantages of the invention will be set forth in part in the description which follows, and in part will be obvious from the description, or may be learned by practice of the invention. Attached Figure Description

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

[0030] Figure 2 This is a flowchart of the CEST image main magnetic field non-uniformity correction according to an embodiment of the present invention;

[0031] Figure 3This is a flowchart of coarse registration of CEST image structure items according to an embodiment of the present invention;

[0032] Figure 4 This is a flowchart of an embodiment of the present invention for obtaining the feature map and deformation field corresponding to the input image;

[0033] Figure 5 This is a flowchart of a structural term fine registration based on a deformation field according to an embodiment of the present invention;

[0034] Figure 6 This is a flowchart of CEST image preprocessing according to an embodiment of the present invention;

[0035] Figure 7 This is a flowchart of a CEST image motion correction based on Z-spectrum similarity according to an embodiment of the present invention;

[0036] Figure 8 This is a flowchart of calculating the Z-spectral similarity map according to an embodiment of the present invention;

[0037] Figure 9 This is an example diagram showing the motion artifact correction and structure registration results according to an embodiment of the present invention;

[0038] Figure 10 This is an example diagram showing the results of main magnetic field non-uniformity correction and structural term registration according to an embodiment of the present invention. Detailed Implementation

[0039] Embodiments of the present invention are described in detail below, examples of which are illustrated in the accompanying drawings, wherein the same or similar reference numerals denote 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 intended to explain the present invention, and should not be construed as limiting the present invention.

[0040] The preprocessing method for CEST images according to an embodiment of the present invention is described below with reference to the accompanying drawings.

[0041] Figure 1 This is a flowchart of the CEST image preprocessing method according to an embodiment of the present invention.

[0042] like Figure 1 As shown, the preprocessing methods for CEST images include:

[0043] S1, obtain the CEST sequence, obtain the structural item image corresponding to each CEST image in the CEST sequence, and determine the magnetic field displacement map corresponding to the CEST sequence.

[0044] Taking the acquisition of pelvic magnetic resonance imaging (MRI) data as an example, a 3T MRI scanner (Ingenia, Philips Healthcare) can be used to perform MRI on the pelvis of patients or volunteers. This scanner has a 32-channel phase array coil. The scanning sequence is a three-dimensional amide proton transfer (CEST) sequence and a T2-weighted sequence, and the readout sequence is a fast spin echo sequence, resulting in CEST and T2-weighted sequences. The CEST images in the CEST sequence correspond one-to-one with the structural term images in the T2-weighted sequence. The structural term images are used to perform structured correction on the corresponding CEST images.

[0045] The CEST imaging parameters can be set as follows: saturation time is 2.5s, saturation field strength is 0.7ut, 30 frequency offsets are acquired within the range of -10ppm to 10ppm, repetition time (TR) is 5600ms, and echo time (TE) is 6.85m. The unsaturated image (S0) is obtained at a frequency offset of -1560ppm. Water saturation imaging (WASSR) is used for magnetic field correction. The T2 weighted imaging parameters are as follows: flip angle (Flip Angle) = 90°, TR = 3200ms, and TE = 80ms.

[0046] In some embodiments, the magnetic field displacement map is obtained by calculating the water signal frequency difference between the water saturation image and the reference image corresponding to the CEST sequence. The reference image is the CEST image acquired when there is no water saturation pulse in the CEST sequence, and the water saturation image is the CEST image acquired after water molecules are saturated by applying a water saturation pulse in the CEST sequence.

[0047] Specifically, in CEST image acquisition, water molecules have a specific resonant frequency, called the water signal frequency. This frequency changes under conditions of inhomogeneity in the main magnetic field. The basic principle of WASSR is to acquire multiple images with different water signal frequencies and estimate the B0 field shift based on the signals in these images. In practice, two types of images need to be acquired: a reference image, acquired without a water saturation pulse, used as a benchmark for estimating the water signal frequency; and a water-saturated image, acquired after water molecules have been saturated with a water saturation pulse. By calculating the signal difference between these two images, the change in water signal frequency can be estimated, thus obtaining the magnetic field displacement map B0map. This B0map characterizes the main shift and is used to guide the correction of magnetic field inhomogeneity in CEST images, reducing the impact of main magnetic field inhomogeneity on the quality and accuracy of CEST images.

[0048] S2, based on the magnetic field displacement map, the CEST images in the CEST sequence are corrected for magnetic field inhomogeneity using cubic spline interpolation, and / or, based on the structural term images, the corresponding CEST images are registered with structural terms.

[0049] In some embodiments of the present invention, the CEST image in the CEST sequence is corrected for magnetic field inhomogeneity using cubic spline interpolation based on the magnetic field displacement map, including: for each voxel x in the CEST sequence, for each pair of adjacent data points (ω) of voxel x. i y i ) and (ω i+1 y i+1 Construct a cubic polynomial S between ) i (w):

[0050] S i (ω)=a i (ω-ω i ) 3 +b i (ω-ω i ) 2 +c i (ω-ω i )+d i

[0051] Where i = 1, 2, ..., n-1, n represents the number of CEST images in the CEST sequence, a i b i c i d i The coefficients are represented by equation S. i (ω i )=y i S i (ω i+1 )=y i+1 S′ i (ω i+1 )=S′ i+1 (ω i+1 ), S″ i (ω i+1 )=S″ i+1 (ω i+1 S″1(ω1)=0, S″ n-1 (ω n ) = 0 is determined, ω i y represents the acquisition frequency of the i-th CEST image in the CEST sequence. i This indicates that the i-th CEST image in the CEST sequence is at frequency ω. i The grayscale value of voxel X; using the cubic polynomial S i (ω) The interpolation sequence is obtained, and the B0 field offset sequence is obtained based on the interpolation sequence and the magnetic field displacement map; the Z-spectrum correction is performed on the corresponding CEST image in the CEST sequence based on each B0 field offset in the B0 field offset sequence.

[0052] Specifically, the CEST image acquisition process is discrete, collecting data only at set pulse frequency points. However, due to magnetic field inhomogeneity, the actual acquired data may not perfectly correspond to the preset frequency. Therefore, interpolation processing can be performed on the discrete data to generate a smooth Z-spectrum, and the Z-spectrum can be frequency-shifted according to the guidance of BOmap. In this invention, a cubic spline interpolation method is used, the specific process of which is as follows: Figure 2 As shown. For a voxel X on the CEST image, a set of discrete Z-spectral data can be obtained (ω1, y1), (ω2, y2), (ω3, y3)......(ω n y n ), where ω1, ω2......ω n Representing the sampling frequency, y1, y2...y n This is the grayscale value of the voxel at the corresponding frequency. In each pair of adjacent data points (ω... i y i ) and (ω i+1 y i+1 Construct a cubic polynomial S between ) i (ω):

[0053] S i (ω)=a i (ω-ω i ) 3 +b i (ω-ω i ) 2 +c i (ω-ω i )+d i

[0054] Among them, i=1, 2...n-1, a i b i c i d i It is a coefficient, determined by the following conditions:

[0055] (1) Interpolation condition: S i (ω i )=y i S i (ω i+1 )=y i+1

[0056] (2) Smoothing condition: S′ i (ω i+1 )=S′ i+1 (ω i+1 ), S″ i (ω i+1 )=S″ i+1 (ω i+1 )

[0057] (3) Boundary conditions: S″1(ω1)=0, S″ n-1 (ω n ) = 0

[0058] By iterating through all voxels, we obtain the cubic polynomial S corresponding to each voxel. i (ω), and then substitute the gray values ​​of each voxel in each CEST image of the CEST sequence into the corresponding cubic polynomial S. i (ω) yields the interpolation sequence. The grayscale value offset of each voxel in each interpolated image relative to the corresponding voxel in the BOmap image is calculated to obtain the magnetic field offset value Δ, i.e., the B0 field offset. The corresponding interpolated image is shifted by Z spectrum Δ units in the frequency direction using the B0 field offset to complete the magnetic field inhomogeneity correction of the corresponding CEST image.

[0059] In some embodiments of the present invention, structural term registration is performed on the corresponding CEST image based on each structural term image, including: coarse structural term registration of the CEST image and its corresponding structural term image based on normalized mutual information and Powell optimization algorithm to obtain the aligned structural term map corresponding to the CEST image; floating deformation field and reference deformation field are obtained based on the CEST image and the aligned structural term map based on a pre-trained ResNet network, and target deformation field is obtained based on the floating deformation field and reference deformation field; and fine structural term registration is performed on the CEST image using the target deformation field.

[0060] Specifically, while CEST images display rich quantitative information, they also exhibit unclear tissue structures. Therefore, structural item information can be used for localization. In this invention, precise localization is achieved using structural item images acquired at the same anatomical level as the CEST images. First, coarse registration of structural items is performed on the CEST images and their corresponding structural item images based on normalized mutual information and the Powell optimization algorithm to eliminate rigid deformations. Then, a deep network and deformation field are introduced. Based on a pre-trained ResNet network, a floating deformation field and a reference deformation field are obtained from the CEST images and the aligned structural item images. A target deformation field is then obtained from the floating deformation field and the reference deformation field. Finally, the target deformation field is used to perform fine registration of structural items on the CEST images, i.e., non-rigid fine registration of the images, such as organ deformation.

[0061] In some embodiments of the present invention, coarse registration of the CEST image and its corresponding structural item image is performed based on normalized mutual information and the Powell optimization algorithm, including: transforming the structural item image corresponding to the CEST image using a second deformation matrix; calculating the current normalized mutual information between the CEST image and the transformed structural item image; if the current normalized mutual information satisfies a second preset condition, then the transformed structural item image is used as the aligned structural item image of the CEST image, wherein the second preset condition includes: the current normalized mutual information is greater than or equal to a preset normalized mutual information threshold, and / or, the difference between the current normalized mutual information and the normalized mutual information of the previous iteration is less than a second difference threshold; if the current normalized mutual information does not satisfy the second preset condition, then the second deformation matrix is ​​optimized using the Powell algorithm, and the process proceeds to the step of transforming the structural item image corresponding to the CEST image using the second deformation matrix.

[0062] The second deformation matrix can be expressed by the following formula:

[0063]

[0064] Where T2 represents the second deformation matrix, θ represents the rotation angle, and t x t represents displacement along the X-axis. y represents the displacement along the Y-axis, and s represents the scaling factor.

[0065] Specifically, the coarse registration process for structural items is as follows: Figure 3 As shown. In medical imaging, the observed target may move between different acquisition sequences; this movement is usually considered rigid. Furthermore, different acquisition sequences may result in differences in the scaling of the target size in the image due to variations in the field of view (FOV). Therefore, coarse alignment can be performed first for the three types of rigid deformation: translation, rotation, and scaling. The deformation model (i.e., the second deformation matrix) is defined as follows:

[0066]

[0067] In this process, the original coordinates (x, y) in the image are transformed into (x′, y′) after being transformed by the second deformation matrix.

[0068] In this invention, normalized mutual information can be used as a loss function to evaluate the statistical correlation between two images, particularly in the field of medical imaging. Mutual information (MI) is a statistic that measures the amount of information shared between two variables. For two images A and B, mutual information can be defined as:

[0069]

[0070] Where P(a, b) is the joint probability distribution function, and P(a) and P(b) are the marginal probability distribution functions of images A and B, respectively. Normalized mutual information (NMI) is defined as:

[0071]

[0072] Here, H(A) and H(B) are the entropies of A and B, respectively.

[0073] See Figure 3 The structure term image is transformed using the second deformation matrix to obtain the transformed structure term image. Then, the normalized mutual information (NMI) between the transformed image and its corresponding CEST image is calculated. Subsequently, the parameters of the second deformation matrix are optimized based on the normalized mutual information NMI.

[0074] The Powell algorithm can be used as an optimization strategy to find the optimal deformation matrix in each iteration, guiding the parameter changes in the deformation matrix after each iteration. The selected transformation parameters are initially P = [t]. x , t y Let θ, s] define an initial transformation direction D = {d1, d2, d3, d4} for each transformation parameter. For each transformation direction d i Find a scalar λ such that NMI(ref, P+λ*d) i Maximize, where ref is the transformed structural term image, and the NMI calculation method is the same as described above:

[0075]

[0076] Update the current parameter using the found value:

[0077] p new ←p old +λ * d i

[0078] After each iteration, a new set of transformation parameters p is obtained. new Next, update the search direction D. First, calculate the new direction:

[0079] d new =p new -p old

[0080] Then, the worst-performing direction (i.e., the one that contributes the least to NMI) is replaced with the value from the new direction, and the next iteration is started until the iteration termination condition (i.e., the second preset condition) is met: the normalized mutual information is large enough (greater than the preset normalized mutual information threshold) or the change in normalized mutual information calculated in each iteration is small enough (the difference is less than the second difference threshold), and the alignment structure item map corresponding to the CEST image is obtained.

[0081] In some embodiments of the present invention, a floating deformation field and a reference deformation field are obtained based on a pre-trained ResNet network according to the CEST image and the aligned structure term map, and a target deformation field is obtained based on the floating deformation field and the reference deformation field. This includes: extracting features from the CEST image using a pre-trained ResNet50 network to obtain a first feature map; dividing the CEST image into 2*2 grids to obtain (m / 2)*(n / 2) grids; and inputting the information of each voxel point on the first feature map into the corresponding grid. The midpoint of each grid cell is used as a control point to construct a floating deformation field. A pre-trained ResNet50 network is used to extract features from the alignment structure termmap corresponding to the CEST image, resulting in a second feature map. This second feature map is then divided into 2x2 grids, resulting in (m / 2)*(n / 2) grids. Information from each voxel point in the second feature map is input into the corresponding grid, and the midpoint of each grid cell is used as a control point to construct a reference deformation field. For each control point D in the deformation field, multiple control points D' around that control point D are selected. i A set of points to be matched is constructed, and the feature vector T(D) of the corresponding control point D is obtained on the CEST image. The feature vector T′(D) of each control point in the set of points to be matched is obtained on the alignment structure term map. i ), calculate T(D) and T′(D) i The Euclidean distance between the control points D in the floating deformation field and the optimal matching point in the reference deformation field are determined based on the Euclidean distance; the target deformation field is obtained based on all the optimal matching points.

[0082] Specifically, after rigid coarse alignment, precise registration is required for non-rigid deformations such as organ movement within the pelvic cavity. Non-rigid deformation registration is typically based on a suitable deformation field. The distortion of the deformation field is manipulated by moving control points on the field, and then the deformation field is applied to the floating image to achieve registration. The structural fine registration process is as follows: Figure 4 and Figure 5 As shown.

[0083] Because of the significant differences between the CEST image and the structure item image, relying solely on grayscale or texture calculations is insufficient to achieve good registration results. Therefore, a ResNet50 network pre-trained on the ImageNet dataset is introduced. This network possesses powerful image recognition capabilities, and its intermediate layers can output features at different levels, gradually transitioning from basic texture to high-level semantic information. See [link to related documentation] Figure 4 An m*n image (CEST image or its corresponding alignment structure term map) is input into the ResNet50 network. The outputs of the first four layers of the network are summed to obtain the final 3D feature map D (the one corresponding to the CEST image is denoted as the first feature map, and the one corresponding to the alignment structure term map is denoted as the second feature map). This feature map is halved in the row and column directions, adding a third dimension of length t, meaning each pixel corresponds to a 1*t one-dimensional feature vector. The original input image is divided into a 2*2 grid, resulting in (m / 2)*(n / 2) grids, corresponding to voxels on the feature map. The midpoint of each small grid is taken as a control point to form a deformation field (the one corresponding to the CEST image is denoted as the floating deformation field, and the one corresponding to the alignment structure term map is denoted as the reference deformation field). Each control point can then obtain a corresponding feature vector.

[0084] Subsequently, the deformation field is guided by the feature description. Currently popular deformation field models include spline surfaces and differential surfaces. The shape of a spline surface is determined by the control points and their surrounding basis functions, while a differential surface simulates the deformation behavior under different materials and stress states by solving partial differential equations. In this invention, both of these deformation field options are incorporated to find the optimal registration result.

[0085] In the previous step, we linked the deformation field control points and feature maps. For a detailed algorithm flowchart using a single control point as an example, please refer to [link to flowchart documentation]. Figure 5 The control points within the floating deformation field correspond one-to-one with the control points in the reference deformation field. For a control point position D, several surrounding control points D'' are also selected. i (i = 1 to n) form the set of points to be matched. The feature vector T(D) of the corresponding control point D is obtained on the floating feature map. The feature vector T′(D) of each control point in the set of points to be matched is obtained on the reference feature map. i (i = 1 to n), calculate T(D) and each T′(D) respectively. i The Euclidean distance between ) is calculated, and the minimum Euclidean distance corresponding to D is calculated. i The optimal matching point of D in the floating deformation field is determined. By traversing all control points, the optimal matching point is calculated. Based on all optimal matching points, the movement of the floating control points is guided, and finally the transformed deformation field, i.e., the target deformation field, is obtained.

[0086] After obtaining the target deformation field, accurate registration is performed based on it. Assuming a voxel point P on the CEST image has a displacement (u, v) in the target deformation field, the position P′ after applying the target deformation field can be calculated using the following formula:

[0087] P′=(x+u,y+v)

[0088] Based on the information of the target deformation field, this formula is applied to each voxel point on the CEST image to achieve non-rigid deformation of the entire CEST image, making it correspond to the target shape. Thus, precise registration of structural terms is achieved.

[0089] In some embodiments of the present invention, before performing main magnetic field inhomogeneity correction and structural term registration on the CEST images, the method further includes: performing motion correction on the CEST images in the CEST sequence based on Z-spectrum similarity.

[0090] Specifically, in some examples, such as Figure 6 As shown, after acquiring the CEST image, the preprocessing workflow can sequentially include run-through correction, spatial magnetic field inhomogeneity correction, and multimodal registration (i.e., structure term registration), outputting the preprocessed CEST image and its corresponding aligned structure term map. The positions of spatial magnetic field inhomogeneity correction and multimodal registration can be interchanged.

[0091] In some embodiments of the present invention, motion correction is performed on each CEST image in a CEST sequence based on Z-spectrum similarity, including: expanding each CEST image in the current CEST sequence into a one-dimensional column vector with m*n rows and one column, where m*n is the size of the CEST image; merging the one-dimensional column vectors in the column direction to obtain a Z-spectrum set matrix of all voxels in the CEST sequence, where the column direction of the Z-spectrum set matrix represents different acquisition frequencies and the row direction represents different voxels; using the entropy weight topsis algorithm to obtain a Z-spectrum similarity map based on the Z-spectrum set matrix, and calculating the current sharpness of the Z-spectrum similarity map; if the current sharpness or the current iteration number meets a first preset condition... If the current CEST sequence is used as the motion-corrected CEST sequence, the first preset condition includes at least one of the following: the current sharpness is greater than or equal to a preset sharpness threshold, the difference between the current sharpness and the sharpness of the previous iteration is less than or equal to a first difference threshold, and the current iteration number is greater than or equal to a preset iteration number threshold. If the current sharpness and the current iteration number do not meet the preset conditions, the parameters of the current first deformation matrix are updated using the gradient descent method, and the updated first deformation matrix is ​​used to transform each CEST image in the current CEST sequence, and then proceeds to the step of expanding each CEST image in the current CEST sequence into an m*n row and column one-dimensional column vector.

[0092] Calculating the current sharpness of the Z-spectral similarity map may include: using the Sobel operator to calculate the gradient G in the x-direction of the Z-spectral similarity map. x The gradients G in the (x, y) and Y directions y (x, y); The gradient value G(x, y) of each voxel point (x, y) in the Z-spectrum similarity map is obtained by the following formula:

[0093]

[0094] The current sharpness of the Z-spectrum similarity map is obtained by averaging the gradient values ​​of all voxel points.

[0095] The first deformation matrix can be expressed by the following equation:

[0096]

[0097] Where T1 represents the first deformation matrix, θ represents the rotation angle, and d x d represents displacement in the X-axis direction. y This indicates displacement along the Y-axis.

[0098] Specifically, motion correction of CEST sequences based on Z-spectrum similarity mainly addresses the misalignment problem of target bodies during CEST sequence acquisition. The specific algorithm flow is as follows: Figure 7 As shown, the CEST sequence image has dimensions of m*n*w, where m*n represents the size of each CEST image and w represents the total number of frequencies collected.

[0099] When entering the correction process, since there is no guiding information in the initial iteration, the parameter matrix is ​​usually assigned an initial value of 0 or a random value. In subsequent iterations, it is continuously optimized based on gradient guidance. For the rigid deformation matrix (i.e., the first deformation matrix), three parameters are selected: θ (rotation angle), d... x (X-axis displacement), d y (Y-axis displacement) can be transformed using the following formula:

[0100]

[0101] Where (x0, y0) represents the coordinates in the original image, which are obtained after transformation as (x0, y0).

[0102] The optimization iterative process requires setting appropriate termination conditions and gradient descent parameters. This invention selects three judgment conditions (i.e., the first preset condition): when the image sharpness reaches a preset high threshold (i.e., preset sharpness threshold); when the difference between the sharpness calculated in two consecutive iterations drops to a low level (i.e., the first difference threshold); and when the number of iterations exceeds the preset maximum number of iterations (i.e., the preset iteration number threshold). The gradient descent parameter, i.e., the learning rate, can be set to 0.01.

[0103] See Figure 7 During motion correction, motion artifacts in the CEST sequence are first quantified to obtain a Z-spectrum similarity map. For example... Figure 8 As shown, each m*n CEST image in the CEST sequence is expanded into an m*n row, one-column one-dimensional column vector. Then, w one-dimensional column vectors are merged along the column direction to obtain the z-spectrum set matrix of all voxels. This matrix represents different acquisition frequencies along the column direction and different voxels along the row direction. Treating the z-spectrum of each voxel as an object and each acquisition frequency as an evaluation dimension, the entropy-weighted TOPSI algorithm can be used to score each z-spectrum. After normalizing the score, the score only represents the similarity between z-spectrums, thus obtaining a z-spectrum similarity image.

[0104] Rigid motion in a sequence of images can produce edge motion artifacts in the Z-spectrum similarity image. The sharpness of the image is related to the edge sharpness of objects in the image, and the sharpness of the edges can be quantified by calculating the image gradient, thus quantifying the motion artifacts. The Sobel operator can be used to calculate the gradient G in the X and Y directions of the image. x (x, y) and G y (x, y), the final scalar gradient value at the voxel point is:

[0105]

[0106] The image sharpness can be obtained by averaging the gradient values ​​of all voxels.

[0107] Motion correction for CEST images primarily involves finding the optimal deformation parameters for the first deformation matrix corresponding to each image. After obtaining the sharpness, a gradient descent method is introduced to guide the transformation of the deformation parameters, and a numerical approximation method is selected to solve for the gradient. Taking the solution of the rotation angle gradient as an example:

[0108]

[0109] Δθ represents a very small deflection angle, F(·) represents the image sharpness, img represents the Z-spectral similarity map before deflection, and img′ represents the Z-spectral similarity map after the original sequence image undergoes a Δθ transformation. The updated θ is as follows:

[0110] θ = θ + gradient d Parameters * d θ

[0111] It should be noted that structural term registration serves as a reference for spectral extraction analysis and is an indispensable step in the preprocessing process. The necessity of magnetic field non-uniformity correction in the preprocessing process of this invention will be illustrated below through experimental data.

[0112] Two subjects were selected in the experiment. The evaluation indicators included asymmetric quantization plots and Z-spectrum curves. Three experiments were conducted: Experiment 1 used a preprocessing procedure of motion correction + main magnetic field correction + structure term registration; Experiment 2 used a preprocessing procedure of main magnetic field correction + structure term registration; and Experiment 3 used a preprocessing procedure of structure term registration.

[0113] For asymmetric quantization maps, the differences before and after asymmetric quantization can be compared. The asymmetric quantization map at a specific frequency signal can be solved using the following formula to obtain information about the solute at the corresponding saturation frequency:

[0114]

[0115] Where S0 is the intensity of the water signal without a saturation pulse. and These represent the distance from the reference water frequency under the action of the saturation pulse. and Signal strength at that location.

[0116] For Z-spectrum curves, the shapes of the Z-spectrum curves before and after processing the same ROI can be compared. The Z-spectrum curves before and after processing, generated based on the selected ROI region, can be observed; this curve reflects the changing trend of the average intensity value of all saturated images within the region of interest. ROIs are typically selected from areas of detail in the image, such as edge regions or lesion areas.

[0117] Based on the above evaluation indicators, the following conclusions can be drawn through experiments:

[0118] The pretreatment effects of Experiment 1 and Experiment 2 are compared, for example... Figure 9 As shown, motion artifacts are artificially added interference, representing data without motion correction. The white box circles the lesion area. Figure 9 This shows that the analysis results of the corrected lesion area are close to the analysis results without added motion damage. It also concludes that motion correction is unnecessary when the original acquired image does not contain motion artifacts, i.e., the target body does not exhibit significant motion.

[0119] The pretreatment effects of Experiments 2 and 3 are compared, for example... Figure 10 As shown, from Figure 10 As can be seen, the asymmetric analysis plot after magnetic field correction is more uniform, and the axis of symmetry of the Z-spectrum (the dashed line in the figure) is closer to 0 ppm. Furthermore, since the magnetic field inhomogeneity is not caused by human factors, this step is indispensable.

[0120] In summary, the CEST image preprocessing method of this invention introduces Z-spectrum similarity as a new metric in the motion correction stage; introduces cubic spline interpolation in the main magnetic field correction process; and applies ResNet neural network to extract feature descriptions in the structural term registration stage, thereby achieving the following beneficial effects:

[0121] 1) Improved the accuracy and reliability of the processing results;

[0122] 2) High Processing Efficiency: In the entire preprocessing workflow, computation time is comprehensively considered as a crucial metric for each step. While ensuring image processing quality, the most time-efficient algorithms are prioritized to guarantee the overall efficiency of the workflow. Furthermore, the standardized processing workflow reduces errors caused by experimental variations, facilitating the evaluation and comparison of future new methods.

[0123] 3) Wide range of applications: This invention is not only applicable to the pelvic region but can also be transferred to similar multi-organ tissue environments, such as the thoracic cavity and kidneys, to improve data quality. Simultaneously, this method is also suitable for processing images of relatively simple targets, such as the human brain and rats. Tests on different datasets show that this invention performs excellently in terms of accuracy and robustness, which is conducive to its widespread application and promotion in the CEST field.

[0124] 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 sequenced list of executable instructions for implementing logical functions, and can be embodied in any computer-readable medium for use by, or in conjunction with, an instruction execution system, apparatus, or device (such as a computer-based system, a processor-included system, or other system that can fetch and execute instructions from, an instruction execution system, apparatus, or device). For the purposes of this specification, "computer-readable medium" can be any means that can contain, store, communicate, propagate, or transmit programs for use by, or in conjunction with, an instruction execution system, apparatus, or device. More specific examples (a non-exhaustive list) of computer-readable media include: an electrical connection having one or more wires (electronic device), a portable computer disk drive (magnetic device), random access memory (RAM), read-only memory (ROM), erasable and editable read-only memory (EPROM or flash memory), fiber optic devices, and portable optical disc read-only memory (CDROM). Alternatively, the computer-readable medium may be paper or other suitable media on which the program can be printed, since the program can be obtained electronically, for example, by optically scanning the paper or other medium, followed by editing, interpreting, or otherwise processing as necessary, and then stored in a computer memory.

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

[0126] In the description of this specification, references to terms such as "one embodiment," "some embodiments," "example," "specific example," or "some examples," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of the invention. In this specification, the illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples.

[0127] In the description of this invention, it should be understood that the terms "center," "longitudinal," "lateral," "length," "width," "thickness," "upper," "lower," "front," "rear," "left," "right," "vertical," "horizontal," "top," "bottom," "inner," "outer," "clockwise," "counterclockwise," "axial," "radial," and "circumferential" indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings. They are used only for the convenience of describing this invention and simplifying the description, and are not intended to indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limitations on this invention.

[0128] Furthermore, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of technical features indicated. Thus, a feature defined as "first" or "second" may explicitly or implicitly include at least one of that feature. In the description of this invention, "a plurality of" means at least two, such as two, three, etc., unless otherwise explicitly specified.

[0129] In this invention, unless otherwise explicitly specified and limited, the terms "installation," "connection," "linking," and "fixing," etc., should be interpreted broadly. For example, they can refer to a fixed connection, a detachable connection, or an integral part; they can refer to a mechanical connection or an electrical connection; they can refer to a direct connection or an indirect connection through an intermediate medium; they can refer to the internal communication of two components or the interaction between two components, unless otherwise explicitly limited. Those skilled in the art can understand the specific meaning of the above terms in this invention according to the specific circumstances.

[0130] In this invention, unless otherwise explicitly specified and limited, "above" or "below" the second feature can mean that the first feature is in direct contact with the second feature, or that the first feature is in indirect contact with the second feature through an intermediate medium. Furthermore, "above," "over," and "on top" of the second feature can mean that the first feature is directly above or diagonally above the second feature, or simply that the first feature is at a higher horizontal level than the second feature. "Below," "below," and "under" the second feature can mean that the first feature is directly below or diagonally below the second feature, or simply that the first feature is at a lower horizontal level than the second feature.

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

Claims

1. A preprocessing method for CEST images, characterized in that, include: Obtain the CEST sequence, and obtain the structural item image corresponding to each CEST image in the CEST sequence, and determine the magnetic field displacement map corresponding to the CEST sequence; Based on the magnetic field displacement map, cubic spline interpolation is used to correct the main magnetic field inhomogeneity of the CEST images in the CEST sequence, and structural term registration is performed on the corresponding CEST images based on each of the structural term images. Before performing main magnetic field inhomogeneity correction and structure term registration on the CEST image, the method further includes: Motion correction is performed on each CEST image in the CEST sequence based on Z-spectrum similarity; The motion correction of each CEST image in the CEST sequence based on Z-spectrum similarity includes: Expand each CEST image in the current CEST sequence into m. A one-dimensional column vector with n rows and one column, where m n is the size of the CEST image; the one-dimensional column vectors are merged along the column direction to obtain the Z-spectrum set matrix of all voxels in the CEST sequence, where the column direction of the Z-spectrum set matrix represents different acquisition frequencies and the row direction represents different voxels; the entropy weight topsis algorithm is used to obtain the Z-spectrum similarity map based on the Z-spectrum set matrix, and motion correction is performed on each CEST image in the CEST sequence based on the Z-spectrum similarity map; The structural item registration of the corresponding CEST image based on each of the structural item images includes: Based on normalized mutual information and the Powell optimization algorithm, coarse registration of structural terms is performed on the CEST image and its corresponding structural term image to obtain the aligned structural term map corresponding to the CEST image. Based on the pre-trained ResNet network, floating deformation field and reference deformation field are obtained from the CEST image and the aligned structural term map, and target deformation field is obtained from the floating deformation field and the reference deformation field. Finally, fine registration of structural terms is performed on the CEST image using the target deformation field.

2. The preprocessing method for CEST images according to claim 1, characterized in that, The motion correction of each CEST image in the CEST sequence based on the Z-spectrum similarity map includes: Calculate the current sharpness of the Z-spectral similarity map; If the current sharpness or the current iteration number meets the first preset condition, then the current CEST sequence is taken as the motion-corrected CEST sequence. The first preset condition is that at least one of the following three conditions is met: the current sharpness is greater than or equal to a preset sharpness threshold, the difference between the current sharpness and the sharpness of the previous iteration is less than or equal to a first difference threshold, and the current iteration number is greater than or equal to a preset iteration number threshold. If none of the three conditions are met, the parameters of the current first deformation matrix are updated using gradient descent, and the updated first deformation matrix is ​​used to transform each CEST image in the current CEST sequence, proceeding to the step of unfolding each CEST image in the current CEST sequence into m... The steps to create an n-row, one-column one-dimensional column vector.

3. The preprocessing method for CEST images according to claim 2, characterized in that, The calculation of the current sharpness of the Z-spectral similarity map includes: The gradient of the Z-spectral similarity map in the X direction is calculated using the Sobel operator. gradient in the Y direction ; The voxel points in the Z-spectrum similarity map are obtained by the following formula. gradient value : ； The current sharpness of the Z-spectrum similarity map is obtained by averaging the gradient values ​​of all voxel points.

4. The preprocessing method for CEST images according to claim 2, characterized in that, The first deformation matrix is ​​represented by the following formula: in, Represents the first deformation matrix. Indicates the rotation angle. Indicates displacement in the X-axis direction. This indicates displacement along the Y-axis.

5. The preprocessing method for CEST images according to claim 1, characterized in that, The magnetic field displacement map is obtained by calculating the water signal frequency difference between the water saturation image and the reference image corresponding to the CEST sequence. The reference image is the CEST image acquired without a water saturation pulse in the CEST sequence, and the water saturation image is the CEST image acquired after water molecules are saturated with a water saturation pulse in the CEST sequence. The correction of the main magnetic field inhomogeneity of the CEST images in the CEST sequence using cubic spline interpolation based on the magnetic field displacement map includes: For each voxel X in the CEST sequence, in each pair of adjacent data points of voxel X and Construct a cubic polynomial between them : in, , n This indicates the number of CEST images in the CEST sequence. Representing the coefficients, from the formula , , Sure, This represents the acquisition frequency of the i-th CEST image in the CEST sequence. This indicates that the i-th CEST image in the CEST sequence has a frequency... The grayscale value of voxel X; Using the cubic polynomial The interpolation sequence is obtained, and the B0 field offset sequence is obtained based on the interpolation sequence and the magnetic field displacement map; Z-spectral correction is performed on the corresponding CEST image in the CEST sequence based on each B0 field offset in the B0 field offset sequence.

6. The preprocessing method for CEST images according to claim 1, characterized in that, The coarse registration of structural terms between the CEST image and its corresponding structural term image based on normalized mutual information and Powell optimization algorithm includes: The structural item image corresponding to the CEST image is transformed using the second deformation matrix; Calculate the current normalized mutual information between the CEST image and the transformed structure term image; If the current normalized mutual information satisfies the second preset condition, then the transformed structure item image is used as the aligned structure item image of the CEST image. The second preset condition includes: the current normalized mutual information is greater than or equal to a preset normalized mutual information threshold, and / or the difference between the current normalized mutual information and the normalized mutual information of the previous iteration is less than the second difference threshold. If the current normalized mutual information does not meet the second preset condition, then the second deformation matrix is ​​optimized using the Powell algorithm, and the process proceeds to the step of transforming the structural item image corresponding to the CEST image using the second deformation matrix.

7. The preprocessing method for CEST images according to claim 6, characterized in that, The second deformation matrix is ​​expressed by the following formula: T2 ， in, This represents the second deformation matrix. Indicates the rotation angle. Indicates displacement in the X-axis direction. represents the displacement along the Y-axis, and s represents the scaling factor.

8. The preprocessing method for CEST images according to claim 1, characterized in that, The method, based on a pre-trained ResNet network, obtains a floating deformation field and a reference deformation field from the CEST image and the aligned structure term map, and obtains a target deformation field from the floating deformation field and the reference deformation field, including: The CEST image is used to extract features using a pre-trained ResNet50 network to obtain a first feature map, and the CEST image is then divided into 2... A grid of size 2 yields (m / 2). (n / 2) grids, and input the information of each voxel point on the first feature map into the corresponding grid, and take the midpoint of each grid as the control point to form the floating deformation field; The alignment structure term map corresponding to the CEST image is extracted using a pre-trained ResNet50 network to obtain a second feature map, and the alignment structure term map is divided into 2... A grid of size 2 yields (m / 2). (n / 2) grids, and input the information of each voxel point on the second feature map into the corresponding grid, and take the midpoint of each grid as the control point to form the reference deformation field; For each control point D in the deformation field, select multiple control points around control point D. A set of points to be matched is constructed, and the feature vectors of the corresponding control points D are obtained on the first feature map. The feature vector of each control point in the set of points to be matched is obtained on the second feature map. ,calculate and The Euclidean distance between them is used to determine the optimal matching point of the control point D on the floating deformation field on the reference deformation field based on the Euclidean distance. The target deformation field is obtained based on all the best matching points.

Images