CEST image denoising method and storage medium

Abstract

The invention discloses a CEST image denoising method and a storage medium, and relates to the technical field of image processing. The method comprises the following steps: acquiring a CEST image; variance stabilization transformation is carried out on the CEST image, and the transformed CEST image is used as a denoising input image; performing denoising processing on the denoised input image based on the spatial information channel and the frequency spectrum information channel to obtain a first denoised image and a second denoised image; fusing the first de-noised image and the second de-noised image to obtain a fused image; and performing reverse variance stabilization transformation on the fused image to obtain a target de-noised image. The method can retain important details and signal components in the CEST image by suppressing noise from two dimensions of space and frequency spectrum, so that the denoising effect can be remarkably improved, and the accuracy and credibility of a processing result are ensured.

Description

Technical Field

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

[0002] Chemical Exchange Saturation Transfer (CEST) Magnetic Resonance Imaging (MRI) is a new type of MRI molecular imaging method that has attracted much attention and can reveal pathological changes at the molecular level. In biological tissues, there are various exchangeable solute protons with different resonance frequencies, such as amide groups and hydroxyl groups. By using a saturation pulse with a specific frequency to label the target solute molecules, the exchangeable solute molecules are induced to undergo chemical exchange with water molecules and reach a steady state, which results in the attenuation of the water signal and can then be detected. Although the concentrations of these small solute molecules in biological tissues are usually in the millimolar (mM) or even micromolar (μM) range, as long as appropriate experimental parameters are selected, the cumulative effect of the chemical exchange between the exchangeable protons and water molecules can be used to amplify the signal. Therefore, CEST has a high sensitivity (compared with traditional MRI, it can reach 102 - 106).

[0003] Although CEST shows great application potential in revealing pathological changes at the molecular level, due to the complexity of the experimental environment and the numerous experimental parameters, it is vulnerable to noise interference in practical applications. The presence of noise not only reduces the Signal-to-Noise Ratio (SNR) of CEST images but also seriously affects the quantitative accuracy and reliability of the imaging results. Therefore, one of the research focuses in the CEST field is how to effectively solve the noise problem and improve the image quality. Developing more efficient denoising algorithms and imaging technologies to reduce the negative impact of noise on CEST images is of great significance for improving the lesion detection ability, enhancing the diagnostic accuracy, and promoting its application in clinical medicine.

[0004] In recent years, in order to improve the signal-to-noise ratio of CEST images, researchers have begun to combine prior knowledge for image processing. For example, the Principal Component Analysis (PCA) method is used to identify the main features of the Z-spectrum to enhance the denoising performance; the Multilinear Singular Value Decomposition (MLSVD) method combined with spatio-temporal low-rank prior, which integrates spatio-temporal low-rank prior and enhances the image quality by utilizing spatial self-similarity; the pixel-level denoising filter based on non-local self-similarity (NlmCED), which restores pixels by weighted averaging of similar pixels within a 3D search window; the image denoising method (BOOST) that combines non-local low-rank constraint and spectral local smoothing regularization, etc.

[0005] Although the above denoising techniques show certain effects as a whole, there are still some challenges in specific applications. For example, when rearranging CEST images into a two-dimensional matrix and applying the PCA method, the inherent spatio-temporal correlation in the images may be lost; MLSVD has difficulties in ranking threshold design and does not fully consider local self-similarity; the implementation of NlmCED at the pixel level may affect the preservation of structural information; the BOOST algorithm has a high computational complexity and is sensitive to parameter adjustment, and sometimes may lead to over-smoothing problems. Therefore, how to effectively overcome the limitation of noise on CEST imaging quality is still one of the key directions for future research and development in this field. Summary of the Invention

[0006] 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 CEST image denoising method and a storage medium to improve the denoising effect of CEST images and ensure the accuracy and reliability of the denoising results.

[0007] In a first aspect, an embodiment of the present invention provides a CEST image denoising method, including: acquiring a CEST image; performing variance stabilization transformation on the CEST image, and using the transformed CEST image as a denoising input image; respectively performing denoising processing on the denoising input image based on a spatial information channel and a spectral information channel to obtain a first denoised image and a second denoised image; fusing the first denoised image and the second denoised image to obtain a fused image; performing inverse variance stabilization transformation on the fused image to obtain a target denoised image.

[0008] According to an embodiment of the present invention, the method further includes: determining whether the fused image meets a preset termination condition; if not, updating the denoising input image to the fused image and returning to the step of respectively performing denoising processing on the denoising input image based on the spatial information channel and the spectral information channel; if so, performing the step of performing inverse variance stabilization transformation on the fused image.

