System and method for quantitative measurement of physical properties

Abstract

Exemplary methods for quantitative mapping of physical properties, systems, and computer-accessible media may be provided to generate images of tissue magnetic susceptibility, transport parameters, and oxygen consumption from magnetic resonance imaging data using a Bayesian inference approach, minimizing data fidelity terms based on constraints imposed by structural prior knowledge. The data fidelity terms are constructed directly from the magnetic resonance imaging data. Structural prior knowledge may be characterized from known anatomical images using image feature extraction operations or artificial neural networks. Thus, according to exemplary embodiments, systems, methods, and computer-accessible media may be provided for determining physical properties associated with at least one structure.

Description

[0001] STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH AND DEVELOPMENT

[0002] “This invention was made with Government support under Grant Nos. R01 CA181566, DK116126, NS072370, NS090464, NS095562, NS105144 and R21EB024366 awarded by the National Institutes of Health. The Government reserves all rights in this invention.” This statement is provided solely for compliance with 37 C.P.R. §401.14(f)(4) and should not be construed as an admission that this application discloses and / or claims only one invention.

[0003] CROSS-REFERENCE TO RELATED APPLICATIONS

[0004] This application claims priority to U.S. Provisional Application No. 62 / 853,290, filed May 28, 2019, which is incorporated herein by reference in its entirety. Technical Field

[0005] The present invention relates to magnetic resonance imaging, and in particular to a method for collecting signals and performing quantitative measurement imaging of inherent physical properties of materials in magnetic resonance imaging. Background Art

[0006] Quantitative susceptibility imaging (QSM) in magnetic resonance imaging (MRI) has received increasing clinical and scientific interest. QSM shows promise in characterizing and quantifying chemical constituents such as iron, calcium, and contrast agents including gadolinium and superparamagnetic iron oxide nanoparticles. The tissue composition of these compounds may be altered in various neurological diseases such as Parkinson's disease, Alzheimer's disease, stroke, multiple sclerosis, hemochromatosis, and tumors, as well as other diseases throughout the body. Deoxyheme iron reflects tissue metabolic oxygen extracted from the circulation. QSM can reveal new information related to magnetic susceptibility, a physical property of the underlying tissue. Due to their ubiquitous presence in living organisms, iron and calcium actively participate in important cellular processes. Due to their important roles in the musculoskeletal system, QSM is often very useful for studying the molecular biology of iron and calcium by tracking circulatory iron, as well as metabolic activities using iron and calcium as surrogate markers. QSM can also be used to quantify contrast agents, capturing contrast agent transmission in tissues in time-resolved imaging. This can be fitted with physical transport equations to generate quantitative maps of transport parameters. Therefore, accurately mapping iron, calcium, and contrast agent-induced magnetic susceptibility will greatly aid clinical researchers in exploring human structure and function, and clinicians in better diagnosing and providing relevant treatments for various diseases. Summary of the Invention

[0007] Implementations of systems and methods for collecting and processing MRI signals from an object, reconstructing maps of the object's inherent physical properties (e.g., magnetic susceptibility, transmission parameters, and deoxyhemoglobin concentration), and reconstructing multiple contrast images are described below. In some embodiments, MR signal data corresponding to the object can be converted into multi-echo or time-resolved images that quantitatively depict the object's structure and / or composition and / or function. Using this quantitative magnetic susceptibility or transmission map, one or more magnetic susceptibility- or transmission-based images of the object can be generated and displayed to a user. The user can then use these images for diagnostic or therapeutic purposes, such as to study the object's structure and composition and function, and to diagnose or treat various conditions or diseases. Because one or more of the described implementations can result in magnetic susceptibility- and / or transmission-based images with higher quality and / or accuracy than other magnetic susceptibility and / or transmission mapping techniques, at least some of these implementations can be used to improve a user's understanding of a subject's structure and / or composition and / or function, and can be used to improve the accuracy of any resulting medical diagnostic or therapeutic analysis.

[0008] Some implementations described below can be used to perform magnetic dipole inversion from multi-echo images, while allowing the incorporation of preprocessing and artificial neural networks to improve the quality, accuracy, and speed of magnetic susceptibility measurements. Some implementations can be used to accelerate imaging of various tissue contrasts by reconstructing images from undersampled data, while allowing the incorporation of numerical optimization or artificial neural networks to enforce structural consistency. Some implementations can be used to extract physical parameters including velocity, diffusion, and pressure gradients from time-resolved images. Some implementations can be used to extract tissue oxygen extraction fraction from echo time-resolved amplitude and phase images.

[0009] In general, one aspect of the invention disclosed herein improves the accuracy and speed of reconstructing the magnetic susceptibility distribution by using an adaptive preconditioner estimated from the amplitude decay rate R2*.This aspect of the invention solves the practical problem of automatically determining a preconditioner in QSM reconstruction.

[0010] In general, another aspect of the invention disclosed herein improves the accuracy of deep learning (DL) solutions for image reconstruction from noisy, incomplete data, including the ill-posed inverse problem of quantitative susceptibility mapping, quantitative tissue mapping, and multi-contrast image reconstruction, by updating the DL network or DL ​​output with test data. This aspect of the invention improves the fidelity of the artificial neural network to the test data. By updating the trained DL network weights based on the test data, a significant amount of error in the DL reconstruction can be reduced.

[0011] Generally speaking, another aspect of the invention disclosed herein enables extraction of transport parameters, including velocity, diffusion, and pressure gradients, from time-resolved image data capturing contrast agent transmission in tissue by solving the inverse problem of fitting the image data to the transport equation. This aspect of the invention enables automated post-processing of time-resolved imaging of contrast agent transmission in tissue and eliminates the inherent errors of manually selected arterial input functions used in the prior art.

[0012] In general, another aspect of the invention disclosed herein enables extraction of tissue oxygen extraction fraction (OEF) from echo time-resolved (multi-echo) amplitude and phase images, including complex data acquired in a multi-echo gradient echo sequence. This aspect of the invention enables robust OEF estimation from multi-echo complex image data without the use of vascular challenges.

[0013] Implementations of these aspects may include one or more of the following features.

[0014] In some embodiments, cerebrospinal fluid in the ventricles of the brain is segmented based on low R2* values ​​by applying an R2* threshold and connectivity.

[0015] In some embodiments, the preconditioner is automatically determined based on the tissue R2* value.

[0016] In some implementations, the prior information about the magnetic susceptibility distribution is determined based on object structure information estimated from acquired images of the object.

[0017] In some embodiments, the transport parameter is an estimate of the convective velocity or diffusion of blood flowing through the tissue.

[0018] In some implementations, the optimization is performed using a quasi-Newton method, an alternating direction method of multipliers, or deep learning of artificial neural networks.

[0019] In some implementations, the artificial neural network includes a U-Net structure, each convolutional block of which consists of a 3D convolution layer, a ReLU activation, and a batch normalization layer.

[0020] In some implementations, the object (subject) includes the cerebral cortex, putamen, globus pallidus, red nucleus, substantia nigra, or subthalamic nucleus, or the liver in the abdomen, and generating one or more images of the object includes generating one or more images depicting the cerebral cortex, putamen, globus pallidus, red nucleus, substantia nigra, or subthalamic nucleus, or the liver in the abdomen.

[0021] In some implementations, the object includes at least one of a multiple sclerosis lesion, a cancerous lesion, a calcified lesion, an ischemic lesion, or a hemorrhage, and generating one or more images of the object includes generating one or more images depicting at least one of the one or more images of the multiple sclerosis lesion, the cancerous lesion, the calcified lesion, or the hemorrhage.

[0022] In some implementations, the operations further include quantifying a distribution of a contrast agent introduced into the object during contrast-enhanced magnetic resonance imaging based on the one or more images.

[0023] The details of one or more embodiments are set forth in the accompanying drawings and the description below. Other features and advantages will be apparent from the description and drawings, and from the claims. BRIEF DESCRIPTION OF THE DRAWINGS

[0024] Figure 1. Preprocessor P auto Flowchart A. Susceptibility estimation χ generated from the full magnetic field f and soft tissue ROI M. est The magnetic resonance sources inside and outside M are obtained from PDF and LBV+TKD processing. The voxels outside M are based on the distance to the boundary of M. Cluster, median of absolute values ​​of magnetic susceptibility within the cluster χ med Fitted to The cubic decay function of M is represented by the value of R2*. Cluster, median of absolute values ​​of magnetic susceptibility within the cluster χ med Fitted to The median fitting value of the absolute value of magnetic susceptibility in the cluster is taken as P auto The weight of .

[0025] Figure 2. Simulation results for non-ICH cases. Actual magnetic susceptibility χ T (a), Mean square error (RMSE) of the real and reconstructed QSM (b), CN ROI measurement (c) and GP ROI measurement (d) with different CG times. auto The lowest RMSE was achieved among the three methods. We found Produces consistently low estimates compared to P emp and P auto·

[0026] Figure 3. ICH simulation results. Actual magnetic susceptibility χ T (a), Mean square error (RMSE) of the real and reconstructed QSM (b), CN ROI measurement (c) and GP ROI measurement (d) with different CG times. auto The lowest RMSE was achieved among the three methods. After 1000 CG cycles, P emp Produces a significantly larger RMSE compared to Pemp+R2* and P auto· CN low estimate in P emp+R2* was consistently observed in the emp and P auto· For GP measurement, P emp Shown in 1000 th -3% low estimate in CG, and P auto The offset can be ignored.

[0027] Figure 4. ICH simulation results: using P emp , and P auto True QSM and brain TFI reconstruction at 1000 CG cycles (top). The bottom row shows the difference from the true QSM. emp The results show a strong underestimation artifact around the ICH, which leads to a larger RMSE relative to the true brain susceptibility shown in Fig. 3b. This artifact obscures part of the nearby GP.

[0028] Figure 5. Example of preprocessor for a healthy subject using COSMOS as reference. Top row: P emp , and P auto Preprocessor diagram. Bottom row: COSMOS and P emp , and P auto Reconstructed QSM.P emp and P auto The results of whole brain and COSMOS were similar. The contrast between CN and adjacent white matter was It is not obvious compared with other methods.

[0029] Figure 6. QSM.P of 5 ICH patients emp Overestimation artifacts appear around the ICH in the output, and P auto Artifacts are suppressed in the output. The average magnetic susceptibility values ​​within the millimeter-wide layers are indicated in the lower right corner of each map. and P auto Compared with P emp , the increased mean susceptibility value reflects the reduction of overestimation artifact.

[0030] Figure 7. 5 healthy subjects (left) and ROI measurements of CN and GP in 18 ICH patients (right). was found to be under-predicted in CN measurements, compared to P emp and P auto· For GP, P empThere was a significant underestimation in ICH patients compared with P auto and

[0031] Figure 8. Cardiac QSM results of 3 healthy subjects. From top to bottom: Amplitude images in the short axis view, using P emp , and P auto The difference in oxygenation levels between the LV and RV calculated based on the sensitivity difference between the LV and RV is indicated in the lower right corner of each map. The low estimate of P emp and P auto·

[0032] Figure 9. Comparison of non-ICH cases (top) and ICH cases (bottom). QSM (1st column) and R2* Figure (2 The |χ| scatter plot is shown on the right. Note the different vertical ranges. The pink cross shows the absolute median sensitivity of all voxels in the 1-Hz wave R2*χ med The sigmoid function (red curve) was fitted to R2* and χ med relationship.

[0033] Figure 10. Non-ICH TFI Use Simulation results. Using global increase ∈ = 10 -6 ×30 2 Numerical conditions have been improved The RMSE of the loop is increased at the expense of the final RMSE and the increased offset of the CN measurement (the true value is marked by the red line). As a comparison, the preprocessor uses adaptive spatial scaling ∈ = 10 -6 ×P 2 , relative to the default ∈=10 -6 This achieves low RMSE and less ROI measurement error without introducing CN measurement offset.

[0034] Figure 11 The U-Net architecture used in this work. Each convolutional block consists of a 3D convolutional layer (kernel 3×3×3), a ReLU activation layer, and a batch normalization layer. The number of output features is indicated above each block. Max pooling is used for downsampling.

[0035] Figure 12 L-curve for selecting λ2. λ2 = 32 is selected (blue triangle) as the maximum value, which has a lower fidelity cost than MEDI (red dashed line). Note that λ2 = ∞ corresponds to an output χ that is identical to the original output of the network. φ ·

[0036] Figure 13An example axial slice from a healthy subject (top) and a patient with MS (bottom). From left to right: local field, QSM reconstructed using MEDI, MEDI+DL, and DL, respectively. MEDI+DL shows better contrast of venous structures, such as in the zoomed region, compared to either MEDI or DL. These structures are also present in the local field.

[0037] Figure 14 MEDI, DL, and MEDI+DL are ROI measurements of the GP, PU, ​​and CN on QSM generated for healthy subjects (a) and patients with MS (b), respectively. COSMOS measurements are also provided for healthy subjects. The values ​​shown are the averages across all subjects within each group. MEDI and MEDI+DL produce similar results, while DL estimates are underestimated compared to MEDI and MEDI+DL.

[0038] Figure 15 Linear regression and Bland-Altman plots of MS lesion measurements between MEDI and DL (a) and between MEDI and MEDI+DL (b). MEDI and MEDI+DL showed higher correlation, smaller bias, and narrower limits of agreement. Figure 16 .a): Weights of each layer of 3D U-Net before FINE. From left to right: U-Net layers from downsampling to upsampling, each square represents Figure 11 Weights of the layers in the U-Net. b) Absolute weight changes after FINE. The encoder and decoder parts of the U-Net undergo significant changes after FINE. c) Weights in each layer that experience a relative change of more than 10% after FINE. d) Reconstructed QSM after FINE. e) Randomized initialization in the deep image prior and f) the corresponding weight changes, showing that all layers of the U-Net experience significant weight changes. g) Weights in each layer that experience a relative change of more than 10% after the deep prior network update. h) The QSM reconstructed by the deep prior network fails to invert the field into a magnetic susceptibility source map.

[0039] Figure 17 a) Fidelity values ​​and SSIM of various QSM reconstructions in healthy subjects, with COSMOS as the gold standard. b) Fidelity values ​​and SSIM of T2w and T2FLAIR in multi-contrast reconstructions.

[0040] Figure 18.a) Comparison of the QSM of a healthy subject reconstructed by (from left to right) COSMOS, MEDI, DL, DLL2, and FINE. More detailed structures are recovered after performing the fidelity term. Structures in the occipital lobe are more clearly delineated in FINE and DLL2 than in MEDI and DL. b) Axial image from a representative MS patient. From left to right: Local field, QSM reconstructed by MEDI, DL, DLL2, and FINE, respectively. The central vein is better delineated in FINE and DLL2 than in MEDI or DL. The same structures can also be discerned in the local field. c) Axial image from a representative ICH patient. From left to right: Local field, QSM reconstructed by MEDI, DL, DLL2, and FINE. Low-intensity artifacts surrounding the ICH are observed in DL and DLL2 but are suppressed in MEDI and FINE.

[0041] Figure 19 Comparison of T2w reconstructions. From left to right: Fully sampled original image, undersampled k-space reconstructions using DTV, DL, DLL2, and FINE, respectively. First row: Reconstructed images. Second row: Reconstruction error relative to the original image. Third row: Zoomed-in image. The FINE-reconstructed image has a sharp edge at the white matter / gray matter boundary, while the DTV and DL-reconstructed images are overly smoothed, and the DLL2-reconstructed image is noisy.

[0042] Figure 20. Comparison of T2FLAIR reconstructions. From left to right: Fully sampled original image, undersampled k-space reconstructions using DTV, DL, DLL2, and FINE, respectively. First row: Reconstructed images. Second row: Reconstruction error relative to the original image. Third row: Zoomed-in images. The FINE reconstruction has a sharp edge at the white matter / gray matter boundary, while the DTV and DL reconstructions are oversmoothed, and the DLL2 reconstructions have excessive noise.

[0043] Figure 21 a) Post-Gd T1-weighted imaging shows meningioma. The vector diagram of the red box area is as follows Figure 22 b) L-curve plot of Equation 37. The curves are annotated with the corresponding regularization parameter values. λ = 0.003 was chosen for this and the rest of the reconstructions. c) Velocity plots for the three spatial directions x (front to back), y (right to left), and z (back to back).

[0044] Figure 22 . Figure 21 Vector diagram of the x and y velocity components of the ROI in a. The QTM method captures the flow into the tumor (circles on the left side of the image) and out of the tumor (circles on the right side of the image).

[0045] Figure 23: Brain tumors imaged by DCE. Comparison of Kety flow and QTM velocity u in axial slices of brains with high-grade gliomas: a) Post-Gd T1-weighted imaging showing the tumor (red ROI), b) tumors in all patients f Kety and the regression between the velocity |u|, c) Kety's blood flow f Kety , AIF placed close to the tumor, and d) QTM velocity |u|. Linear regression (b) showed good agreement with a coefficient of determination of 0.60.

[0046] Figure 24 : DCE imaging of liver tumors (colon cancer metastases). Comparison of Kety flow and QTM velocity |u| in axial sections of liver with metastases: a) Post-Gd T1-weighted imaging shows tumor edge enhancement (white arrows), b) in tumors of all patients, f Kety and the regression between speed |u|c)Kety's blood flow f Kety and d) QTM velocity |u|, corresponding to a). Linear regression (b) showed good agreement with a coefficient of determination of 0.68.

[0047] Figure 25 B DCE imaging of breast tumors. Comparison of Kety flow and QTM |u| velocity in sections showing breast cancer: a) post-Gd T1-weighted imaging showing the tumor (red ROI), b) in tumors of all patients, f Kety and the regression between the velocity |u|. c) Kety's blood flow f Kety and d) QTM velocity |u|, corresponding to a). Linear regression (b) shows the coefficient of determination Good consistency.

[0048] Figure 26 Differentiation of invasive breast cancer from benign breast lesions: a) QTM velocity showed significant differences (p = .04), but b) Kety blood flow did not show significant differences (p = .12).

[0049] Figure 27: a) The effect of SNR on the sensitivity of calculating Y values ​​based on the initial values. Shown is the relative error between the estimated Y and the true value. Y0 and v0 are the initial values ​​of Y and v, respectively. As the SNR decreases, Y becomes more and more sensitive to the initial assumptions, which will lead to larger errors when the initial assumptions are far away from the true values. The gray box represents the true value. b) Hypothetical comparison between OEF maps obtained using previous QSM+qBOLD and CAT QSM+qBOLD at different SNRs. Across all SNRs, CAT QSM+qBOLD captures low OEF values, while previous QSM+qBOLD is insensitive to low OEF values ​​at low SNRs. The white numbers represent the OEF mean and standard deviation for the entire brain, and the black numbers represent the root mean square error (RMSE)

[0050] Figure 28 : OEF, CMRO2, v, R2, and χ based on QSM+qBOLD (previously QSM+qBOLD) and CATQSM+qBOLD with constant OEF initial value in healthy subjects nb Comparison of images. CAT QSM+qBOLD demonstrates lower noise and more uniform OEF, as well as good CMRO2 contrast in cortical gray and white matter without extremes. Anatomical structures at corresponding locations are shown in T1-weighted images, CBF maps, and susceptibility maps for reference.

[0051] Figure 29 :In the onset OEF, CMRO2, v, R2, and χ in cluster analysis of QSM+qBOLD based on constant OEF initial value (previously QSM+qBOLD) and QSM+qBOLD based on time evolution (CAT) in stroke patients at 24 hours. nb Comparison of Figures. In both the CMRO2 and OEF maps, lesions were more clearly distinguished using CAT QSM+qBOLD. With CATQSM+qBOLD, regions of low OEF were clearly visible and included in the lesion region defined on DWI, but the low OEF regions obtained using the previous QSM+qBOLD method were not visualized and were not included in the lesion visualized on DWI. CAT QSM+qBOLD generally showed lower v values ​​in DWI-defined lesions. The contrast of v in the previous QSM+qBOLD results was similar in appearance to CBF. CAT QSM+qBOLD showed higher R2 and χ than the previous QSM+qBOLD. nb picture.

[0052] Figure 30: Histogram of OEF values ​​for the lesion side and contralateral side of a second stroke patient imaged 12 days after stroke onset. CAT QSM+qBOLD shows a different distribution in the lesion compared to the contralateral side. The lesion shows eight peaks, with the two strongest peaks at 0 and 17.5%, while the contralateral side has six peaks, with the main peak at 35-45%. However, the previous QSM+qBOLD did not show a specific distribution of low OEF values ​​in the lesion, but both the lesion and contralateral sides had bell-shaped distributions (wider on the contralateral side), with peaks at 47% and 49%, respectively.

