A Phase Preprocessing Method and System Based on a Priori-Guided 3D Convolutional Network

Abstract

A phase preprocessing method and system based on a prior-guided three-dimensional convolutional network, which can perform one-stop processing of integrated phase unwrapping and background field removal on the original phase map, including an encoder module, a decoder module, and an amplitude map correction module. At the network input end, the phase map is cascaded with its three-dimensional gradient map to quickly locate the wrapped boundaries of phase jumps. In the encoder module, multiple groups of Res2Net modules are connected in series to enhance the multi-scale representation ability of the network. The amplitude correction module helps the network accurately distinguish the field map perturbations caused by susceptibility artifacts. Tests are carried out on simulated brain data, healthy brain data, and clinical brain data. Qualitative and quantitative analyses prove that the proposed method can quickly and accurately perform phase preprocessing, especially in areas with complex wrapped distributions, uneven background fields, severe lesions, etc., and accurate field map reconstructions can be obtained.

Description

Technical Field

[0001] The present invention relates to the technical field of magnetic resonance image processing based on convolutional neural networks, and particularly to a phase preprocessing system and method based on a prior-guided three-dimensional convolutional network. Background Technique

[0002] Magnetic Resonance Imaging (MRI) is a medical imaging technique. MRI technology captures and analyzes the responses of different tissues inside the human body or an object to magnetic fields and high-frequency electromagnetic waves to generate high-resolution images. Through magnetic resonance technology, clear magnitude images and phase images can be obtained. Based on these images with different contrasts, the tissue structure information can be deeply analyzed, thereby providing effective auxiliary information for various diseases. Quantitative susceptibility mapping is a post-processing technique of MRI for calculating tissue susceptibility values. Susceptibility is an inherent property of a substance used to measure the degree of magnetization of the substance under the action of an external magnetic field. Quantitative susceptibility technology uses the relationship between the dipole field and susceptibility to inversely solve the susceptibility using the phase image and obtain the susceptibility distribution of tissues, providing auxiliary information for the diagnosis of diseases in the central nervous system, heart, abdominal cavity, etc. Since the susceptibility information of tissues is contained in the phase data of MRI, and due to factors such as phase wrapping and background field interference in the phase data, a series of preprocessing is required to obtain the local field map. Therefore, susceptibility imaging is a complex process involving multiple steps. The basic processing flow of QSM includes phase unwrapping, background field removal, and susceptibility reconstruction. Among them, phase unwrapping and background field removal are called phase preprocessing.

[0003] Under the action of the echo time TE, the dipole field B Δ and the phase relationship can be expressed as

[0004]

[0005] where γ is the gyromagnetic ratio, is the phase value at the spatial position r, is the phase when the echo time is 0, which can generally be ignored. In actual imaging, the acquired MRI signal is in complex form, and the original phase image is obtained by taking the arctangent of the complex MRI signal. The phase information of the signal is obtained using the four-quadrant arctangent function. However, due to the periodicity of the four-quadrant arctangent function, its value range is (-π, π]. The phase value will increase with the increase of the echo time, and the phase exceeding this interval range will be restricted to this interval, resulting in the phase wrapping phenomenon. Therefore, the unwrapped phase value can be expressed as

[0006]

[0007] where is the winding phase and N(r) is the winding number.

[0008] In human brain susceptibility imaging, the perturbation magnetic field caused by the influence of magnetic field inhomogeneity and many strong susceptibility sources outside the brain is called the background field, while the magnetic field generated by the internal tissues of the brain is called the local field. Since the order of magnitude of the background field is two to three orders of magnitude larger than that of the local field, the background field will penetrate deep into the brain and interfere with the local field. That is, the total perturbation magnetic field B Δ in the brain can be regarded as the superposition of the local field B loc (r) generated by the internal tissues of the brain and the background field B bkg (r) generated by external factors:

[0009] B Δ = loc (r) + bkg ()

[0010] Therefore, it is necessary to remove the background field from the unwrapped phase to eliminate the interference of external factors. Since at the junction of tissues and air, such as strong interference regions like the skull and sinuses, the background field has a serious impact. If the background field is not removed cleanly, phase residues will be left in these strong interference regions, thus affecting the effect of susceptibility reconstruction.

[0011] At present, many domestic and foreign teams have proposed many algorithms for phase preprocessing. The commonly used phase unwrapping algorithms are divided into methods based on path tracking and methods based on the Laplacian operator. Methods based on path tracking include Best Path, PRELUDE, SEGUA, etc. Common methods for removing the background field include the dipole field projection method and the method based on the Laplace equation. However, these traditional methods usually introduce time and space limitations for unwrapping, and a trade-off between time and accuracy is required; while for the background field, boundary conditions usually need to be assumed, and suitable optimization parameters need to be selected according to different data characteristics, and the processing time is generally relatively long.