[0009] According to an embodiment of the present invention, denoising the denoising input image based on a spatial information channel to obtain the first denoised image, including: dividing the denoising input image into a plurality of first three-dimensional data blocks, superimposing similar first three-dimensional data blocks, respectively processing each superimposed result by one-dimensional decorrelation linear transformation to obtain first spectral coefficients corresponding to each superimposed result, and respectively performing hard threshold filtering processing on each of the first three-dimensional data blocks based on the first spectral coefficients, and aggregating the first three-dimensional data blocks after hard threshold filtering back to the original position by an adaptive weighted average method to obtain a first denoised sub-image; dividing the first denoised sub-image into a plurality of second three-dimensional data blocks, superimposing similar second three-dimensional data blocks, respectively processing each superimposed result by one-dimensional decorrelation linear transformation to obtain second spectral coefficients corresponding to each superimposed result, and respectively performing Wiener filtering processing on each of the second three-dimensional data blocks based on the second spectral coefficients, and aggregating the second three-dimensional data blocks after Wiener filtering back to the original position by an adaptive weighted average method to obtain the first denoised image.

[0010] According to an embodiment of the present invention, similar three-dimensional data blocks are determined by the following method:

[0011] Calculate the distance between any two three-dimensional data blocks by the following formula:

[0012]

[0013] where, represents the i-th three-dimensional data block and the j-th three-dimensional data block The distance between them, z represents the CEST image, x represents the three-dimensional coordinates in the signal domain X, L represents the side length of the three-dimensional data block, and ‖·‖ 2 represents the 2-norm; if the distance is less than a preset distance threshold, it is determined that these two three-dimensional data blocks are similar.

[0014] According to an embodiment of the present invention, the adaptive weight used in the adaptive weighted average method is obtained by the following formula:

[0015]

[0016] where, represents the image after filtering processing, x R ∈X represents traversing each three-dimensional data block in the image X before filtering processing, represents other three-dimensional data blocks in X that are similar to the reference three-dimensional data block x R similar, represents x RFor each x i the weight of represents x after filtering i the estimated value of, y represents the data before filtering, represents the characteristic function for indicating x i whether it is within the valid range, represents the number of non-zero coefficients after filtering, and σ represents the standard deviation of the noise.

[0017] According to an embodiment of the present invention, denoising the denoising input image based on the spectral information channel includes: constructing the following objective function;

[0018]

[0019] wherein, X represents the denoising input image, W represents the first weight tensor, Y represents the denoised image, represents the square of the Frobenius norm, λ 1 and λ 2 respectively represent regularization parameters, ‖‖ 1 represents the L1 norm, S represents the second weight tensor, D represents the three-dimensional first-order forward difference operator, ⊙ represents the element-wise product, P represents a local region, and P i represents the index of the i-th local block, ‖‖ w,* represents the weighted nuclear norm, s.t. X = U × 3 V represents the constraint condition: X is reconstructed by the three-dimensional tensor product of the projection data U and the spectral subspace basis V;

[0020] Denoise the denoising input image based on the objective function.

[0021] According to an embodiment of the present invention, after determining that the fused image does not meet the preset termination condition, the method further includes: calculating the first energy of the first denoised image and the second energy of the second denoised image; obtaining the spatial channel weight and the spectral channel weight according to the first energy and the second energy for obtaining the fused image next time; wherein, the fused image is obtained by weighting the first denoised image and the second denoised image using the spatial channel weight and the spectral channel weight.

[0022] According to an embodiment of the present invention, the spatial channel weight and the spectral channel weight are obtained by the following formula:

[0023]

[0024] wherein, λ 3 and λ 4respectively represent the spatial channel weight and the spectral channel weight, energy 1 and energy 2 respectively represent the first energy and the second energy, and ∈ represents a constant.

[0025] According to an embodiment of the present invention, determining whether the fused image meets a preset termination condition includes: calculating a ratio of a noise level of the fused image to a noise level of the CEST image; if the ratio is less than a preset ratio threshold, determining that the fused image meets the preset termination condition.

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

[0027] The CEST image denoising method and storage medium according to the embodiments of the present invention perform variance stabilization transformation on the acquired CEST image, and use the transformed CEST image as a denoising input image; then perform denoising processing on the denoising input image based on a spatial information channel and a spectral information channel respectively to obtain a first denoised image and a second denoised image; then fuse the first denoised image and the second denoised image to obtain a fused image; and then perform inverse variance stabilization transformation on the fused image to obtain a target denoised image. Thus, by suppressing noise from two dimensions of space and spectrum, important details and signal components in the CEST image can be retained, thereby significantly improving the denoising effect and ensuring the accuracy and reliability of the processing result.

[0028] Additional aspects and advantages of the present invention will be given in part in the following description, become apparent in part from the following description, or be understood through the practice of the present invention. Description of the Drawings

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

[0030] Figure 2 is a flowchart of a CEST image denoising method according to another embodiment of the present invention;

[0031] Figure 3 is a comparison diagram of quantization effects obtained by using different denoising methods according to an example of the present invention;

[0032] Figure 4 is a Bland-Altman diagram of the Z-spectrum difference between a reference value and results obtained by using different denoising methods according to an example of the invention. Detailed Embodiments

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

[0034] The CEST image denoising method and storage medium according to an embodiment of the present invention will be described below with reference to the accompanying drawings.

[0035] Figure 1 It is a flowchart of the CEST image denoising method according to an embodiment of the present invention.

[0036] As Figure 1 shown, the CEST image denoising method includes:

[0037] S11, obtaining a CEST image.

[0038] Taking the magnetic resonance image data of a rat brain as an example, a 9.4T nuclear magnetic resonance scanner (Bruker) can be used to perform magnetic resonance imaging on the rat brain. The scanning sequences are water saturation imaging (WASSR) sequence, CEST sequence (RARE), and T 2 sequence (TurboRARE_sat), and a CEST sequence containing a CEST image can be obtained. Among them, the water saturation imaging sequence is used to correct the magnetic field of the CEST image, and the T 2 sequence is used to perform structural correction on the corresponding CEST image.

[0039] Due to the complexity of the experimental environment and diverse experimental parameters, the CEST imaging process is prone to be contaminated by Rician distribution noise, that is, the acquired CEST image usually contains a certain degree of Rician distribution noise. Exemplarily, the CEST image for denoising can be the one after magnetic field correction and / or structural correction.

[0040] S12, performing a variance stabilization transformation on the CEST image, and using the transformed CEST image as the denoising input image.

[0041] Specifically, since CEST images are vulnerable to Rice distribution noise, denoising the data with Rice distribution noise is a highly challenging task. The characteristic of this noise is that the calculation of its standard deviation is not only related to the unknown signal amplitude, but also there is a non-linear relationship between the expected value of the noise and the noise-free signal amplitude. The existence of this non-linear relationship makes it difficult to directly denoise CEST images. Therefore, the present invention first performs a variance stabilizing transformation (VST) on the CEST image, and uses the transformed CEST image as the denoising input image for subsequent denoising processing.

[0042] After applying the variance stabilizing transformation, the noise distribution can be converted from Rice distribution to Gaussian distribution, and the noisy image obtained can be expressed as the sum of a clean image and Gaussian noise, as shown in the following formula (1):

[0043] z(x) = y(x) + η(x), x ∈ X (1)

[0044] where z(·) represents the noisy image, y(·) represents the clean image, x represents the three-dimensional coordinates belonging to the signal domain and η(·) represents independent and identically distributed Gaussian noise with a mean of zero and a variance of σ 2 .

[0045] The purpose of subsequent denoising is to remove η(·) in the above formula (1). Compared with Rice distribution, Gaussian distribution is more convenient for analysis and calculation. Thus, through the variance stabilizing transformation, the quality and accuracy of subsequent CEST image denoising can be improved.

[0046] S13. Denoise the denoising input image based on the spatial information channel and the spectral information channel respectively to obtain a first denoised image and a second denoised image.

[0047] Specifically, a block matching filter can be used to denoise the denoising input image in the spatial dimension, and a method combining low-rank constraint and frequency-domain total variation regularization can be used to denoise the denoising input image in the spectral dimension.

[0048] In some embodiments of the present invention, denoising the denoising input image based on the spatial information channel to obtain a first denoised image, including: dividing the denoising input image into a plurality of first three-dimensional data blocks, superimposing similar first three-dimensional data blocks, respectively processing each superimposed result by one-dimensional decorrelation linear transformation to obtain first spectral coefficients corresponding to each superimposed result, and performing hard threshold filtering on each first three-dimensional data block based on the first spectral coefficients, and aggregating the first three-dimensional data blocks after hard threshold filtering back to the original position in an adaptive weighted average manner to obtain a first denoised sub-image; dividing the first denoised sub-image into a plurality of second three-dimensional data blocks, superimposing similar second three-dimensional data blocks, respectively processing each superimposed result by one-dimensional decorrelation linear transformation to obtain second spectral coefficients corresponding to each superimposed result, and performing Wiener filtering on each second three-dimensional data block based on the second spectral coefficients, and aggregating the second three-dimensional data blocks after Wiener filtering back to the original position in an adaptive weighted average manner to obtain a first denoised image.

[0049] Exemplarily, when determining whether two three-dimensional data blocks are similar, the distance between these two three-dimensional data blocks can be calculated first by the following formula (2). If the distance is less than a preset distance threshold, it is determined that these two three-dimensional data blocks are similar.

[0050]

[0051] Wherein, represents the i-th three-dimensional data block and the j-th three-dimensional data block The distance between them, L represents the side length of the three-dimensional data block, and ‖·‖ 2 represents the 2-norm.

[0052] Specifically, for spatial dimension denoising, a block matching filter can be used to process the denoising input image. By searching for and matching similar image blocks (i.e., three-dimensional data blocks) in the denoising input image, and using the local correlation of voxels within each image block and the non-local correlation of corresponding voxels between different image blocks, the spatial information can be effectively enhanced and the noise can be suppressed.

[0053] Specifically, the block matching filter can be implemented based on two stages: hard threshold filtering and Wiener filtering. Both stages can perform grouping, collaborative filtering, and aggregation steps on each three-dimensional data block of the denoising input image (which can be composed of multiple three-dimensional data blocks) to effectively separate the clean signal and the noise signal in the transform domain, and achieve efficient denoising through coefficient reduction. In actual implementation, hard threshold filtering can be performed first, and then Wiener filtering; or Wiener filtering can be performed first, and then hard threshold filtering. Here, taking hard threshold filtering first and then Wiener filtering as an example for illustration.