[0053] Figure 31 :In healthy subjects and stroke patients (after stroke onset Figure 3 Segmentation and OEF maps using different numbers of clusters (K = 1, 5, 10, 15, 20, and x-mean results) in the image segmentation dataset. Different colors represent different clusters in the segmentation. The obtained OEF observations are almost constant when . For healthy subjects, the x-means method selects K = 13, and for stroke patients, K = 16 (indicated in red).

[0054] Figure 32 : OEF, CMRO2, v, R2, and χ between QSM+qBOLD based on constant OEF initial value (previous QSM+qBOLD) and cluster analysis (CAT) QSM+qBOLD in gray matter of healthy subjects (N=11). nb CAT QSM+qBOLD showed smaller mean CMRO2, OEF, and v than previous QSM+qBOLD, but CAT QSM+qBOLD showed higher mean R2 and χ nb *p<0.01 (paired t-test).

[0055] Figure 33 Shown are example processes for a) mapping tissue susceptibility, b) mapping tissue transport parameters, and c) generating multiple contrast images.

[0056] Figure 34 is a block diagram of an example computer system.

[0057] Like reference symbols in the various drawings denote like elements. DETAILED DESCRIPTION

[0058] The implementation of a system and method for collecting and processing MRI signals of an object, reconstructing a map of the object's inherent physical properties (e.g., magnetic susceptibility, transport parameters, deoxyhemoglobin concentration), and reconstructing multiple contrast images is described below. Physical processes in tissue, including magnetism, transport, and relaxation, affect the temporal evolution of MRI signals. Phase data acquired at different echo times in a gradient echo sequence allows the determination of the magnetic field generated by the tissue's magnetic susceptibility, and thus the magnetic susceptibility. Image data acquired in a time-resolved manner during the passage of a contrast agent through the tissue allows the determination of tissue transport properties, including convective velocity, diffusion, pressure gradients, and permeability. Images of various contrasts, including T1, T2, T2*, and diffusion weights, are structurally consistent. Accurate mapping of physical parameters and reconstruction of multiple contrast images are achieved through the profound or incisive use of image data in formulating and solving inverse problems.

[0059] Some disclosed applications overcome limitations in quantitative susceptibility mapping (QSM), where computational speed and accuracy, as well as Bayesian prior information, impact QSM image quality and practical use. In some applications, to reconstruct a magnetic susceptibility map of an object, automatic adaptive preprocessing is used based on the signal amplitude decay rate R2* to accelerate the iterative computational process of finding the magnetic susceptibility distribution that best fits the available field data and tissue structure prior information. In some implementations, an artificial neural network is trained to reconstruct the magnetic susceptibility map of the object. In some embodiments, an artificial neural network is used to improve the tissue structure prior information in the Bayesian method to estimate magnetic susceptibility. The implementation of these improvements makes the QSM technique fast and accurate, extending the use of QSM from the brain to include the skull, bones, heart, and liver, as well as the quantification of deoxyhemoglobin associated with aerobic energy metabolism. The QSM for processing magnetic susceptibility χ(r) is based on minimizing a cost function that fits the measured tissue magnetic field b(r) (unit B0) to the dipole field b=d*χ, and Gaussian noise b is assumed. The sensitivity solution can be estimated under regularization associated with prior information. For the mathematical analysis of the magnetic susceptibility solution, the equivalent partial differential equation can be obtained directly from the Maxwell equations with the Lorentz sphere correction, which is a wave equation for the magnetic susceptibility χ (the z axis is time) and a Laplace equation for the field b. Then, any field data that is incompatible with the dipole field will lead to artifacts based on wave propagation with large values ​​(divergence) at the magic angle cone. The dipole incompatible part, including noise, discretization errors and anisotropy sources, produces artifacts defined by the wave propagator, with the B0 field direction as the time axis. Granular noise will cause stripes to appear, while continuous errors will cause stripes and shadow artifacts to appear. The invention described here uses this mathematical fact to improve the speed and accuracy of QSM. The stripe artifact is characterized by the presence of a wave propagator along the magic angle in k space. The edges of the magnetic field or the complementary magic angle in image space are almost completely different from the tissue edges. Therefore, they can be minimized during numerical optimization by a regularization term based on tissue structural information or an artificial neural network to identify a solution with minimal artifacts or minimal streaking and shadowing. This is the Morphologically Enhanced Dipole Inversion (MEDI) method, which imposes structural consistency between the magnetic susceptibility map and prior knowledge (e.g., derived from amplitude, phase, and other structural images). The tissue structural information is numerically evaluated using edge gradients and the L1 norm in some implementations or using artificial neural networks in others.

[0060] Some disclosed applications overcome limitations in quantitative transport imaging (QTM), where the appropriate physical model for the underlying transport process to fit time-resolved image data affects the performance and accuracy of QTM. Time-resolved imaging of tracers through tissue allows the determination of the tracer concentration time course. For example, contrast agent concentration in time-resolved MRI can be mapped using QSM or other methods. Time-resolved mapping of tracer concentration can be modeled by transport equations that govern conservation of mass and momentum. The conservation of mass or continuity equation describes the temporal variation of local concentration based on the local divergence of mass flux consisting of a convection or drift velocity term and a diffusion term. The conservation of momentum or Navier-Stokes equations describes the flow in tissue based on pressure gradients and other forces. In some implementations, an artificial neural network is used to estimate transport parameters (velocity, diffusion, and pressure gradient) from time-resolved images. In some embodiments, flow in a porous medium is modeled by a velocity profile, and the data fit of the velocity map is solved by using quasi-Newton, the alternating direction method of multipliers, or optimization of an artificial neural network.

[0061] Some disclosed implementations overcome limitations of tissue oxygen extraction fraction (OEF) mapping, where noise in the amplitude and phase data of multi-echo gradient echo (mGRE) sequences challenges the separation of deoxyheme iron in venules and veins from diffusely stored iron in tissue (primarily present in ferritin and its pathological form, hemosiderin). Physical modeling of the amplitude and phase data can separate deoxyheme iron from stored iron, and OEF can be estimated from mGRE data without any vascular challenges. The problem that the inverse solution of this physical model can easily generate noise in the data can be addressed by grouping voxels of similar tissue together. In some embodiments, this grouping is achieved by using automatic cluster analysis of the signal evolution with echo time, or by optimizing quasi-Newton, alternating direction methods of multipliers, or artificial neural networks to impose sparsity in modeling the signal time course.

[0062] Some disclosed implementations overcome limitations in numerical optimization that are important in solving complex and often ill-posed inverse problems in determining physical parameters from image data and reconstructing multiple contrast images from undersampled data. Numerical optimization can be performed using a quasi-Newton method, an alternating direction method of multipliers, or an artificial neural network. Recently, powerful computing power has made artificial neural networks or deep learning (DL) very effective in extracting complex image features that can be recognized by the human eye but are often difficult to express using specific mathematical expressions. A large number of weights in a DL network can be trained to model the complexity of desired but unexpressible image features. However, these DL network weights may not be able to recognize new image features in the test data. In some implementations, the DL output is updated by numerical optimization or by imposing data fidelity by updating the network weights.

[0063] Precise mapping of physical properties improves numerous applications: including brain QSM with reduced shadowing and streak artifacts, using uniform ventricular CSF as an automatic and consistent zero reference without skull stripping; mapping the oxygen extraction fraction of brain tissue; mapping magnetic susceptibilities in the heart, including oxygenation levels and intramyocardial hemorrhage, and mapping blood flow velocity in tumors via dynamic imaging of contrast agents.

[0064] 1. Automatic Adaptive Preprocessor for QSM

[0065] In this section, we describe an automated adaptive preprocessing method that can obtain high-quality QSM from the total field for imaging a wide range of anatomical structures over a dynamic range of magnetic susceptibilities.

[0066] 1.1. Introduction

[0067] Quantitative susceptibility imaging (QSM) has become an increasingly active area of ​​MRI for studying tissue susceptibility properties, as summarized in recent review articles by numerous research groups. Susceptibility is altered in numerous pathological processes, including neurodegeneration, inflammation, hemorrhage, oxygenation, and calcification. QSM has proven to be a powerful tool for lesion characterization, iron deposition quantification, oxygenation measurement, and deep brain stimulation surgical guidance. The fundamental approach to reconstructing tissue susceptibility distribution from measured magnetic fields uses Bayesian inference, leveraging additional structural knowledge to optimally solve the unstable field-susceptibility inverse problem. Current numerical optimization solvers are iterative and require methods to accelerate iterative convergence to achieve target accuracy. These iterative solvers are not sensitive to image content, use generic termination criteria, and may terminate before sufficient convergence at small volumetric lesions with high susceptibility values, thereby introducing severe streak artifacts.

[0068] A preprocessing method has been proposed to accelerate iterative convergence in QSM. Preprocessed QSM can reconstruct the QSM of the entire image region (including bone and air cavities) directly from the total field, without explicitly or implicitly removing the background field. By incorporating R2* information into the preprocessor, preprocessed QSM can reduce the hypointense artifacts associated with intracerebral hemorrhage (ICH) with extremely high magnetic susceptibility. However, previous preprocessors in the literature consist of simple binary images whose values ​​correspond to given anatomical structures or are manually selected.

[0069] In this work, we propose a general automated framework for constructing an adaptive preprocessor from total field and R2* derived from input gradient echo data, eliminating the need for assumptions about the underlying image content or manual selection of preprocessing weights. We evaluate the automated adaptive preprocessor by its reconstruction performance in total field inversion (TFI) in different imaging scenarios, including healthy brain, ICH brain, and cardiac MRI, and compare it with previous empirical methods in the literature.

[0070] 1.2. Theory

[0071] Based on previous automatic zero-reference work on TFI and QSM, we propose the following model for QSM:

[0072] χ * =Py * ,

[0073]

[0074] where χ is the magnetic susceptibility distribution in image space, * is the convolution operation with the magnetic dipole kernel d, w is the noise weight, and f is the total field. is the gradient operator, M G is the edge mask derived from the magnitude image. The L2 regularization in the last term is used for smoothness in a predefined M2 region, such as the QSM in the cerebrospinal fluid (CSF). Calculate the mean value in M2. P is a preprocessor, currently set to a binary matrix:

[0075]

[0076] Here r indicates the position of each voxel, and M is the tissue mask ROI. S >1 is determined empirically for each anatomical structure and application scenario. Based on R2* information, P can be constructed as:

[0077]

[0078] here This is the area with low R2* The difference in these weights (1 vs. PS ) represents the contrast between the weak magnetic susceptibility of soft tissue and sources of relatively strong magnetic susceptibility, including air, bone, and hemorrhage.

[0079] To fully exploit the power of preprocessing to accelerate the reconstruction of QSM for various image anatomy, P needs to be extended from binary values ​​to a continuous weight range. We propose to automatically generate an adaptive preprocessor from the total field f and R2*, as shown in Figure 1. First, an approximate susceptibility map is quickly estimated from the field input f. Strong susceptibility sources outside the ROI are estimated using the projection to dipole field method (PDF), while the susceptibility inside is calculated using the Laplace boundary value method (LBV) and truncated K-space partitioning (TKD). These two susceptibility components are combined into a single map χ est .

[0080] The noise and artifacts are then addressed by averaging and fitting as follows to generate the weights for the preprocessor P: We group voxels with similar spatial or R2* properties based on their location relative to the ROI and calculate χ est The median of the absolute susceptibility values ​​within each group, χ med Within M, voxels are sorted according to the R2* value. The width of each group is 1 Hz, and χ med By fitting |χ est | to the sigmoid function to determine:

[0081]

[0082] Here [σ1, σ2] controls the output range and (s1, s2) controls the shape. The selection of R2* groups and the determination of the sigmoid function are inspired by the positive correlation between R2* and magnetic susceptibility.

[0083] Outside of M, due to the lack of R2* information, voxels are based on the distance to the ROI to group them, each group has a width of 1 mm, and χ med By fitting |χ est | to the inverse cubic function to determine:

[0084]

[0085] Here σ0 and r0 control the amplitude and decay rate respectively.

[0086] The final P structure is:

[0087]

[0088] Note that, consistent with previous work, we divide P by σ1 so that the weight of soft tissue (low R2*) is close to 1. The resulting preconditioner is fed into TFI (Eq. 1) and solved using either the quasi-Newton method combined with a conjugate gradient (CG) solver or the alternating direction method of multipliers (ADMM).

[0089] Materials and Methods

[0090] We design a numerical simulation to analyze the proposed automatic preconditioner P auto (equation ) for healthy brain (non-ICH) and intracerebral hemorrhage (ICH), and its empirical corresponding method P emp (Equation 2) and (Eq. 3) for comparison. We then compared their performance on in vivo brain QSM in healthy subjects and ICH patients. Finally, they were tested on cardiac QSM in healthy subjects.

[0091] Implementation details: For P emp (Equation 2) and (Equation 3), for the brain / head QSM, the weight P S Select 30, for cardiac QSM, weight P S Choose 20. In estimating the approximate solution χ est When PDF is used CG iterations. For internal sources, LBV was performed with a 3-voxel erosion of the ROI mask to exclude the influence of noise at the boundaries, followed by TKD with a threshold of 0.2. In the weight generation stage, nonlinear least squares fitting of the parameters (σ0, σ1, σ2, r0, s1, s2) was performed using Matlab's fmincon. The distance map D was calculated using the Euclidean distance transform with variable data aspect ratio. For brain / head QSM, λ １ = 0.001 and λ2 = 0.1 are used in (Equation 1). M2 is chosen to be the ventricular cerebrospinal fluid, which can be automatically generated from the multi-echo amplitude image. The reference magnetic susceptibility value is set to the average magnetic susceptibility within M2. For cardiac QSM, λ1 = 0.001 and λ2 = 0 are used. This paper implements the Quasi-Newton Conjugate Gradient (GNCG) solver. When calculating the derivative of the L1 norm in Equation 1, the numerical condition parameter ∈ is replaced by

[0092]

[0093] To avoid zero division, ∈ = P 2 ×10 -6 is chosen so that the numerical condition weight for weak susceptibility is close to 10 -6 (P≈1), while the weight increases at strong magnetic susceptibility (P>1), thereby improving the GNCG convergence speed.

[0094] Simulation experiment. This paper constructs a digital head model with a size of 180×220×128 from the Zubal model. The magnetic susceptibility of different brain tissues is known: white matter (WM) Gray matter (GM) Thalamus (TH) 0.073ppm, caudate nucleus (CN) 0.093ppm, putamen (PU) 0.093ppm, globus pallidus (GP) 0.193ppm, superior sagittal sinus (SSS) 0.27ppm, cerebrospinal fluid (CSF) 0ppm. Air (9ppm), muscle (0ppm) and skull (-2ppm) are distributed outside the brain (Figure 2a). The model was then used to simulate two different scenarios: (A) non-ICH subjects and (B) intracerebral hemorrhage (ICH) patients, characterized by a spherical magnetic susceptibility source of 2ppm located inside the brain mask (Figure 3a). In both cases, the total field was calculated from the true magnetic susceptibility map using the forward model: f=d*χ. Gaussian white noise (SNR=100) was added to the field. R2* was simulated for each tissue type: WM20Hz, GM 20Hz, TH 20Hz, CN 30Hz, PU 30Hz, CSF 4Hz and ICH 100Hz. For each iteration i, the estimated brain magnetic susceptibility Mχ i and the true value Mχ T The root mean square error (RMSE) between the two and the average magnetic susceptibility of CN and GP were used as measures of reconstruction accuracy.

[0095] In vivo experiments: Healthy brains. Five healthy subjects were scanned at 3T (GE, Waukesha, WI) using multi-echo GRE with monopolar readout. The imaging parameters were: TR=39ms, ΔTE = 4.6 ms, acquisition matrix = 512 × 512 × 144, voxel size = 0.5 × 0.5 × 1 mm 3 , BW = ± 62.5kHz, total scanning time 7min. For COSMOS reconstruction of brain QSM, The scans were repeated in different directions, with brain masks generated by BET. In one direction, the total field was estimated from the multi-echo GRE image, followed by phase unwrapping SPURS based on graph cutting. The ROI mask for the entire head was determined by thresholding the amplitude image: M = I > 0.15 × I max , here I max The maximum value of the amplitude map I is used to estimate the R2* map using ARLO. With different preprocessor selections P emp (Equation 2), (Equation 3) and P auto (equation ) was used to reconstruct a whole-head QSM. During QSM reconstruction, CN and GP susceptibility was measured within manually drawn ROIs. Similar measurements were also performed for COSMOS. For each preprocessor choice, a linear regression between TFI and COSMOS susceptibility values ​​was performed for all voxels within the brain.