[0012] Nowadays, with the rapid development of hardware devices and artificial intelligence technologies, deep neural networks have also developed and been applied rapidly. Among them, convolutional neural networks are widely used in the fields of image processing and computer vision and have demonstrated powerful performance. Convolutional neural networks have achieved good results in three-dimensional image reconstruction problems. Deep learning and convolutional neural networks can automatically learn the features required in the image processing process without the need for manual design and selection of features, and have achieved better results than traditional methods and machine learning methods in many image reconstruction tasks. With sufficient support of training data, a convolutional neural network structure suitable for the task is designed, and a suitable loss function is selected to train the network parameters, enabling the network to learn the mapping relationship from the original phase image to the local field map and automatically complete the phase preprocessing of three-dimensional images. In recent years, many deep learning algorithms for image reconstruction have been proposed. UNet designed based on FCN connects the feature maps between the encoder and the decoder to reuse features, improves the performance of the network and enhances the efficiency of network training, and is widely used in the field of medical image processing. For the phase preprocessing problem in the QSM field, in 2018, Johnson et al. first proposed a phase unwrapping method based on convolutional neural networks. In 2019, Bollmann et al. applied UNet to the step of removing the background field. However, these networks all need to separately process phase unwrapping and removing the background field, have simple network structures, and do not utilize the amplitude map information of the GRE sequence.

[0013] In summary, phase preprocessing mainly includes two major steps: phase unwrapping and removing the background field. The distribution of phase wrapping is variable, the interference of the background field is complex and is also restricted by the unwrapping result. The strongly interfering regions such as the air-tissue interface regions of the sinuses and skulls are difficult to process. At the same time, traditional algorithms also have disadvantages such as numerous parameters and long processing time. In these aspects, deep learning methods have significant performance advantages and application prospects, but a reasonable convolutional neural network structure needs to be designed to effectively learn phase wrapping and background field information simultaneously. Summary of the Invention

[0014] The main purpose of the present invention is to overcome the above-mentioned defects in the prior art and propose a phase preprocessing method and system based on a prior-guided three-dimensional convolutional network.

[0015] The present invention adopts the following technical solutions:

[0016] A phase preprocessing method based on a prior-guided three-dimensional convolutional network, characterized in that a prior-guided three-dimensional convolutional network is pre-constructed, and the prior-guided three-dimensional convolutional network includes an encoder module, a decoder module, and an amplitude map correction module. The specific processing method is as follows

[0017] S1 Obtain the original phase map, input the original phase map into the encoder module to extract multi-scale information, and obtain the first feature map;

[0018] S2 inputs the first feature map into the decoder module for preliminary reconstruction to obtain a second feature map;

[0019] S3 obtains an amplitude map, inputs the second feature map and the amplitude map into the amplitude map correction module for fusion and feature screening, and outputs a local field map.

[0020] Preferably, the encoder module includes a convolutional layer, a gradient cascading module, and four Res2Net modules arranged in sequence, and downsampling is performed once between every two of the Res2Net modules.

[0021] Preferably, the gradient cascading module first calculates the gradient distribution maps in the x, y, and z directions of the input image, adds them up, and then cascades them with the input image to obtain the output

[0022]

[0023] where represents the input image, Y represents the output of the gradient cascading module, Concat represents cascading, respectively represent performing gradient operations in three directions.

[0024] Preferably, the Res2Net module is used to pass the input through a Conv1 convolutional layer, reduce its dimension to 64 when the number of channels is greater than 64, then evenly divide its channels into four groups, marked as x = {x1, x2, x3, x4}, and the operations performed after grouping are as follows:

[0025]

[0026] where f i represents different Conv3 convolutional layers, and the output corresponding to x i after convolution is marked as y i , after cascading the four groups of outputs y i , and then passing through a Conv1 convolutional layer and the activation function LeakyReLU in sequence to obtain the output of the Res2Net module.

[0027] Preferably, in the Res2Net module, there are three Conv3 convolutional layers that respectively perform convolution on three of the grouped channels. The output of the first Conv3 convolutional layer also skip-connects to the input of the second Conv3 convolutional layer, the output of the second Conv3 convolutional layer skip-connects to the input of the third Conv3 convolutional layer, and the input of the Res2Net module directly skip-connects to the output.

[0028] Preferably, the decoder module includes three residual modules, and through two skip connection structures, the feature maps output by two of the Res2Net modules in the encoder module are respectively cascaded with the feature maps output by two of the residual modules in the decoder module one by one, and the first feature map is restored to the original size through the three residual modules.

[0029] Preferably, in the amplitude map correction module, first, the normalized amplitude map M, the second feature map F p and the feature map C after cascading the two are respectively passed through a Conv3 convolutional layer to obtain the corresponding feature maps of the three and Then, through two sets of processing units composed of a Conv3 convolutional layer and a Sigmoid activation function, the weight distribution matrices W m and W p are obtained respectively:

[0030]

[0031]

[0032] The weight distribution matrix W m is multiplied pointwise with the feature map to obtain the filtered feature map and label it as Similarly, the weight distribution matrix W p is multiplied pointwise with the feature map to obtain the filtered feature map and label it as On this basis, the feature maps and are cascaded, and then passed through a Conv3 convolutional layer and LeakyReLU to obtain the output of the amplitude map correction module

[0033]

[0034]

[0035] The output end of the amplitude map correction module is connected to a Conv3 convolutional layer, so as to compress the number of channels to 1.

[0036] Preferably, the mean square error L mse (θ) is used as the main loss function to guide the three-dimensional convolutional network with prior knowledge,

[0037]

[0038] where T is the batch size of each batch of data, i represents the index of the batch data, G θdenotes a prior-guided three-dimensional convolutional network, θ denotes the trainable variables in the network, and M i and respectively denote the amplitude map and the phase map input to the network, denotes the local field map output by the network, is the label image for local field information used during training; in order to retain the texture structure information in the image as much as possible, the gradient loss L gdl (θ) is added as an additional term,

