High Dynamic Range Processing for Magnetic Resonance Imaging

Purpose To minimize feature loss in T1- and T2-weighted MRI by merging multiple MR images acquired at different TR and TE to generate an image with increased dynamic range. Materials and Methods High Dynamic Range (HDR) processing techniques from the field of photography were applied to a series of acquired MR images. Specifically, a method to parameterize the algorithm for MRI data was developed and tested. T1- and T2-weighted images of a number of contrast agent phantoms and a live mouse were acquired with varying TR and TE parameters. The images were computationally merged to produce HDR-MR images. All acquisitions were performed on a 7.05 T Bruker PharmaScan with a multi-echo spin echo pulse sequence. Results HDR-MRI delineated bright and dark features that were either saturated or indistinguishable from background in standard T1- and T2-weighted MRI. The increased dynamic range preserved intensity gradation over a larger range of T1 and T2 in phantoms and revealed more anatomical features in vivo. Conclusions We have developed and tested a method to apply HDR processing to MR images. The increased dynamic range of HDR-MR images as compared to standard T1- and T2-weighted images minimizes feature loss caused by magnetization recovery or low SNR.

Introduction T 1 -and T 2 -weighted imaging are ubiquitous in clinical and research MRI studies [1]. In these imaging methods, contrast can be diminished due to complete magnetization recovery or low signal-to-noise ratio (SNR) [1][2][3][4][5]. For MR images with a limited range of T 1 and T 2 , optimizing the acquisition parameters T R (repetition time) and T E (echo time) circumvents the issue. However, for MR images with a wide range of T 1 and T 2 originating from lesions, tissues, or contrast media, feature loss cannot be avoided regardless of the choice of T R and T E [6]. This scenario occurs, for example, when imaging the whole-animal biodistribution of novel contrast agents [7] or labeled biomolecules [8,9], and can be thought of as a limitation in the ''dynamic range'' of the imaging method.
The formal definition of dynamic range is the ratio between the maximum and the minimum measurable value above noise in an acquisition system [10]. In MRI, T 1 and T 2 are not limited by the hardware dynamic range per se because they are not directly acquired. However, contrast from a limited range of T 1 and T 2 occupies a disproportionately large portion of the signal dynamic range because magnetization recovers and decays nonlinearly according to the two parameters [1]. In other words, voxels with T 1 and T 2 values outside of a limited range either vanish into the background noise or saturate the intensity scale ( Figure 1) [1]. Therefore, a useful definition of ''dynamic range'' in T 1 -and T 2weighted MRI is the range of T 1 and T 2 values that result in detectable and differentiable signal. In an ideal image with unlimited dynamic range as defined here, T 1 or T 2 is monotonically represented by intensity. In some other fields of research, dynamic range is defined in a similar way [11,12].
Dynamic range limitation is a well-recognized issue in the field of photography [10,[13][14][15]. When a scene has a large range of illumination, the resulting photograph invariably has regions that are either saturated or hidden in darkness regardless of the camera exposure time setting. High Dynamic Range (HDR) photography is a technique invented to overcome this challenge by capturing a scene with multiple exposure times and then computationally merging the Low Dynamic Range (LDR) photographs together [16][17][18][19]. Since the over-or under-exposed regions in any one LDR picture are properly exposed in the other pictures in the set, the merging process produces an HDR image that appears properly exposed throughout. HDR processing is highly automated and has been applied widely to films and photography in the last decade [10,[20][21][22][23][24][25].
We applied HDR processing to merge MR images of different T R and T E combinations in a technique we call High Dynamic Range MRI (HDR-MRI). MR images acquired at different T R and T E are analogous to photographs captured at different exposure times because both produce a set of images with increasing brightness. Similar to HDR photography, HDR-MRI displays more bright and dark features from specimens with a large range of T 1 or T 2 compared to standard MRI ( Figure 1). This is possible because each T R and T E combination produces an image with complementary information that is then merged to minimize feature loss.
HDR-MRI follows a long history of post-processing algorithms that have been proposed since the 1980s to transform or combine multiple MR images. These algorithms include principle component analysis [26], eigen-images [27], weighted averaging [28], image synthesis [6,29], and other linear filters [30][31][32][33]. The processed images aid in diagnostic interpretation by increasing contrast-to-noise ratio (CNR), suppressing interfering features, or allowing dynamic viewing of images synthesized with arbitrary T R and T E . HDR-MRI is most similar to image synthesis in that the outputs are readily interpretable by experienced readers, as the hyper-and hypo-intensity relationships of standard acquisitions are preserved. Images produced by other more sophisticated techniques are often more enhanced, but require additional training to read [31]. Compared to image synthesis, HDR-MRI exaggerates or diminishes contrast, is more resistant to uncertainties in T 1 and T 2 fitting, and in theory, displays more information simultaneously when shown on an HDR monitor to potentially reduce reading time [34].
To apply HDR processing to MRI, we derived a mathematical expression to transform the MR parameters T R and T E into the HDR parameter exposure value (EV). The performance of HDR-MRI on solution phantoms doped with contrast agents and a live mouse was investigated. Compared to standard MRI, HDR-MRI differentiated signals from a larger range of contrast agent concentrations in phantoms and revealed more anatomical features in vivo.

Ethics Statement
Animal studies were conducted in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals. The protocol was approved by the Northwestern Institutional Animal Care and Use Committee (Permit Number: 2010-2027).

Theory
A limitation in optical image capture technology is the low dynamic range of films and sensors relative to the human eye Figure 1. Simulation of intensity scale of T 1 -and T 2 -weighted imaging compared to HDR-MRI. T 1 and T 2 values cover the physiological range. For comparative purposes, dynamic range was defined to span from 10% to 90% of each intensity scale. Regions of an image with T 1 and T 2 outside this range would appear featureless as either completely white or black. HDR-MRI improved the dynamic range compared to T 1 -and T 2weighted imaging. Furthermore, in conventional imaging, dynamic range was sacrificed to visualize short T 1 and T 2 features whereas HDR-MRI did not suffer from this limitation. doi:10.1371/journal.pone.0077883.g001 [14,15]. As a result, some features are obscured in darkness while others appear saturated. HDR processing overcomes this limitation by merging multiple LDR photos with varying exposure times to capture the full range of features observed naturally by the human eye. There are a number of algorithms for HDR processing [16,17,[35][36][37]. The commercial package utilized here (HDR Pro built into Adobe Photoshop CS5) implements the algorithm developed by Debevec and Malik [16].
The Debevec algorithm estimates the physical illumination (I) of an object based on the principles of digital image capture systems. An image is captured through light sensor exposures (E) that scale with object illumination and exposure time (t) (Equation 1).
The light sensor converts the exposure (E) to a finite range of pixel brightness through a digital conversion function f, otherwise known as the characteristic curve ( Figure 2A) (Equation 2).
Z ij denotes the brightness of pixel i at exposure time t j . With the reasonable assumption of f as a monotonically increasing function, an inverse function, f 21 , exists (Equation 3).
Taking the natural logarithm defines the function g (Equation 4) In Equation 4, g and I i are the unknowns and can be estimated by minimizing a quadratic error function (Equation 5).
N is the number of pixels per picture and P is the number of pictures to be merged. When N and P are sufficiently large, the error function can be seen as an over-determined system of equations that is solved readily using available algorithms [16]. The first term of the error function ensures fidelity of Equation 4. The second term ensures smoothness of g. The strength of the second term is controlled by the choice of l. The second derivative g0(Z) is defined in discrete form (Equation 6).
The weighing function, w, emphasizes pixels closer to the center of the light sensor dynamic range so that pixels without signal or that are saturated weigh less in the determination of g (Equations 7 and 8).
Once g has been computed, the estimation of physical illumination I i can be further improved by a weighted average (Equation 9).
In this equation, the determination of I i is weighted toward pixels with proper exposure because the weighing function approaches zero when the pixel is under-or over-exposed (i.e. Z ij approaches Z min or Z max ). The re-constructed HDR image presents the calculated physical illumination (I) of the scene and consequently spans over a wider dynamic range than any of the source images ( Figure 2B). To utilize HDR processing for MR images, it is necessary to transform the MRI parameters T R and T E into the photography parameter exposure time (t). T R and T E in MRI are analogous to exposure time in photography because they are acquisition parameters that affect voxel intensity just as exposure time affects pixel intensity. For a Multi-Echo Spin Echo (MESE) pulse sequence, the voxel signal intensity is governed by Equation 10 [1,38].
M o is the equilibrium magnetization, T R is the repetition time, T E is the echo time, and T 1 and T 2 are the longitudinal and transverse relaxation time, respectively. In this equation, T R and T E have monotonic relationships with the MR signal (S), analogous to the relationship between exposure time (t) and film exposure (E) in photography (Equation 1) [39]. By setting equations 1 and 10 equal to each other, Equation 11 is obtained: x is a constant consisting of M o and I. The calculated t can be used to solve Equation 5, with the index for t j running from 1 to 4 to represent images acquired at 4 different T R s or T E s used in the experiments. In practice, t is further transformed into exposure values (EV) by Equation 12 for input into the HDR software [16,39].
F refers to the relative aperture. Since EV is a relative value in HDR processing, F and, to a first order approximation of constant M o , x can be arbitrarily chosen. We set F to 4 because it is a commonly used setting in photography and x to 1 for simplicity. It is important to note that the T 1 , T 2 , and M o chosen to calculate t affect processing. This complication should not significantly impact the qualitative interpretation of HDR-MR images because it does not alter the hyper-and hypo-intensity relationships in the image. The effect of non-uniform T 1 , T 2 , and M o on HDR processing is examined more closely in the discussion section.

High Dynamic Range Processing
For T 1 -weighted image series with varying T R , the RecoSeries macro in ParaVision software (version 5.1, Bruker BioSpin, Billerica, MA) was used to reconstruct all 4 images on the same intensity scale. This pre-processing was not required for images in (5) T E series because they were acquired in a single scan. All images were converted from signed 16-bit integer encoding to unsigned 16-bit integer encoding in ImageJ (NIH, Bethesda, MD).
T 1 and T 2 of solution phantoms were fitted using the built-in Image Sequence Analysis tool in ParaVision or a custom program written in Matlab (ver. R2010b, MathWorks, Natick, MA). HDR processing was performed using HDR Pro built into Adobe Photoshop CS5. The only input parameters for processing were EV and the choice of white point. The EV of each source image was calculated using Equations 11 and 12 based on T 1 and T 2 values specified in Table S1. The T 1 and T 2 used represent intermediate values within the estimated range for each figure.
The white point of the intensity scale was always set to the brightest pixel of the image. No additional adjustments, such as tone mapping or detail enhancement, were performed. The output was saved in 32-bit TIFF format. For the phantom results, each HDR image was normalized using ImageJ such that the intensity of the water signal is equivalent to the water signal of the brightest input image. All images were displayed on the same intensity scale within each figure except for the in vivo images. For pseudocoloring, Fire LUT in ImageJ was applied.
Simulation of Dynamic Range of T 1 -weighted, T 2weighted, and High Dynamic Range MRI The T 1 -and T 2 -weighted intensity scales were created by simulation in Matlab using equation 10. T E = 0 and T R = ' was used for T 1 -and T 2 -weighted simulation, respectively. Each scale is displayed between pixel values 0 and 65535. The HDR intensity scale was generated by merging the LDR scales using HDR processing as described previously.

Quantification
Image intensity quantification was done using ImageJ with the Measure function on ROIs. Normalized contrast was defined by Equation 13.
I x,y is the normalized contrast between samples x and y. It is calculated by subtracting the intensity of sample y (I y ) from that of sample x (I x ), followed by normalization with the intensity of water (I water ). To calculate SNR, water signal was divided by noise either in the void region or within the water signal as indicated in the figures. Noise in the void is calculated as s~s m =0:65, where s is the physical noise and s m is the measured standard deviation, to account for the Rician distribution of noise in MR magnitude images [42]. Noise within water is approximated as s = s m and requires no correction because the noise distribution is largely Gaussian at high SNR. Error bars represent s that were propagated from the source images according to the arithmetic operations performed. For example, the s of I x -I y is calculated as  In vivo Imaging A female Balb/C athymic nude mouse was acquired from Harlan (Indianapolis, IN) and housed under pathogen-free conditions. The mouse was maintained under anesthesia (1-3% isoflurane) and respiration monitor. Tubing containing heated water was positioned under the mouse to maintain constant body temperature. Imaging was performed on an 89 mm bore PharmaScan 7.05 T MR imager fitted with shielded gradient coils (Bruker BioSpin, Billerica, MA, USA) using a RF RES 300 1H 089/038 quadrature transceiver volume coil (Bruker BioSpin, Billerica, MA, USA). Respiration-gated MESE scans with fat suppression were used: T R = 5632 ms, T E = 10.635…212.7 ms in 10.635 ms increments, NEX = 1, Bandwidth = 391 Hz/pixel, FOV = 35635 mm 2 , slice thickness = 1 mm, and matrix size = 1286128. HDR processing was performed as described previously. The T 2 map was generated by fitting every voxel in the T E series to A?exp(2T E /T 2 )+C using Matlab. Masking was done by manually thresholding the T E 21 LDR image at 1500 such that the voxels below and above the threshold were set to 0 and 1, respectively. The intensity scale of each image was adjusted individually for visibility.

Results
To assess the applicability of HDR processing to MR images, we prepared solution phantoms consisting of T 1 (Gd(III)HPN 3 -DO3A) [41] and T 2 (NF) [40] contrast agents (Figure 3). This assortment of phantoms simulated scenarios when many different tissues or a large range of contrast agent concentrations are imaged in the same view. The concentrations, T 1 , and T 2 values of each phantom are listed in Figure 3. For samples with high concentrations of contrast agents (samples 1, 2, 7, and 8), the T 2 was short compared to the shortest T E available for the pulse sequence used. As a result, these samples were invisible in most of the MR images acquired [1,38].

HDR Processing of MR Images
HDR-MR images were generated and compared to conventional images. For T 1 -weighted HDR-MRI, images were taken at constant T E and varying T R . The result showed that a T R of 200 ms gave the highest contrast among samples 3-5, while longer T R s better distinguished the lower concentration samples from background ( Figure 4A-D). By merging the LDR images, HDR processing combined the advantage of pronounced contrast at short T R and improved delineation of weak signals at longer T R ( Figure 4E, F).
Similarly, for T 2 -weighted HDR-MRI using images with varying T E , each image spanned a different dynamic range and provided the optimal contrast between different pairs of samples ( Figure 5A-D). For example, the T E used in Figure 5A and 5D are best at differentiating sample 9 from 10 and sample 10 from 11, respectively. The HDR image had an extended dynamic range and was capable of differentiating samples 9, 10, and 11, on the same image ( Figure 5E, F). However, this wider dynamic range sacrificed quantitative contrast between samples since the absolute range of pixel intensities is fixed by the display device and file encoding capability.

EV Interval Selection
The EV interval of the source LDR images has a significant impact on the quality of the final HDR image. In the extreme, a selection of LDR images with identical EV produces a HDR image without any enhancement. In photography, the optimal EV interval is typically between 1 and 2 to ensure overlap of the characteristic curves at different exposure times (Equation 2, Figure 2A) and a wide dynamic range [16,43]. To investigate whether MR images have similar optimal EV gaps, a series of images with constant T R and varying T E were acquired; selections of 4 single LDR images with EV intervals of 0, 0.5, 1, and 1.4 were utilized for HDR processing. Each selection had a different T E range that centered around T E = 270 ms, with the largest spanning from 11 ms to 521 ms (Table S1). In addition, simple averaging of 4 LDR images with T E 11 ms and 521 ms were performed to compare with HDR processing.
As previously discussed, HDR improves the dynamic range of an image to include dark and bright features by the sacrifice of relative contrast, and this is further illustrated in Figure 6A. The curves shown can be interpreted as rough representations of characteristic curves (pixel intensity vs. film exposure) because iron concentration is correlated to T 2 . The two LDR curves had similar slopes, but were shifted horizontally relative to each other, indicating similar dynamic ranges that cover different regions. With HDR processing using EV interval 0, the slope was similar to the LDR images, consistent with expectation; as the EV interval increased, the curve flattened to span through a wider range of concentrations, indicating an increase in dynamic range with a concomitant decrease in relative contrast. The same conclusion regarding EV intervals could be drawn visually from the acquired images ( Figure 6B). Wider EV intervals generated images with greater dynamic range that better delineated low signal features from background while retaining differential contrast in the high signal features. It is evident that an EV gap of 1-2 yielded the best combination of range and contrast, similar to the guidelines used in photography. Functionally, the increased dynamic range of HDR-MRI as compared to conventional imaging allowed for the capture of contrast enhancement from a larger range of contrast agent concentrations.

In Vivo Imaging and HDR
Based on the in vitro findings, a live mouse was scanned using a carefully selected T E series that was merged into a T 2 -weighted HDR image (Figure 7). The HDR image contains features that one or more of the individual LDR images lack. For example, the red highlight shows the low signal features present when T E = 21 ms but not at longer T E s, and the yellow highlight outlines the contrast between features that were only evident when T E $74 ms. For comparison, a T 2 map was generated using the same dataset. The T 2 map generally captures the same features, but required manual masking and is noisier in the regions with low SNR due to insufficient signal for fitting at long T E . A similar result was obtained when T 2 was mapped using images acquired at T E = 21, 43, 64, and 85 ms to improve accuracy ( Figure S1). In contrast, HDR processing is able to capture low signal features accurately even in the case when they are only visible in a single source image.

Discussion
Our results demonstrate a method of applying HDR processing to MRI. HDR-MRI exchanges relative contrast for dynamic range to minimize feature loss caused by complete magnetization recovery and low SNR. The improved dynamic range captured contrast agent enhancements over a larger concentration range and preserved more anatomical details in vivo. These properties make HDR-MRI suitable for imaging a large number of tissue types and contrast agent concentrations simultaneously, such as when performing small animal biodistribution studies by MRI. The Effect of T 1

, T 2 , and M o on HDR Processing
When EV is calculated using Equations 11 and 12, the analogy between photography and MR imaging is complicated by the fact that each voxel has a different T 1 , T 2 , and M o . One way to understand the effect of this complication is to think of the inputted EV as under-and over-estimating the EV of different voxels. Consequently, their physical illuminations become miscalculated according to Equations 4 and 9. In the case when T 1,voxel ,T 1,input , T 2,voxel .T 2,input , or M o,voxel .M o,median , the HDR algorithm overestimates the illumination of the voxel ( Figure S2). Incidentally, shorter T 1 , longer T 2 , and larger M o all lead to higher signals on an MR image to coincide with the overestimated illumination (Equation 10). Similarly, when T 1,voxel .T 1,input , T 2,voxel ,T 2,input , or M o,voxel ,M o,median , the HDR algorithm underestimates the voxel illumination. In other words, the effect of voxel-dependent T 1 , T 2 , and M o can be understood as a brightening of high signal features and a darkening of low signal features. The result is the exaggeration and diminishment of quantitative contrast without alteration to the hyper-and hypointensity relationships between the different features.
To validate the theoretical expectations and quantify the contrast modifications, HDR processing was performed on the same set of source images using EVs calculated from different T 1 and T 2 combinations (Figures S3 and S4). Only the choice of T 1 affected the processing of T 1 -weighted HDR-MRI because T E %T 2 and was fixed across the T R series (Equation 11 and Figure S3). Similarly, only the choice of T 2 was important in T 2weighted HDR-MRI ( Figure S4). The result showed that in both cases, the quantitative contrast varied, but the relative feature brightness remained the same regardless of the EV used, consistent with the theoretical predications.
A further assessment showed that the degree of intensity variation depended on the choice of T 1 or T 2 in calculating the input EV. When extreme values of T 1 or T 2 were used, the contrast modification can be significant. For example, when T 2 = 43 ms was used to calculate EV in the T 2 -weighted HDR-MR image, the water illumination (T 2 = 400 ms) was overestimat- Images A-D were computationally merged by the HDR algorithm to generate the (E) HDR image. HDR-MRI accurately captured the intensity difference between samples 9, 10, 11, and 12, which was not achieved in any of the source images. This was shown quantitatively by contrasts normalized against water (F). Image A was poor at differentiating samples 10-12, B and C were poor at differentiating samples 11 and 12, and D was poor at differentiating samples 9 and 10. Error bars represent standard deviation. doi:10.1371/journal.pone.0077883.g005 ed to such a degree that the darker features became invisible in comparison ( Figure S4). However, when intermediate T 1 or T 2 was chosen, the variations were mitigated. For example, in the T 2weighted HDR-MR image (T 2 range = 43-400 ms), EVs calculated using a T 2 anywhere between 170-370 ms resulted in intensity variations of no more than 17% relative to water ( Figure  S4). This degree of variation is commonly experienced in conventional T 2 -weighted imaging as a result of operator differences in choosing T E ( Figure S4). A similar conclusion can be drawn for T 1 -weighted HDR-MRI ( Figure S3).
In HDR-MRI, voxel-dependent T 1 , T 2 , and M o results in the exaggeration or diminishment of quantitative contrast while conserving hyper-and hypo-intensity relationships in accordance with the source images. When intermediate T 1 or T 2 is used to calculate the input EV, the contrast modifications introduced are no greater than the variations introduced by operator differences in conventional imaging due to T R and T E choices. Similar outputs can be obtained using a range of T 1 or T 2 in the calculation of EV, making the method somewhat robust. Based on these observations, we expect the effect of voxel-dependent EV to only minimally interfere with qualitative interpretation, especially on viewing software with dynamic brightness and contrast adjustments. If quantitative information is desired, a more sophisticated approach could be taken in the future to modify the HDR algorithm based on the MR signal equation.

The Effect of T R and T E Choice on HDR Processing
In HDR photography, the choice of exposure times affects the resulting HDR image. Similarly, the choice of T R s or T E s is an important consideration in HDR-MRI. The choice of T R s in T 1weighted HDR-MRI and T E s in T 2 -weighted HDR-MRI determines the EVs and, in turn, the dynamic range of the output ( Figure 6). For an MR image set that covers a fixed dynamic range, variations in T R or T E change the calculated EV for each image (Equation 11); as a result, the feature intensities in the HDR-MR image would vary in a manner similar to when the choice of T 1 or T 2 is varied in the calculation of EV. As discussed previously, this degree of variation is not expected to hinder qualitative interpretation. To obtain unique HDR-MR images for quantitative purposes, the underlying HDR algorithm would need to be modified.

Recommendations
Based on the findings presented, a recommended protocol for HDR-MRI is summarized in a flowchart (Figure 8). Choosing an intermediate T 1 or T 2 value for EV calculation and a T R /T E Figure 6. EV interval analysis. (A) Signal intensities normalized against water at various NF concentrations. Sigmoidal curves illustrate trends in the data and can be thought of as characteristic curves. As the EV interval increased, the characteristic curve slope flattened to span over a wider dynamic range. The characteristic curves of the two averaged LDR images (T E 11 and T E 521) had the steepest slopes and were analogous to the EV interval 0 condition. (B) Comparison between averaged LDR images and HDR images. The HDR image with EV interval 1.4 showed contrast between every sample from 8 to 12, indicating its large dynamic range. N = 4 was used for the averaged images. Error bars represent standard deviation. doi:10.1371/journal.pone.0077883.g006 series with EV intervals between 1-2 have been shown to work well with the commercial package utilized here. Although no precise method is prescribed for choosing T 1 or T 2 in calculating EV and some operator differences cannot be avoided, similar results are obtained within a range. Consequently, the choice should not significantly impact qualitative interpretation based on hyper-and hypo-intensity relationships. Once the EV interval has been determined, the number of images to acquire depends on the estimated range of T 1 or T 2 in the image. For best results, the shortest and the longest T R or T E in the sequence should approximate the shortest and the longest estimated T 1 or T 2 , respectively. In our experience, four images comfortably cover many scenarios.

Cost Benefit Analysis
HDR-MRI preserves features with extreme T 1 or T 2 values at the expense of relative contrast and imaging time compared to standard T 1 -and T 2 -weighted imaging. Relative contrast in HDR-MRI is decreased due to flattening of the characteristic curve such that an expanded portion of the intensity scale is used to represent T 1 and T 2 contrast. With the advent of HDR displays, this limitation of HDR-MRI may be alleviated. HDR displays exhibit a contrast ratio of 50000:1 and a maximum intensity of 8500 cd/ m 2 , compared to a typical desktop display with a contrast ratio of 300:1 and a maximum intensity of 300 cd/m 2 [34]. Trading contrast for feature visibility using HDR-MRI would be advantageous on an HDR display due to the ample availability of contrast.
Compared to single image acquisitions, T 1 -weighted HDR-MRI lengthens imaging time by 7 to 8 fold while T 2 -weighted HDR-MRI has no additional cost. The reason is that T 1 -weighted HDR-MRI uses long T R s while T 2 -weighted HDR-MRI uses a multi-echo sequence that acquires in a single T R . Therefore, the time cost of HDR-MRI is the same as in T 1 and T 2 mapping.
For the same amount of imaging time or with the same set of source images, many algorithms can be applied to generate synthetic images of higher quality. A simple way to improve quality is by averaging multiple images to obtain increased SNR and CNR (contrast-to-noise ratio). A direct comparison between averaging and HDR-MRI in terms of SNR is difficult because of the complex SNR behavior produced by the HDR algorithm ( Figure S5). However, HDR-MRI likely depicts features at the extremes of the T 1 /T 2 scale with higher contrast due to the use of different T R and T E combinations.  Compared to T 1 and T 2 mapping, HDR-MRI is more resistant to noise but does not provide quantitative parameters. Uncertainty in the fitting of voxel T 1 or T 2 does not affect HDR processing because the individual voxel relaxation times are not used as parameters. In addition, Equation 4 guarantees that the low SNR voxels that are problematic for relaxation time fitting would always remain dark in the HDR-MR image.
Transformations based on component analysis are another method of combining MR images [26,27,30,31,33]. These algorithms enhance the CNR of the features of interest and can additionally suppress the intensities of other interfering features. HDR-MRI does not achieve the CNR improvement of these methods. Instead, HDR-MRI produces outputs similar to T 1 -and T 2 -weighted images that are familiar to experienced readers and is suitable when the appearance of the feature of interest is not known a priori. In this regard, HDR-MRI is closely related to image synthesis [6,29]. Image synthesis is a method to dynamically generate synthetic T 1 or T 2 -weighted images at arbitrary T R and T E using a map of relaxation times. In theory, HDR-MRI compresses the information of all possible synthetic images into one image at the expense of relative contrast. On an HDR display, the higher information content of HDR-MR images may reduce reading time without sacrificing accuracy.
HDR processing can be easily built into viewing software to provide an alternative way of visualizing multi-image data alongside other algorithms. If the images have already been acquired for a different purpose such as relaxation time mapping, the cost of HDR-MRI is purely computational. In addition to the use of HDR processing with the MESE pulse sequence used here, it is possible to apply the same principle to other contrast parameters or pulse sequences such as flip angle, inversion time, phase-based flow imaging, or diffusion-weighted imaging.

Conclusions
HDR processing is a powerful technique in photography used to simultaneously capture dark and bright features. In MRI, having a large range of T 1 or T 2 results in diminished feature visibility, analogous to the challenge faced in photography. We have applied the fundamental principles of HDR photography to MR imaging and provided both phantom and in vivo examples showing the comparisons between standard T 1 -or T 2 -weighted images and HDR-MR images. HDR-MRI provides an alternative to standard imaging by merging multiple non-optimal images highlighting different features into a single image that displays all features simultaneously. This technique may increase contrast agent visibility in whole-animal biodistribution studies and reduce diagnostic workload by maximizing the information content of an image. Figure S1 Additional comparison between T 2 mapping and T 2 -weighted HDR-MRI. T 2 mapping and HDR-MRI based on the same source images display similar features, with the T 2 map appearing noisier in the low signal regions. A similar result was obtained when a more accurate T 2 map was produced using shorter T E s. (TIF) Figure S2 Qualitative picture for the effect of T 1 , T 2 , and M o on HDR Processing. In the presented scenario, there are only two voxels. The inputted exposure times t j , j = 1-5 (or more strictly, EVs) are accurate for voxel A. Voxel B is physically brighter. (A-C) The g function before HDR processing when illumination (I) is arbitrarily assumed to be 1 by the algorithm at both voxel A and B. (A) In photography, voxels A and B share the same set of five exposure times. Therefore, the inputted t j are accurate for both A and B. Voxel B is physically brighter, resulting in its larger voxel intensities. I B can be accurately estimated by the algorithm. (B) If the input overestimates the true exposure times of voxel B, as when T 1,voxel .T 1,input , T 2,voxel ,T 2,input , or M o,voxel ,M o,median (Equation 11), I B is underestimated (Equations 4 and 9). (C) Conversely, if the input underestimates the true exposure times of voxel B, as when T 1,voxel ,T 1,input , T 2,voxel .T 2,input , or M o,voxel .M o,median , I B is overestimated. (D) In all three cases, HDR processing calculates I B to obtain the identical final monotonically increasing g function (Equations 4 and 5). (TIF) Figure S3 Effects of the choice of T 1 and T 2 in calculating EV on T 1 -weighted HDR processing. The values tested are the T 1 and T 2 of the samples imaged. The choice of T 2 has almost no effect while the choice of T 1 has moderate effect. In general, the quantitative contrast changes, but the relative order of the feature intensities is preserved regardless of the T 1 chosen for EV calculation. The variations caused by the different choices of T 1 are similar to variations seen in conventional T 1 -weighted imaging due to different T R settings. Water signal has been normalized to 100 across all images. (TIF) Figure S4 Effects of the choice of T 1 and T 2 in calculating EV on T 2 -weighted HDR processing. The values tested are the T 1 and T 2 of the samples imaged. The choice of T 1 has almost no effect while the choice of T 2 has moderate effect. In general, the conclusion is the same as for T 1 -weighted HDR processing. The quantitative contrast changes, but the relative order of the feature intensities is preserved regardless of the T 2 chosen for EV calculation. The variations caused by the different choices of T 2 are similar to variations seen in conventional T 2 -weighted imaging due to different T E settings. Water signal has been normalized to 100 across all images. (TIF) Figure S5 SNR comparison between image averaging and HDR-MRI. Data from Figure 6 was used for the analysis. Dashed line (----) and dotted line (????) represent the SNR of the averaged LDR images with T E 11 ms and T E 521 ms, respectively. (A) When calculated against void noise, HDR-MRI SNR improved with increasing EV interval and outperformed averaging. (B) When calculated against noise in the water signal, HDR-MRI SNR remained constant across EV intervals and underperformed averaging. SNR is inhomogeneous in HDR-MRI because the characteristic curve is nonlinear and each voxel is processed by the algorithm independently. N = 4 was used for the averaged images. Error bars represent standard deviation. (TIF)