[0096] In vivo Experiment: ICH Brain. Five patients with intracerebral hemorrhage (ICH) were scanned at 3T (GE, Waukesha, WI) using multi-echo GRE with monopolar readout. The imaging parameters were: FOV = 24cm, TR=49ms, ΔTE = 5 ms, acquisition matrix = 512 × 512 × 64, voxel size = 0.47 × 0.47 × 2 mm 3 , Scan time 4 min. Total field estimation from multi-echo GRE images, followed by graph cut-based unwrapping algorithm SPURS. R2* maps were calculated using ARLO. emp (Equation 2), (Equation 3) and P auto (equation The TFI of 20 μg / cm2 was used to reconstruct the QSM of the brain, where the brain ROI was determined by BET. During the QSM reconstruction, the magnetic susceptibility of CN and GP was measured using manually drawn ROI. The magnetic susceptibility around the ICH site was calculated for each QSM. The mean magnetic susceptibility within wide slices (segmented by an experienced radiologist) was used to quantify low-intensity artifacts surrounding ICH.

[0097] In vivo: 2D multi-echo GRE in the short axis view using breath-hold electrocardiogram triggering Cardiac MRI was performed on three healthy subjects using a GE Healthcare MRI scanner (GE Healthcare, Waukesha, WI). The imaging parameters were: TR=23ms, ΔTE = 2.2 ms, in-plane voxel size = 1.25 × 1.25 mm 3 , Number of slices = 20. Flow correction is applied in readout and slice directions. Acceleration factor Each time you hold your breath, you scan one slice. The average scan time was approximately 12 minutes. The total field was estimated from the multi-echo GRE images using the graph cut-based unwrapping algorithm SPURS. The R2* map was calculated using ARLO. emp(Equation 2), (Equation 3) and P auto The TFI (Equation 7) was used to reconstruct the QSM. Using manually drawn ROIs, the magnetic susceptibility difference between the right ventricle (RV) and left ventricle (LV) was measured on the QSM for each method. The blood oxygenation difference, ΔSaO2, between arterial blood (LV) and venous blood (RV) was then calculated from the magnetic susceptibility difference, Δχ (in ppm), using the following relationship:

[0098]

[0099] 1.4. Experimental results

[0100] Simulation experiment

[0101] The results from the non-ICH scenario are shown in Figure 2. All three preconditioner selections achieved RMSE < 0.005 ppm at 1000 CG iterations, while P auto The number of CG iterations required for each preprocessor to converge to within ±5% of the true value of each ROI is: for CN, 170 (P emp ), and 180(P auto )( Figure 2 c) For GP, (P emp ), and 170(P auto )( Figure 2 d).

[0102] The results of the ICH scenario are shown in Figure 3. For the ICH TFI, and P auto The RMSE was less than 0.005ppm in 1000 CG iterations, and P auto The RMSE of P is low, while emp It failed to fall below 0.010 ppm (Figure 3b). This can be seen from the low-intensity artifacts around the ICH, namely P emp The results of 1000 CG iterations are shown in Figure 4. The number of CG iterations required for each preprocessor to converge to the true value of each ROI within ±5% is: CN, 90 (P emp ), and 90(P auto )(Figure 3c); GP, (P emp ), and 110(P auto )( Figure 3 d).

[0103] In vivo brain imaging

[0104] picture An example of whole-head TFI reconstruction of a healthy subject is shown. Note that, When using the R2* threshold (30Hz), regions with higher R2* (such as GP) are given the same weight as the background region (P S =30). On the contrary, P auto Weights are adaptively assigned to different regions according to their R2* values ​​(Equation 7). Also shown is the QSM reconstructed by COSMOS using P emp , and Paut o Reconstructed TFI. Using P auto The TFI full-head QSM results are visually consistent with COSMOS and P emp At the same time, the gray / white matter contrast The results are not very clear, especially at the CN boundary. CN magnetic susceptibility measurements also confirm this: and The ROI measurements of CN and GP for all five healthy subjects are summarized in Figure 7a. The mean difference relative to COSMOS is: and GP, -9% (P emp ), and -4% (P auto The slope of the linear fitting of brain tissue magnetic susceptibility between TFI and COSMOS for the five subjects is P emp 0.77, and

[0105] The QSM in ICH patients is shown in Figure As shown. emp Low-intensity artifacts around the ICH site are and P auto Suppressed. Around the ICH site The average magnetic susceptibility within the layer within the range is shown in the lower right corner of each QSM, showing the relative P emp of promote and Improvement (P auto ). However, compared to P emp and P auto , The CN susceptibility was significantly underestimated, as shown by ROI measurements (Fig. 7b): Among patients, the mean CN measurement was and The average GP measurement is and 0.211ppm(P auto ).

[0106] Heart or

[0107] Figure 8 shows how this is done by using P emp , and P auto Cardiac QSM in 3 healthy subjects reconstructed using TFI. RV-LV magnetic susceptibility measurements compared to P emp and P auto , The difference in estimated blood oxygen levels was low: 14.7% to 23.4% (P emp ), and The reference value of blood oxygen difference in healthy subjects is 18.8%.

[0108] 1.5. Conclusion

[0109] This paper describes a novel algorithm for automatically constructing a general adaptive preconditioner for QSM, using total fields and R2* derived from input gradient-echo data. This automated adaptive preconditioner overcomes the limitations of binary values ​​and manual selection in previous implementations of preprocessed QSM. Our data demonstrate that the automated adaptive preconditioner achieves the lowest error metric in numerical simulations compared to previous R2*-based preconditioners. The automated adaptive preconditioner suppresses hemorrhage-related artifacts while preserving surrounding brain tissue contrast.

[0110] An important benefit of preprocessing QSM is the elimination of errors in conventional QSM reconstruction, which involves two separate steps: background field removal and local field inversion. Background field removal aims to resolve the strong contrast between background magnetic susceptibility sources (such as air and bone) and local sources including parenchymal tissue. The fitting process is inherently separate and therefore susceptible to error propagation from the unresolved background field to the local field. To address this error propagation, Laplace-based methods have been proposed in the literature to facilitate implicit background field removal. They avoid the separate fitting of background and local fields by applying the Laplace operator to the forward signal model, which essentially eliminates the harmonic background field components. However, the Laplace implementation requires eroding the region of interest (ROI) due to unreliable phase measurements at the ROI boundaries, thus preventing the complete delineation of the brain parenchyma.

[0111] In the field of numerical optimization, preconditioning methods for accelerating conjugate gradient (CG) solvers have been well developed. Traditional methods involve spectral analysis of the eigenvalues ​​of the system matrix, but it is challenging to directly apply this approach to the TFI reconstruction problem due to the huge problem scale and the inherent convolution operation. Although recent work has shown that a preconditioner can be constructed in k-space to approximate the inverse of the dipole kernel and the gradient operator, its implementation in TFI is not intuitive due to the SNR weighting and morphological constraints of Equation 1. At the same time, another aspect of research has focused on developing preconditioners by exploiting prior information about the unknown solution, resulting in a branch called prior enhanced preconditioning. In this regard, the behavior of the preconditioned CG is related to the generalized Tikhonov regularization:

[0112]

[0113] The principle idea is that we can simulate the inverse covariance matrix Γ H F≈∑ -1 To construct the whitening operator Γ, if the prior distribution Satisfies Gaussian distribution. Studies have shown that if the preprocessor is selected as P -1 ≈Γ The following problems will converge faster in a limited number of iterations:

[0114]

[0115] This provides a guideline for choosing a preprocessor P: the transformed variable P -1 χ should have unit variance: In other words, P should be proportional to the estimated susceptibility magnitude of each voxel. Following this guideline, previous TFI work used a binary preprocessor in which an empirical weight of 30 was assigned to sources of strong susceptibility, such as air, bone, and hemorrhage, while a weight of 1 was assigned to soft tissue. This work improves the preprocessor to allow greater flexibility in depicting the spatial distribution and dynamic contrast of susceptibility maps. It uses PDF and LBV+TKD to efficiently generate an approximate solution χ from the field input. est , and use its absolute susceptibility value to construct the preprocessing weight. It should be noted that compared to directly using χ est Within each voxel, we grouped voxels by their spatial (distance to the object boundary) or relaxation properties (R2*) and calculated the median of the absolute susceptibility values ​​observed within each group. This approach aims to mitigate the effects of noise and artifacts within individual voxels, and the median was chosen because it is more robust to outliers.

[0116] The R2* information is essential in our preprocessing technique because it describes the location of strongly magnetized tissues in the subject, such as ICH. ), R2*( and P auto ) significantly suppresses low-intensity artifacts around the ICH site. Otherwise, QSM may suffer from qualitative and quantitative losses, especially for structures close to the ICH, as observed in Figure 3d where GP is underestimated. Previous work has adopted a simple binary thresholding framework to distinguish weak / strong sources by their R2* and give the same weight as the background to high R2* regions. The method proposed in this paper incorporates R2* contrast into the construction of the preprocessor in a more adaptive manner. After we obtain a rough estimate of the tissue magnetic susceptibility, the sigmoid function (Eq. ) to simulate R2* and the observed absolute susceptibility value χ med This monotonic relationship incorporates the prior knowledge that higher R2* values ​​generally imply larger absolute susceptibility values, as is the case with hemorrhage (strong positive values) or calcification (strong negative values). Another view is to consider the susceptibility of each voxel as a random variable that follows a conditional probability distribution based on the value of R2*. use As the preprocessing weight, P -1 χ can have unit variance. It is worth noting that if σ1=1,σ2=30,s1=30,s2<<1,and the sigmoid function (Formula ) will be simplified to The hard threshold used. Therefore, the proposed preprocessor can be considered a generalized version of its empirical counterpart. Figure 9 shows how the fitted sigmoid function explains the different dynamic contrasts of non-ICH and ICH subjects: ICH data produce a higher sigmoid curve compared to non-ICH brain, which corresponds to an increased population of ICH voxels. High R2* - high susceptibility regions. This built-in scalability preserves the ability to suppress ICH-related artifacts (Figures 4 and ), and solved the previous The problem of underestimation of CN structure in the results.

[0117] In this work, the numerical tuning parameter ∈ is scaled using a preconditioner that assigns higher weights to sources with strong magnetization. To evaluate its impact on the reconstruction performance, we use three choices: ∈ = 10 -6 ,∈=P 2 ×10 -6 (proposed) and ∈=30 2 ×10 -6 The TFI experiment was repeated. Figure 10 shows: ∈ = 10 -6 ×30 2 Achieved ratio ∈ = 10 -6 forward We observe a trade-off between convergence speed and accuracy for spatially averaged ∈ 10. Meanwhile, ∈ = 10 -6 ×P 2 The convergence of TFI continues to improve up to the 4000th CG in terms of lower RMSE and smaller CN measurement deviation. To understand this, we can observe the performance of different stages of GNCG in the presence of non-uniform ∈. In the early iterations where the solution is close to 0, we have Then the concept of numerical regulation ∈ is closely related to L2 regularization This suggests that the ∈ scaling using the preconditioner essentially applies different regularization levels to voxels with different susceptibilities in the early GNCG. On the other hand, since this scaling is most significant in strong susceptibility sources (P>>1), which only consist of a small part of the object, it does not significantly bias the final susceptibility estimate (CN: 92.34 ppb (∈=10 -6 ), (∈=P 2 ×10 -6 ) and 89.94ppb(ε=30 2 ×10 -6 ).

[0118] The proposed automatic preprocessor construction method is easily applicable to QSM problems outside the head. Previous empirical selections required a new preprocessor weight search for each anatomical structure and application. This is avoided in the proposed method, which results in a low reconstruction RMSE (Figures 3a and 3b) and reliable ROI measurements (Figure 7a). When applied to cardiac QSM, the proposed LV-RV oxygenation level difference estimation is reasonable. With the Experience Preprocessor The measured results were 18.8% closer to the reported levels.

[0119] The current work is still limited in the following aspects, and possible solutions are proposed here for further research:

[0120] 1) In Eq. Constructing the preprocessor in

[15] requires R2* information, which is not available in single-echo acquisitions. One solution is to change the sigmoid function in Equation 7 to a simple uniform value, which can be the median or average of the absolute magnetic susceptibility of all tissues within region M, but the ability to suppress ICH-related artifacts may be reduced. Other alternative solutions may be to use image contrast such as T1w or T2w, as Equation 7 2) The preconditioner structure depends on the approximate map of the magnetic susceptibility map χ est Currently, we use PDF and LBV+TKD to estimate the susceptibility sources outside and inside the region, respectively. This crude estimate of the susceptibility still contains a lot of noise or artifacts (e.g., due to the low signal-to-noise ratio within the ICH). One can improve the susceptibility by providing a χ that is more robust to noise and artifacts. est For example, the final reconstructed QSM of TFI can be used as χ est A new round of TFI reconstruction is initiated with the updated preprocessor, at the cost of doubling the overall reconstruction time. 3) The skull mask used in this work is obtained using a fixed thresholding method. Automatic skull stripping has been extensively studied in neuroimaging and can be introduced here. 4) The performance of the proposed preprocessing method will be evaluated on more image content, such as carotid artery, musculoskeletal, and animal QSM.

[0121] In summary, this article describes an automated adaptive preprocessor for QSM reconstruction that allows full-field inversion without background field removal. The automated adaptive preprocessor uses the total field and R2* derived from the input gradient-echo data to adapt to the actual susceptibility content. The automated preprocessor improves the convergence speed and accuracy of reconstruction iterations and produces reliable susceptibility measurements in a variety of pathologies, particularly suppressing hemorrhage-related artifacts while preserving surrounding normal tissue contrast.

[0122] 2. MEDI for QSM based on deep learning priors

[0123] In this section, we propose a deep learning-based Bayesian QSM regularization term. This approach can provide excellent anatomical detail while maintaining susceptibility quantification accuracy. Artificial neural networks, including convolutional neural networks, are trained to identify target susceptibility structures by imposing structural consistency between target susceptibility maps and known gradient echo image data.

[0124] 2.1. Introduction

[0125] Quantitative susceptibility reconstruction (QSM) has been an active research area due to its ability to describe magnetic susceptibility, an intrinsic tissue property that is closely related to various pathologies, including demyelination, hemorrhage, and calcification. The core problem of QSM reconstruction is a Bayesian optimization problem to solve for a magnetic susceptibility distribution consistent with inhomogeneous field measurements after dipole convolution. At the same time, due to the instability of dipole inversion, prior structural information is required to suppress streaks and other artifacts. However, explicit regularization solidifies the prior knowledge in image features. For example, fine venous structures may be suppressed due to the smooth or blocky uniformity assumption of the L2 or L1 regularization term in the prior knowledge.

[0126] Recently, researchers have proposed using convolutional neural networks (CNNs) in deep learning for image feature representation for QSM reconstruction. CNNs directly map local fields to susceptibility distributions, bypassing traditional iterative optimization and reducing computational cost to a single network forward pass. The model can be trained on labeled datasets, such as healthy subjects with COSMOS or digital phantoms, both of which demonstrate the feasibility of neural network QSM solutions. However, this approach does not ensure consistency between the measured tissue field and its corresponding CNN susceptibility output, and this deficiency can become significant when transferred to different populations.

[0127] In this work, we propose to incorporate a data fidelity term into CNN-based QSM, thereby combining traditional optimized reconstruction with CNN reconstruction: the CNN output is introduced as Bayesian prior knowledge into the cost function of the QSM reconstruction problem. We obtain preliminary results with this combined fidelity and CNN approach, showing improvements over QSM.

[0128] 2.2. Theory

[0129] The fundamental inverse problem in QSM reconstruction is to obtain the magnetic susceptibility distribution χ from the measured tissue field f. Given a sample with known magnetic susceptibility, a neural network φ(·) can be trained to simulate the inversion process, i.e., to learn the mapping from field to magnetic susceptibility:

[0130] χ φ =φ(f)

[11]

[0131] However, this network does not consider the fidelity between the output and the measured field. In this work, we adopt a Bayesian reconstruction framework to integrate the results of the neural network into the optimization problem:

[0132]

[0133] The first term is a data fidelity penalty in conventional QSM reconstruction with noise weight w and dipole kernel d. The second term is an L2 regularization to enforce the desired structural similarity or penalize the estimated graph χ given by the network φ. φ The difference between χ and χ. Here, E is a linear transformation that can be a unit operation or a spectral filter to enhance specific spectra in the solution. The design details of E are described in the Methods section.

[0134] The Bayesian principle of this L2 regularization is E χ The prior on the iid Gaussian at each voxel A similar Bayesian framework can be found in Morphologically Enabled Dipole Inversion (MEDI):

[0135]

[0136] where the L1 norm of the image gradient weighted by a binary edge mask is used for regularization. This consists of a Considering the necessity of zero reference for susceptibility quantization, MEDI+0 imposes an additional L2 regularization to ensure the consistency of CSF:

[0137]

[0138] Here MCSF is the mask of cerebrospinal fluid CSF, and is the average magnetic susceptibility. Similarly, the CSF regularization term can be added to Equation 12, leading to the following final cost function:

[0139]

[0140] In this work, the proposed method (Formula ) is called MEDI+DL. It will be combined with the original result of deep learning network (called DL) φ , and compared with MEDI+0 (Equation 14), which is referred to as MEDI throughout this work for brevity.

[0141] Materials and methods

[0142] All data in this work were acquired using a 3D multi-echo GRE sequence on a 3T GE scanner (voxel size 0.5 × 0.5 × 3 mm). 3 , field of view 24cm, dTE 0.48ms, number of layers Matrix size 512 × 512 × 50-60. Twenty healthy subjects and eight patients with multiple sclerosis were scanned according to a protocol approved by the Institutional Review Board. Local tissue fields were estimated by multi-echo phase fitting, followed by phase unwrapping and background field removal. Different orientations and 1×1×1mm 3 Voxel size. Four of the 20 healthy subjects were scanned repeatedly for COSMOS reconstruction and quantitative comparison.

[0143] Neural network settings. U-Net, a fully convolutional neural network, was chosen as the target network architecture. We used a 3×3×3 kernel, 1 voxel stride, and reflection padding in each convolutional layer. Detailed parameters regarding the number of layers, up / down sampling factors, and feature size can be found in Figure 11 The network input / output is designed for 128×128×2 local image patches, while the original 3D volume field map uses the local image patches between adjacent local patches. Overlapping solutions are split into small blocks. Each local block of the network output is recompiled to restore the full volume size.

[0144] The network was trained on data from 20 healthy subjects, with tissue fields as input and QSM as output. 4199 local patches were extracted from all 20 cases, and the target QSM was reconstructed using MEDI. The training / validation / test data ratio was 0.7 / 0.1 / 0.2. The L1 norm of the difference map between the estimated QSM and the true MEDI result was chosen as the training cost function. The Adam optimizer (learning rate 0.001, maximum iteration 80) was used to optimize the trainer. This work was implemented on a platform with a Core i7 CPU, 128GB of memory, and 4 GTX TITAN XP GPUs (12GB each).

[0145] Iterative reconstruction: MEDI+DL. The DL result χ estimated by the network φ It is inserted into the formula in the form of L2 regularization We propose a high-pass filter E as a linear transformation before the penalty, which is implemented by convolution with a spherical kernel: E(χ) = χ - S * χ. This penalty emphasizes the difference χ - χ φ The reason for this is that since our neural network is patch-based, it cannot capture spatial variations beyond the range of a single patch. Therefore, regularization allows for smooth variations across the entire range by high-pass filtering the difference map. The radius of S is chosen to be 10 mm.

[0146] MEDI+DL(formula ) was optimized using the quasi-Newton conjugate gradient method. The L-curve (see Figure 12 ) was used to determine the regularization strength λ2 = 32, so as to achieve a fidelity similar to MEDI. MEDI+DL was applied to 4 healthy subjects (untrained) and 8 MS patients. For comparison, MEDI (Equation 14)λ １ = 0.001 and λ CSF = 0.1 was also used for the same data.

[0147] Quantitative Analysis. For healthy and patient subjects, regions of the three major gray matter structures (SGM): the globus pallidus (GP), putamen (PU), and caudate nucleus (CN) were manually segmented, and the mean magnetic susceptibility was measured within each segmented region. In addition, the magnetic susceptibility of multiple sclerotic lesions was measured in manually drawn ROIs in each MS patient, referenced to normal-appearing white matter on a cerebroscopic image. Linear regression and Bland-Altman analysis were used to quantify the agreement between the different methods.

[0148] Results

[0149] Healthy subjects. The QSM reconstruction of one subject is shown in Figure 13 (Top row). Compared with MEDI or DL, MEDI+DL shows the venous structure more clearly, e.g. Figure 13 The yellow box highlights the part. The same structure is also reflected in the local field ( Figure 13 center, left column).

[0150] The ROI measurement results of SGM can be found in Figure 14 The average susceptibility measurements compared to MEDI show that the difference in MEDI+DL is The relative susceptibility of DL is -13ppb (GP), -4ppb (PU), and -12ppb (CN). Compared with COSMOS, the relative susceptibility differences are: MEDI, -14% (GP), -30% (PU), and 3% (CN); DL, -23% (GP), -20% (CN); MEDI+DL, -17% (GP), -29% (PU), -7% (CN).

[0151] The average reconstruction time for MEDI is seconds, DL is 2 seconds, MEDI+DL is 30 seconds (excluding network computing time).

[0152] The average fidelity term (normalized) is and

[0153] MS patients. A QSM of an MS patient is shown in Figure 13 (Bottom row) Compared with results using either MEDI or DL, the venous structures close to the corpus callosum were better delineated using MEDI+DL.

[0154] SGM ROI measurement can be done at Figure 14 b. It shows the mean magnetic susceptibility measurements of 8 patients, indicating that the difference between DL and MEDI is 3ppb (PU), -2ppb (CN), the difference between MEDI+DL and MEDI is 23ppb (GP), -7ppb (PU),

[0155] Linear regression and Bland-Altman plots for MS lesion measurements are shown in Figure The measurement between MEDI and MEDI+DL is shown in 2 =0.81) showed a better correlation than that between MEDI and DL (R 2 =0.62). In addition, Bland-Altman analysis showed that the agreement interval between MEDI and MEDI+DL The ppb range is narrower, compared to [-27, 13] ppb between MEDI and DL.

[0156] The average reconstruction time is 112 seconds for MEDI, 3 seconds for DL, and 41 seconds for MEDI+DL (excluding network inference time).

[0157] The fidelity cost of MEDI (normalized) is DL is 1.84%, MEDI+DL is

[0158] in conclusion

[0159] Our data demonstrate that QSM reconstruction, which combines the output of a deep learning (DL) neural network with Bayesian optimization reconstruction, provides superior contrast and maintains fidelity compared to both DL and traditional Bayesian optimization alone. DL alone suffers from significant fidelity deviations, which we mitigate by enhancing data fidelity through Bayesian optimization.

[0160] Conventional Bayesian QSM reconstruction consists of a fidelity term that incorporates the measured data and one or more regularization terms that penalize any deviation of the magnetic susceptibility distribution from prior knowledge. Typically, regularization is expressed in terms of L1 and L2 norms, and these terms have known limitations. Applying L2-regularization to the gradient promotes a substantially smooth solution, while L1-regularization of the gradient promotes a piecewise smooth solution. Therefore, the total variation term in MEDI (Equation 14) suppresses the fine structures close to the corpus callosum (small veins, Figure 13 ) tissue structure. These limitations of traditional regularized Bayesian MRI reconstruction can be addressed by regularization based on deep convolutional neural networks, as deep learning outperforms traditional methods in defining subtle image features. We introduce deep learning into the Bayesian optimization framework to solve the inverse QSM problem, utilizing the approximate magnetic susceptibility distribution generated by the neural network and constructing a penalty term relative to its difference to guide the optimization process. Compared with MEDI, our proposed MEDI+DL can better display small veins. At the same time, MEDI+DL only requires about 30% of the time cost of MEDI, because the calculation of the total variation in Equation 14 is replaced by Equation L2 regularization calculation in .