[0039]

[0040] wherein, denotes the gradient operator for calculating gradients in three directions respectively.

[0041] Preferably, the total loss function of the prior-guided three-dimensional convolutional network is

[0042] L total (θ) = L mse (θ) + λL gdl (θ)

[0043] where λ is a hyperparameter for adjusting the two loss terms.

[0044] Preferably, a phase preprocessing system based on a prior-guided three-dimensional convolutional network is characterized by comprising

[0045] a network construction module for constructing a prior-guided three-dimensional convolutional network, and the prior-guided three-dimensional convolutional network includes an encoder module, a decoder module, and an amplitude map correction module.

[0046] The encoder module is used to extract multi-scale information from the input original phase map to obtain a first feature map;

[0047] The decoder module is used to perform preliminary reconstruction on the input second feature map to obtain a second feature map;

[0048] The amplitude map correction module is used to fuse the input second feature map and the amplitude map and perform feature screening, and output to obtain a local field map.

[0049] From the above description of the present invention, compared with the prior art, the present invention has the following beneficial effects:

[0050] In the present invention, a prior-guided three-dimensional convolutional network with phase preprocessing is utilized to complete a one-step process of integrated phase unwrapping and background field removal. The prior-guided three-dimensional convolutional network is based on a multi-scale convolutional structure and is mainly divided into an encoder module, a decoder module, and an amplitude map correction module. The encoder module includes a gradient cascade module, a Res2Net module, and downsampling, etc. The Res2Net module is used to extract multi-scale information and distinguish different ranges of wrapped distributions and background fields without introducing too many parameters. The gradient cascade module is used to pre-extract wrapped boundaries to help the network quickly locate strong interference regions with severe phase wrapping. The decoder module is composed of residual modules and upsampling; the amplitude map correction module uses amplitude map information to guide field map reconstruction, improving the accuracy and robustness of the network.

[0051] The present invention is compared with other traditional methods and tested on simulated brain data, healthy brain data, and clinical brain data. Through qualitative and quantitative analysis, it is proved that phase preprocessing can be quickly and accurately achieved, especially in regions with complex wrapped distributions, uneven background fields, severe lesions, etc., and accurate field map reconstruction can be obtained. BRIEF DESCRIPTION OF THE DRAWINGS

[0052] To more clearly illustrate the technical solutions in the embodiments of the present invention, the following briefly introduces the drawings required for the description of the embodiments. Obviously, the following drawings are only some embodiments of the present invention. For those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative efforts.

[0053] Figure 1 It is a flowchart of the method of the present invention.

[0054] Figure 2 It is a schematic diagram of the composition of the prior-guided three-dimensional convolutional network of the present invention, where (a) is the overall network structure, (b) is the residual module, (c) is the gradient cascade module, (d) is the Res2Net module, and (e) is the amplitude map correction module.

[0055] Figure 3 It is a comparison of the convergence curves of ablation experiments on different functional modules of the prior-guided three-dimensional convolutional network of the present invention, including PNet_NoGCB with the gradient cascade module removed, PNet_NoR2B with the Res2Net module removed, and PNet_NoMCB with the amplitude map correction module removed.

[0056] Figure 4The following is a comparison of the reconstructed results of the local field maps of simulated brain data from the ablation experiments of different functional modules by the prior-guided three-dimensional convolutional network of the present invention, including PNet, PNet_NoGCB with the gradient concatenation module removed, PNet_NoR2B with the Res2Net module removed, and PNet_NoMCB with the amplitude map correction module removed.

[0057] Figure 5 The local field results of different methods on simulated human brain data; the first column is the original phase map and amplitude map; the upper two rows are the local field map results of the cross-section and coronal plane; the numbers in the figure are the PSNR and SSIM values between this figure and the reference image; the lower two rows are the difference maps between the results of different methods and the reference image.

[0058] Figure 6 The local field reconstruction results of different methods on healthy brain data, with the difference map referring to the iRSHARP method. The first two rows are the processing results of 23 layers, and the last two rows are the processing results of 38 layers.

[0059] Figure 7 The reconstruction results of different methods on the calcified data of the 2019 QSM Challenge, with the first two rows being the coronal views and the last two rows being the sagittal views.

[0060] The present invention will be further described in detail below in conjunction with the accompanying drawings and specific embodiments. Specific Embodiments

[0061] The present invention will be further described below through specific embodiments.

[0062] See Figure 1 , a phase preprocessing method based on a prior-guided three-dimensional convolutional network, which pre-constructs a prior-guided three-dimensional convolutional network. The prior-guided three-dimensional convolutional network includes an encoder module, a decoder module, an amplitude map correction module, etc.

[0063] Among them, the encoder module includes a convolutional layer, a gradient concatenation module (Gradient Concatenation Block, GCB), and four Res2Net modules (Res2Net Block, R2B) arranged in sequence. The convolutional layer is a 3×3×3 convolution, and downsampling is performed once between every two Res2Net modules.

[0064] The decoder module includes three Residual Blocks (RB), and through two skip connection structures, the feature maps output by two of the Res2Net modules in the encoder module are cascaded one-to-one with the feature maps output by two of the residual modules in the decoder module respectively. Through the three residual modules, the first feature map is restored to the original size. Through the two skip connection structures, the feature maps of two of the Res2Net modules in the encoder module are cascaded with the feature maps of two of the residual modules in the decoder module, realizing the effective fusion of the deep semantic features of the decoder and the shallow features of the encoder, strengthening the reuse of features and accelerating the convergence of the network.