[0054] In the grouping step of the hard threshold filtering stage, cube (i.e., three-dimensional data block) matching grouping is performed, and the denoising input image z is divided into multiple three-dimensional data blocks of size L×L×L. (R represents the selected reference cube). The similarity between two three-dimensional data blocks is measured by calculating the distance (see Equation (2) above), and similar three-dimensional data blocks are grouped into a set. After that, the similar three-dimensional data blocks are superimposed, and a one-dimensional decorrelating linear transformation is applied to each dimension of each superimposed result to obtain the spectrum of a four-dimensional group, forming a four-dimensional array. This four-dimensional array includes the first spectral coefficients corresponding to each superimposed result.

[0055] In the collaborative filtering step of the hard threshold filtering stage, in the transform domain, the following Equation (3) is used to perform hard threshold processing on each first three-dimensional data block based on the first spectral coefficients:

[0056]

[0057] where, represents the result after hard threshold processing, represents the four-dimensional transform of the four-dimensional array ht represents that hard threshold filtering transforms the signal from the spatial domain to the transform domain, represents the four-dimensional transformation process, threshold() represents performing hard threshold processing on the data after four-dimensional transformation, and σλ 4D represents the threshold of hard threshold processing, which is related to the noise standard deviation σ, and λ 4D represents the control parameter of hard threshold processing, which determines how much signal to retain and how much noise to remove during the denoising process.

[0058] In the aggregation step of the hard threshold filtering stage, all the first three-dimensional data blocks after hard threshold processing are aggregated using adaptive weighted averaging through the following Equation (4) to obtain the denoised data estimate (i.e., the above-mentioned first denoised sub-image). The specific weighted aggregation process can be described as follows: a reference voxel block x R and some voxel blocks x i similar to it (they are selected from different positions, and voxel block means three-dimensional data block). These similar voxel blocks will be separately subjected to hard threshold denoising to obtain the denoising estimate value of each voxel block. Then, a weighting coefficient is calculated according to the quality of each voxel block and its similarity to the reference voxel block. The results of all these denoised similar voxel blocks are weighted and averaged to obtain the final denoising estimate value of the reference voxel block x R . After the weighted aggregation of all voxel blocks at all positions, the final denoised data estimate is obtained. (i.e., the first denoised image mentioned above).

[0059]

[0060] Among them, represents the data estimate after denoising (i.e., the first denoised image mentioned above), and \(x\) R \(\in X\) means traversing each data point (voxel block) in the entire image \(X\). represents other voxel blocks similar to the reference voxel block \(x\) R , and these similar voxel blocks will be used for weighting. represents the weight of the reference voxel block \(x\) R for each similar voxel block \(x\) i . This weight determines the contribution degree of the voxel block in the final reconstruction result. represents the estimated value of the voxel block \(x\) i after hard thresholding. \(y\) represents the data before denoising. The hard thresholding process will remove noise and retain important signal components. represents a characteristic function used to indicate whether a voxel block \(x\) i is within the effective range, that is, whether it participates in weighted averaging.

[0061] In addition, is the number of non-zero coefficients after hard thresholding, and \(\sigma\) represents the standard deviation of the noise.

[0062] Similarly, the Wiener filtering stage also includes three steps: data grouping, collaborative filtering, and aggregation. In the data grouping step, based on the result of the hard threshold filtering stage (i.e., the first denoised image), cube matching grouping is performed again. Subsequently, in the collaborative filtering step, in the transform domain, the Wiener filtering method is used based on the second spectral coefficient to further improve the separation effect of the signal and noise. The corresponding noise group is extracted from the observed data \(z\) and Wiener filtering is applied. The Wiener filtering coefficient formula is shown in Equation (5) below:

[0063]

[0064] Among them, represents the Wiener filtering coefficient, represents the four-dimensional transform applied to the three-dimensional data block (wie represents Wiener filtering), represents the voxel block data after hard thresholding, represents the region of similar voxel blocks related to the reference data block.

[0065] In the aggregation step of the Wiener filtering stage, similar to the hard threshold filtering stage, the estimates of all cubes (the second three-dimensional data block) are aggregated using adaptive weights to obtain the final reconstruction result.

[0066] As described above, hard threshold filtering is a fast denoising method. On this basis, Wiener filtering is further refined and optimized, so that the final reconstruction result can achieve a better balance between retaining details and reducing noise. The two filtering stages complement each other, focusing on analyzing the spatial characteristics of CEST images, which can significantly improve the adaptability and robustness of the algorithm to various noise levels and different types of data.

[0067] In some embodiments of the present invention, denoising the denoising input image based on the spectral information channel includes:

[0068] Constructing the objective function as follows (6);

[0069]