[0161] DL neural networks are fast for image processing and can handle complex nonlinear mappings, which has been widely demonstrated in denoising, super-resolution, segmentation, and reconstruction problems. In this paper, a feed-forward 3D convolutional neural network is used to perform a direct transformation from tissue fields to magnetic susceptibility maps. However, previous work has not considered fidelity to the measured data. In this work, we evaluate the necessity of fidelity by pairing neural networks with a fidelity cost in a numerical optimization framework. As shown here, DL has a much higher fidelity cost than MEDI, while MEDI+DL achieves a lower fidelity cost than MEDI.

[0162] Because COSMOS does not require regularization, it serves as a reasonable reference standard for this study. DL generally underestimates susceptibility compared to COSMOS. MEDI+DL improves this underestimation, as demonstrated in measurements of deep gray matter ROIs. MEDI+DL and MEDI provide comparable susceptibility measurements in deep gray matter and multiple sclerosis lesions. However, some questions remain to be answered: while MEDI+DL underestimates Putamen better than MEDI, it underestimates the globus pallidus and caudate nucleus more than MEDI.

[0163] We use a pre-trained network to provide χ φ Instead of training a cascade network to simulate a dedicated iterative projection process, we use an iterative quasi-Newton algorithm to solve the target problem (Eq. 1). The uncertainty of the Gaussian assumption (Eq. 12 and ) can be replaced by statistical probability distribution estimates. However, fundamentally, this framework still involves an explicit formulation of regularization, which undermines the advantages brought by deep learning. Future work should include exploring more effective ways to combine deep learning neural networks with traditional Bayesian optimization.

[0164] In summary, we proposed a Bayesian approach to the inverse QSM problem using neural networks for regularization, which showed consistency with methods using regularization while providing superior anatomical detail at less than half the computational cost.

[0165] 3. Fidelity-Preserving Imposed Network Editing (FINE) for Pathological Image Reconstruction

[0166] In this section, we describe Fidelity-Preserving Imposed Network Editing (FINE), where networks are artificial neural networks in deep learning, including convolutional neural networks. The pre-trained weights of a prior network in FINE can be modified based on a physical model of a test case. We experimentally demonstrate that FINE achieves superior performance on two important inverse problems in neuroimaging: quantitative susceptibility mapping (QSM) and undersampled multi-contrast reconstruction in MRI.

[0167] 3.1 Introduction

[0168] Image reconstruction from noisy and / or incomplete data is often addressed through various forms of regularization, which can be expressed as maximum a posteriori (MAP) estimates in Bayesian inference. Traditionally, these regularizations promote the desired properties of explicitly extracted image features, such as image gradients or wavelet coefficients. Compared to explicit feature extraction, deep learning (DL) using multi-layer convolutional neural networks has demonstrated a superior ability to capture all desired image features and has achieved great success in a wide range of computer vision applications. Recently, DL has been applied to image reconstruction.

[0169] The most popular approach formulates the problem as supervised learning, where a DL model is first trained on input data and target image pairs, and then reconstructs images directly from the test data. However, the performance of this supervised DL approach depends heavily on the similarity between the test and training data. Any even slight deviation can lead to significant errors in the reconstruction. For example, if the test case exhibits a certain pathology that is not present in the training data, the DL model may fail to capture the pathology. This is partly because DL-based techniques are agnostic to the underlying physics and do not incorporate a known physical model of the imaging system.

[0170] A common approach to combining DL with physical models of imaging systems is to use the DL model to define explicit regularization in the classic Bayesian MAP framework, typically via an L1 or L2 penalty. However, traditional explicit regularization terms in Bayesian reconstruction may provide imperfect feature descriptions and limit image quality.

[0171] The advantage of DL over explicit feature extraction may come from the fact that the explicit feature representations used during training are buried deep in many convolutional layers via backpropagation. Therefore, we propose to incorporate into the DL layers a physical model of the test data, or a data fidelity term, defined by the difference between the forward model of the measured data and the target image. One way to achieve this is to edit the DL network weights via backpropagation according to the data fidelity of the given test data, and we call this method Fidelity-Imposed Network Editing (FINE). We report preliminary FINE results on two neuroimaging problems, quantitative susceptibility mapping (QSM) and MRI reconstruction from undersampled k-space data.

[0172] 3.2 Theory

[0173] Consider the linear forward problem:

[0174] y=Ax+n,

[16]

[0175] Where x is the desired image, y is the measured data, A is the imaging system matrix defined by the known physical model, and n is the noise in the measured data. Bayesian reconstruction under Gaussian noise is:

[0176]

[0177] where W is the square root of the inverse of the noise covariance matrix, and R(x) is a regularization term that represents prior knowledge. The L2 term in Equation 17 is referred to as the data fidelity. The method of Equation 17 can be solved using numerical optimization methods such as the quasi-Newton method or the alternating direction iterative multiplier method. Common choices for R(x) include the total variation (TV) or the sparsity of the L1 norm of the wavelet coefficients in an appropriate wavelet domain. These types of priors are crucial for solving ill-posed inverse problems but can also limit the quality of the reconstructed image, for example by exhibiting blocking constants.

[0178] Fundamentally, regularization promotes desirable image features that can be advantageously implemented using deep learning (DL) rather than traditional explicit feature extraction. The convolutional neural network φ(; Θ0) with convolutional weights Θ0 from data to image can be obtained by {α i , β i} data pairs for supervised training, where α i is the real image, β i is the input data. The weights of each convolutional layer and the nonlinear activation function can be viewed as a set of feature extractors for the desired image reconstruction. The large number of weights in DL can explain the advantages of DL over explicit feature extraction using a single or a few weights. Although the training data β may usually be of a different type (size, contrast, etc.) than the test data y, the test data y can be treated as the same type as the training data β to generate DL reconstructions:

[0179]

[0180] An example of such an approach for solving ill-posed inverse problems is QSMnet, which aims to solve the field-to-susceptibility dipole inversion in QSM.

[0181] If there is a mismatch between the test and training data, this supervised DL strategy may perform poorly because it is agnostic to the forward physical model generated by the data defined for the imaging system (Eq. ). In particular, when the measured data differs significantly from the training data, this lack of data fidelity can lead to substantial errors in the network output. To address this lack of fidelity, it has been proposed to treat the network output in Equation 18 as a regularization of Equation 17, using an L2 form to penalize the difference between the network output and the final optimized solution:

[0182]

[0183] We refer to this reconstruction as DL with L2 regularization (DLL2). The main drawback of this DLL2 approach is the use of an explicit L2 norm, which is known to be flawed and may not be effective in reducing the deviation of the final solution from the fidelity error of the fixed network result φ(y; Θ0).

[0184] To leverage DL to achieve performance superior to explicit feature extraction, we propose embedding data fidelity terms deeply into all layers via backpropagation in the DL network. One approach to implementing this approach to reconstruct the desired image x is to edit the weights in a pre-trained DL network guided by the data fidelity of the given test data y. The weights of the network φ(.; Θ) are initialized with Θ0 and edited based on the physical fidelity of the imaging system for the test data:

[0185]

[0186] The output of the updated network is the reconstructed image x with data fidelity and deep learning regularization:

[0187]

[0188] We call this approach "Fidelity-Preserving Imposed Network Editing (FINE)" for solving ill-posed inverse problems using deep learning and imaging physics.

[0189] 3.3 Materials and methods

[0190] We applied the proposed FINE to two inverse problems in MRI: QSM and multi-contrast reconstruction. The human subjects studies followed IRB-approved protocols.

[0191] QSM

[0192] First, we apply FINE to QSM, which is ill-posed due to the zeros at the magic angle of the forward dipole kernel. Therefore, stripe artifacts appear in the image domain after the non-regular dipole inversion. Bayesian methods have been widely used to solve this problem. An example is the Morphologically Enabled Dipole Inversion (MEDI) method, which adopts the following objective function:

[0193]

[0194] where χ is the susceptibility distribution to be solved, f is the measured field, and d is the dipole kernel. Regularization is a weighted total variation where is the gradient operator, M G is a binary edge mask determined from the magnitude image that enforces morphological consistency between magnitude and susceptibility.

[0195] Data acquisition and preprocessing. A 3T system (GE, Waukesha, WI) with a multi-echo 3D gradient echo (GRE) sequence was used to acquire the data. MRI was performed on healthy subjects. Detailed imaging parameters included TR=39ms, Voxel size = 1 × 1 × 3 mm3, The local tissue field is estimated by using nonlinear fitting across the multi-echo phase data, followed by phase unwrapping and background field removal based on a graph cut algorithm. COSMOS reconstruction is based on the GRE imaging was repeated in 10 different orientations, and the reconstructions were used as the gold standard for brain QSM. In addition, GRE MRI was performed on 8 patients with multiple sclerosis (MS) and 8 patients with intracerebral hemorrhage (ICH), using the same 3T system and the same sequence only in the standard supine position.

[0196] Dipole inversion network. We implemented 3D U-Net, a fully convolutional neural network architecture, for mapping from the local tissue field f to the COSMOS QSM. The convolution kernel size is 3×3×3. The detailed network structure is shown in Figure 11 shown. Among healthy subjects Name used for training, in-plane rotation For data enhancement. Each 3D volume data is divided into The total number of blocks is We randomly select 20% of these blocks as validation sets. We use the same loss function combination as used in training the QSMnet network with the Adam optimizer (learning rate 10 -3 , epoch 40), and obtain 3D U-Netφ(;Θ0).

[0197] Fidelity Imposed Network Edit (FINE). For given test data, the network weights Θ0 from training are used to initialize weights Θ in the following minimization:

[0198]