[0065] The Magnitude Correction Block (MCB) effectively combines the amplitude map structure information to guide the reconstruction process of the field map. The prior-guided three-dimensional convolutional network finally also includes a 3×3×3 convolution for channel normalization.

[0066] In the embodiments of the present invention, all are three-dimensional convolutional kernels. The input of the network is a phase map and a normalized amplitude map, and the output is a local field map. The number of channels in the first layer of the network encoder module is 32, and the number of channels doubles after each downsampling. The maximum number of channels is 128, and the changes in the number of channels are all shown by numbers in the figure. Since the image reconstruction task is sensitive to the data distribution of the image, batch normalization is not used in the network. The downsampling in the network is implemented by max pooling, and the upsampling is implemented by trilinear interpolation. The network uses LeakyReLU as the activation function.

[0067] The specific processing method of the present invention is as follows:

[0068] S1 Obtain the original phase map, input the original phase map into the encoder module to extract multi-scale information, and obtain the first feature map.

[0069] Among them, the input phase map passes through the gradient cascading module to calculate the phase gradient distribution map, which is cascaded with the phase map and then input into the first Res2Net module to map the image to the feature space. Then, through three Res2Net modules and three downsamplings, while learning the features, the size of the feature map becomes 1 / 8 of the original, and the first feature map is obtained.

[0070] Specifically, refer to Figure 2In (c), the designed gradient cascade module GCB pre - extracts the wrapped boundary of the phase map. Since there are obvious features in phase wrapping, which is caused by the boundary jump of the phase reaching (-π, π], the gradient of the wrapped boundary is usually larger than that of the normal area. However, the network model may not be able to quickly learn this feature. In this application, the wrapped boundary where phase jumps occur can be quickly extracted through the gradient operator, helping the network locate the key area. Then the gradient cascade module first calculates the gradient distribution maps in the x, y, and z directions of the input phase image, and after adding them, cascades with the input image to obtain the output

[0071]

[0072] Where represents the input image, Y represents the output of the gradient cascade module, Concat represents cascading, respectively represent performing gradient operations in three directions.

[0073] Refer to Figure 2 In (d), the wrapped distribution of the MRI phase map is complex and changeable. In the area with a small sinus range, the number of wraps is large, and there is a large - area wrapping phenomenon at the skull position. Therefore, the extraction of multi - scale features is very crucial for the phase processing problem. However, if multiple stacks of 3×3×3 convolutions are used, it will lead to too many parameters or channel redundancy. Therefore, the Res2Net module is introduced. In the Res2Net module, there are three Conv3 convolutional layers that respectively perform convolutions on three of the grouped channels. The output of the first Conv3 convolutional layer also has a skip connection to the input of the second Conv3 convolutional layer, the output of the second Conv3 convolutional layer has a skip connection to the input of the third Conv3 convolutional layer, and the input of the Res2Net module is directly skip - connected to the output.

[0074] Specifically, the Res2Net module is used to pass the input through a Conv1 convolutional layer, reduce its dimension to 64 when the number of channels is greater than 64, and then evenly divide its channels into four groups, marked as x = {x1, x2, x3, x4}. The operations performed after grouping are as follows:

[0075]

[0076] Where f i represents different Conv3 convolutional layers, x i The corresponding output after convolution is marked as y i , after cascading the four groups of outputs y i , and then passing through a Conv1 convolutional layer and the activation function LeakyReLU in sequence to obtain the output of the Res2Net module. S2 inputs the first feature map into the decoder module for preliminary reconstruction to obtain the second feature map.

[0077] Specifically, the first feature map passes through three residual modules in sequence, and an upsampling is performed after each residual module. The inputs of the second and third residual modules are the concatenation of the output of the previous residual module and the feature map of the same size in the corresponding encoder module. After the third upsampling, a second feature map with the same size as the original phase map is obtained.

[0078] In a specific embodiment, referring to Figure 2 (b) therein, the residual module includes two layers of Conv3 and a LeakyReLU activation function in the middle. The result after two convolutional layers and a LeakyReLU activation function is subjected to a residual connection with the input of the residual module, directly connecting the input information to the output, and only the residual texture between the two needs to be learned. The gradient vanishing problem can be effectively avoided through the residual module, improving the learning efficiency of the network.

[0079] S3 obtains the magnitude map, inputs the second feature map and the magnitude map into the magnitude map correction module for fusion and feature screening, and outputs to obtain the local field map. The magnitude map correction module (Magnitude Correction Block, MCB) can effectively combine the structural information of the magnitude map to guide the reconstruction process of the field map.

[0080] Referring to Figure 2 (e) therein, the deep layer of the decoder of the network uses the residual module to perform preliminary reconstruction in combination with the encoder features. After the upsampling is completed, the magnitude map information is introduced to construct the magnitude map correction module MCB to guide the reconstruction process of the local field map. Due to the superposition of background field interference and wrapping phenomena in the phase map, the texture structure features are severely lost. The GRE sequence can obtain rich magnitude map information of tissue without multiple repeated acquisitions, and the texture structure features in the magnitude map are basically the same as those in the local field map. Introducing the magnitude map can supplement the information lost in the phase map due to artifacts or lesion interference.