[0070] Where X represents the denoising input image, W represents the first weight tensor, Y represents the denoised image, represents the square of the Frobenius norm, λ 1 and λ 2 respectively represent regularization parameters for balancing the fidelity term and the regularization term; ‖‖ 1 represents the L1 norm, S represents the second weight tensor, D represents the three-dimensional first-order forward difference operator, ⊙ represents the element-wise product, P represents a local region, and P i represents the index of the i-th local block, ‖‖ w,* represents the weighted nuclear norm, s.t. X = U × 3 V represents the constraint condition: X is reconstructed by the three-dimensional tensor product of the projection data U and the spectral subspace basis V; denoising the denoising input image based on the objective function.

[0071] Specifically, the channel that focuses on extracting spectral information integrates the weighted non-local low-rank model and the adaptive total variation regularization. At the same time, the prior knowledge and noise characteristics of the image are modeled, which can effectively capture the local smoothness along the frequency spectrum dimension. The denoising process based on the spectral information channel mainly includes the following four parts:

[0072] First, use the Gaussian mixture model of non-independent and identically distributed to describe the complex noise (i.e., the denoising input image), corresponding to optimizing the weight tensor W involved in a weighted fidelity function, and its probability density function p(E ijb ) is as follows (7):

[0073]

[0074] Where E ijb represents the element at the i-th row, j-th column, and b-th frequency offset in the denoising input image, and π bkrepresents the mixing ratio of the k-th Gaussian component at the b-th frequency offset, γ bk represents the precision (reciprocal of variance) of the k-th Gaussian component at the b-th frequency offset, and K is the number of components of the Gaussian mixture model. represents the probability density function of the k-th Gaussian component at the b-th frequency offset.

[0075] This Gaussian mixture model estimates the noise parameters through variational inference to obtain the weighted fidelity function, as shown in Equation (8) below:

[0076]

[0077] where W = ∑∑∑W ijb , represents the first weight tensor, used to adjust the noise weights at different positions and bands; Loss(Y,X) represents the loss function; represents the square of the Frobenius norm, used to calculate the sum of the squares of the differences of matrix elements.

[0078] In addition, represents E ijb corresponding weight tensor, Z ijbk represents E ijb corresponding latent variable of the k-th Gaussian component, characterizing the distribution of noise.

[0079] Second, use adaptive total regularization to solve the weighted non-local low-rank model, that is, update the low-rank components U and V by means of the non-local similarity and spatial-spectral correlation priors of the data, as shown in Equation (9) below:

[0080]