[0199] This minimization fine-tunes the pre-trained dipole inversion network φ(f;Θ) to produce outputs that are consistent with the forward dipole model for the given test field data f. Equation 23 is minimized using Adam (learning rate 10 -3 , the iterations stop when the relative reduction of the loss function between two consecutive iterations reaches 0.01. The final reconstruction of the fine-tuned network is

[0200] FINE was applied to one healthy subject (excluded from the training set), eight MS patients, and eight ICH patients. The MEDI reconstruction algorithm (λ = 0.001) was used as one of the baseline algorithms. As another baseline, we also implemented the following equation based on Equation 19:

[0201]

[0202] where λ2 is chosen so that the fidelity term ||W(fd*χ)||2 is similar to the fidelity cost in MEDI (Eq. 22).

[0203] Quantitative Analysis. The fidelity term ||W(fd*χ)||2 was calculated for DL, MEDI, DLL2, and FINE. For healthy subjects, the reconstructed QSM was compared with COSMOS (a widely used (but expensive) QSM reference standard) in terms of fidelity and the structural similarity index (SSIM), a metric that quantifies image intensity similarity, structural similarity, and contrast similarity between pairs of image patches.

[0204] Multi-contrast MRI reconstruction

[0205] Secondly, we apply FINE to multi-contrast MRI reconstruction with undersampled data. To speed up the time-consuming acquisition of certain contrasts, such as T2-weighted (T2w) or T2 fluid-attenuated inversion recovery (T2FLAIR) images, k-space is undersampled, so a regularization algorithm is needed to recover images with minimal artifacts. To help solve this ill-posed problem, fully sampled images of another contrast, such as T1-weighted (T1W) images, are combined to exploit shared structural information in the contrasts. Bayesian inference with Directed Total Variation (DTV) regularization can be used for image reconstruction:

[0206]

[0207] where b is the measured undersampled k-space data, U is the binary k-space undersampled mask, is the Fourier transform operator, u is the image to be solved (T2w / T2FLAIR), and γ is the measured reference image (Tlw). b and γ form the test data. Anisotropic regularization is a penalty to encourage parallel image gradients in the T2w / T2FLAIR and Tlw images at each voxel.

[0208] Data acquisition and preprocessing. We acquired and registered T1w, T2w, and T2FLAIR axial images of 237 patients with multiple sclerosis. Matrix size and 1mm 3 Isotropic resolution. For each contrast, we extract from each volume axial 2D image slices, resulting in a total of The grayscale value of each slice is normalized to the range [0, 1].

[0209] Contrast Conversion Network. We adopt 2D U-Net as the network architecture for converting fully sampled T1w images to fully sampled T2w images, assuming that T1w and T2w images share common tissue structures but have different contrasts. The network is trained using a 3×3 convolution kernel ( Figure 11 We use the L1 difference between the network output and the target image as the loss function for training the network with the Adam optimizer (the learning rate is initialized to 10 -3 , epoch 40). There are samples. This training builds a U-net φ(; Θ0). Similarly, we build another 2D U-Net for mapping from fully sampled T1w images to fully sampled T2FLAIR images.

[0210] Fidelity-imposed network editing. The test data b of the test object are obtained by undersampling the subject's axial T2w image by a factor of 4.74 using a sampling pattern of a Poisson disk in k-space. Similar to Equation 23, we initialize the network weights Θ with Θ0 and update them using the following minimization:

[0211]

[0212] Adam optimizer (learning rate 10 -3 , the iteration stops when the relative reduction of the loss function between two consecutive epochs reaches 0.01). The reconstructed 2w image is the final result of the editing network:

[0213]

[0214] Similarly, T2 FLAIR images were reconstructed.

[0215] The FINE reconstruction is compared with DTV (λ = 0.001) and DLL2, where Equation 19 for DLL2 has the following form:

[0216]

[0217] Here, λ2 is selected so that the value of the fidelity term is similar to that of DTV.

[0218] Quantitative analysis. By calculating the fidelity terms of DTV, DL, and DLL2

[0219]

[0220] And SSIM is used to quantify the similarity between the reconstructed image and the real image.

[0221] 3.4 Results

[0222] Our experimental results include network weight editing, QSM, and multi-contrast reconstruction.

[0223] Editing network weights

[0224] For the example case of applying FINE in reconstructing the QSM of an MS patient, the difference between Θ0 and Θ is shown in Fig. During pre-training, the weights Θ0 of 3D U-Net are estimated almost uniformly along all layers (Fig. ). FINE significantly changes the weights in the encoder and decoder parts of the network, but changes the weights in the intermediate encoding vector layer. (picture and c). Compared to FINE, the random Θ in Equation 20 is initialized using a truncated normal distribution centered at 0 (Figure ),

[0225] n is the number of input units in the weight tensor (depth image prior), which leads to a significant change in weights in all layers (Figure yang) and resulted in significantly poorer QSM (Fig. and h).

[0226] QSM

[0227] Healthy subjects. QSM reconstructed by COSMOS, MEDI, DL, DLL2 and FINE are shown in Figure 17 and 18 Structures in the occipital lobe ( Figure 18 The enlarged area in a) is clearly delineated in the FINE and DLL2 reconstructions but is more blurred in the MEDI and DL. In this case, the fidelity term and SSIM are as follows Figure 17 As shown in a, FINE shows the best performance (minimum fidelity value and maximum SSIM).

[0228] MS patients. QSMs of two representative MS patients reconstructed by MEDI, DL, DLL2, and FINE are shown in Figure 18 b. Compared with MEDI or DL, FINE and DLL2 showed the fine structure of the central vein near the ventricle ( Figure 18 Magnified area in b).

[0229] ICH patients. QSM reconstructed by MEDI, DL, DLL2 and FINE is shown in Figure 18In c, severe shadowing artifacts were observed in DL and DLL2. These artifacts were significantly suppressed in MEDI and FINE.

[0230] Multi-contrast MRI reconstruction

[0231] T2w and T2 FLAIR images reconstructed by DTV, DL, DLL2 and FINE are shown in Figure 19 and Figure 20 In DTV and DL, structural details such as the white matter / gray matter boundary are blurred. They are clearly delineated in FINE and DLL2, with DLL2 being visually noisier.

[0232] The fidelity value and SSIM measurement results are as follows Figure 17 As shown in b, FINE shows the best performance.

[0233] in conclusion

[0234] Our results demonstrate that the proposed Fidelity-Preserving Imperative Network Editing (FINE) approach can be highly effective in solving ill-posed inverse problems in medical image reconstruction. FINE embeds the desired physical model of the test data into a multi-layer deep learning (DL) network via backpropagation. Consequently, FINE enables the robust use of DL as an implicit regularizer when constraining the output manifold of the DL network and proposes an approach that outperforms traditional explicit regularizers for solving ill-posed inverse problems in Bayesian reconstruction. Compared to traditional total variation (TV) regularization, DL, and DL-based L2 regularization, FINE is superior in recovering subtle anatomical structures missed in TV and resolving pathologies not encountered in the DL training data.

[0235] DL has recently been used to solve inverse problems in medical image reconstruction. Deep neural networks (DNNs) are often trained to map directly from the data domain to the image domain. For example, DL can be used to map tissue fields into QSMs. This approach bypasses the time-consuming iterative reconstruction typical of traditional numerical optimization and significantly reduces the reconstruction time (forward pass through the network). However, this approach does not take into account the deviation in fidelity between the reconstructed image and the actual measured data. This problem of lack of data fidelity has recently been recognized in DNN image reconstruction. Data fidelity can be approximately enforced by using iterative projections of many convolutional networks. A more precise way to achieve data fidelity is to use a Bayesian framework with explicit regularization, typically the L2 norm of the difference between the desired image and the network output (DLL2). However, the use of the L2 norm or other forms of explicit regularization may introduce artifacts in the final reconstructed image and is considered inferior to DL for image feature characterization. In the case of measured data containing pathological features that deviate significantly from the training data, the deviation in the network output may not be effectively compensated. This is in Figure 18 As an example, Figure c illustrates that DL and DLL2 fail to correctly capture hemorrhage features not encountered in the training dataset of healthy subjects. This issue is addressed in the proposed FINE method, where the bias of the pre-trained network is effectively reduced by updating the network weights guided by the measured data. Compared to DLL2, which only optimizes the image itself, FINE updates more network weights, which can make the proposed method more flexible and effective than DLL2.

[0236] Because QSM requires prior information to perform ill-posed dipole inversion, the search for better image features to regularize the reconstruction has been a major research effort in QSM development. Mathematically, regularization should suppress both streaking artifacts associated with granular noise and shading artifacts related to smoothing model errors. L1-type regularization has been effective in reducing streaking artifacts, but these techniques introduce stair-step artifacts. Shading artifacts have not been effectively suppressed. These challenges in QSM reconstruction can be more effectively addressed using complex image features. DL provides the required but difficult-to-describe complex image features. The FINE implementation of DL reported here realizes the potential of DL for QSM reconstruction.

[0237] Related previous work is deep image priors

[15] , which trains DL networks from scratch on a single dataset for the inverse problems of denoising, super-resolution, and inpainting. We The work in this paper shows that QSM is related to Figure 11 The network structure in

[15] was consistent, deep image priors failed to produce satisfactory results, and the use of pre-trained weights or FINE was necessary. In FINE, the network was initialized as a pre-trained network instead of training from scratch. Our empirical analysis showed that for the case of QSM reconstruction in MS patients, FINE mainly changed the weights of the encoder and decoder parts of U-Net (Figure 5). ), which reflects high spatial frequency content or patient-specific details. We expect FINE to be effective in improving other ill-posed inverse problems, such as in image reconstruction, using various DL networks, where noisy and / or incomplete data are available and the physical model of the data generation is known.

[0238] FINE generates the desired image using the network weights trained on the corresponding fully sampled T1-weighted data, and the sampled T2-weighted and T2 FLAIR data, along with structural priors. Similarly, FINE can be extended to generate images from a large amount of undersampled data of other tissue contrasts, including perfusion / diffusion / susceptibility weighting.

[0239] Future work may involve evaluating FINE in a wide range of applications including super-resolution and denoising. Due to the additional network updates based on iterative optimization, the computational cost of FINE is much higher than a single pass through the DL network. The computational cost can be reduced by updating a subset of layers instead of the entire network.

[0240] In summary, a data fidelity term can be used to update a deep learning network on a single test data set to produce high-quality image reconstructions. This fidelity-imposed network editing (FINE) strategy is promising for solving ill-posed inverse problems in medical imaging.

[0241] 4. Quantitative Transport Imaging

[0242] Contrast-enhanced MRI is a routine sequence in clinical practice. The contrast agents used in MRI are highly paramagnetic, allowing for quantitative measurement using QSM. Indeed, QSM originated from the need to address the issue of contrast quantification in time-series contrast-enhanced MR angiography. Images of contrast agent concentration versus time, quantified by QSM, can be used to study contrast agent or tracer transport in tissues.

[0243] In this section, we describe the inversion of transport field quantities from time-lapse imaging of tracers, which is formulated as an optimization problem and named quantitative transport mapping (QTM). QTM is clinically feasible in the context of porous vascular structures, and its velocity maps computed from dynamic contrast-enhanced (DCE) MRI data can be used to characterize transport processes in tumors.

[0244] 4.1. Introduction

[0245] Currently, the quantitative interpretation of time-lapse imaging data containing tracers is based on the Kety equation, which requires a global arterial input equation (AIF) for all local tissues. The Kety equation was originally based on the conservation of tracer mass in an organ, with the AIF defined at the aorta supplying the organ. In the first attempts to apply the Kety equation to time-lapse tomography, the difficulty of measuring the AIF at the voxel level was recognized. The arterial supply of tissues in a voxel is multiple, and voxel tissue deviates from the global AIF due to delays and dispersion in the transport process. Therefore, flow quantification based on the Kety equation is subject to error. In addition, identifying the AIF can be difficult and is related to the operator's experience and the patient's condition.

[0246] When considering voxel-level mass conservation in tomographic imaging, including computed tomography (CT), magnetic resonance imaging (MRI), and positron emission tomography (PET), the 3D distribution of the contrast agent during arterial transport to the voxel of interest should be considered. Naturally, the transport of the tracer in tissue should be described by transport equations derived from first principles of physics. These transport equations reflect local mass conservation using scalar, vector, and tensor field quantities, eliminating the unphysical assumption of global arterial input used in Kety's equations. Recently, this transport field equation has been employed to describe the mass conservation observed in tomographic imaging.

[0247] We are developing a quantitative interpretation of tracer time-series imaging based on transport equations. As described below, we report the calculation of drift velocities using a weak microstructure approximation. In this approach, the complex structure affecting the transport process is approximated as a homogeneous porous medium, allowing for a simple forward problem based on the continuity equation. The inverse problem of estimating flow from imaging data can then be solved using the alternating direction multiplier method. We report MRI-based results, where tracer particle concentrations are estimated directly from time-resolved contrast-enhanced (DCE) MRI.

[0248] 4.2. Theory

[0249] The transport of molecules is described by the Fokker-Planck transport equation, which describes the time evolution of the probability density function of a particle under the action of water pressure, drag, and random forces. In continuous space with infinite resolution, the Fokker-Planck transport equation can be expressed as:

[0250]

[0251] where c(r, t) is a probability density function that can be interpreted as the particle concentration, is the derivative with respect to time, The derivative of three-dimensional space, u(r)=[u x ,u y ,u z ] is the velocity vector of the particle flux, which is related to the particle concentration, and D(r) is the diffusion coefficient matrix, which reflects the particle flux proportional to the concentration gradient. u(r) and D(r) are assumed to be constant over the experimental observation time period. The partial differential equation (PDF) in Eq. 31 reflects local mass transport and is therefore called the continuity equation. For tomographic imaging experiments, Eq. 31 can be modified to account for tracer decay; for example, arterial spin labels in MRI decay at a rate of 1 / T1, and radiolabels in PET decay at a rate of ln2 / half-life.

[0252] Time series tomographic data of tracer concentrations can be fitted to Eq. 31 to estimate u(r) and D(r) for each voxel, but their calculation and interpretation may require microstructural information. The tissue in a voxel can be divided into vascular compartments defined by vascular structures (v, microvascular structures Ω v , blood volume v v and the contrast agent concentration c v ) and the extravascular compartment (e, extravascular space Ω e , volume v e and the contrast agent concentration c e ), the average concentration in the voxel is expressed as:

[0253] <c>=v v c v +v e c e

[32]

[0254] The motion of the tracer is assumed to be the same as that of the carrier: water. In the v-chamber (water density is ρ0), blood flow is mainly generated by the drift velocity, which is determined by external forces, including the viscosity μ = vρ0 and the water pressure p. According to the Navier-Stokes equation

[0255]

[0256] Assume blood is incompressible For healthy tissue, contrast agents stay in the v-chamber, and for diseased tissue, they can pass through the vessel wall into the e-chamber. The movement of contrast agents in the e-chamber can be approximated as follows:

[0257]

[0258] The contrast agent concentration in the two chambers v and e can be determined by the difference between the two chambers at the boundary ( Normal vector ) is related to the permeability (ψ) of the vascular wall, the flow rate through the vascular wall is:

[0259]

[0260] For known vascular structures Ω v Image space, Eqs.32, 34 and can be solved to obtain the voxel average <c>The contrast agent concentration c v , c e , and Eq.33 can be solved to obtain These solutions can be substituted into Eq. 31 and integrated over the voxels to obtain a discretized equation. Subsequently, the tissue transport parameters p, D, and ψ can be solved by fitting the difference equation to the time series images. We call this process voxelization of the transport equation.

[0261] Typically, transport voxelization is very complex and depends on detailed prior information about tissue composition and vascular structure within the voxel. It can be used to generate time series images corresponding to the parameters p, D, ψ. These corresponding data can be used to train an artificial neural network, which can then be used to estimate the transport parameters p, D, ψ from the acquired time-resolved images.

[0262] In our first attempt to voxelize transport processes, hoping to obtain solvable and reasonable equations, we assumed that the tissue is either uniformly porous or has a weak microstructure with little compartmentalization. As a zero-order approximation, this porous tissue model roughly describes the tumor, where the vessel walls are porous, the blood flow is high, and the tracer can diffuse freely. We can approximately describe the dynamics of transport in porous media by the average drift velocity over the voxels, without any detailed microstructural information. The estimation of the kinetic parameters of the average drift velocity only requires fitting the data to the continuity equation. This greatly simplifies the estimation of the tissue parameters p, D, ψ, which requires detailed information on the vascular structure. We define the following discrete spatial derivative, ξfor c i (ξ)=c(ξ,t i ): where Δx, Δy, and Δz are the voxel sizes. Subsequently, the calculation of the average drift velocity can be interpreted as fitting the simplified Eq. 32:

[0263]

[0264] here, It is c i The discrete spatial gradient of (x, y, z), u = [u x ,u y ,u z For L1 regularization It is a LASSO (least absolute convergence and selection operator) regression problem that can be solved efficiently using the alternating direction method of multipliers (ADMM). Afterwards The augmented Lagrange multiplier can be written as

[0265]

[0266] here and The ADMM algorithm is updated according to the following process

[0267]

[0268]

[0269] in Defined as Where p, q = 0, x, y, z, t, and In Eq. 38, all multiplications are performed point by point, Δ is the Laplace operator, and · is the divergence operator.

[0270]

[0271]

[0272] in Otherwise 0 represents the soft threshold. The update of the dual variable is as follows:

[0273]

[0274] 4.3. Methods and Experiments

[0275] Our proposed QTM method was applied to 4D (3D time series) dynamic contrast-enhanced (DCE) data using gadolinium to enhance the MRI signal of tissue. This study was conducted under a HIPAA-compliant IRB-approved protocol for retrospective analysis of MRI data. DCE MRI data from patients with tumors in various organs were individually identified and processed.

[0276] Data collection

[0277] DCE of brain tumor patients (n = 10; 2 meningiomas, 7 high-grade and 1 low-grade glioma): T1-weighted 3D DCE data before and after contrast injection (0.1 kg gadobutrol / kg body weight) ( Patients) and a 3T MRI scanner (GE Healthcare, Milwaukee, USA). Acquisition parameters included: voxel size, matrix, Flip angle, Temporal resolution and 30-40 time points.

[0278] DCE of patients with liver tumors (n=7, 1 metastasis, 1 hemangioma, 2 focal nodular hyperplasias, and 3 hepatocellular carcinomas). T1-weighted images were acquired using a 3D stacked spiral sequence from a Tesla (GE Healthcare, Milwaukee, USA). Whole-liver 4D images were reconstructed using a sliding window. Imaging parameters included: scan time Frame rate Flip angle, Voxel size, TR / TE The slice with the most prominent tumor was selected as the ROI, and the aorta and portal vein were selected as the AIF.

[0279] DCE in breast cancer patients Benign tumors and All patients underwent MRI examinations on a 3T GEMRI system. Ultrafast DCE-MRI using differential sampling of Cartesian coordinates (DISCO) was previously used. Continuous acquisition within seconds Stages ( The contrast agent injection (0.1 mmol gadobutrol per kg body weight) was started simultaneously with 10 phases (8-channel breast coil) or 10 phases (8-channel breast coil). Additional acquisition parameters included TR / TE = 3.8 / 1.7 ms, flip angle = 10°, Seconds, axial.

[0280] Image analysis

[0281] QTM processing is based on Eq.38 and is fully automatic. We process all data in the same way. Patient 4D image data c i (x, y, z), where x, y, z are coordinate names, i = 1, ..., N t is the number of frames, N t is the total number of time frames. We use time series DCE images to generate contrast agent concentration data, assuming that the contrast agent concentration change and image signal intensity are linearly related. We use the gradient descent method to solve the system of Eq. 38, and the relative error is 10 -3 Maximum number of cycles The regularization coefficient λ is determined based on L-curve analysis.

[0282] The voxel drift velocity value is calculated according to the following equation:

[0283]

[0284] Regional blood flow map of traditional Kety method (f Kety )(ml / 100g / min) is also calculated to provide a comparison with the QTM drift velocity graph |u|:

[0285]

[0286] where c4(t) is the global arterial input equation (AIF) and v(r) is the regional blood volume. For the DCE brain tumor data, the AIF was selected from enhancing arteries close to the tumor volume. For the DCE liver data, the AIF was selected from the combination of the aorta and portal vein in the same axial slice at the tumor center. For the DCE breast data, the AIF was selected from the ipsilateral internal mammary artery. We evaluated the dependence of Kety's method on AIF selection in the tumor cases by comparing blood flow calculated with another AIF selected from the internal carotid artery, as well as by calculating blood flow from another AIF selected from an axial slice away from the tumor center.

[0287] We drew regions of interest (ROIs) within the tumor bulk and other tissues of interest for linear regression analysis between QTM and Kety's method in all subjects.

[0288] Results

[0289] brain tumors

[0290] Figure 21 Gd-enhanced T1-weighted axial slice ( Figure 21 a) shows the anatomical information. QTM velocity component map ( Figure 21 c) and the velocity vector u image ( Figure 22 ) shows blood flow in and within the vessels supporting the tumor

[0291] Figure 23 Glioblastoma (WHO grade IV) is shown. Gd-enhanced T1-weighted axial slice ( Figure 23 a) Shows anatomical tumor information. Kety flow f Kety ( Figure 23 c) Calculated from the AIF close to the tumor, the QTM velocity |u|( Figure 23 d) Depicts similar blood flow in tumors, blood vessels, and other tissues. Figure 23 b shows the average f in the tumor ROI for all 10 brain tumor patients. Kety The linear regression of the correlation between |u| and |u| shows the correlation coefficient Good consistency.

[0292] When AIF is closer to the tumor, the f in the tumor ROI is larger than that in the internal carotid artery when AIF is taken. Kety Higher: Averaged across all patients, the relative difference within the tumor ROI is

[0293] Liver tumors

[0294] Figure 24 Showing metastatic lung cancer. Kety flow diagram f Kety ( Figure 24 c) and QTM speed |u|( Figure 24 d) shows similar contrast. Kety flow map f Kety The velocity map |u| is higher (~4 times) in the tumor area (white arrow) compared with normal liver tissue, which corresponds to the tumor enhancement on post-Gad T1-weighted images ( Figure 24 The linear regression showed that the average velocity of QTM and the flow rate of Kety's method were well correlated (R 2 =0.68).

[0295] Kety's method was sensitive to the selection of AIFs: in all cases, tumor f Kety The average relative difference of normal liver tissue was 11.9%.

[0296] Breast tumors

[0297] DISCO images, QTM velocity maps, and f of a malignant tumor Kety The picture is in the picture. The diseased area is shown in the QTM velocity map and f Kety Linear regression showed a good correlation between the QTM average velocity and the Kety's method flow rate (R2 = 0.62).

[0298] In the figure. In Figure 2, box plots show the mean QTM velocity of all diseased areas and the flow rate comparison of Kety method. The t-test showed that the QTM velocity was significantly different between benign and malignant tumors (p = 0.04), but f Kety (p = 0.12) were not significantly different (Fig. ).

[0299] in conclusion

[0300] Our preliminary results demonstrate the feasibility of calculating drift velocity maps from time-series image data of tracer concentrations based on a quantitative description of tissue transport processes. This quantitative transport imaging (QTM) approach eliminates a major obstacle in the conventional Kety method, namely the need to use an arterial input function (AIF) to quantitatively interpret tracer transport in time-series imaging. The elimination of the AIF makes QTM fully automated and user-friendly in routine clinical practice, and more stable than Kety's method. Preliminary results suggest that drift velocity maps from QTM have the potential to be used to characterize tumor blood flow. The transport partial differential equation (Eqs. 31- ) reflects the physical laws of mass conservation and dynamics: the temporal variation of contrast agent concentration in a voxel is related only to the contrast agent concentrations in nearby voxels. Voxelization, performed by integrating the solutions of transport equations over the in vivo tissue microstructure, provides a systematic way to construct tissue compartments based on fundamental tissue properties, which may be an ideal alternative to traditional phenomenological compartments. Fitting time-series tracing image data to mass conservation equations has been used to quantify arterial blood flow in high-resolution X-ray projection angiography, but without voxelization. Voxelizing transport equations using microstructure allows understanding of tissue transport, including regional blood flow, blood pressure, vessel wall permeability, and diffusion in tissues. In the first report of QTM, an approximation for large-flow porous media based on the physical law of mass conservation is proposed, which eliminates the need for the AIF used in the Kety method, used in all current perfusion quantification techniques in medical imaging.

[0301] Historically, the Kety equation describes the exchange of diffusible tracers in tissues by accounting for the blood partition coefficient or regional blood volume, an extension of the principle of mass conservation for blood input and output to organs. A fundamental flaw in current perfusion quantification methods is the global extension of the Kety model to characterize the entire organ as a voxel-based perfusion, which must be based on local quantities or fields. The Kety equation treats the entire organ as a lumped element, in a manner similar to how Ohm's law treats a conductor as a lumped element without considering the field distribution within the conductor. As in electrodynamics, local tissue perfusion requires consideration of distributed quantities characterized by vector and tensor fields and their corresponding partial differential equations (PDEs) that govern their local variations.

[0302] Results from preliminary QTM reported drift velocities that characterize blood flow in a variety of tumors with high flow and porous vessel walls. Although a rough correlation between Kety flow and QSM velocities was observed, QTM differs significantly from Kety’s method. The first difference is that manual ROI definition of the AIF is used in Kety’s method, whereas QTM is fully automated and does not use an AIF. It is well known that direct or indirect estimation of the AIF leads to poor reproducibility in perfusion imaging, as observed here in tumor data from the brain and liver. Attempts to estimate the local AIF or set the delay as an unknown variable cannot determine the important local spread of the AIF. Perhaps more importantly, it is often difficult to manually draw an arterial ROI, making Kety’s method difficult to use and difficult to reproduce in practice.

[0303] A second difference is that the quantitative correlation between QTM velocity and Kety flow is variable from tumor to tumor. In particular, QTM velocity showed statistically significant differences between malignant and benign breast lesions, whereas Kety flow did not. The superior sensitivity of QTM over the Kety method is likely due to 1) appropriate biophysical modeling of contrast agent transport in the tumor and 2) elimination of variance typically associated with manually entered AIFs in data processing. This very encouraging preliminary result warrants further investigation in a larger number of cases to determine the relationship between flow velocity and tumor malignancy and benignity.

[0304] A third difference is that the drift velocity maps from QTM are directional, whereas the flow maps from Kety's method are scalar. Tissue perfusion, as a transport process, is directional and should be considered directional when interpreting the data. This QTM work brings the physics of transport equations to bear on previous attempts to introduce directionality into perfusion imaging. The directional information from QTM could be very useful in guiding lesion interventions. For example, when delivering therapeutics via intravascular catheters or intratissue injections, it is important to know the direction of blood flow.

[0305] The mean velocity magnitudes observed here within the tumor ROI ranged from 0.1 mm / s to 1 mm / s. This velocity magnitude is the volumetric average of the tracer drift velocity within the voxel. Although tumor blood flow has been measured and discussed in numerous studies, tumor blood flow velocity has not been well quantitatively explored. Recent studies of tumor blood flow velocity measurements using flow MRI, capillaroscopy, and finite element simulations have reported similar velocity ranges within tumor tissue (from 0.1 mm / s to 0.8 mm / s in vessels ranging from 20 μm to 70 μm in diameter in rodent tumors and 0.2 mm / s in porous liver tumor tissue). This suggests that the QTM method provides reasonable estimates of tumor blood flow velocity. We have verified that QTM can accurately estimate flow velocities within straight vessels, consistent with X-ray angiography results. Our results require further experiments to verify the accuracy of QTM-derived velocities under various imaging conditions.

[0306] First attempts to solve the QTM were limited to modeling tissue as a homogeneous porous medium with weak microstructure, allowing for a simple solution of the forward problem. Important future work to further develop the QTM includes adding vascular microstructural information to the voxelization of the transport equations, which would allow voxel-specific averaging of drift velocity, tissue diffusion, and vascular permeability, as well as calculation of flux into and out of a voxel. Vascular wall permeability deserves special consideration because of its clinical importance in distinguishing pathological tissue from normal tissue. Vascular permeability is defined in the MRI literature as the flux through the vessel wall proportional to the contrast agent concentration gradient or, as in the Starling equation, as the flux through the vessel wall proportional to the pressure gradient; both definitions are equivalent because osmotic pressure can dominate the pressure gradient, and the relationship can be approximated by the van't Hoff equation. Radial velocity in the vessel can also be used to account for tracer leakage from the vascular compartment, but in many cases this radial velocity can be small.

[0307] Determine the transport parameters p, D, ψ and Deep learning methods can be used to solve the inverse problem of . Artificial neural networks (ANNs), such as the U-Net used above, can be trained by using simulated data covering all possible values ​​of p, D, and ψ. This simulation can be based on The ANN generates images based on tissue information including vasculature details and appropriate noise added to simulated image signals. This simulation is applicable to MRI, CT, PET, and other imaging modalities. The simulated images can be simulated from other patient data and used as input to the ANN, while the image labels p, D, and ψ are the output of the ANN. The trained ANN can then be applied to in vivo experimental data to estimate transport parameters.

[0308] Other important future work includes extending the physics-based QTM to perfusion quantification to other imaging modalities and validating the QTM against reference standards. The QTM applied here to DCE MRI data can be extended to other imaging modalities, including CT, PET, and single-photon emission computed tomography (SPECT), as well as arterial spin labeling data in MRI. A major problem in quantitative perfusion imaging is the lack of experimental validation with known true values, and the Kety method lacks experimental validation relative to true values ​​of flow calculated using methods independent of the Kety method. Critical validation of perfusion quantification by any method, including the Kety method and the QTM, must be performed in the future using experimental setups in which perfusion parameters are known.

[0309] In summary, quantitative transport imaging (QTM) overcomes the AIF problem in Kety's method by fitting local transport equations to time-lapse imaging of tracer transport in tissue (4D image data). QTM, under the weak microstructure approximation, can process time-lapse imaging data to automatically generate drift velocity maps that are promising for characterizing blood flow in brain, breast, and liver tumors.

[0310] Brain metabolic oxygen uptake fraction imaging (OEF)

[0311] Implementations of the above techniques can be used in a variety of applications. As examples, the following sections describe the application of one or more of these techniques in numerical simulations, phantom experiments, and human brain imaging. In particular, the temporally evolving cluster analysis (CAT) method is used to determine the iron concentration of deoxyhemoglobin and the oxygen extraction fraction of tissue, which significantly improves the signal-to-noise ratio of images and the accuracy of OEF measurements.

[0312] Introduction

[0313] Cerebral metabolic rate of oxygen (CMRO2) and oxygen extraction fraction (OEF) are important markers of brain tissue viability and function, such as in stroke, and their measurement using MRI has attracted great interest. Quantitative models have been proposed to investigate the effect of deoxyhemoglobin (a strongly paramagnetic substance) in the blood on the MRI signal: 1) amplitude signal modeling methods, such as Quantitative Imaging of Oxygen Uptake Fraction and Tissue Consumption (QUIXOTIC), calibrated fMRI, and quantitative BOLD (qBOLD), and 2) phase signal modeling methods for whole-brain CMRO2 imaging and voxel-based quantitative susceptibility imaging (QSM) for CMRO2 imaging.