[0081] In the magnitude map correction module, first, the normalized magnitude map M, the second feature map F p and the feature map C after concatenating the two are respectively passed through a Conv3 convolutional layer to obtain the corresponding feature maps of the three and Then, through two sets of processing units composed of a Conv3 convolutional layer and a Sigmoid activation function, the weight distribution matrices W m and W p .:

[0082]

[0083]

[0084] Weight distribution matrix W m With feature map Perform point multiplication to obtain the filtered feature map and mark it as Similarly, the weight distribution matrix W p With feature map Perform point multiplication to obtain the filtered feature map and mark it as On this basis, the feature map and Cascade, then pass through a Conv3 convolution layer and LeakyReLU to get the output of the amplitude map correction module

[0085]

[0086]

[0087] The output end of the amplitude map correction module is connected to a Conv3 convolutional layer, thereby compressing the number of channels to 1. The present invention introduces the amplitude map, and the network learns and screens the features extracted from the amplitude map, which not only supplements information in areas with severe structural loss, but also enhances the generalization performance of the network.

[0088] In addition, the present invention adopts the mean square error L mse (θ) is used as the main loss function of the prior-guided 3D convolutional network,

[0089]

[0090] Where T is the batch size of each batch of data, i represents the index of the batch data, G θ represents a priori guided 3D convolutional network, θ represents the trainable variables in the network, M i and Respectively represent the amplitude and phase diagrams of the input network, represents the local field map output by the network, is the label image about the local field information used in training; in order to retain the texture structure information in the image as much as possible, the gradient loss L is added gdl (θ) as an additional term,

[0091]

[0092] in, Represents the gradient operator that performs gradient calculations in three directions.

[0093] Preferably, the total loss function of the prior-guided three-dimensional convolutional network is

[0094] L total (θ) = L mse(θ) + λL gdl (θ)

[0095] where λ is a hyperparameter for adjusting the two loss terms.

[0096] Based on this, the present invention also proposes a phase preprocessing system based on a prior-guided three-dimensional convolutional network, including:

[0097] A network construction module for constructing a prior-guided three-dimensional convolutional network, which includes an encoder module, a decoder module, and an amplitude map correction module.

[0098] The encoder module is used to extract multi-scale information from the input original phase map to obtain a first feature map;

[0099] The decoder module is used to perform preliminary reconstruction on the input second feature map to obtain a second feature map;

[0100] The amplitude map correction module is used to fuse and perform feature screening on the input second feature map and the amplitude map, and output a local field map.

[0101] The prior-guided three-dimensional convolutional network is denoted as PNet. PNet takes the original phase map as the input and the local field map as the output, and completes the phase processing process end-to-end, so there is no need for step-by-step training. During the network training process, the batch data size is 8, the network is built using the Pytorch deep learning framework, ADAM is used as the optimizer, and the initial learning rate is 10 -4 , and a total of 30 epochs are trained. After every 10 epochs, the learning rate is reduced to 0.5 of the original. The network parameters are initialized using the Xavier method. The computing platform used is a 12GB NVIDIA GTX Titan X graphics processor, and the time required to train PNet using this platform is approximately 16 hours.

[0102] The healthy human brain data is collected from healthy volunteers using a 7T Philips magnetic resonance imaging device. Among them, some volunteers are scanned with 5 echoes in 5 directions using the GRE sequence, the repetition time TR is 28ms, the first echo time is 5ms, the echo time interval ΔTE is 5ms, the matrix size is 224×224×126, and the imaging field of view is 220×220×126mm 3 . For the other part of the volunteers, 9 echoes are scanned in 4 directions, the repetition time TR = 45ms, TE / ΔTE = 2 / 2ms, the matrix size is 224×224×110, and the imaging field of view is 220×220×110mm 3, a total of 216 three-dimensional brain data. The phase data is unwrapped using the Laplacian-based method, and the background field is removed using the iRSHARP method to obtain the field map. Then, the corresponding COSMOS results are calculated based on the multi-directional data collected, and the gray level range of the amplitude map is normalized to [0,1].

[0103] Since there are no completely correct results as labels for the two steps of phase unwrapping and background field removal in phase preprocessing, healthy human brain data is preprocessed to create a training dataset. First, the local field map calculated by the forward model using the COSMOS multi-directional susceptibility results of each sample and the dipole kernel is used as the label image. Then, a susceptibility source with a random size is set in the sinus region to simulate a strong interference background field area. Additionally, 3 - 5 ellipsoids are randomly placed outside the brain tissue to simulate susceptibility sources, with the radius size randomly distributed in the range [10, 40] and the susceptibility value randomly distributed in the range [0, 1] to generate background fields with different distribution characteristics. The background field is calculated by the forward model using the susceptibility source and the dipole kernel and is superimposed on the label local field map to form the total field map. Referring to the multi-echo imaging time used in actual imaging, four echo times TE = {0.01s, 0.015s, 0.02s, 0.025s} are set for the simulated data, and the corresponding phase maps are derived. The numerical distribution range of all phase images is calculated to the interval (-π, π] through phase wrapping, and together with the amplitude Figure 1 after normalization is used as the network input.