[0081] where ‖‖ ASSTV represents the adaptive spatial-spectral total variation (ASSTV) regularization term, X represents the input noisy data, ‖‖ 1 represents the L1 norm, used to calculate the sum of the absolute values of matrix elements, and ⊙ represents element-wise multiplication, that is, the corresponding elements of the matrices are multiplied. D = [D h , D v , D s , represents the three-dimensional first-order forward difference operator, along the horizontal spatial, vertical spatial, and spectral directions respectively; S = [S h , S v , S s , represents the second weight tensor, used to adjust the smoothing intensity in different directions.

[0082] The calculation method of the second weight tensor S is as follows:

[0083]

[0084] Among them, X 0 represents the denoised data obtained by low-rank approximation, which is used to reduce the influence of noise. δ represents the threshold to avoid excessive or too small weights. The non-local low-rank model characterizes the non-local similarity and spatial-spectral correlation of hyperspectral images through low-rank matrix approximation. The formula is as follows in Equation (10):

[0085]

[0086] Among them, ‖X‖ NSS represents the non-local spatial similarity (NSS) regularization term. P represents a local region, and Pi represents the index of the i-th local block, which is used for low-rank approximation of the local block. ‖‖ w,* represents the weighted nuclear norm, which is used to enhance the low-rank property. s.t. X = U × 3 V represents the constraint condition, that is, the noisy image X can be reconstructed by the three-dimensional tensor product (× 3 ) of the projection data U and the spectral subspace basis V.

[0087] Third, using the low-rank decomposition results (U and V) obtained in the above steps, recover the data X through inverse transformation.

[0088] Fourth, update the first weight tensor W with the recovered data X and the original noise model (Gaussian mixture model), and adjust the second weight tensor S in the spatial-spectral total variation regularization. Integrate the above parts to obtain the final denoising model (i.e., the objective function shown in Equation (6) above), and the denoising model can be solved by the ADMM (Alternating Direction Method of Multipliers) algorithm.

[0089] S14, fuse the first denoised image and the second denoised image to obtain a fused image.

[0090] Specifically, the spatial channel weight and the spectral channel weight can be used to weight the first denoised image and the second denoised image to obtain a fused image. The following formula can be used:

[0091] Ir = λ 3 *I 3 + λ 4 *I 4

[0092] Among them, Ir represents the fused image, λ 3 and λ 4 represent the spatial channel weight and the spectral channel weight respectively, and I 3 and I 4respectively represent a first denoised image and a second denoised image.

[0093] S15. Perform an inverse variance stabilization transform on the fused image to obtain a target denoised image.

[0094] Specifically, apply the inverse transform of the variance stabilization transform to the data after denoising processing to return the noise-free amplitude under the maximum likelihood estimate. This process is a key step in extracting the signal amplitude closest to the true noise-free state from the denoised data.

[0095] In some embodiments of the present invention, as Figure 2 shown, after performing step S14, the CEST image denoising method further includes:

[0096] S21. Determine whether the fused image meets a preset termination condition.

[0097] If not, perform step S22 and return to step S13; otherwise, perform step S15.

[0098] Specifically, determining whether the fused image meets a preset termination condition may include: calculating the ratio of the noise level of the fused image to the noise level of the CEST image; if the ratio is less than a preset ratio threshold, it is determined that the fused image meets the preset termination condition.

[0099] S22. Update the denoised input image to the fused image.

[0100] Specifically, during the denoising operation, the noise level of the denoising result needs to be monitored in real time. The variance stabilization transform can be used to robustly estimate the noise level to obtain the estimated values of the noise levels before and after denoising, where the noise level is the estimated value of the σ parameter in the Rice distribution. When the noise level after denoising is less than a preset proportional threshold (such as 20%) of the noise level before denoising, it can be determined that the denoising effect has reached the expected standard, and the iterative process can be terminated. This judgment criterion provides an objective basis for the iterative optimization of the denoising algorithm, helps to improve the execution efficiency of the algorithm while ensuring the denoising quality, and avoids unnecessary waste of computing resources. If the termination condition is not met, this denoising result is used as the input data for the next denoising.

[0101] In some embodiments of the present invention, after determining that the fused image does not meet the preset termination condition, the method further includes: calculating the first energy of the first denoised image and the second energy of the second denoised image; obtaining a spatial channel weight and a spectral channel weight based on the first energy and the second energy for obtaining the fused image next time.

[0102] Among them, the energy of the image can reflect features such as image pixel intensity, texture variation, or information density. The image energy can be calculated using the square of grayscale, the sum of squared gradients, and statistical models such as entropy or PCA. Exemplarily, the calculation of the first energy can be based on the spatial channel, combining the spatial domain pixel intensity of the first denoised image with the inter-channel correlation, and quantifying the image information density and variation features through mathematical statistics or transform domain analysis; the calculation of the second energy can be based on the spectral channel, combining time domain and frequency domain analysis methods.

[0103] Specifically, after determining that the termination condition is not met, the dual-channel fusion weights can be updated based on the denoising effect of this time. That is, in each iteration process, after the denoising image is updated, the energies of the two channels are recalculated, and the fusion weights of the two channels are updated accordingly. This dynamic adjustment mechanism enables the continuous and adaptive allocation of the contribution ratios of different denoising methods in the final result based on the actual performance of the current denoising result throughout the iteration process.

[0104] Exemplarily, the spatial channel weight and the spectral channel weight are obtained through the following formula:

[0105]

[0106] where λ 3 and λ 4 represent the spatial channel weight and the spectral channel weight respectively, energy 1 and energy 2 represent the first energy and the second energy respectively, ∈ represents a constant with a small value close to 0, which is used to prevent division-by-zero errors.

[0107] Specifically, the denoising results of the two channels (i.e., the first denoised image and the second denoised image) are fused by weighted averaging. At the first fusion, the fusion ratio (i.e., the spatial channel weight and the spectral channel weight) can be set based on experience. In subsequent fusion processes, the weights can be dynamically adjusted based on the energies of the denoising results of the two channels. This dynamic adjustment mechanism can make the advantages of the two denoising methods complement each other, thereby optimizing the overall denoising performance. The larger the weight, the stronger the denoising effect of the method. By adopting an energy-based weighting strategy, the contributions of the two channels can be adaptively allocated according to the characteristics of the denoising results, effectively combining the strengths of different denoising techniques, and thus a more refined and accurate denoising effect can be achieved.

[0108] The beneficial effects of the CEST image denoising method proposed by the present invention are illustrated through experiments as follows:

[0109] (I) Experimental data

[0110] The brain data of rats were collected using a 9.4T MRI scanner. Single-layer CEST imaging was performed using a rapid acquisition and relaxation enhancement (RARE) sequence. The parameters of CEST imaging can be set as follows: the saturation field strength is 0.7 μT, the repetition time (TR) is 5500 ms, the echo time (TE) is 3.5 ms, and the unsaturated image (S 0 ) was obtained at a frequency offset of -200 ppm. The image slice thickness was 1.5 mm, and the acquired matrix size was 72×96. Under these conditions, 31 frequency offsets in the range from -10 ppm to 10 ppm were acquired. Meanwhile, 21 frequency offsets were collected in the frequency offset range of ±1 ppm using a water saturation imaging sequence. The acquisition parameters and geometry were the same as those of other CEST experiments. T 2 The imaging parameters of the sequence were set as follows: TR was 4500 ms, TE was 36 ms, the flip angle was 180°, the slice thickness was 0.5 mm, and a total of 41 slices were acquired, obtaining 31 frequency offsets in the range from -10 ppm to 10 ppm.

[0111] (II) Evaluation Metrics

[0112] To comprehensively evaluate the denoising effect, three comparison methods were set in the experiment, specifically as follows: the first is the non-local transform domain filter (BM4D) for volume data denoising and reconstruction, the second is the hyperspectral image denoising method (WNLRATV) based on a weighted non-local low-rank model and adaptive total variation regularization, and the third is the denoising method (BOOST) based on spatial spectral redundancy.

[0113] In terms of evaluation metrics, it mainly covers the following contents: one is the quantitative images of three important frequency offsets, specifically including the amide proton transfer (APT) image, the nuclear Overhauser effect (NOE) image, and the phosphocreatine (PCr) quantification result map; the second is the difference between Z-spectrum curves; the third is the structural similarity index (SSIM) used to measure the image similarity.

[0114] Regarding the acquisition methods of the quantitative result images of the three important frequency offsets: the quantification result of the APT image was obtained by averaging the offset images at 3.4 ppm to 3.7 ppm; the NOE quantification result map was obtained by averaging the offset images at -3.5 ppm to -3.2 ppm; the PCr quantification result map was calculated by averaging the offset images at 2.4 ppm to 2.8 ppm.

[0115] (3) Experiment and evaluation

[0116] In the experiment, a representative live rat was selected for data acquisition to obtain CEST images. The CEST images may inevitably be disturbed by noise. Artificial Rice noise with an intensity of 0.2% was added to the reference image to generate the noise image to be processed.

[0117] Figure 3 The contrast diagram of the quantization effects obtained by applying different denoising methods in an example of the present invention is shown. As Figure 3 shown, after artificially adding 0.2% additional Rice noise, the quantization result of the CEST image has decreased significantly. All methods have improved the quantization quality of the CEST image to a certain extent. Specifically, although the WNLRATV method and the BOOST method have achieved certain results in improving the quantization result, the quantization result diagrams show an over-smoothing phenomenon, resulting in the loss of a large amount of detailed information. Among them, the BOOST method is particularly prominent in this regard. It is worth emphasizing that the present method is significantly superior to other comparison methods in terms of noise reduction performance. This advantage can not only be intuitively reflected from the quantization result diagram, but also the higher average structural similarity index value provides strong quantitative support for this conclusion, further confirming the effectiveness and superiority of the present method in improving the quality of CEST images.

[0118] Figure 4 The Bland-Altman diagram of the Z-spectrum difference between the reference value in an example of the present invention and the results obtained by using different denoising methods is shown. In this experiment, four representative regions of interest (ROIs) were carefully selected by experienced researchers in this professional field. On this basis, a detailed Z-spectrum analysis was performed on these four ROIs, and the differences in the Z-spectrum of the reference values of these four ROIs were further explored, and the accuracy and consistency of different denoising methods were evaluated from a statistical perspective.

[0119] By comparing with other denoising methods, the present method produces the smallest deviation among all comparison methods, providing a quantitative result closer to the true value. This result strongly proves the reliability and superiority of the present method in terms of noise reduction performance.

[0120] Based on the above experiments, the CEST image denoising method proposed by the present invention can achieve the following beneficial effects:

[0121] 1) Enhances the accuracy and credibility of processing results: This method introduces non-local similarity matching and low-rank adaptive total variation regularization techniques to achieve efficient denoising and retention of image details. In addition, this method adopts a dynamic weighted fusion strategy based on energy information, which can dynamically adjust the fusion weights according to the contribution degree of each denoising channel, making the final denoising result more balanced and effectively avoiding the deviation that may be caused by a single denoising channel;

[0122] 2) Improves processing efficiency: Through the iterative optimization process, this method achieves the optimal denoising effect within a small number of iterations, significantly improving the processing rate and meeting the requirements of real-time imaging;