[0314] Recently, a combined model based on QSM and qBOLD (QSM+qBOLD) was introduced to eliminate the assumptions in the individual models and simultaneously consider the effects of OEF on the amplitude and phase signals from the same gradient echo (GRE) data. Although it solves the problems of methods based on either QSM or qBOLD alone, it still has drawbacks. For example, the qBOLD model is nonlinear and has a strong coupling between venous oxygen saturation (Y) and venous blood volume (v), so its inversion is a difficult non-convex optimization and is highly sensitive to noise in data acquisition. This can lead to critical errors in parameter estimation for a specific signal-to-noise ratio (SNR) level. Since the QSM+qBOLD model partially relies on the qBOLD model, it is also sensitive to noise.

[0315] In this study, we introduced cluster analysis of temporal evolution (CAT) to overcome the sensitivity of QSM+qBOLD to noise and obtain more accurate parameters, such as OEF, by improving the effective SNR. Voxels with similar GRE signal evolution should have similar parameter values. Using machine learning methods, the whole brain can be divided into multiple clusters much smaller than the number of voxels based on the similarity in GRE signal evolution. For each cluster, a set of model parameters is assumed, so the signal-to-noise ratio is effectively increased by calculating the average value for each cluster. This improvement in SNR is expected to lead to more accurate OEF calculation. The results of the CAT method are then used as the initial assumptions for voxel-based optimization. The proposed CATQSM+qBOLD method is compared with the previous QSM+qBOLD method in healthy subjects and patients with ischemic stroke, using the OEF with a constant initial assumption for the whole brain.

[0316] theory

[0317] CMRO2 (μ mol / 100 g / min) and OEF (%) can be expressed as CMRO2 = CBF·OEF·[H] a , CBF is cerebral blood flow (ml / 100g / min), [H] a is the molar concentration of oxyhemoglobin in the artery (7.377 μmol / m1), based on is the hematocrit, Y is the venous oxygen content (unitless), and Y a is arterial oxygenation (set to 0.98). Assuming that deoxyheme iron is distributed in the draining veins and veins of a randomly oriented cylindrical geometry, and other magnetic substances including non-heme iron are diffusely distributed in the tissue, the qBOLD method simulates the gradient echo amplitude signal of the voxel as follows:

[0318]

[0319] where G(t) is the macroscopic field inhomogeneity contribution of GRE at time t, and F BOLD The GRE signal is attenuated due to the presence of deoxygenated blood in the vascular network within the voxel: F BOLD (Y, v, χ nb , t) = exp[-v·f s (δω·t)], where f s is the signal of vascular network attenuation, and δω is the characteristic frequency caused by the difference in magnetic susceptibility between deoxygenated blood and surrounding tissue:

[0320]

[0321] γ gyromagnetic ratio ( MHz / T), B0 is the main magnetic field (3T in this study), Hct hematocrit Δχ0 is the difference in magnetic susceptibility between fully oxygenated and fully deoxygenated red blood cells (4π×0.27ppm), χ ba The magnetic susceptibility of fully oxygenated red blood cells is calculated using Estimated to be -108.3ppb.

[0322] The optimization problem of fitting the data to Equation 43 is nonconvex, and F BOLD Display v, Y, and χ nb High coupling between. Accurate calculation of v, Y, and χ nb depends on the initial value and is highly sensitive to the signal-to-noise ratio. The sensitivity of qBOLD to noise propagates to the combined QSM+qBOLD model ( Figure 27 a). Therefore, a high signal-to-noise ratio is required to obtain accurate OEF.

[0323] In order to increase the effective signal-to-noise ratio before fitting, we propose a new cluster analysis method for QSM+qBOLD time evolution (hereinafter referred to as "CAT QSM+qBOLD"). qBOLD Voxels with similar tissue parameters (Y, v, R2) in (t) / G(t) form a cluster. Assuming that a cluster contains a large number of voxels, and the cluster represents the time evolution of all voxel signals, it can be used as a constraint to make the QSM+qBOLD maximum intensity resist noise. K-means can be used to identify all clusters. First, by assuming that the parameters Y, v, and R2 are constant within each cluster, and S0 and χ nb is different between voxels, and the utilization depends mainly on S0 and χ nb The magnitude and phase signals of the QSM+qBOLD problem are solved by cluster-based optimization. Then, a voxel-based optimization is performed by using the solution from the cluster optimization as the initial value.

[0324] The formula for joint QSM+qBOLD optimization is:

[0325]

[0326] where w is the weight of the first QSM term:

[0327]

[0328] χ is the magnetic susceptibility calculated from QSM, χ ba =ψ Hb ·χ oHb +(1-ψ Hb )·χ p The magnetic susceptibility is measured for fully oxygenated blood, α is the ratio of venous blood to total blood volume v / CBV (assumed to be a constant of 0.77), and ψ Hb is the volume fraction of hemoglobin (set to 0.0909, ), χ p is the magnetic susceptibility of plasma (set to -37.7 ppb), Δχ Hb is the difference in magnetic susceptibility between deoxygenated and oxyhemoglobin is arterial oxygen saturation. The third term is the physiological regulation on Y, so that OEF is averaged over the whole brain. Should be used with OEF wb Similarly, the OEF of the whole brain was estimated from the main draining vein [straight sinus (SS)]. ss =1-Y ss / Y a , where Y ss is the magnetic susceptibility of the straight sinus (SS), ψ in Equation 46 Hb =0.1197, v=1, χ nb = 0. In the formula where S(t) is the measured GRE data. w and λ are the optimized intensities determined by the L-curve method.

[0329] Materials and methods

[0330] simulation

[0331] Numerical simulations were performed to investigate the relationship between SNR and the dependence of OEF on the initial value. Using Equation 43 and Equation The GRE signal and QSM value were simulated separately. The following input parameters were used: Y = 60%, v = 3%, χ nb = -0.1 ppm, S0 = 1000 au, R2 = 20 Hz. The same 7 TEs obtained in healthy subjects (see below) were used: TE1 / ΔTE / TE7 = 2.3 / 3.9 / 25.8 ms. Gaussian noise was then added to the simulated GRE signal and QSM values, resulting in SNRs of: ∞ (noise-free), 1000, 100, and The following method is then used to estimate Y and v. In order to test the dependence on the initial value, different initial values ​​are used. and v To optimize the guess and use the true value as S0, χ nb and the initial value of R2. Repeat the example to calculate the simulated noise. w / λ=5×10 -3 / 0. The relative error is calculated as

[0332] In order to compare the accuracy of the proposed CAT QSM+qBOLD with the previous QSM+qBOLD, we performed simulation calculations. First, we used Equation 43 and Equation 44, respectively. The GRE signal and QSM value of each brain voxel are simulated. The input is the CAT QSM+qBOLD result of a stroke patient. Brain imaging was performed on the same day. We used the same eight TEs obtained from stroke patients (see below): Gaussian noise is added to the simulated GRE signals and QSM values, resulting in datasets with SNRs: ∞ (no noise), 1000, 100, and The simulated data were processed using two methods: 1) the previous QSM+qBOLD method with a constant initial OEF value for the whole brain and 2) the CAT QSM+qBOLD method. The same optimization settings were used in healthy subjects and stroke patients. wb Set to the average true OEF of the whole brain (29%). For CAT QSM+qBOLD, w / λ=5×10 -3 / 10 3 (L-curve analysis from 2 healthy subjects and 2 stroke patients.) The root mean square error (RMSE) was calculated to measure the difference from the true value.

[0333] In vivo live imaging:

[0334] Healthy subjects: This study was approved by the local institutional review board. Healthy volunteers (n = 11; 10 males, 1 female, mean age 34 ± 12 years) were recruited and underwent brain magnetic resonance imaging on a 3T scanner (HDxt, GE Healthcare) using an 8-channel head receiving coil. After signing the informed consent, all subjects avoided caffeine or alcohol for 24 hours before MRI. Magnetic resonance imaging was performed in the resting state using a 3D fast spin echo (FSE) arterial spin labeling (ASL) sequence, a 3D multi-echo spoiled gradient echo (GRE) sequence, and an inversion recovery T1w SPGR sequence (BRAVO). The 3DFSE ASL sequence parameters were: FOV 20 cm, in-plane resolution Layer thickness Marking cycle Post-mark delay time bandwidth / pixel, 8 interwoven spiral samples, each Reading points, scanning layers Layer axial data, TE10.1ms, The sum signal average is 3. The 3D GRE sequence parameters are: in-plane resolution 0.78 mm, slice thickness 1.2 mm, the same volume coverage as the 3D FSE ASL sequence, 7 equally spaced echoes, the first TE is 2.3 ms, and the echo interval is 3.9 ms. Bandwidth 488.3Hz / pixel, flip angle The pulse sequence was flow compensated in all three directions. The inversion recovery T1w SPGR sequence parameters were: in-plane resolution 0.78 mm, slice thickness 1.2 mm, volume coverage the same as the 3D FSE ASL sequence, TE 2.92 ms, Preparation time bandwidth / pixel, flip angle

[0335] The QSM reconstruction steps are as follows: First, the total magnetic field is estimated using an adaptive quadratic fit of the GRE phase map. Second, the local magnetic field is obtained using the map to dipole field (PDF) method. Finally, the magnetic susceptibility values ​​are calculated using the MEDI method with an automatic homogeneous cerebrospinal fluid reference (MEDI+0) algorithm. CBF maps (ml / 100 g / min) are calculated from the ASL data using the FuncTool software package (GE Healthcare, Waukesha, WI, USA). All images are simultaneously registered and interpolated into the QSM map using the FSL FLIRT algorithm.

[0336] Stroke patients: Ischemic stroke patients were randomly assigned to a 32-channel brain receiving coil on a clinical 3T scanner (GEMR Discovery ) were scanned with 3D ASL, 3D GRE and DWI sequences. The time interval between stroke onset and MRI examination was Hours to 12 days. All lesions were located in the unilateral cerebral artery supply area. 3D FSE ASL sequence parameters were: FOV 24cm, in-plane resolution 1.9mm, slice thickness 2.0mm, marking cycle Post-mark delay time bandwidth / pixel, number of scanning layers layer, TR 4787ms, signal average value 3. 3D GRE sequence parameters are: in-plane resolution 0.47mm, slice thickness 2mm, same volume coverage as 3D FSE ASL sequence, 8 equally spaced echoes, first TE The echo interval is 4.9ms, TR is 42.8ms, bandwidth is 244.1Hz / pixel, and flip angle is 20°. The DWI sequence parameters are: FOV 24cm, in-plane resolution 0.94mm, slice thickness 3.2mm, bandwidth b value 0, 1000s / mm 2 , TE71ms, TR3000ms, signal average value 4.

[0337] The processing of QSM and CBF in stroke patients was the same as that in healthy subjects, except that a linear fit of the GRE phase was used to estimate the total magnetic field because 3D flow compensation was not used on the scanners used in the stroke patient studies.

[0338] Cluster analysis:

[0339] The amplitude data of the GRE signal S(t) obtained after macroscopic field inhomogeneity processing is used for cluster analysis, and G is removed. The GRE signal of each voxel is the signal average over all echoes. The k-means cluster analysis algorithm is then applied to the standardized amplitude signal. Each voxel can be considered as a point in a vector space whose dimension is the number of echoes. Points close to each other (voxels with similar signal evolution) fall into the same cluster. The cluster center represents the representative signal evolution. The squared Euclidean distance is used to measure the proximity of points in the vector space (how many similar signal evolutions the voxels have). The appropriate number of clusters is automatically determined by the x-means method, which is an improved fast k-means method. First, it performs conventional k-means with a given initial cluster. Secondly, each centroid is replaced by two sub-centroids, and local k-means (k=2) is performed within each cluster using the sub-centroids as initial values. This is done for each cluster obtained in the first step. Third, a decision is made whether each cluster should be replaced with its two subclusters based on the Bayesian Information Criterion (BIC), which is the sum of the clustering log-likelihood and the number of clusters. For example, if the BIC values ​​of the two child nodes are larger than the BIC value of the parent node, then the child cluster should replace the parent cluster. As the number of clusters increases, the goodness of fit (log-likelihood) also increases, but this can lead to overfitting of the noise. The number of clusters in the BIC is adjusted to reduce this possibility. The x-means algorithm terminates when no more clusters can be split or when the number of clusters reaches a maximum value set a priori. Then, the regular k-means is repeated one last time, setting the number of clusters to the value obtained with the x-means method. In this study, x-means uses 1 to represent the initial number of clusters, Indicates the maximum number of clusters.

[0340] For speed, the x-means algorithm to obtain the number of clusters K was performed on 10% of the total voxels selected randomly. This process was repeated 10 times, and the K with the largest BIC value among the 10 trials was selected. The corresponding centroid was then used as the initial centroid for the final k-means on all voxels.

[0341] optimization:

[0342] Optimization of QSM+qBOLD (Formula ) is solved by dividing it into the following optimization problem: Update S0 based on qBOLD (Eq. 47); then χ nb Based on the QSM+qBOLD optimization (Equation 48) and the calculated magnetic susceptibility map χ; and update the Y, v, R2 values ​​according to the QSM+qBOLD optimization (Equation 49). Using k to represent the number of iterations, these steps are:

[0343]

[0344]

[0345]

[0346] Equation 47 is solved using a closed-form expression. Equations 48 and 49 are then solved using an iterative method (see below).

[0347] First, get a rough initial guess of Y, v, and R2 by following these steps: From the formula OEF wb Y0 was estimated. SS was automatically obtained using global and regional thresholds on QSM and position (posterior inferior part of the brain) and spatial (straight sinus) constraints. For the initial value of v (v0), the whole brain was roughly divided into three parts by FSL FAST using T1w (11 healthy subjects and 4 stroke patients) or T2-FLAIR images (1 stroke patient without T1w images), namely gray matter (GM), white matter (WM), and cerebrospinal fluid (CSF). According to the literature, v0 of GM / WM / CSF was set to 3 / / 1%. Set χ nb，0 To satisfy the formula with Y0 and v0 The initial value S is obtained by solving formula 43 0，0 and R 2，0 , and Y, v and χ nb Fit to these rough initial values. Single exponential fitting was performed using ARLO. Before fitting, 3D Gaussian smoothing was performed on S and G to improve SNR (with voxel diagonal length R2>100Hz or Voxels with are considered outliers for all subsequent processing.

[0348] Secondly, a cluster-based optimization is performed, where the unknowns Y, v, R2 are assumed to be constant within each cluster. The average of the voxel-based initial values ​​of Y, v, R2 is used as the initial value for the cluster-based optimization. In order to improve the convergence behavior during the nonlinear fitting, the unknowns Y, v, χ nb , R2 is considered to have approximately the same amplitude signal: where x is the unknown in the original scale, c is Y, v, χ nb , R2 proportional factors: respectively |χ nb，0 |,avg(R2,0)+4·SD(R2,0). avg(R 2，0 ) and SD(R 2，0 ) represent the R 2，0 The mean and standard deviation of . In the scaled Y, 0.4 and 2 initial values ​​v, and Before scaling the R_2, the lower and upper limits are set to 0.0 and 0.98. nb , set the lower and upper limits to the values ​​from the formula Calculated χ nb The values ​​are Y / v = 0.98 / 0.1 and 0.0 / 0.1 respectively. The optimization is performed on all clusters together. At the beginning of each optimization, the qBOLD term in Equations 48 and 49 is given by Normalization is performed to compensate for the metrics of the input MRI data, where is the average value of the first echo amplitude in the whole brain, N voxel is the number of voxels, N TE is the number of echoes. The QSM term in equations 48 and 9 is also expressed by 2 Normalization. The normalization weighting factor (λ) and the weight of QSM (w) are selected by performing L-curve analysis: first λ is selected when w = 0, and then w is defined by the previously determined λ. The limited memory Broyden-Fletcher-Goldfarb-Shanno-Bound (L-BFGS-B) algorithm is used for constrained optimization. When the relative residual Less than 10 in formulas 48 and 49 -5 The optimization is stopped when E n is the energy of the nth iteration. For the outer kth loop, the optimization is done at r k <10 -3 To prevent L-BFGS-B from not updating the Hessian because the residual is lower than the preset threshold, the cost function is multiplied by a factor of 10 before L-BFGS-B starts. 4 .

[0349] Third, the QSM+qBOLD optimization was repeated, but now Y, v, and R2 were allowed to vary across voxels. The cluster-based results were used as initial values. The method was the same as for the cluster-based optimization. The lower and upper bounds of the recalculated unknowns were set to 0.7 and 1.3 of the initial values, except that Y was set to 0.0 and 0.98 before the calculation. The same L-BFGS-B algorithm was used as in the cluster-based optimization. For Equations 8 and 9, when the relative residual was less than 2×10 -4 The optimization stops when the relative residual is less than 10 for the outer loop. -2 , the optimization stops.

[0350] For numerical simulation 1, the optimization for healthy subjects and stroke patients uses the same optimization settings, except that fixed lower and upper limits of v, 0.01 and 0.1, are set before the calculation.

[0351] The CAT QSM+qBOLD method was compared with a previous QSM+qBOLD method with a constant initial OEF value for the whole brain (hereinafter referred to as "previous QSM+qBOLD"). For the QSM+qBOLD method based on the optimized constant initial OEF value (previous QSM+qBOLD), we followed the optimization method described above. Based on the L-curve analysis, the weight of QSM(w) was set to 100 for healthy subjects and stroke patients. When the relative residual of the healthy subjects was less than And for stroke patients, the optimization was stopped when the value was 0.001, and the lower and upper limits were set the same as those of the cluster-based method.

[0352] All algorithms were implemented in Matlab (Mathworks Inc., Natick, MA).

[0353] Statistical analysis:

[0354] ROI analysis (mean and standard deviation) and paired t-tests were performed to compare CMRO2 and OEF values ​​between the previous QSM+qBOLD and CAT QSM+qBOLD. For the ROIs of healthy subjects, the cortical gray matter (CGM) was reconstructed based on T1-weighted images by an experienced neuroradiologist (SZ 7 years of experience). For stroke patients, the area on the lesion side and its corresponding contralateral side was delineated based on DWI by an experienced neuroradiologist (SZ 7 years of experience). To investigate the resulting dependence on the number of clusters, a conventional k-means method was performed, where and x-mean results (K = 13 in healthy subjects and K = 13 in stroke patients) Onset d). CAT QSM+qBOLD was then performed for each K. The same optimized protocol as used for healthy subjects and stroke patients was used, including To investigate the differences in OEF between different K values, repeated measures ANOVA was performed.

[0355] result

[0356] Among the optimal number of clusters determined by the x-means method, there were 11 healthy subjects and 11 stroke patients. For example, the mean of the 10% subsampling plan and the 100% sampling plan is < 1 A 10% subsampling plan (consisting of 10 trials) is faster than a 100% subsampling plan. The optimal number of clusters for the x-means method was 10±2 (N=11) for healthy subjects and 10±2 (N=11) for stroke patients.

[0357] In the L-curve analysis, four randomly selected subjects (two healthy subjects and two stroke patients) were randomly selected at λ = 1000 and w = 5 × 10 -3 Determined under the circumstances. Figure 27 a is the effect of the signal-to-noise ratio on the calculated Y sensitivity to the initial value in the numerical simulation. In the absence of noise, the relative error is small, but as the signal-to-noise ratio decreases from 1000 to When the initial value deviates from the true value, the relative error tends to increase. For example, when the signal-to-noise ratio is When Y0=0.6 and v0=0.03, the relative error is However, the relative error is 29.1% when Y0=0.1 and v0=0.03.

[0358] Figure 27 b shows a comparison of the oxygen uptake fraction maps obtained using the previous QSM+qBOLD and CAT QSM+qBOLD during a stroke simulation. CAT QSM+qBOLD provides more accurate oxygen uptake fraction maps than the previous QSM+qBOLD, especially in low signal-to-noise ratio conditions. For example, in a signal-to-noise ratio of CAT QSM+qBOLD captures regions of low oxygen uptake fraction, whereas the previous QSM+qBOLD is insensitive to low oxygen uptake fraction values. CAT QSM+qBOLD provides lower root mean square error than the previous QSM+qBOLD for all signal-to-noise ratios.