[0104] To save computing resources and expand the training data, during network training, the phase map, amplitude map, and local field map are all made into 64×64×32 three-dimensional chunks, with 10 pixels of overlap reserved between adjacent image chunks. For the three-dimensional imaging data of the brain tissue, the 64×64 image dimension within a layer can usually contain most of the tissue structures in the original data, and 32 layers of data in the slice direction can provide rich inter-layer information while significantly reducing the computing cost. The prepared phase, amplitude, and local field image chunk datasets are used to train the network. 180 three-dimensional brain data are cut into 9000 groups of 64×64×32 data, among which 7500 groups are used as training data and 1500 groups are used as validation data. During network testing, the entire image with a size of 224×224×110 is directly used as the network input.

[0105] The hematoma human brain data comes from 1 volunteer and is collected using a 3T General Electric magnetic resonance imaging device. There are a total of 8 groups of echo data, where the first echo time is 3.6 ms, the echo interval is 5 ms, and the imaging field of view is 256×256×146 mm 3, the matrix size is 240×240×73. Similarly, the phase data of 8 echoes is unwrapped using the Laplacian-based method, and the iRSHARP method is used to remove the background field to obtain the field map. The gray range of the amplitude image is normalized to [0,1].

[0106] The test calcified brain data is the data of the 2019 QSM Challenge, which is a set of susceptibility images with an isotropic resolution of 0.64 mm based on in vivo T1 and the brain map. Among them, the T1 and images are acquired using a 7T instrument, with TR / TEs = 50 / 4 / 12 / 20 / 28 ms and a flip angle of 15°. This data contains two datasets, Sim1 and Sim2. The main differences between them are the presence of calcified regions and the susceptibility contrast between tissues. Both of these datasets contain measurement noise with typical intensities of clinical scans and are simulated using the same regularization method. These two datasets will be used as test data to verify the generalization ability of the network.

[0107] To explore the influence of each functional module in PNet on the reconstruction results, ablation experiments are conducted on different modules, namely PNet_NoGCB marked without the gradient cascade module, PNet_NoR2B without the Res2Net module, and PNet_NoMCB without the amplitude map correction module. The convergence curves of the ablation experiments on the PNet network are as Figure 3 shown. The three curves of the ablation experiments show significant differences from the curve of the original PNet network, proving that each functional module improves the reconstruction effect of the network to a certain extent. Among them, the PNet_NoMCB corresponding to the yellow convergence curve has the highest convergence value, that is, the amplitude map correction module has the most significant improvement on the network. Figure 4 This is the local field map of the 60th layer reconstructed by each ablation experiment network on the simulated brain data. The result of PNet_NoMCB shows a deviation consistent with the tissue structure, indicating the important role of supplementing the structural information of the amplitude map in local field map reconstruction. Secondly, after removing the Res2Net module, in the globus pallidus region of the reconstruction result, PNet_NoR2B shows the most serious error. The multi-scale ability of the Res2Net module can help the network capture the organizational structure information in a large range. In the first five cycles, it can be seen that without using the gradient cascade module, the convergence speed is significantly slower than that of other ablation experiment networks.

[0108] PNet will be compared with two traditional phase preprocessing methods and a CNN-based method, namely Laplacian and VSHARP, Laplacian and iRSHARP, and ResUNet. Since there is currently no one-step neural network method for phase unwrapping and background field removal, and most of the networks used for phase processing are based on the UNet network structure, the ResUNet network that introduces residual learning based on the UNet network is used as a comparison algorithm. The parameter settings of the comparison methods are as follows: the convolution kernel radius of the VSHARP and iRSHARP methods is 6, the truncated singular value decomposition threshold of iRSHARP is 0.1, and the Gaussian kernel weight is 0.25. The initial number of channels of the ResUNet network is 32, and the number of channels doubles after each downsampling. The overall number of parameters is basically the same as that of PNet, and batch normalization is not used. The datasets, loss functions, and environmental settings used for training the two networks are kept consistent.

[0109] In subsequent experiments, PNet and other methods will be tested on healthy brain data, simulated brain data, and patient brain data. The evaluation metrics used include Root Mean Squared Error (RMSE), High Frequency Error Norm (HFEN), Peak Signal Noise Ratio (PSNR), and Structural Similarity (SSIM).

[0110] Figure 5 The local field map results of reconstructing the simulated brain data by different methods and the difference maps with the label images are shown. The first row is the cross-section, and the second row is the coronal section. The numbers in the figure represent the PSNR and SSIM values between the image and the label image. The PNet method achieves the best results, reaching 38.81 and 0.854 respectively. It can be seen from the difference maps that the results of the Laplacian and VSHARP methods have poor contrast and large-scale structural errors. Laplacian and iRSHARP leave obvious phase residues in the sinus region. Compared with the traditional methods, the CNN-based method does not have large-scale errors and is more accurate in processing strong interference regions. However, the ResUNet network has errors in the overall structure and has a deviation consistent with the winding boundary at the skull position, indicating that the ResUNet network cannot effectively identify phase winding and background field information without guidance. In the area pointed by the arrow in the figure, the comparison methods have obvious errors and artifacts in the strong interference region and the organizational structures with large fluctuations in magnetic susceptibility values such as blood vessels. Due to the introduction of the gradient distribution map and amplitude map priors, PNet can obtain a more accurate local field estimation.

[0111] Table 1 Quantitative Comparison of Reconstruction Results of Different Methods on Simulated Brain Data

[0112]