[0123] 3) Has broad application prospects: The denoising method of the present invention is widely applicable to CEST image denoising, which has important value for clinical diagnosis and disease research. Especially, it shows remarkable effects when dealing with medical images with low signal-to-noise ratio and scientific research image data. This method not only performs well in the field of CEST image denoising, but also is applicable to other types of medical image denoising tasks and has clinical application potential. In addition, it can be widely applied in the fields of hyperspectral image processing, remote sensing image denoising, and machine vision, effectively improving the image quality.

[0124] Based on the CEST image denoising method of the above embodiments, the present invention proposes a computer-readable storage medium.

[0125] In this embodiment, a computer program is stored on the computer-readable storage medium. When the computer program is executed by a processor, the CEST image denoising method of the above embodiments is implemented.

[0126] It should be noted that the logic and / or steps represented in the flowchart or described in other ways herein, for example, can be considered as a definite sequence list of executable instructions for implementing logical functions, and can be specifically implemented in any computer-readable medium for use by an instruction execution system, apparatus, or device (such as a computer-based system, a system including a processor, or other systems that can fetch and execute instructions from the instruction execution system, apparatus, or device), or in combination with these instruction execution systems, 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 part with one or more wirings (electronic device), a portable computer diskette (magnetic device), a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or flash memory), an optical fiber device, and a portable compact disc read-only memory (CDROM). Additionally, the computer-readable medium can even be paper or other suitable media on which the program can be printed, because the program can be obtained electronically, for example, by optically scanning the paper or other media, followed by editing, interpretation, or other suitable processing as necessary, and then stored in a computer memory.