[0359] Figure 28 Comparison of the previous QSM+qBOLD and CAT QSM+qBOLD. The oxygen uptake fraction is less noisy and more uniform, whereas the previous QSM+qBOLD is very noisy with extreme values, such as >80% in dark gray matter. CAT QSM+qBOLD shows good contrast in cerebral oxygen metabolism between CGM and WM, without the extreme values ​​seen in the previous method. In the new method, v, which shows the contrast between cortical gray matter and white matter, is generally lower than in the previous QSM+qBOLD values.

[0360] Figure 29 A stroke patient (post-stroke) was shown using the before and CAT QSM+qBOLD methods. day) oxygen uptake fraction, cerebral oxygen metabolic rate, v, R2 and χ nb Figure. CAT QSM+qBOLD can more clearly identify lesions in oxygen uptake fraction and cerebral oxygen metabolic rate maps. CAT QSM+qBOLD low oxygen uptake fraction areas are clearly contained within the lesions shown on diffusion-weighted imaging. However, previous QSM+qBOLD did not show obvious localized low oxygen uptake fraction areas, either within or outside the lesions defined by diffusion-weighted imaging. CAT QSM+qBOLD generally showed low v-specific areas of lesions clearly defined by diffusion-weighted imaging, while previous QSM+qBOLD methods showed similar v values ​​compared with cerebral blood flow methods. CAT QSM+qBOLD showed R2 and χ nb The values ​​are generally higher than those shown previously for QSM+qBOLD.

[0361] Figure 30 shows the histogram of oxygen uptake fraction values ​​of the lesion and its contralateral area before and after CAT QSM+qBOLD in the second stroke patient (12 days after stroke). The distribution of oxygen uptake fraction of CAT QSM+qBOLD in the lesion is different from that in the contralateral area. The lesion has 8 peaks, of which the two strongest peaks are 0 and On the opposite side peaks, of which the dominant peak is However, previous QSM+qBOLD did not have a specific distribution of hypoxia uptake fraction values ​​in the lesion, but rather a bell-shaped distribution in both the lesion and contralateral side (wider in the contralateral side), with peak values ​​at 47% and 49%, respectively.