[0113] As shown in Table 1, a quantitative analysis was performed on a set of simulated brain data using four methods. Compared with the other three comparison methods, all indicators of the PNet of the present invention are significantly better than the comparison methods. The root mean square error (RMSE) of PNet is more than 10 percentage points lower than that of the traditional algorithm. The high-frequency error norm (HFEN) of PNet is nearly 3 percentage points lower than that of the second-ranked Laplacian and iRSHARP methods. At the same time, the structural similarity index measure (SSIM) of the other three methods is 0.80, while PNet reaches 0.84. This shows that the local field maps reconstructed by PNet not only have smaller errors in the low-frequency gray and white matter regions of the brain, but also have more accurate reconstruction effects for high-frequency tissue structures.

[0114] Although healthy volunteer data was used as a reference to produce the training data set for the production of simulated data, there must be certain differences in their data distributions. Therefore, before testing the network with in vivo data, an in vivo data set was used for fine-tuning the network. The in vivo brain data used to produce the training data set was used for transfer learning and testing. There were a total of 16,000 groups of 64×64×32 cut data from 180 three-dimensional brain data, and another 1 sample was used for testing. The Laplacian method was used for unwrapping, and the iRSHARP method was used to remove the background field. When performing transfer learning based on the CNN method, the learning rate was set to 10 -5 , and the number of training epochs was 10, to learn the data distribution of the real human brain to a certain extent.

[0115] Figure 6 In [reference], the local field maps calculated by Laplacian, VSHARP, two deep learning-based methods, Laplacian and iRSHARP on healthy brain data were compared. The first two rows had 23 layers, and the last two rows had 38 layers. As can be seen from the figure, in the same individual data, the number of phase wraps and their complexity at different layers also showed significant differences. As indicated by the black arrows in the figure, compared with other methods, PNet had the smallest difference between the local field in the red nucleus and substantia nigra regions and the reference map, resulting in accurate reconstruction. The difference map further showed the area indicated by the red arrow. Conventional convolutional networks had errors at the edges of brain tissue and regions with large susceptibility fluctuations. The distribution of the number of phase wraps in the 38-layer phase map containing the globus pallidus was relatively simple, and its unwrapping and background field removal results also showed the same results as those of layer 23. Due to the large susceptibility fluctuations at the location of blood vessels, there were large residuals in the results of the comparison network in this area. The method of the present invention can better calculate the local field values of blood vessels and the surrounding tissues. The above results show that the PNet network has good reconstruction effects and generalization capabilities.

[0116] Finally, two sets of data from the 2019 QSM Challenge are used to further verify the generalization ability of the network. Since there are certain differences in the susceptibility contrast and numerical distribution between this data and the dataset used to train the network, one set of non-calcified dataset Sim1 is first used for transfer learning of the two networks respectively, and then the data Sim2 containing calcifications is applied for testing.

[0117] As Figure 7 shown, the results are presented using the coronal view and sagittal view that are convenient for observing calcifications. As can be seen from the figure, the calcifications present high-brightness circular artifacts due to the large difference in the susceptibility values from the surrounding tissues. In the results of the three comparison algorithms, the accurate contour of the calcifications is difficult to recover, and there are still large artifacts in the surrounding area, and serious structured errors can be seen in many places in the corresponding difference map. The method of the present invention is slightly superior to other results in terms of the coherence of the calcifications, and there are obvious deviations only at the calcification sites in the difference map. Given that PNet integrates the amplitude map information and can achieve more accurate reconstruction of the calcified area, it proves that PNet has a certain robustness. Although one set of data from the 2019 QSM Challenge is used to fine-tune the network parameters, limited by the fact that this data does not contain calcifications, the network fails to effectively learn the characteristics of the calcifications. If a certain amount of clinical calcification data is used for sufficient transfer learning of the network, the results of phase processing are expected to be further improved.

[0118] The above is only the specific implementation manner of the present invention, but the design concept of the present invention is not limited thereto. Any non-substantial modification made to the present invention using this concept shall fall within the scope of infringement of the protection of the present invention.

Claims

1. A phase preprocessing method based on a prior-guided three-dimensional convolutional network, characterized in that Pre-construct a prior-guided three-dimensional convolutional network, where the prior-guided three-dimensional convolutional network includes an encoder module, a decoder module, and an amplitude map correction module. The encoder module includes a convolutional layer, a gradient concatenation module, and four Res2Net modules arranged in sequence. Downsampling is performed once between every two Res2Net modules. The gradient concatenation module first calculates the gradient distribution maps in the x, y, and z directions of the input image, and after adding them, cascades them with the input image to obtain the output. The Res2Net module is used to pass the input through a Conv1 convolutional layer, reduce the number of channels to 64 when the number of channels is greater than 64, and then evenly divide its channels into four groups, labeled as x = {x1, x2, x3, x4}. In the Res2Net module, there are three Conv3 convolutional layers that perform convolution on three of the grouped channels respectively. The output of the first Conv3 convolutional layer also skip-connects to the input of the second Conv3 convolutional layer, the output of the second Conv3 convolutional layer skip-connects to the input of the third Conv3 convolutional layer, and the input of the Res2Net module directly skip-connects to the output; The decoder module includes three residual modules, and through two skip-connection structures, the feature maps output by two of the Res2Net modules in the encoder module are respectively and correspondingly cascaded with the feature maps output by two of the residual modules in the decoder module, and the first feature map is restored to the original size through three residual modules; In the amplitude map correction module, first, the normalized amplitude map M, the second feature map F p and the feature map C after cascading the two are respectively passed through a Conv3 convolutional layer to obtain the corresponding feature maps of the three and Then, through two sets of processing units composed of a Conv3 convolutional layer and a Sigmoid activation function, the weight distribution matrices W m and W p are obtained respectively: Weight distribution matrix W m is multiplied element-wise with the feature map to obtain the filtered feature map and label it as Similarly, the weight distribution matrix W p is multiplied element-wise with the feature map to obtain the filtered feature map and label it as On this basis, the feature maps and are concatenated, and then passed through a Conv3 convolutional layer and LeakyReLU to obtain the output of the amplitude map correction module The output end of the amplitude map correction module is connected to a Conv3 convolutional layer, so as to compress the number of channels to 1; The specific processing method includes the following: S1 Obtain the original phase map, input the original phase map into the encoder module to extract multi-scale information, and obtain the first feature map; S2 Input the first feature map into the decoder module for preliminary reconstruction to obtain the second feature map; S3 Obtain the amplitude map, input the second feature map and the amplitude map into the amplitude map correction module for fusion and feature screening, and output to obtain the local field map.