[0127] 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), etc.

[0128] 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 a suitable manner in any one or more embodiments or examples.

[0129] 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. is based on the orientation or positional relationship shown in the drawings, and is only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and thus should not be construed as a limitation on the present invention.

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

[0131] In the present invention, unless otherwise clearly specified 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.

[0132] In the present invention, unless otherwise clearly specified 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 indirectly in 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.

[0133] 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 limitations 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 denoising method, characterized in that: include: Acquire CEST images; Performing a variance stabilization transformation on the CEST image, and using the transformed CEST image as a denoising input image; Performing denoising processing on the denoised input image based on the spatial information channel and the spectral information channel respectively to obtain a first denoised image and a second denoised image; Fusing the first denoised image and the second denoised image to obtain a fused image; An inverse variance stabilization transform is performed on the fused image to obtain a target denoised image.

2. The method according to claim 1, characterized in that: The method further comprises: Determining whether the fused image meets a preset termination condition; If not, the denoised input image is updated to the fused image, and the process returns to the step of performing denoising processing on the denoised input image based on the spatial information channel and the spectral information channel respectively; If the conditions are met, the step of performing an inverse variance stabilization transform on the fused image is executed.

3. The method according to claim 2, characterized in that Performing denoising processing on the denoised input image based on the spatial information channel to obtain the first denoised image includes: Dividing the denoised input image into a plurality of first three-dimensional data blocks, superimposing similar first three-dimensional data blocks, processing each superposition result by using a one-dimensional decorrelation linear transformation, obtaining a first spectrum coefficient corresponding to each superposition result, performing hard threshold filtering processing on each of the first three-dimensional data blocks based on the first spectrum coefficient, and aggregating the first three-dimensional data blocks after the hard threshold filtering processing back to the original position by means of adaptive weighted averaging, to obtain a first denoised sub-image; The first denoised sub-image is divided into a plurality of second three-dimensional data blocks, and similar second three-dimensional data blocks are superimposed, and each superposition result is processed by using a one-dimensional decorrelation linear transformation to obtain a second spectral coefficient corresponding to each superposition result, and each second three-dimensional data block is subjected to Wiener filtering based on the second spectral coefficient, and the second three-dimensional data blocks subjected to Wiener filtering are aggregated back to their original positions by means of adaptive weighted averaging to obtain the first denoised image.

4. The method according to claim 3, characterized in that Similar 3D data blocks are determined as follows: The distance between any two 3D data blocks is calculated by the following formula: in, Represents the i-th three-dimensional data block With the jth three-dimensional data block , z represents the CEST image, x represents the three-dimensional coordinates in the signal domain X, L represents the side length of the three-dimensional data block, and ‖·‖2 represents the 2-norm; If the distance is less than a preset distance threshold, it is determined that the two three-dimensional data blocks are similar.

5. The method according to claim 3, characterized in that: The adaptive weight used in the adaptive weighted average method is obtained by the following formula: in, represents the image after filtering, x R ∈X represents each three-dimensional data block in the image X before traversal filtering. Represents and references a 3D data block x R Similar to other three-dimensional data blocks in X, Represents x R For each x i The weight of represents the filtered x i The estimated value of , y represents the data before filtering, Represents the characteristic function, which is used to indicate x i Is it within the valid range? represents the number of non-zero coefficients after filtering, and σ represents the standard deviation of the noise.

6. The method according to claim 2, characterized in that Performing denoising processing on the denoised input image based on the spectrum information channel includes: Construct the following objective function; Wherein, X represents the denoised input image, W represents the first weight tensor, and Y represents the denoised image. represents the square of the Frobenius norm, λ1 and λ2 represent regularization parameters, ‖‖1 represents the L1 norm, S represents the second weight tensor, D represents the three-dimensional first-order forward difference operator, ⊙ represents the element-wise product, P represents a local region, and P i represents the index of the i-th local block, ‖‖ w,* represents the weighted nuclear norm, stX = U × 3V represents the constraint: X is reconstructed by the three-dimensional tensor product of the projected data U and the spectral subspace basis V; The denoising input image is denoised based on the objective function.

7. The method according to claim 6, characterized in that After determining that the fused image does not meet the preset termination condition, the method further includes: calculating a first energy of the first denoised image and a second energy of the second denoised image; Obtaining a spatial channel weight and a spectral channel weight according to the first energy and the second energy, so as to obtain a fused image next time; The fused image is obtained by weighting the first denoised image and the second denoised image using the spatial channel weight and the spectral channel weight.

8. The method according to claim 7, characterized in that The spatial channel weight and spectral channel weight are obtained by the following formula: Among them, λ3 and λ4 represent the spatial channel weight and the spectral channel weight respectively, energy1 and energy2 represent the first energy and the second energy respectively, and ∈ represents a constant.

9. The CEST image denoising method according to claim 2, characterized in that: The determining whether the fused image satisfies a preset termination condition includes: Calculating a ratio of a noise level of the fused image to a noise level of the CEST image; If the ratio is less than a preset ratio threshold, it is determined that the fused image meets the preset termination condition.

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