[0362] exist Figure 31 In the figure, a series of clustering results are shown ( x means the difference between healthy subjects and stroke patients (post-stroke Brain segmentation and final oxygen uptake fraction map of the patients with onset of disease on the first day of the first month. The obtained oxygen uptake fraction map is for all patients with disease greater than The K values ​​of the healthy subjects showed similar performance, and there was no significant difference between the K values: p = 1.0000, Stroke patients p = 0.9999, (Repeated measures ANOVA). The x-mean clusters were automatically selected with K = 13 for healthy subjects and K = 13 for stroke patients.

[0363] Figure 32 The results were compared between the previous QSM+qBOLD and CAT QSM+qBOLD analysis of cortical gray matter (CGM) in healthy subjects. CAT QSM+qBOLD showed a smaller oxygen uptake fraction, cerebral oxygen metabolic rate, and v value compared to the previous QSM+qBOLD: oxygen uptake fraction 32.7±4.0% and 37.9±3.1% (p<0.01), cerebral oxygen metabolic rate 148.4±23.8, 171.4±22.4μmol / 100g / min (p<0.01), and v value 1.00±0.2% and (p<0.01). CAT QSM+qBOLD values ​​were higher, respectively 13.1±0.7Hz (p<0.01), -20.2±8.1ppb, -33.8±9.0ppb (p<0.01).

[0364] in conclusion

[0365] Our results demonstrate that the QSM+qBOLD method, based on cluster analysis of temporal evolution (CAT), provides more homogeneous and less noisy oxygen uptake fraction maps in healthy subjects than previous QSM+qBOLD methods. Regions of low oxygen uptake fraction were contained within lesions defined on diffusion-weighted imaging in stroke patients, whereas this clear representation was not observed with previous QSM+qBOLD methods. By significantly improving the effective signal-to-noise ratio using a clustering approach, the previous QSM+qBOLD model achieved more accurate oxygen uptake fraction maps in numerical simulations. Finally, CAT QSM+qBOLD estimates oxygen uptake fraction solely from gradient echo data, whereas previous QSM+qBOLD methods rely on oxygen uptake fraction measurements.

[0366] Compared with the previous QSM+qBOLD method, the oxygen uptake fraction maps of healthy subjects obtained by the new method are more uniform and have smaller extreme values ​​( Figure 28 ), which is consistent with previous positron emission tomography (PET) studies. The reason for the noise suppression may be the use of cluster optimization, in which the unknowns are assumed to be constant across the cluster, thus creating an effective signal average. Numerical simulation results show that the high signal-to-noise ratio makes the estimated parameters insensitive to measurement errors. Compared with the previous QSM+qBOLD method, the overall noise reduction effect propagated to Y and v is comparable to that mapped to χ nb Same noise reduction effect as R2 ( Figure 28 ).

[0367] CAT QSM+qBOLD showed a smaller oxygen uptake fraction value ( Figure 27 c and 32). CAT and previous QSM+qBOLD in cortical gray matter were 32.7±4.0% and 37.9±3.1%, respectively (p<0.01). Both oxygen uptake fraction values ​​are within the range of oxygen uptake fractions previously reported using PET: and 40 ± 9%, and other oxygen uptake fraction values ​​obtained by MRI-based techniques: The smaller the oxygen uptake score value of the CAT QSM+qBOLD method, the higher the corresponding oxygen uptake score value ( vs. 13.1 ± 0.7 Hz). This is expected, since to obtain the same measured amplitude signal attenuation, the oxygen uptake fraction decreases with increasing θ (Equations 43 and 4).

[0368] The values ​​of the CAT QSM+qBOLD method in healthy subjects were smaller than those of the previous QSM+qBOLD method: 1.00±0.2% for cortical gray matter and (p < 0.01). Compared with PET and other magnetic resonance techniques, the v previously calculated from CAT and QSM + qBOLD are slightly smaller and slightly larger, respectively: the v values ​​obtained using PET, e.g. 2.0 ± 0.2%, other magnetic resonance techniques: The new method also has larger cortical gray matter values, which are and 13.1±0.7Hz, compared with the values ​​calculated by other magnetic resonance techniques of 14.9±0.2Hz, The consistency of 17.1±2Hz is good.

[0369] In patients with ischemic stroke, the areas of low oxygen uptake scores obtained by CAT QSM+qBOLD were mainly contained within the lesions defined by diffusion-weighted imaging ( Figure 29 and 30), however, the oxygen uptake fraction maps obtained by the previous QSM+qBOLD method cannot identify lesions. The low area of ​​CAT QSM+qBOLD roughly corresponds to the location of the lesion; this observation is consistent with the decrease in blood volume that occurs in ischemic stroke lesions. While the previous QSM+qBOLD method used a constant initial value of oxygen uptake fraction, the contrast of the v map was similar to that of the cerebral blood flow map ( Figure 29 ). This shows that the v value results are not much different from the original hypothesis, which was based on the phenomenological relationship based on cerebral blood flow.

[0370] CAT QSM+qBOLD showed lesions in the oxygen uptake fraction histogram 12 days after symptom onset in stroke patients, which were different from the contralateral area (Figure 30): In the new method, there were two strongest oxygen uptake fraction peaks, 0% and and several high oxygen uptake fraction peaks This low oxygen uptake fraction area may represent an area of ​​tissue necrosis (oxygen uptake fraction <10%), while another area of ​​the brain may be salvageable. Previous QSM+qBOLD results showed no specific distribution of low oxygen uptake fraction values ​​within lesions. These results suggest that the use of CAT makes QSM+qBOLD increasingly sensitive to expected low oxygen uptake fraction values ​​within stroke lesions.

[0371] The CAT QSM+qBOLD method described can be further improved. The appropriate number of clusters is automatically selected using the k-means method. This is based on a well-known algorithm, BIC. However, the optimal number of clusters may differ when using different algorithms (e.g. AIC). In general, clustering results may differ when using different clustering methods, such as hierarchical clustering. This may affect the generated oxygen uptake fraction maps. In our study, the oxygen uptake fraction maps were not sensitive to the number of clusters selected by k-means clustering ( Figure 31 Performing voxel-based optimization after cluster-based optimization can alleviate the consequences of incomplete cluster analysis.

[0372] Optimization problem formulation It can be solved by deep learning. Artificial neural networks (ANN), such as the U-Net mentioned above, can simulate data in a network covering Y, v, R2, S0, χ nb All these parameters are trained within the range of possible pathophysiological values. The ANN simulates complex data by adding Gaussian noise to the simulated data, using equations related to MRI signals such as MRI and MRI-related signals. Simulated complex multi-echo MRI data obtained from a patient is approximated and fed into the ANN, serving as the ANN's MRI data "label" output. Deep learning studies show that the numerous weights in a neural network can learn a sparse representation of the MRI signal amplitude versus echo time, along with other features, resulting in optimal denoising or robustness to noise. The trained neural network can generate oxygen uptake scores using m-gradient echo datasets acquired from subjects.

[0373] To improve the accuracy of calculations of cerebral metabolic rate of oxygen, it is also necessary to improve the accuracy of measurements of cerebral blood flow. The arterial spin labeling used in this study to measure cerebral blood flow has lower resolution (compared to gradient echo acquisition) and is less accurate in cortical white matter. The OEF and estimates of large veins may be inaccurate because they are treated the same as normal brain tissue, and the oxygen uptake fraction and estimates of large veins can be alleviated by setting v = 1 for large veins. However, this will require delineating the large veins before calculation and using a thresholding method on the magnetic susceptibility map, as was done in this study to delineate the transverse sinus (SS). The optimization of CAT QSM+qBOLD is still nonlinear, which means that the fusion may affect the solver implementation, parameter calibration, and stopping criteria. Finally, no true values ​​or reference measurements are available, so it will need to be performed on the PET-MR scanner. Studies have shown its accuracy in vivo.

[0374] In summary, this study demonstrated the feasibility of temporally evolving cluster analysis (CAT) for QSM+qBOLD in healthy subjects and patients with ischemic stroke by effectively improving the signal-to-noise ratio. Numerical simulation results showed that this method achieved higher accuracy than existing QSM+qBOLD methods. CAT QSM+qBOLD provided a lower-noise and more uniform oxygen uptake fraction in healthy subjects. In patients with ischemic stroke, regions of low oxygen uptake fraction were contained within the stroke lesions. CAT QSM+qBOLD can be used to study tissue viability in a variety of diseases, such as Alzheimer's disease, multiple sclerosis, tumors, and ischemic stroke.

[0375] Equipment system application

[0376] In some cases, one or more of the quantitative mapping techniques described above may be implemented using process 800 shown in FIG. 8 . Figure 33 As an example, process 800 can be used to map the tissue magnetic susceptibility of an object, such as a patient (or a portion of a patient) or an experimental sample (e.g., an imaging phantom, a material or tissue sample, etc.). In some embodiments, process 800 can be used to convert magnetic resonance (MR) signal data corresponding to the object into a plurality of images that quantitatively depict the structure and / or composition and / or function of the object. For example, in some cases, process 800 can be used to obtain multi-echo MR data corresponding to the object and process the multi-echo MR data to generate a quantitative magnetic susceptibility map of the object. As another example, in some cases, process 800 can be used to obtain fast undersampled multi-contrast MR data corresponding to the object and process the data to generate a multi-contrast MR image. As another example, in some cases, process 800 can be used to obtain time-resolved MR data corresponding to a subject as a contrast agent is passed through tissue in an organ of the subject and process the time-resolved MR data to generate a quantitative subject traffic map. Using these physical properties and the plurality of contrast quantitative maps, one or more images of the object can be generated and displayed to a user. The user can then use these images for diagnostic, therapeutic, or experimental purposes, such as to study the structure and / or composition and / or function of an object, and / or diagnose various conditions or diseases, and / or treat various conditions or diseases based at least in part on the images. Because process 800 can produce mappings of physical quantities with higher quality and / or accuracy than other mapping techniques, implementation of process 800 can be used to improve a user's understanding of the structure and / or composition and / or function of an object, and can be used to improve the accuracy of any resulting medical diagnosis, treatment, or experimental analysis.

[0377] Process 800 begins by acquiring MR data corresponding to an object (step 810). In some cases, the MR data may correspond to a patient (e.g., the entire patient or a portion of the patient, such as a specific portion of the patient's body). In some cases, the MR data may correspond to one or more samples (e.g., imaging phantoms, samples of one or more materials, samples of one or more types of tissue, and / or one or more other objects).

[0378] exist Figure 33 In some embodiments shown in a, an MRI scanner can be used to acquire MR data using one or more suitable pulse sequences. For example, in some cases, a gradient echo sequence can be used to acquire MR data, which acquires MR data at a single echo time or at multiple echo times (e.g., two, three, four, five, etc.). Various scan parameters can be used. As another example, different echo times (TE) (e.g., 4TE, such as 3.8, 4.3, and ) in an interleaved manner with the following imaging parameters: echo time (TR) = 22 ms; voxel size = 0.5 × 0.5 × 0.5 mm 3 , Flip angle = 15°. As another example, MR data can be acquired using a 3D multi-gradient echo sequence on a 3T scanner with imaging parameters Although example sequences and example parameters are described above, these are merely illustrative. In practice, other sequences and parameters may be used, depending on various factors (e.g., the size of the area to be inspected, the known or assumed sensitivity range of the object, scan time considerations, equipment limitations, etc.).

[0379] After acquiring MR data corresponding to the object, process 800 continues by determining a magnetic field based on the MR data (step 820). As described above, in some cases, MRI signal phase is affected by a magnetic field corresponding to the object's magnetic susceptibility distribution and chemical shifts of tissue components.

[0380] The magnetic field can be determined from complex MRI data by fitting the detected signal as a sum of tissue spectral components, where each component signal is characterized by an exponential, with a negative real part representing signal attenuation and an imaginary part representing the magnetic field-dependent phase field and chemical shift (step 830). This fitting of such a complex signal model can be performed iteratively using numerical optimization. A robust initialization of the numerical optimization can be estimated using a simplified single-species model and graph cuts to separate the smoothed field from the chemical shift and phase unwrapping.

[0381] After determining the magnetic field based on the MR data, process 800 continues by determining the relationship between the magnetic field and the magnetic susceptibility (step 840). As described above, in some embodiments, the relationship between the magnetic field and the magnetic susceptibility can be expressed as a relationship between the magnetic field at a given location and the magnetic susceptibility at that location. In some embodiments, this field-susceptibility relationship can be expressed in integral form or differential form. In differential form, the relationship between the magnetic field at a given location and the magnetic susceptibility at that location can include an equation where the Laplace operator of the magnetic field is equal to one-third of the Laplace operator of the magnetic susceptibility minus the second derivative of the magnetic susceptibility, the magnetic susceptibility. In some implementations, this field-susceptibility relationship can be represented by weights of an artificial neural network.

[0382] After determining the relationship between magnetic field and magnetic susceptibility, process 800 continues by determining a priori knowledge about the tissue magnetic susceptibility distribution (step ). One estimate of the tissue susceptibility distribution is to use R2* values ​​derived from the amplitude signal. High R2* values ​​can be used to identify high sensitivity regions, such as hemorrhage, for preprocessing to accelerate convergence of the numerical optimization. An example of preprocessing is to divide the entire image volume into regions with high sensitivity (including hemorrhage, air regions, and background) and regions with normal tissue sensitivity. This preprocessing reduces the sensitivity search range, thereby accelerating convergence. Another estimate of the tissue sensitivity distribution is to use low (near zero) R2* regions to identify pure water regions, such as cerebrospinal fluid in the ventricles of the brain, oxygenated arterial blood in the aorta, or pure fat regions in the abdominal wall. Since water and fat have known susceptibility values, they can be used to provide a zero reference (to water) to generate absolute susceptibility values ​​using minimum variance regularization. After obtaining the magnetic susceptibility prior information and the signal-data noise characteristics, process 800 continues by estimating the magnetic susceptibility distribution of the object based at least in part on the prior information and the data noise characteristics (step ). In some embodiments as described above, estimating the magnetic susceptibility distribution of an object by optimization may include determining a cost function corresponding to the magnetic susceptibility distribution, the magnetic field, and a mask corresponding to the region of interest. The cost function includes a data fidelity term based on noise characteristics of the data, a relationship between the magnetic field and the tissue susceptibility in integral or differential form, and a regularization term representing prior information. The estimated magnetic susceptibility distribution of the object can be determined by identifying a specific magnetic susceptibility distribution that minimizes one or more of these cost functions. As described above, in some embodiments, this can be determined numerically using a quasi-Newton, an alternating direction method of multipliers, or an artificial neural network.

[0383] After estimating the magnetic susceptibility distribution of the object, process 800 continues by generating one or more images of the object based on the estimated magnetic susceptibility distribution of the object (step 870). These images can be displayed electronically on a suitable display device (e.g., an LCD display device, an LED display device, a CRT display device, or other display device that can be configured to display images) and / or physically displayed on a suitable medium (e.g., printed, etched, painted, embossed, or otherwise physically presented on paper, plastic, or other material). In some cases, a color scale can be used to visualize the estimated magnetic susceptibility distribution of the object, wherein each of several colors is mapped to a specific magnetic susceptibility value or a range of magnetic susceptibility values. Thus, a two-dimensional or three-dimensional image can be generated in which each pixel or voxel of the image corresponds to a specific spatial location of the object, and the color of the pixel of the voxel depicts the magnetic susceptibility value of the object at that location. In some cases, the color scale can include a color gradient and / or a grayscale gradient (e.g., grayscale) to depict a range of magnetic susceptibility values. In some cases, for a color scale comprising a color gradient or a grayscale gradient, one end of the gradient can correspond to the lowest susceptibility value within a particular window of susceptibility values ​​(e.g., an arbitrarily selected window of values), while the other end of the gradient can correspond to the highest susceptibility value within that window of values. For example, for a grayscale comprising a grayscale gradient between pure white and pure black, pure white can be used to indicate the highest susceptibility value within that particular window of arbitrary values, while pure black can be used to indicate the lowest susceptibility value within that window of values. Other relationships between color / grayscale levels and susceptibility values ​​are also possible. For example, in some cases, for a grayscale comprising a grayscale gradient between pure white and pure black, pure white can be used to indicate the lowest susceptibility value within that window of arbitrary values, while pure black can be used to indicate the highest susceptibility value within that window of values. While the examples are described in the context of grayscale, similar relationships between color scales and susceptibility values ​​are also possible. Susceptibility values ​​and colors can be mapped linearly (e.g., each absolute change in susceptibility value corresponds to a proportional change in color), logarithmically (e.g., each exponential change in susceptibility value corresponds to a linear change in color), or according to any other mapping. Although examples of mappings between color scales and magnetic susceptibility values ​​are described above, these are merely illustrative examples. In practice, other mappings are possible, depending on the specific implementation.

[0384] exist Figure 33 In some embodiments shown in FIG. 2 , an MRI scanner can be used to repeatedly use one or more suitable pulse sequences in a time-resolved manner to acquire MR data during the passage of the contrast agent through the tissue. For example, in some cases, MR data can be acquired using a gradient echo sequence that acquires MR data at a single echo time with multiple time frames. Various scan parameters can be used. As another example, MR data can be acquired using a 3D spoiled gradient echo sequence on a 3T scanner (e.g., a GE Excite HD MR scanner) using the following imaging parameters: voxel size Matrix size Rotation angle Temporal resolution Time points 30-40. Although example sequences and example parameters are described above, these are merely illustrative. In practice, other sequences and parameters may be used, depending on various factors (e.g., the size of the area to be examined, the known or assumed sensitivity range of the object, scan time considerations, equipment limitations, etc.).

[0385] After acquiring MR data corresponding to the object, process 800 continues by determining a contrast agent concentration based on the MR data (step 821). As described above, in some cases, MRI signal amplitude and phase are affected by a magnetic field corresponding to a magnetic susceptibility distribution of a high paramagnetic contrast agent.

[0386] The contrast agent concentration for each time frame can be determined from the complex MRI data: QSM processing of the phase data, amplitude signal model fitting, or a combination of both (step 831). This fitting of the data to a complex signal model to extract agent concentration can be performed using numerical optimization.

[0387] After determining the concentration of the agent or tracer based on the MR data, process 800 continues by determining the relationship between the tracer concentration and the potential transport parameter (step 841). As described above, in some embodiments, the concentration-transport relationship can be expressed as the relationship between the concentration at a given location and the transport parameter at that location, including the temporal variation of the concentration and the divergence of the mass flux. In some cases, the mass flux is the concentration multiplied by the convection or drift velocity.

[0388] After determining the relationship between tracer concentration and transport parameters (including velocity, diffusion, permeability, and pressure gradient), process 800 continues by determining a priori knowledge about the distribution of tissue transport parameters (step One estimate of the distribution of tissue transport parameters in space is tissue morphology information. In some cases, the velocity distribution can be sparse, characterized by a cost function of the L1 norm of the gradient of the velocity spatial distribution. In some cases, the tissue transport parameter distribution can be characterized by an artificial neural network.

[0389] After obtaining the a priori information on the spatial distribution of transmission, process 800 continues by estimating the transmission distribution of the object based at least in part on the a priori information and the concentration data (step As described above, estimating the transmission distribution of an object can include determining a cost function corresponding to the distribution of transmission parameters (or transmission for simplicity) and the tracer concentration. The cost function includes a data fidelity term based on the noise characteristics of the lumped data and a regularization term representing prior information. The estimated transmission distribution of the object can be determined by identifying a specific transmission distribution that minimizes one or more of these cost functions. As described above, in some cases, this can be determined numerically using a quasi-Newton method, an alternating direction method of multipliers, or an artificial neural network.

[0390] After estimating the transmission distribution of the object, process 800 continues by generating one or more images of the object based on the estimated transmission distribution of the object (step 871). These images can be displayed electronically on a suitable display device (e.g., an LCD display device, an LED display device, a CRT display device, or other display device that can be configured to display images) and / or physically displayed on a suitable medium (e.g., printed, etched, painted, embossed, or otherwise physically presented on paper, plastic, or other material). In some cases, the estimated transmission distribution of the object can be visualized using a color scale, where each of several colors is mapped to a specific transmission value or range of transmission values. Thus, a two-dimensional or three-dimensional image can be generated, where each pixel or voxel of the image corresponds to a specific spatial location of the object, and the color of the pixel of the voxel depicts the transmission value of the object at that location. In some cases, the color scale can include a color gradient and / or a grayscale gradient (e.g., grayscale) to depict a range of transmission values. In some cases, for a color scale that includes a color gradient or a grayscale gradient, one end of the gradient may correspond to the lowest transmission value within a particular window of transmission values ​​(e.g., an arbitrarily selected window of values), while the other end of the gradient may correspond to the highest transmission value within that window of values. For example, for a grayscale that includes a grayscale gradient between pure white and pure black, pure white may be used to indicate the highest transmission value within that window of values, while pure black may be used to indicate the lowest transmission value within that window of values. Other relationships between color / grayscale levels and transmission values ​​are also possible. For example, in some cases, for a grayscale that includes a grayscale gradient between pure white and pure black, pure white may be used to indicate the lowest transmission value within that window of values, while pure black may be used to indicate the highest transmission value within that window of values. While the examples are described in the context of grayscale, similar relationships between color scales and transmission values ​​are also possible. Transmission values ​​and colors can be mapped linearly (e.g., each absolute change in transmission value corresponds to a proportional change in color), logarithmically (e.g., each exponential change in transmission value corresponds to a linear change in color), or according to any other mapping. While examples of mappings between color scales and transmission values ​​are described above, these are merely illustrative examples. In practice, other mappings are possible, depending on the implementation.

[0391] exist Figure 33 In some embodiments shown in FIG. 8 , MR data can be acquired using an MRI scanner using one or more suitable pulse sequences that sensitize tissue to multiple contrasts (step 822). For example, in some cases, a T1-weighted image dataset can be used to acquire all of the MR data, but a small portion of the T2-weighted, T2 FLAIR, and diffusion-weighted image datasets can be used. Various scan parameters can be used. As another example, MR data can be acquired from T1-weighted imaging with a matrix size of 256×176 and a resolution of 1 mm. 3 , isotropic. Although example sequences and example parameters are described above, these are merely illustrative. In practice, other sequences and parameters may be used, depending on various factors (e.g., the size of the area to be examined, the known or assumed range of contrast values ​​for the object, scan time considerations, equipment limitations, etc.).

[0392] After acquiring MR data corresponding to the object, process 800 continues by determining prior knowledge about tissue structural information or morphology (step One estimate of tissue morphology consistency between various tissue contrast images is sparsity, characterized by a cost function of the L1 norm of the gradient of the velocity space distribution. In some cases, morphology consistency can be characterized by artificial neural networks.

[0393] After obtaining the tissue morphology prior information, process 800 continues by reconstructing images of the object at various contrasts based at least in part on the prior information and the undersampled data (step As described above, reconstructing an image of an object from undersampled noisy data can include determining a cost function that includes a data fidelity term based on noise characteristics of the concentration data and a regularization term that expresses prior information. An estimated transport distribution for the object can be determined by identifying a particular image that minimizes one or more of these cost functions. As described above, in some cases, this can be determined numerically using a quasi-Newton method, an alternating direction method of multipliers, or an artificial neural network.

[0394] After reconstructing images of the object at various weights, process 800 continues by displaying one or more images of the object (step 872). These images can be displayed electronically on a suitable display device (e.g., an LCD display device, an LED display device, a CRT display device, or other display device configured to display images) and / or physically displayed on a suitable medium (e.g., printed, etched, painted, embossed, or otherwise physically presented on paper, plastic, or other material). In some cases, the reconstructed image of the object can be visualized using a color scale, wherein each of several colors is mapped to a specific intensity value or range of intensity values. Thus, a two-dimensional or three-dimensional image can be generated, wherein each pixel or voxel of the image corresponds to a specific spatial location of the object, and the color of the voxel depicts the intensity value of the object at that location. In some cases, the color scale can include a color gradient and / or a grayscale gradient (e.g., grayscale) to depict a range of intensity values. In some cases, for a color scale including a color gradient or grayscale gradient, one end of the gradient can correspond to the lowest intensity value in a specific intensity value window (e.g., an arbitrarily selected value window), while the other end of the gradient can correspond to the highest intensity value in the value window. For example, for a grayscale that includes a grayscale gradient between pure white and pure black, pure white can be used to indicate the highest intensity value in a particular arbitrary value window, while pure black can be used to indicate the value in the lowest intensity value window. Other relationships between color / grayscale and intensity values ​​are also possible. For example, in some cases, for a grayscale that includes a grayscale gradient between pure white and pure black, pure white can be used to indicate the lower intensity value in a particular arbitrary value window, while pure black can be used to indicate the highest intensity value in the indicated value window. Although the examples are described in the context of grayscale, in practice, similar relationships between color scales and intensity values ​​are also possible. Intensity values ​​and colors can be mapped linearly (e.g., each absolute change in intensity value corresponds to a proportional change in color), logarithmically (e.g., each exponential change in intensity value corresponds to a linear change in color), or according to any other mapping. Although examples of mappings between color scales and intensity values ​​are described above, these are merely illustrative examples. In practice, other mappings are also possible, depending on the implementation.

[0395] Implementation of the above-described techniques may be performed using a computer system. Figure 34 is a block diagram of an example computer system 900 that can be used, for example, to perform an implementation of process 800. In some implementations, computer system 900 can be communicatively connected to another computer system (e.g., another computer system 900) so that it receives data (e.g., an MRI dataset) and analyzes the received data using one or more of the techniques described above.

[0396] System 900 includes a processor 910, a memory 920, a storage device 930, and an input / output device 940. Each of components 910, 920, 930, and 940 may be connected to a system bus, for example. The processor 910 is capable of processing instructions for execution within the system 900. In some implementations, the processor 910 is a single-threaded processor. In some embodiments, the processor 910 is a multi-threaded processor. In some embodiments, the processor 910 includes a graphics processing unit. The processor 910 is capable of processing instructions stored in the memory 920 or on the storage device 930. The processor 910 may perform operations such as executing one or more of the above-described techniques.

[0397] Memory 920 stores information within system 900. In some embodiments, memory 920 is a computer-readable medium. In some embodiments, memory 920 is a volatile memory unit. In some embodiments, memory 920 is a non-volatile memory unit.

[0398] The storage device 930 can provide mass storage for the system 900. In some embodiments, the storage device 930 is a non-transitory computer-readable medium. In various embodiments, the storage device 930 may include, for example, a hard disk device, an optical disk device, a solid-state drive, a flash drive, a magnetic tape, or some other mass storage device. In some embodiments, the storage device 930 may be a cloud storage device, for example, a logical storage device comprising multiple physical storage devices distributed across a network and accessed using the network. In some examples, the storage device may store long-term data. The input / output device 940 provides input / output operations for the system 900. In some embodiments, the input / output device 940 may include one or more network interface devices, such as an Ethernet card, a serial communication device, such as an RS-232 port, and / or a wireless interface device, such as an 802.11 card, a 3G wireless modem, a 4G wireless modem, etc. The network interface devices allow the system 900 to communicate, such as sending and receiving data. In some embodiments, the input / output devices may include driver devices configured to receive input data and send output data to other input / output devices, such as a keyboard, a mouse, a printer, sensors (e.g., sensors that measure component or system-related properties, sensors that measure environment-related properties, or other types of sensors) and display devices. In some implementations, mobile computing devices, mobile communication devices, and other devices may be used.

[0399] The computing system may be implemented by instructions that, when executed, cause one or more processing devices to perform the processes and functions described above, for example, to store, maintain, and display artifacts. Such instructions may include, for example, interpreted instructions, such as script instructions, or executable code, or other instructions stored in a computer-readable medium. The computing system may be implemented distributed over a network, such as a server farm or a group of widely distributed servers, or may be implemented in a single virtual device that includes multiple distributed devices operating in coordination with each other. For example, one of the devices may control the other devices, or the devices may operate under a coordinated set of rules or protocols, or the devices may be coordinated in another manner. Multiple distributed devices operate in coordination to give the appearance of a single device operating. Although in Figure 33 An example processing system has been described in the specification. Implementation of the subject matter and functional operations described above may be implemented in other types of digital electronic circuitry, or in computer software, firmware, or hardware, including the structures disclosed in this specification and their structural equivalents, or in a combination thereof to implement one or more of them. Implementation of the subject matter described in this specification, such as performing one or more of the processes described above, may be implemented as one or more computer program products, such as one or more modules of computer program instructions encoded on a tangible program carrier, such as a computer-readable medium for execution by a processing system or to control the operation of a processing system. The computer-readable medium may be a machine-readable storage device, a machine-readable storage substrate, a memory device, a composition of matter effecting a machine-readable propagated signal, or a combination of one or more of them.

[0400] The term "processing module" may encompass all devices, apparatuses, and machines for processing data, including, for example, a programmable processor, a computer, or multiple processors or computers. In addition to hardware, a processing module may also include code that creates an execution environment for the computer program in question, such as code constituting processor firmware, a protocol stack, a database management system, an operating system, or a combination thereof, or more.

[0401] A computer program (also referred to as a program, software, software application, script, executable logic, or code) can be written in any programming language, including compiled or interpreted languages, or declarative or procedural languages, and can be deployed in any form, including as a stand-alone program or as a module, component, subroutine, or other unit suitable for use in a computing environment. A computer program does not necessarily correspond to a file in a file system. A program can be stored as part of a file containing other programs or data (e.g., one or more scripts stored in a markup language document), as a single file dedicated to a related program, or as multiple coordinated files (e.g., files storing one or more modules, subroutines, or portions of code). A computer program can be deployed and executed on a single computer or on multiple computers located at one site or distributed across multiple sites and interconnected by a communications network. Computer-readable media suitable for storing computer program instructions and data include all forms of nonvolatile or volatile memory, media, and storage devices, including, for example, semiconductor memory devices such as EPROM, EEPROM, and flash memory devices; magnetic disks, such as internal hard disks or removable disks or tapes; magneto-optical disks; and CD-ROM and DVD-ROM disks. The processor and memory can be supplemented by, or incorporated into, special-purpose logic circuitry. Sometimes the server is a general-purpose computer, sometimes it is a customized, special-purpose electronic device, and sometimes it is a combination of these. An implementation may include a back-end component, such as a data server, or a middleware component, such as an application server, or a front-end component, such as a client computer with a graphical user interface or a web browser, through which a user can interact with an implementation of the subject matter described in this specification, or any combination of one or more such back-end, middleware, or front-end components. The components of the system may be interconnected by any form or medium of digital data communication, such as a communication network. Examples of communication networks include local area networks ("LANs") and wide area networks ("WANs"), such as the Internet.

[0402] Certain features described above in the context of separate implementations can also be implemented in combination in a single implementation. Conversely, features described in the context of a single implementation can be implemented in multiple implementations separately or in any subcombination.

[0403] The order in which the above operations are performed may be changed. In some cases, multitasking and parallel processing may be advantageous. The separation of system components in the above embodiments should not be understood as requiring such separation.

[0404] A number of embodiments of the present invention have been described. However, it will be appreciated that various modifications can be made without departing from the spirit and scope of the present disclosure. Accordingly, other embodiments are within the scope of the following claims.

[0405] For example, in one example, a method for mapping tissue oxygen extraction fraction (OEF) of an object includes the acts of receiving a magnetic resonance imaging (MRI) signal obtained by a magnetic resonance scanner, wherein the MRI signal includes multi-echo amplitude and phase data; performing physical modeling of the amplitude and phase data, wherein a QSM-based method and a combined qBOLD (QSM+qBOLD) model are used to account for the effects of OEF on amplitude and phase signals from the same underlying gradient echo data, using automatic cluster analysis of signal evolution over echo time, wherein a time-evolved cluster analysis method is used to overcome noise sensitivity of QSM+qBOLD in obtaining model parameters; and presenting one or more images of the object generated based on the model parameters including the OEF on a display device.

[0406] In some examples of the above methods, clustering is performed using the x-means method.

[0407] In some examples of the above methods, optimization is used to obtain model parameters.

[0408] In some examples of the above methods, deep learning is used to solve and obtain model parameters.

[0409] In some examples of the aforementioned methods, an artificial neural network for deep learning is trained on simulated complex multi-echo MRI data.

[0410] In some examples of the aforementioned methods, the object includes one or more of a cortex, a white matter region, a deep gray matter region, an ischemic brain region, or a diseased tissue lesion.

[0411] In another example, a system for mapping tissue oxygen extraction fraction (OEF) of an object includes a processor; a graphics output module communicatively coupled to the processor; an input module communicatively coupled to the processor; and a non-transitory computer storage medium encoded with a computer program comprising instructions that, when executed by the processor, cause the processor to perform operations comprising: receiving magnetic resonance imaging (MRI) signals acquired by a magnetic resonance scanner, wherein the MRI signals comprise multi-echo amplitude and phase data; physically modeling the amplitude and phase data, wherein a combined QSM-based approach and qBOLD (QSM+qBOLD) model are used to account for the effects of the OEF on amplitude and phase signals from the same underlying gradient echo data; using automatic cluster analysis of signal evolution over echo time, wherein a time-evolved cluster analysis approach is used to overcome noise sensitivity of the QSM+qBOLD in deriving model parameters; and presenting one or more images of the object generated based on the model parameters including the OEF on a display device.< / c> < / c>

Claims

1. A method for generating one or more images of a tissue, the method comprising: capturing time-resolved image data of contrast agent transmission in tissue; Solving the inverse problem involves fitting the time-resolved image data to a transport equation that expresses the relationship between the concentration at a given location and the transport parameter at that location, namely: where c(r,t) is the probability density function of concentration, is the derivative with respect to time, is the derivative of three-dimensional space, u(r)=[u x ,u y ,u z ] is the velocity vector, D(r) is the diffusion coefficient matrix and belongs to the transport parameter, Solving the inverse problem includes: estimating the distribution of the transport parameters of the tissue based on determining prior knowledge about the distribution of the transport parameters, wherein the prior knowledge includes velocity distribution and is sparse; One or more images of the tissue are generated based on the estimated distribution of the transport parameter of the tissue. 2 . The method of claim 1 , wherein solving the inverse problem comprises determining a contrast agent concentration time course by time-resolved imaging of the contrast agent in the tissue. The method of claim 1 , wherein sparse is a cost function characterized by the L1 norm of the gradient of the velocity spatial distribution.

4. The method of claim 1, wherein solving the inverse problem comprises using a quasi-Newton method, an alternating direction method of multipliers, or an artificial neural network. 5 . The method according to claim 1 , wherein the image data is image data of an imaging modality comprising at least one of MRI, CT, PET, SPECT and / or arterial spin labeling data in MRI.

6. A method for generating one or more images of a tissue, the method comprising: capturing time-resolved image data of contrast agent transmission in tissue; Solving the inverse problem fits the image data to the transport equation expressed as the relationship between the concentration at a given location and the transport parameter at that location, namely: where c(r,t) is the probability density function of concentration, is the derivative with respect to time, is the derivative of three-dimensional space, u(r)=[u x ,u y ,u z ] is the velocity vector, D(r) is the diffusion coefficient matrix and belongs to the transport parameter, Deep learning methods are used to solve inverse problems, including determining the distribution of transport parameters of tissues, where an artificial neural network is trained by simulating data based on transport equations; Based on the estimated distribution of the transport parameter, one or more images of the tissue are generated.

7. The method of claim 6, wherein solving the inverse problem comprises determining a contrast agent concentration time course by time-resolved imaging of the contrast agent in the tissue.

8. The method according to claim 6, wherein the image data is image data of at least one of the imaging modalities MRI, CT, PET, SPECT and / or arterial spin labeling data in MRI.

9. A computer-readable storage medium storing a computer program, wherein the computer program enables a processor to execute the method according to any one of claims 1 to 8.

10. A system for generating one or more images of a tissue, the system comprising: processor; a graphics output module communicatively coupled to the processor; an input module communicatively coupled to the processor; A non-transitory computer storage medium encoded with a computer program, the program comprising instructions When executed by a processor, causes the processor to perform operations comprising: capturing time-resolved image data of contrast agent transmission in tissue; Solving the inverse problem involves fitting the time-resolved image data to a transport equation that expresses the relationship between the concentration at a given location and the transport parameter at that location, namely: where c(r,t) is the probability density function of concentration, is the derivative with respect to time, is the derivative of three-dimensional space, u(r)=[u x ,u y ,u z ] is the velocity vector, D(r) is the diffusion coefficient matrix and belongs to the transport parameter, Solving the inverse problem includes: estimating the distribution of the transport parameter of the tissue based on determining prior knowledge about the distribution of the transport parameter of the tissue, wherein the prior knowledge includes velocity distribution and is sparse; One or more images of the tissue are generated based on the estimated distribution of the transport parameter of the tissue, and the one or more images of the tissue generated based on the determined distribution of the transport parameter are presented on a display device.

11. The system of claim 10, wherein solving the inverse problem comprises determining a contrast agent concentration time course by time-resolved imaging of the contrast agent in the tissue.

12. The system of claim 10, wherein sparse is a cost function characterized by the L1 norm of the gradient of the spatial distribution of velocities.

13. The system of claim 10, wherein solving the inverse problem comprises using a quasi-Newton method, an alternating direction method of multipliers, or an artificial neural network.

14. The system of claim 10, wherein the image data is image data of an imaging modality comprising at least one of MRI, CT, PET, SPECT and / or arterial spin labeling data in MRI.

15. A system for generating one or more images of a tissue, the system comprising: processor; a graphics output module communicatively coupled to the processor; an input module communicatively coupled to the processor; A non-transitory computer storage medium encoded with a computer program, the program comprising instructions When executed by a processor, causes the processor to perform operations comprising: capturing time-resolved image data of contrast agent transmission in tissue; Solve the inverse problem, which involves fitting the time-resolved image data to a transport equation that expresses the relationship between the concentration at a given location and the transport parameter at that location, namely: where c(r,t) is the probability density function of concentration, is the derivative with respect to time, is the derivative of three-dimensional space, u(r)=[u x ,u y ,u z ] is the velocity vector, D(r) is the diffusion coefficient matrix and belongs to the transport parameter Using deep learning methods to solve inverse problems, including determining the distribution of transport parameters of tissues, where artificial neural networks are trained on simulated data based on transport equations; One or more images of the tissue are generated based on the estimated distribution of transport parameters of the tissue, and the one or more images of the tissue generated based on the determined distribution of transport parameters are presented on a display device.

16. The system of claim 15, wherein solving the inverse problem comprises determining a contrast agent concentration time course by time-resolved imaging of the contrast agent in the tissue.

17. The system of claim 15, wherein the image data is image data of an imaging modality comprising at least one of MRI, CT, PET, SPECT and / or arterial spin labeling data in MRI.

Images