2. The phase preprocessing method based on a prior-guided three-dimensional convolutional network according to claim 1, characterized in that, The output of the gradient concatenation module is: Among them represents the input image, Y represents the output of the gradient cascade module, and Concat represents concatenation, respectively representing gradient operations in three directions.

3. The phase preprocessing method based on a priori-guided three-dimensional convolutional network according to claim 1, wherein The operations performed after the Res2Net module is grouped are expressed as follows: where f i represents different Conv3 convolutional layers, and x i is marked as y after convolution. The corresponding output i After cascading the four groups of outputs y i After output cascading, it passes through a Conv1 convolutional layer and the activation function LeakyReLU in sequence to obtain the output of the Res2Net module.

4. The phase preprocessing method based on a prior-guided three-dimensional convolutional network according to claim 1, wherein Adopt the mean square error L mse (θ) as the main loss function to guide the 3D convolutional network a priori, where T is the batch size of each batch of data, i represents the index of the batch data, and G θ represents a prior-guided three-dimensional convolutional network, θ represents the trainable variables in the network, and M i and represent the amplitude map and the phase map input to the network respectively, represents the local field map output by the network, is the label image of the local field information used during training; in order to retain as much texture structure information in the image as possible, the gradient loss L gdl (θ) is added as an additional term. Among them, denotes the gradient operator for calculating gradients in three directions respectively.

5. The phase preprocessing method based on a priori-guided three-dimensional convolutional network according to claim 4, characterized in that The total loss function of the prior-guided 3D convolutional network is \(L\) total (\(\theta\)) = \(L\) mse (\(\theta\)) + \(\lambda L\) gdl (\(\theta\)) where λ is a hyperparameter that adjusts the two loss terms.

6. A phase preprocessing system based on a prior-guided three-dimensional convolutional network, characterized in that Include A network construction module for constructing a prior-guided three-dimensional convolutional network, where the prior-guided three-dimensional convolutional network includes an encoder module, a decoder module, and an amplitude map correction module; An encoder module for extracting multi-scale information from the input original phase diagram to obtain a first feature map; it includes a convolutional layer, a gradient cascade module, and four Res2Net modules arranged in sequence, and downsampling is performed once between every two of the Res2Net modules. The gradient cascade module first calculates the gradient distribution maps in the x, y, and z directions of the input image, and after adding them, cascades them with the input image to obtain the output. The Res2Net module is used to pass the input through a Conv1 convolutional layer, reduce its dimension to 64 when the number of channels is greater than 64, and then evenly divide its channels into four groups. In the Res2Net module, there are three Conv3 convolutional layers that respectively perform convolution on three of the grouped channels. The output of the first Conv3 convolutional layer is also skip-connected to the input of the second Conv3 convolutional layer, the output of the second Conv3 convolutional layer is skip-connected to the input of the third Conv3 convolutional layer, and the input of the Res2Net module is directly skip-connected to the output; A decoder module for preliminarily reconstructing the input second feature map to obtain a second feature map; it includes three residual modules, and through two skip connection structures, the feature maps output by two of the Res2Net modules in the encoder module are respectively cascaded with the feature maps output by two of the residual modules in the decoder module in a one-to-one correspondence. The first feature map is restored to the original size through the three residual modules; The amplitude map correction module is used to fuse the input second feature map and the amplitude map and perform feature screening to output a local field map. In the amplitude map correction module, first, the normalized amplitude map M, the second feature map F p and the feature map C after concatenating the two are respectively passed through a Conv3 convolutional layer to obtain the corresponding feature maps of the three and Then, through two sets of processing units composed of a Conv3 convolutional layer and a Sigmoid activation function, the weight distribution matrices W m and W p are obtained respectively: Weight distribution matrix W m Perform a dot product with the feature map to obtain the filtered feature map and label it as Similarly, the weight distribution matrix W p Perform a dot product with the feature map to obtain the filtered feature map and label it as On this basis, cascade the feature maps and and then pass through a Conv3 convolutional layer and LeakyReLU to obtain the output of the amplitude map correction module The output end of the amplitude map correction module is connected to a Conv3 convolutional layer, so as to compress the number of channels to 1.

Images