Introduction
Motivation for This Work
Magnetic resonance imaging (MRI) can not only provide static anatomical information but also capture dynamic processes such as the uptake of an intravenously injected contrast agent with dynamic contrastenhanced MRI (DCEMRI). Rather than using a fixed temporal resolution for such dynamic scans, emerging methods based on goldenangle sampling such as XDGRASP (1) or CIRCUS (2) allow retrospective temporal resampling at the time of image reconstruction. This framework allows retrospective control of the temporal “footprint” by varying the number of kspace rays or spokes per undersampled image, without necessarily requiring the use of view sharing, which can introduce temporal blurring of fast dynamics (2, 3). Doing so allows one to flexibly optimize the tradeoff between better spatial quality (ie, low undersampling factor) and better temporal quality (ie, high undersampling factor).
With additional reconstruction flexibility, however, comes the need to optimize image reconstruction of that specific data set and its associated spatiotemporal features. This is particularly true when one considers the large and intertwined joint optimization space associated with these reconstructions, in which the undersampling factor (ie, temporal footprint), the form of regularization enforced, and the degree of regularization must all be simultaneously examined.
Recently, it has become common to use objective image quality metrics (IQMs) (4, 5) as an outcome measure or within the cost function of compressed sensing (CS) or with deep learning–based image reconstruction (6–11). A recent study by our group examined a broad range of IQMs to determine the varying ability of these metrics to correlate with expert radiologists’ rating of diagnostic image quality (12). Although such an approach is sufficient in the reconstruction of a single image data set with fixed undersampling factor, in a data set with goldenangle sampling that allows multiple potential reconstructions, most IQMs will be maximal for the lowest undersampling factor possible, that is, with no ability to resolve temporal dynamics. Therefore, additional research is required to make use of IQMs in such a scenario, and the approach using expert radiologist opinion may be impractical.
Furthermore, the kspace sampling and image reconstruction factors that maximize the quality of diagnostic images from a dynamic data set may not be the same ones that minimize the parametric error when a model is subsequently fit to the data set in a quantitative dynamic application, such as the Tofts pharmacokinetic model in the case of quantitative DCEMRI (13, 14). For example, while a radiologist may visually prefer a diagnostic image with less noise, a reconstruction with increased image noise but improved temporal frame rate may potentially result in more accurate parameter fits. It is also not feasible to present radiologists with thousands of images resulting from the numerous permutations of kspace sampling and reconstruction choices. The ability to automate the joint tuning of goldenanglesampled dynamic data reconstruction, for example, through the use of an IQMbased approach, would, therefore, potentially improve the accuracy of recovered parameters in quantitative dynamic applications.
Goals of This Study
We determined the relationship between 2 wellestablished IQMs—root mean square error (RMSE) and structural similarity index (SSIM) (5)—and the subsequent fidelity of recovered dynamic model parameters when applied to DCEMRI with CS reconstruction of goldenanglesampled data.
We began by assessing the properties of RMSE and SSIM in static CSMRI reconstructions using retrospectively undersampled data from pelvic scans as well as a numerical phantom. This verified that each IQM was sensitive to the CS reconstruction parameters on a patientbypatient basis.
With verification of IQM response in basic CSMRI application, we then evaluated their ability to predict the quality of goldenanglesampled dynamic MRI data reconstructed with CS. To this end, we developed a simplified simulation framework for a dynamically evolving numerical phantom as a surrogate for the assessment of CSMRI in dynamic MRI procedures. The use of a numerical phantom allowed for groundtruth knowledge of the instantaneous temporal evolution. To validate our simulation framework, we first analyzed cylindrical features that evolved according to a monoexponential decay within a digital phantom. Once validated on a simple closedform signal evolution model, we analyzed the more physiologically relevant extended Tofts model of pharmacokinetics, which models the uptake of a contrast media in tissue.
Methodology
Modified CIRCUS kSpace Acquisition
Goldenangle sampling of kspace was simulated using a modified CIRCUS acquisition (15), with a repetition time of 3.7 milliseconds. Although multiple goldenangle sampling techniques are possible, we chose CIRCUS because it is a Cartesian trajectory, making for easy implementation and analysis. CIRCUS subdivides kspace into nested squares on a Cartesian grid, with rudimentary sampling patterns formed by selecting one point from each square according to a goldenangle rotation (2). We refer to these rudimentary sampling patterns as “quanta.” The 2D CIRCUS patterns are extended for use in 3D imaging by extending the sampling mask along the third dimension (ie, the frequency encoded dimension). The shapes of each quanta are controlled by parameters
b and
c, which apply linear and nonlinear shifts to the position of successive samples and may loosely be thought of as contributing to the “shear” and “twist” of the quanta. The value of
b must be positive, and
c should be between 1 and 2; using
b=0 and
c=1 yields quanta that are straight lines, whereas
b≫1 and
c=2 maximizes randomness. For this work, we set
b=40 and
c=1.5, which are consistent with the values recommended by the authors of CIRCUS.
Individual quanta are too highly undersampled to yield useful images, so they were combined into groups of 10 to form “packets.” Quantities of packets ranging from 3 to 21 in steps of 2 were then interleaved in various numbers when forming image reconstructions. The number of packets interleaved then dictated the temporal resolution of the resulting reconstructions, where more packets meant a larger temporal footprint.
The kspace data were sampled by taking the Fourier transform of the image and selecting the coefficients indicated by the CIRCUS packets. In addition, we opted to perform the studies presented herein without the addition of simulated noise to the signal. We acknowledge this to be a simplification of the underlying signal acquisition physics; however, it allows a preliminary investigation of the application of IQMs to the analysis of dynamically evolving data and to determine their fundamental response to such data in the absence of external factors, which is the goal of the present work.
Retrospective Static Pelvic CSMRI
To show the performance of IQMs in static MRI data as a function of reconstruction parameters, 15 static pelvic images were retrospectively undersampled. The acquisition of patient data for this research was approved by the Nova Scotia Health Authority Research Ethics Board and complied with all ethical obligations. Data were acquired on a GE MR750 3 T system using a 32channel abdominal radiofrequency coil array and a 3D spoiled gradient echo sequence (DISCO) (16) with fat separation. Data were acquired with autocalibrated parallel imaging (acceleration factor 2) in the phaseencoding direction, with all omitted lines reconstructed before retrospective undersampling.
Retrospective undersampling of the static pelvic images was performed using the modified CIRCUS acquisition described previously. Combinations of packets were chosen resulting in undersampling factors ranging from 1 to 12 in steps of 0.5. CS reconstructions used the Berkeley Advanced Reconstruction Toolbox (17) to solve the following system:
where
F is the Fourier transform operator,
ρ is the estimated image reconstruction,
k are the measured kspace data, and
λ is the regularization weight given to enforcing sparsity in the transform domain
Ψ at the expense of data consistency. The transform domain
Ψ should be a basis within which the image can be represented by relatively few elements; said another way, the image must be compressible within
Ψ (
18). A limiting factor in the choice of
Ψ is the ability of the basis to represent the image in this way. A regularization weight
λ of 0 represents no enforcement of sparsity in the domain
Ψ, and increasing values of
λ represent a compromise wherein data consistency requirements are relaxed to allow for a solution to be found that is sparse in
Ψ. Choosing
λ=0, therefore, results in no application of CS, and very large values of
λ will result in reconstructions that are of poor quality owing to a lack of consistency with the acquired data. An appropriate choice of
λ is typically relatively small, and the final value of
λ to use in a CS reconstruction has frequently been an empirical choice.
For this reconstruction, we chose a wavelet basis for
Ψ and used values of
λ ranging from 0 to 0.08 in steps of 0.002. The quality of the resulting image reconstructions was assessed by RMSE and SSIM.
Dynamic CSMRI Simulations
Simulations were written in Matlab 2016b (The MathWorks, Natick, MA). A digital phantom of size 200 × 200 × 32 voxels was designed to contain features that evolve over time. The embedded features of the phantom are shown in Figure 1. The plane of the embedded cylindrical features at the very center of the phantom was chosen for assessment in this work, with the remaining features in the phantom intended to provide flexibility for future work. Within the plane of the embedded cylindrical features used, we placed regions of interest (ROIs) inside the inner cylindrical features of a row of cylinders, shown in yellow in Figure 1, for quantitative parameter reconstructions. The ROIs for each embedded cylindrical feature were combined and averaged for analysis. In total, 596 voxels were used in the combined ROIs.
Figure 1.
Internal features of the digital phantom embedded in a homogeneous static background. Within each plane of pingrids (shown in brown), all grids can be made to evolve together according to a userdesigned function. Within each plane of embedded cylinders, the outer cylinders (shown in orange) can be made to evolve together according to a userdesigned function, and the inner cylinders (shown in yellow) can be made to evolve independently along each row. The interior cylinders in one plane of embedded cylindrical features, indicated by a red arrow, were used for this work; a 2D view of this slice is given in Figure 2F.
For this work, we simulated 2 scenarios. To begin, a simple monoexponential decay model was simulated. This model consisted of 2 free parameters—the initial amplitude and the decay rate of the signal amplitude. In total, 200 packets were acquired during the simulated evolution of the phantom, with individual kspace samples acquired at a simulated interval of 3.7 milliseconds. This result corresponded to 2 minutes and 32 seconds of simulated data acquisition. From these simulations, we fit the model to the reconstructed timecourses extracted from each voxel to recover both the initial magnitude and decay rate parameters. The use of a simple model such as this one validated that the simulations were providing rational results before extending it to a more realistic and complex one such as the extended Tofts model. For validation, we expected the following:
A tradeoff between temporal resolution and image quality would be evident.
The temporal footprint (ie, the undersampling factor) used would have a larger effect than the degree of regularization applied.
Following the simple simulations described previously, features in an extended Tofts model were then simulated to evolve according to equation [2]:
The concentrations of contrast agent in the tissue and in the blood plasma are represented in equation [2] by
Ct(t) and
Cp(t), respectively. The terms
ve and
vp represent the volume fraction in tissue of the extravascular extracellular space and blood plasma, respectively. The term
Ktran in equation [2] is a transfer constant, which is indicative of the permeability of the vasculature. From these simulations, we fit equation [2] to the reconstructed timecourses extracted from each voxel to recover
Ktran. To implement equation [2], we used a populationaveraged arterial input function (AIF) to define
Cp(t) at each point in time (19). The Tofts model–defined timecourses evolved for 2 minutes and 30 seconds, with the contrast agent introduced after 30 seconds. This allowed for the collection of 197 packets.
We interleaved differing numbers of packets before CSMRI image reconstruction when making reconstructions of both the monoexponential decay and the Tofts model scenarios. Combinations of 3–21 packets, in steps of 2 packets, were supplied to Berkeley Advanced Reconstruction Toolbox’s pics command, which performs an eSPIRiTbased parallel imaging and CS image reconstruction (20), enforcing a wavelet sparsity constraint with regularization weights of 0–0.05 in steps of 0.005. This corresponded to undersampling factors ranging from 12.5 to 2, where the undersampling factor is defined as the reciprocal of the fraction of the kspace matrix that is sampled, for example, sampling 50% of the kspace coefficients results in an undersampling factor of 2. We chose to present the results of this portion of the study as a function of the number of packets, rather than of the undersampling factor, because it is the number of packets that the user has control over. These packet combinations corresponded to simulated temporal windows of ∼2–13 seconds.
For the simple dynamic simulations consisting of a monoexponential decay, the reference image was taken to be the instantaneous value of the numerical phantom at the average time of all constituent packets; this is referred to as the instantaneous reference image. For the dynamic simulations using the Tofts model of signal evolution in equation [2], 3 additional reference images were generated from the simulated data. In practice, an instantaneous reference image would not be available, so this step allowed us to test the feasibility of using a datagenerated reference image. The following were the 3 datagenerated references:
A precontrast reference image made by acquiring an extra 30 seconds of data before contrast enhancement and combining all packets up until 5 seconds before contrast enhancement.
A peakAIF reference image, made by combining all packets in a 20second window starting at the peak of the AIF used to get
Cp(t) in equation [2].
A latecontrast reference image, made by combining all packets in the last 20 seconds of data acquisition where the contrast dynamics are relatively slow.
In testing these 3 additional datagenerated references, reconstructions were made using 3–21 packets as described previously, but only a single regularization weight of 0.05 was used.
Image Quality Assessment
RMSE and SSIM were implemented as 3D IQMs. RMSE was implemented by calculating the voxelwise errors between the reconstructed image
ID and the reference image
IR, which is given by equation [3]:
The terms
Nx,
Ny, and
Nz in equation [3] represent the number of voxels in the images along each dimension. We then normalized the value given by equation [3] by dividing by the average intensity of the reference image. A perfect RMSE score is 0, indicating no difference at all between the voxel contents in
ID and
IR.
The image processing toolbox in Matlab includes a 3D implementation of SSIM (21), which we used for this work. SSIM was calculated by equation [4]:
where
μR and
μD represent the average of the reference image and the reconstructed image in an 11 × 11 × 11 voxel window, , respectively;
σR and
σD represent the standard deviations of voxel values within the window in each image; and
σRD represents the crosscovariance for the images in the window. The constants were chosen as
C1=0.01L and
C2=0.03L, where
L is the maximum value of the image. The window size and the parameters specified were chosen for calculating SSIM because they represented the default values suggested by Wang (
5) and implemented within Matlab (
21). SSIM ranges between −1 and 1, with the score of 1 indicating optimal structural similarity between the
ID and
IR.
For IQM calculations during the retrospective static pelvic MRI investigation, IQM scores were calculated by comparing the CSMRI reconstruction with the original reference image.
For IQM calculations during the dynamic simulations, RMSE and SSIM were calculated for each reconstructed image by comparing the reconstruction with the references outlined previously. This resulted in many IQM scores over time for a single simulation, which were then averaged to output a single IQM score that was representative of the entire simulation. Timecourses were extracted voxel by voxel from the expected feature locations in the reconstructed images. The exponential model parameters were recovered using the lsqcurvefit function in the Matlab optimization toolbox (22). The Tofts model parameters were recovered using a modified version of the ROCKETSHIP toolbox (23). The accuracy of recovered model parameters was compared with the known input parameters to assess the quantitative accuracy of the reconstructed temporal data.
Results
Retrospective Static Pelvic CSMRI
Figure 2 shows the effect of varying the regularization weight for a wavelet sparsity constraint and the undersampling factor for static CSMRI reconstructions of retrospectively undersampled data, from a representative patient and from static numerical phantom data. In general, higher SSIM score corresponded to lower RMSE score. The regularization weight resulting in preferential IQM score (defined as the maximal SSIM or the minimal RMSE) tended to differ between the IQMs at each undersampling factor. This is shown in Figure 3, which displays the regularization weight resulting in preferential IQM score for static CSMRI reconstructions at each undersampling factor for all 15 patients in the study and our numerical phantom. It was additionally observed that the preferential regularization weight varied among patients in the study. The preferential regularization weight increased monotonically with an increase in the undersampling of kspace for most patients, but not for all.
Figure 2.
Examples of the effect of undersampling factor and regularization weight on the root mean square error (RMSE) and structural similarity index (SSIM) scores for static compressed sensing–magnetic resonance imaging reconstructions. RMSE and SSIM scores for representative retrospectively undersampled patient data are shown in (A) and (C), respectively, whereas RMSE and SSIM scores for numerical phantom data are shown in (B) and (D), respectively. Unprocessed images of the patient and of the numerical phantom are shown in (E) and (F), respectively.
Figure 3.
Preferential regularization weights with increasing undersampling factor are shown for the RMSE and the SSIM in (A) and (B), respectively. The preferential score was defined as the maximal SSIM or the minimal RMSE at each undersampling factor.
Simulations of Dynamic Digital Phantom With Monoexponential Signal Evolution
Figure 4 shows the effect of changing the regularization weight and the number of packets used in a simple CSDCEMRI scenario, where features evolve in intensity according to a monoexponential decay. While choosing a regularization weight of 0 consistently resulted in poorer performance, nonzero regularization weights were observed to improve the timecourse fidelity. Using too few packets resulted in poor timecourse fidelity owing to image degradation, and it was generally observed that proper choice of the number of packets included in a reconstruction had a larger effect than tuning of nonzero regularization weights.
Figure 4.
Examples of the fidelity of the average reconstructed time course for an exponentially decaying feature. There were 11 packets used in (A) as the regularization weight was varied, whereas (B) used a regularization weight of 0.01 and varied the number of packets included in image reconstruction.
Figure 5 shows the relationship between the recovery accuracy of quantitative parameters describing a monoexponentially decaying feature and the image quality as assessed via RMSE and SSIM. Within Figure 5, reconstructions made using different numbers of packets have been color coded, and each colorcoded group contains reconstructions made with regularization weights ranging from 0 to 0.05. Within each colorcoded group, an outlier data point may be observed at a worse IQM score. These outliers corresponded to reconstructions made with a regularization weight of 0. Increasing the number of packets used in image reconstruction was only beneficial to a point, beyond which the error in each parameter again grew in magnitude. This outcome shows the expected tradeoff between temporal resolution and image quality.
Figure 5.
Relationship between RMSE scores with the recovery accuracy of quantitative parameters in simulated data of a monoexponentially decaying feature are shown in (A) and (C), whereas the relationship with the SSIM are shown in (B) and (D).
For both the amplitude and the decay parameter, a correlation between the error of the parameter recovery and each IQM score was observed. These correlations were quantified using the Pearson linear correlation coefficient (PLCC), which measures linear correlations, and the Spearman rank correlation coefficient (SRCC), which measures the extent of a more general monotonic relationship. The values of the PLCC and SRCC for the data in Figure 5 are shown in Table 1.
Table 1.
Correlations between IQM Scores and Error in the Recovered Monoexponential Decay Parameters for the Data in Figure 5


RMSE

SSIM

Initial Amplitude 
PLCC 
−0.973 
0.982 
SRCC 
−0.986 
0.949 
Decay Rate 
PLCC 
0.959 
−0.963 
SRCC 
0.937 
−0.882 
Simulations of Dynamic Digital Phantom With Tofts Model Signal Evolution
Figure 6 shows the effect of changing the regularization weight and number of packets in a CSMRI reconstruction on the recovered timecourse fidelity of features whose intensity evolves according to equation [2]. The same observations seen in the monoexponentially evolving phantom simulation were observed in this case. In addition, the choice of temporal resolution, made through the choice of the number of packets to include in an image reconstruction, most affected the ability to resolve the period of rapid contrast enhancement.
Figure 6.
Examples of the fidelity of the average reconstructed time course for a voxel evolving according to the Tofts model described by equation [2], where the number of packets was held constant in (A) and the regularization weight was held constant in (B).
Figure 7 shows the relationship between the recovery accuracy of
Ktran and the image quality as assessed via RMSE and SSIM. Color coding follows the same convention as described for Figure 5. Both RMSE and SSIM show the ability to define a range of IQM scores within which
Ktran error is minimized, similar to the results for the simple simulation with a monoexponentially evolving phantom. The correlations between the recovery accuracy of
Ktrans and RMSE and SSIM scores are shown in Table 2.
Figure 7.
Relationship between RMSE scores with the recovery error in quantitative parameters in simulated data of a feature evolving according to the Tofts model described by equation [2] are shown in (A), whereas the relationship with the SSIM is shown in (B).
Table 2.
Correlations between IQM Scores and the Error in the Recovered Tofts Model Parameter for the Data in Figure 7


RMSE

SSIM

Ktran

PLCC 
−0.986 
0.968 
SRCC 
−0.975 
0.929 
Figure 8 shows the relationship between IQM scores and the recovery accuracy of
Ktran for the 4 different reference image scenarios described earlier, reconstructed with a regularization weight of 0.05. The choice of reference image was found to change the absolute values of the IQM score corresponding to a given parameter recovery accuracy and, in particular, the score corresponding to minimum parameter error. However, a clear relationship between IQM score and parameter accuracy was observed for all choices, and the observed correlations are shown in Table 3.
Figure 8.
Effect of choosing each of the different reference images, and the relationship that results between each image quality metric and the recovery accuracy of
Ktran. Results for the RMSE are shown in (A), whereas results for the SSIM are shown in (B). Images were reconstructed with a regularization weight of 0.05.
Table 3.
Correlations between IQM Scores and Error in the Recovered
Ktran for the Data Shown in Figure 8


RMSE

SSIM

Instantaneous Truth 
PLCC 
−0.997 
0.998 
SRCC 
−1.000 
1.000 
Precontrast 
PLCC 
−0.976 
0.988 
SRCC 
−1.000 
1.000 
Peak AIF 
PLCC 
−1.000 
0.999 
SRCC 
−1.000 
1.000 
Latecontrast 
PLCC 
−0.995 
0.998 
SRCC 
−1.000 
1.000 
Discussion
Retrospective Static Pelvic CSMRI
We assessed the response of RMSE and SSIM to wavelet regularization weight and kspace undersampling factor. Both IQMs show that there is a general decline in image quality as the undersampling factor is increased. Furthermore, at each undersampling factor, the effects of under or overregularizing the image are captured by the IQMs as a variation in their reported scores. There was also a tendency for the preferential regularization weight at each undersampling factor to differ slightly between patients and to not always change monotonically. This was true for both RMSE and SSIM, showing their sensitivity to CSMRI reconstruction process for individual patients. The patient results potentially suggest the presence of 2 subpopulations of preferred CS regularization weights, evident in Figure 3; future research can be conducted to assess which image properties caused differences between the responses of the IQMs to each patient on an individual basis that may lead to such subpopulations in CSMRI reconstruction parameters. The present results serve to confirm that each IQM is sensitive to the CSMRI reconstruction process on a patientbypatient basis when reconstructing goldenanglesampled data.
These IQMs could thus be used on a patientbypatient basis to guide the selection of image reconstruction parameters. These results motivate the assessment of IQMs in selecting image acquisition parameters for CSDCEMRI, where the appropriate choice of temporal resolution (ie, undersampling factor) and CS regularization weight is essential.
We note that the magnitude of each IQM, as well as the preferential regularization weight at each undersampling factor, tended to differ in magnitude between the synthetic phantom and the patient data. This is in part owing to the highly nonanthropomorphic nature of our simple synthetic phantom. Although our present work has established a basic understanding of the expected trends in RMSE and SSIM scores as a function of CS regularization weight and undersampling factor, we stress that our current work does not offer a quantitative rule to select these parameters. Future work may use a more anthropomorphically correct phantom that is expected to more closely mirror the reconstruction parameters desirable for real patient data.
We acknowledge that in practice, RMSE and SSIM may not be suitable for the assessment of prospectively accelerated static images. This is owing to the lack of a suitable reference image and the redundancy in generating one. However, the response of RMSE and SSIM IQMs in our simulated study indicates that alternative referencefree IQMs might be a worthwhile area of investigation for this purpose.
Simulations of Dynamic Digital Phantom With Monoexponential Signal Evolution
To validate our simulation framework in the context of CSDCEMRI, we first simulated a phantom containing features that evolved according to a simple monoexponential function. This simple model allowed us to determine the effects of CSDCEMRI reconstructions on the accuracy of recovery of parameters describing contrast enhancement.
Image reconstruction parameters affected the recovered timecourses as we expected they should have. The reconstructed timecourses were most affected by the number of packets included in the reconstruction. While the choice of regularization weight was also found to influence the accuracy of the reconstructed timecourse, in practice, any reasonable nonzero regularization weight was found to be beneficial.
Choosing too high of a temporal resolution (ie, a lower number of packets per reconstruction) degraded the image quality owing to the increased presence of incoherent aliasing artifacts, whereas choosing too low of a temporal resolution (ie, a larger number of packets per reconstruction) reduced the ability to capture the dynamic signal evolution. This outcome shows the expected compromise between temporal resolution and image quality.
We observed that the scalerelated (ie, amplitude) and raterelated (ie, decay constant) parameters in our monoexponential model responded differently to the adjustment of CS regularization weight. In general, the use of a nonzero regularization weight was found to improve the IQM scores without affecting the recovery accuracy of the rate parameter. When recovering the scale parameter with a nonzero regularization weight, the recovery accuracy was found to have a positive offset compared with recovery with a regularization weight of 0. This shows the sensitivity of various model parameters to aspects of the image reconstruction process. For example, the scale parameter should be sensitive to losses in overall image intensity, whereas the rate parameter should be relatively robust to such an effect.
We also note that in the limit of a perfect IQM score (SSIM = 1 or RMSE = 0), the results would not yield perfect recovery accuracy of quantitative parameters. To achieve a perfect IQM score, assuming a noiseless simulated acquisition, a complete sampling of kspace would be required for a given image reconstruction. The temporal span of all the interleaved data for each image reconstruction would then span a period long enough to introduce blurring of the signal evolution. This tradeoff in image reconstruction quality for the appropriate temporal resolution is captured by each IQM in our results, where continuously increasing the number of packets included in each image reconstruction was found to improve the IQM scores while degrading parameter recovery accuracy past a certain point.
The results from the monoexponentially evolving features thus agree with predictable outcomes, verifying the utility of our simulation framework for CSDCEMRI investigations. In particular, these results allow us to confidently determine a more complex contrastenhancement model.
Simulations of Dynamic Digital Phantom with Tofts Model Signal Evolution
The Tofts model exhibits an underlying signal evolution with more complex and physiologically relevant features. This increased temporal signal complexity further motivates the requirement to properly select the image reconstruction parameters in a physiologically realistic CSDCEMRI application.
Similar to the monoexponentially evolving signal model, the reconstructed timecourse for the Tofts model was most sensitive to the number of packets included in the reconstruction. The effect of choosing too low or too high of a temporal resolution was especially pronounced on the recovered timecourses owing to their more complex structure; the period of rapid signal enhancement between 0 and 10 seconds was especially prone to losing its structure.
If RMSE or SSIM were to be implemented in a clinical setting, there would be no “instantaneous reference image,” unlike in our simulated data, and some appropriate choice of a reference image would need to be made. In addition to IQM scores calculated by using the instantaneous reference image, we determined the possible use of 3 realistic methods of obtaining a datagenerated reference image. Despite the more complex model for signal evolution under the Tofts model, the recovery accuracy of
Ktran was found to correlate very well with both IQMs for all 4 reference images tested. This suggests that RMSE and SSIM may be used to select the appropriate reconstruction parameters in CSDCEMRI to reach a minimum tolerance of parameter recovery accuracy. These results motivate further assessment of the use of these IQMs to guide CSDCEMRI reconstructions for patients in clinical settings. For example, algorithms incorporating RMSE or SSIM could be developed for the automated selection of image reconstruction parameters in goldenanglesampled CSDCEMRI.
We note that we have made direct use of a populationaveraged AIF in our present study. Patientspecific measurement of pharmacokinetic parameters should ideally make use of a patientspecific AIF, which affects the
Cp(t) term in equation [2] through measurement uncertainties introduced by inadequate temporal resolution. Our direct use of a populationaveraged AIF ignores the extra error that would be introduced into parameter reconstructions from a patientspecific AIF measurement; however, the goal of our present study was to establish a baseline from which such factors could be studied. Future work may begin to determine the effect of simulated “in vivo” measurements of AIF on the accuracy of the reconstructed parameters and the resulting correlations with IQMs such as RMSE or SSIM.
Although the extended Tofts model (equation [2]) is parametrized by 3 parameters, namely,
ve,
vp, and
Ktran, we chose to determine the relationship between RMSE and SSIM scores and the accuracy of recovery of
Ktran only. We chose to focus exclusively on
Ktran because previous work has found it to be of primary clinical performance. In the diagnosis of prostate cancer, for example, it was found that the inclusion of pharmacokinetic parameters aside from
Ktran does not improve performance of datadriven cancer detection using DCEMRI techniques (24), and several successful computeraided diagnosis tools exist that make use of only
Ktran (25, 26). A subject for future assessment would be the characterization of the relationship of different pharmacokinetic parameters to objective measures of image quality, such as RMSE or SSIM.
In general, our current work has shown preliminary results showing that objective measures of image quality can serve as predictors of dynamic MRI data quality, which opens new research avenues. Future work may determine the use of other IQMs in CSDCEMRI applications or in alternative dynamic MRI applications. This includes both IQMs requiring reference images and those that do not require a reference image, such as the perceptionbased image quality evaluator IQM (27) or the blind/referenceless image spatial quality evaluator IQM (28).
Conclusion
We have determined the relationship between 2 IQMs, RMSE and SSIM, and their correlation with the recovery accuracy of
Ktran in CSDCEMRI. To accomplish this, we began with a study of the properties of each IQM in retrospectively undersampled static pelvic CSMR images. Each IQM was found to be sensitive to changes in the CSMRI reconstruction process on a patientbypatient basis.
We then studied RMSE and SSIM in CSDCEMRI by designing a simulation framework in Matlab. To validate our simulation framework, we assessed a numerical phantom containing features that evolved according to a simple monoexponential decay. The simulation results agreed with expected outcomes, validating the use of our framework for more complex contrast evolution models such as the extended Tofts model.
Using our simulation framework, we determined the ability of RMSE and of SSIM to predict the recovery accuracy of
Ktran. Both IQMs require the use of a defined reference image. For 4 different methods of selecting a reference image (one based on the instantaneous image of the simulated phantom, one based on precontrast imaging, one based on latecontrast imaging, and one based on the peak of the AIF used in the Tofts model), we found that both RMSE and SSIM consistently correlate well with the recovery accuracy of
Ktran. This motivates further assessment into the use of IQMs to guide CSDCEMRI in clinical settings, such as in the design of automated algorithms for tuning the reconstruction parameters of CSDCEMRI acquisitions. More generally, these results suggest that objective measures of image quality have utility in the guidance of dynamic MRI techniques.
Research Articles
Download PDF (2.92 MB)
TOMOGRAPHY, December 2020, Volume 6, Issue 4:362372
DOI: 10.18383/j.tom.2020.00045
Evaluation of GoldenAngleSampled Dynamic ContrastEnhanced MRI Reconstruction Using Objective Image Quality Measures: A Simulated Phantom Study
Nathan Murtha^{1}, Allister Mason^{2}, Chris Bowen^{3}, Sharon Clarke^{3}, James Rioux^{3}, Steven Beyea^{3}
Abstract
We aim to extend the use of image quality metrics (IQMs) from static magnetic resonance imaging (MRI) applications to dynamic MRI studies. We assessed the use of 2 IQMs, the root mean square error and structural similarity index, in evaluating the reconstruction of quantitative dynamic contrastenhanced (DCE) MRI data acquired using goldenangle sampling and compressed sensing (CS). To address the difficulty of obtaining groundtruth knowledge of parameters describing dynamics in real patient data, we developed a Matlab simulation framework to assess quantitative CSDCEMRI. We began by validating the response of each IQM to the CSMRI reconstruction process using static data and the performance of our simulation framework with simple dynamic data. We then extended the simulations to the more realistic extended Tofts model. When assessing the Tofts model, we tested 4 different methods of selecting a reference image for the IQMs. Results from the retrospective static CSMRI reconstructions showed that each IQM is responsive to the CSMRI reconstruction process. Simulations of a simple contrast evolution model validated the performance of our framework. Despite the complexity of the Tofts model, both IQM scores correlated well with the recovery accuracy of a central model parameter for all reference cases studied. This finding may form the basis of algorithms for automated selection of image reconstruction aspects, such as temporal resolution, in goldenanglesampled CSDCEMRI. These further suggest that objective measures of image quality may find use in general dynamic MRI applications.
Introduction
Motivation for This Work
Magnetic resonance imaging (MRI) can not only provide static anatomical information but also capture dynamic processes such as the uptake of an intravenously injected contrast agent with dynamic contrastenhanced MRI (DCEMRI). Rather than using a fixed temporal resolution for such dynamic scans, emerging methods based on goldenangle sampling such as XDGRASP (1) or CIRCUS (2) allow retrospective temporal resampling at the time of image reconstruction. This framework allows retrospective control of the temporal “footprint” by varying the number of kspace rays or spokes per undersampled image, without necessarily requiring the use of view sharing, which can introduce temporal blurring of fast dynamics (2, 3). Doing so allows one to flexibly optimize the tradeoff between better spatial quality (ie, low undersampling factor) and better temporal quality (ie, high undersampling factor).
With additional reconstruction flexibility, however, comes the need to optimize image reconstruction of that specific data set and its associated spatiotemporal features. This is particularly true when one considers the large and intertwined joint optimization space associated with these reconstructions, in which the undersampling factor (ie, temporal footprint), the form of regularization enforced, and the degree of regularization must all be simultaneously examined.
Recently, it has become common to use objective image quality metrics (IQMs) (4, 5) as an outcome measure or within the cost function of compressed sensing (CS) or with deep learning–based image reconstruction (6–11). A recent study by our group examined a broad range of IQMs to determine the varying ability of these metrics to correlate with expert radiologists’ rating of diagnostic image quality (12). Although such an approach is sufficient in the reconstruction of a single image data set with fixed undersampling factor, in a data set with goldenangle sampling that allows multiple potential reconstructions, most IQMs will be maximal for the lowest undersampling factor possible, that is, with no ability to resolve temporal dynamics. Therefore, additional research is required to make use of IQMs in such a scenario, and the approach using expert radiologist opinion may be impractical.
Furthermore, the kspace sampling and image reconstruction factors that maximize the quality of diagnostic images from a dynamic data set may not be the same ones that minimize the parametric error when a model is subsequently fit to the data set in a quantitative dynamic application, such as the Tofts pharmacokinetic model in the case of quantitative DCEMRI (13, 14). For example, while a radiologist may visually prefer a diagnostic image with less noise, a reconstruction with increased image noise but improved temporal frame rate may potentially result in more accurate parameter fits. It is also not feasible to present radiologists with thousands of images resulting from the numerous permutations of kspace sampling and reconstruction choices. The ability to automate the joint tuning of goldenanglesampled dynamic data reconstruction, for example, through the use of an IQMbased approach, would, therefore, potentially improve the accuracy of recovered parameters in quantitative dynamic applications.
Goals of This Study
We determined the relationship between 2 wellestablished IQMs—root mean square error (RMSE) and structural similarity index (SSIM) (5)—and the subsequent fidelity of recovered dynamic model parameters when applied to DCEMRI with CS reconstruction of goldenanglesampled data.
We began by assessing the properties of RMSE and SSIM in static CSMRI reconstructions using retrospectively undersampled data from pelvic scans as well as a numerical phantom. This verified that each IQM was sensitive to the CS reconstruction parameters on a patientbypatient basis.
With verification of IQM response in basic CSMRI application, we then evaluated their ability to predict the quality of goldenanglesampled dynamic MRI data reconstructed with CS. To this end, we developed a simplified simulation framework for a dynamically evolving numerical phantom as a surrogate for the assessment of CSMRI in dynamic MRI procedures. The use of a numerical phantom allowed for groundtruth knowledge of the instantaneous temporal evolution. To validate our simulation framework, we first analyzed cylindrical features that evolved according to a monoexponential decay within a digital phantom. Once validated on a simple closedform signal evolution model, we analyzed the more physiologically relevant extended Tofts model of pharmacokinetics, which models the uptake of a contrast media in tissue.
Methodology
Modified CIRCUS kSpace Acquisition
Goldenangle sampling of kspace was simulated using a modified CIRCUS acquisition (15), with a repetition time of 3.7 milliseconds. Although multiple goldenangle sampling techniques are possible, we chose CIRCUS because it is a Cartesian trajectory, making for easy implementation and analysis. CIRCUS subdivides kspace into nested squares on a Cartesian grid, with rudimentary sampling patterns formed by selecting one point from each square according to a goldenangle rotation (2). We refer to these rudimentary sampling patterns as “quanta.” The 2D CIRCUS patterns are extended for use in 3D imaging by extending the sampling mask along the third dimension (ie, the frequency encoded dimension). The shapes of each quanta are controlled by parametersb and
c , which apply linear and nonlinear shifts to the position of successive samples and may loosely be thought of as contributing to the “shear” and “twist” of the quanta. The value of
b must be positive, and
c should be between 1 and 2; using
b = 0 and
c = 1 yields quanta that are straight lines, whereas
b ≫ 1 and
c = 2 maximizes randomness. For this work, we set
b = 40 and
c = 1.5 , which are consistent with the values recommended by the authors of CIRCUS.
Individual quanta are too highly undersampled to yield useful images, so they were combined into groups of 10 to form “packets.” Quantities of packets ranging from 3 to 21 in steps of 2 were then interleaved in various numbers when forming image reconstructions. The number of packets interleaved then dictated the temporal resolution of the resulting reconstructions, where more packets meant a larger temporal footprint.
The kspace data were sampled by taking the Fourier transform of the image and selecting the coefficients indicated by the CIRCUS packets. In addition, we opted to perform the studies presented herein without the addition of simulated noise to the signal. We acknowledge this to be a simplification of the underlying signal acquisition physics; however, it allows a preliminary investigation of the application of IQMs to the analysis of dynamically evolving data and to determine their fundamental response to such data in the absence of external factors, which is the goal of the present work.
Retrospective Static Pelvic CSMRI
To show the performance of IQMs in static MRI data as a function of reconstruction parameters, 15 static pelvic images were retrospectively undersampled. The acquisition of patient data for this research was approved by the Nova Scotia Health Authority Research Ethics Board and complied with all ethical obligations. Data were acquired on a GE MR750 3 T system using a 32channel abdominal radiofrequency coil array and a 3D spoiled gradient echo sequence (DISCO) (16) with fat separation. Data were acquired with autocalibrated parallel imaging (acceleration factor 2) in the phaseencoding direction, with all omitted lines reconstructed before retrospective undersampling.
Retrospective undersampling of the static pelvic images was performed using the modified CIRCUS acquisition described previously. Combinations of packets were chosen resulting in undersampling factors ranging from 1 to 12 in steps of 0.5. CS reconstructions used the Berkeley Advanced Reconstruction Toolbox (17) to solve the following system:
(1)
For this reconstruction, we chose a wavelet basis forΨ and used values of
λ ranging from 0 to 0.08 in steps of 0.002. The quality of the resulting image reconstructions was assessed by RMSE and SSIM.
Dynamic CSMRI Simulations
Simulations were written in Matlab 2016b (The MathWorks, Natick, MA). A digital phantom of size 200 × 200 × 32 voxels was designed to contain features that evolve over time. The embedded features of the phantom are shown in Figure 1. The plane of the embedded cylindrical features at the very center of the phantom was chosen for assessment in this work, with the remaining features in the phantom intended to provide flexibility for future work. Within the plane of the embedded cylindrical features used, we placed regions of interest (ROIs) inside the inner cylindrical features of a row of cylinders, shown in yellow in Figure 1, for quantitative parameter reconstructions. The ROIs for each embedded cylindrical feature were combined and averaged for analysis. In total, 596 voxels were used in the combined ROIs.
Figure 1.
Internal features of the digital phantom embedded in a homogeneous static background. Within each plane of pingrids (shown in brown), all grids can be made to evolve together according to a userdesigned function. Within each plane of embedded cylinders, the outer cylinders (shown in orange) can be made to evolve together according to a userdesigned function, and the inner cylinders (shown in yellow) can be made to evolve independently along each row. The interior cylinders in one plane of embedded cylindrical features, indicated by a red arrow, were used for this work; a 2D view of this slice is given in Figure 2F.
For this work, we simulated 2 scenarios. To begin, a simple monoexponential decay model was simulated. This model consisted of 2 free parameters—the initial amplitude and the decay rate of the signal amplitude. In total, 200 packets were acquired during the simulated evolution of the phantom, with individual kspace samples acquired at a simulated interval of 3.7 milliseconds. This result corresponded to 2 minutes and 32 seconds of simulated data acquisition. From these simulations, we fit the model to the reconstructed timecourses extracted from each voxel to recover both the initial magnitude and decay rate parameters. The use of a simple model such as this one validated that the simulations were providing rational results before extending it to a more realistic and complex one such as the extended Tofts model. For validation, we expected the following:
A tradeoff between temporal resolution and image quality would be evident.
The temporal footprint (ie, the undersampling factor) used would have a larger effect than the degree of regularization applied.
Following the simple simulations described previously, features in an extended Tofts model were then simulated to evolve according to equation [2]:
(2)
The concentrations of contrast agent in the tissue and in the blood plasma are represented in equation [2] byC t ( t ) and
C p ( t ) , respectively. The terms
v e and
v p represent the volume fraction in tissue of the extravascular extracellular space and blood plasma, respectively. The term
K t r a n in equation [2] is a transfer constant, which is indicative of the permeability of the vasculature. From these simulations, we fit equation [2] to the reconstructed timecourses extracted from each voxel to recover
K t r a n . To implement equation [2], we used a populationaveraged arterial input function (AIF) to define
C p ( t ) at each point in time (19). The Tofts model–defined timecourses evolved for 2 minutes and 30 seconds, with the contrast agent introduced after 30 seconds. This allowed for the collection of 197 packets.
We interleaved differing numbers of packets before CSMRI image reconstruction when making reconstructions of both the monoexponential decay and the Tofts model scenarios. Combinations of 3–21 packets, in steps of 2 packets, were supplied to Berkeley Advanced Reconstruction Toolbox’s pics command, which performs an eSPIRiTbased parallel imaging and CS image reconstruction (20), enforcing a wavelet sparsity constraint with regularization weights of 0–0.05 in steps of 0.005. This corresponded to undersampling factors ranging from 12.5 to 2, where the undersampling factor is defined as the reciprocal of the fraction of the kspace matrix that is sampled, for example, sampling 50% of the kspace coefficients results in an undersampling factor of 2. We chose to present the results of this portion of the study as a function of the number of packets, rather than of the undersampling factor, because it is the number of packets that the user has control over. These packet combinations corresponded to simulated temporal windows of ∼2–13 seconds.
For the simple dynamic simulations consisting of a monoexponential decay, the reference image was taken to be the instantaneous value of the numerical phantom at the average time of all constituent packets; this is referred to as the instantaneous reference image. For the dynamic simulations using the Tofts model of signal evolution in equation [2], 3 additional reference images were generated from the simulated data. In practice, an instantaneous reference image would not be available, so this step allowed us to test the feasibility of using a datagenerated reference image. The following were the 3 datagenerated references:
A precontrast reference image made by acquiring an extra 30 seconds of data before contrast enhancement and combining all packets up until 5 seconds before contrast enhancement.
A peakAIF reference image, made by combining all packets in a 20second window starting at the peak of the AIF used to getC p ( t ) in equation [2].
A latecontrast reference image, made by combining all packets in the last 20 seconds of data acquisition where the contrast dynamics are relatively slow.
In testing these 3 additional datagenerated references, reconstructions were made using 3–21 packets as described previously, but only a single regularization weight of 0.05 was used.
Image Quality Assessment
RMSE and SSIM were implemented as 3D IQMs. RMSE was implemented by calculating the voxelwise errors between the reconstructed imageI D and the reference image
I R , which is given by equation [3]:
(3)
The termsN x ,
N y , and
N z in equation [3] represent the number of voxels in the images along each dimension. We then normalized the value given by equation [3] by dividing by the average intensity of the reference image. A perfect RMSE score is 0, indicating no difference at all between the voxel contents in
I D and
I R .
The image processing toolbox in Matlab includes a 3D implementation of SSIM (21), which we used for this work. SSIM was calculated by equation [4]:
(4)
For IQM calculations during the retrospective static pelvic MRI investigation, IQM scores were calculated by comparing the CSMRI reconstruction with the original reference image.
For IQM calculations during the dynamic simulations, RMSE and SSIM were calculated for each reconstructed image by comparing the reconstruction with the references outlined previously. This resulted in many IQM scores over time for a single simulation, which were then averaged to output a single IQM score that was representative of the entire simulation. Timecourses were extracted voxel by voxel from the expected feature locations in the reconstructed images. The exponential model parameters were recovered using the lsqcurvefit function in the Matlab optimization toolbox (22). The Tofts model parameters were recovered using a modified version of the ROCKETSHIP toolbox (23). The accuracy of recovered model parameters was compared with the known input parameters to assess the quantitative accuracy of the reconstructed temporal data.
Results
Retrospective Static Pelvic CSMRI
Figure 2 shows the effect of varying the regularization weight for a wavelet sparsity constraint and the undersampling factor for static CSMRI reconstructions of retrospectively undersampled data, from a representative patient and from static numerical phantom data. In general, higher SSIM score corresponded to lower RMSE score. The regularization weight resulting in preferential IQM score (defined as the maximal SSIM or the minimal RMSE) tended to differ between the IQMs at each undersampling factor. This is shown in Figure 3, which displays the regularization weight resulting in preferential IQM score for static CSMRI reconstructions at each undersampling factor for all 15 patients in the study and our numerical phantom. It was additionally observed that the preferential regularization weight varied among patients in the study. The preferential regularization weight increased monotonically with an increase in the undersampling of kspace for most patients, but not for all.
Figure 2.
Examples of the effect of undersampling factor and regularization weight on the root mean square error (RMSE) and structural similarity index (SSIM) scores for static compressed sensing–magnetic resonance imaging reconstructions. RMSE and SSIM scores for representative retrospectively undersampled patient data are shown in (A) and (C), respectively, whereas RMSE and SSIM scores for numerical phantom data are shown in (B) and (D), respectively. Unprocessed images of the patient and of the numerical phantom are shown in (E) and (F), respectively.
Figure 3.
Preferential regularization weights with increasing undersampling factor are shown for the RMSE and the SSIM in (A) and (B), respectively. The preferential score was defined as the maximal SSIM or the minimal RMSE at each undersampling factor.
Simulations of Dynamic Digital Phantom With Monoexponential Signal Evolution
Figure 4 shows the effect of changing the regularization weight and the number of packets used in a simple CSDCEMRI scenario, where features evolve in intensity according to a monoexponential decay. While choosing a regularization weight of 0 consistently resulted in poorer performance, nonzero regularization weights were observed to improve the timecourse fidelity. Using too few packets resulted in poor timecourse fidelity owing to image degradation, and it was generally observed that proper choice of the number of packets included in a reconstruction had a larger effect than tuning of nonzero regularization weights.
Figure 4.
Examples of the fidelity of the average reconstructed time course for an exponentially decaying feature. There were 11 packets used in (A) as the regularization weight was varied, whereas (B) used a regularization weight of 0.01 and varied the number of packets included in image reconstruction.
Figure 5 shows the relationship between the recovery accuracy of quantitative parameters describing a monoexponentially decaying feature and the image quality as assessed via RMSE and SSIM. Within Figure 5, reconstructions made using different numbers of packets have been color coded, and each colorcoded group contains reconstructions made with regularization weights ranging from 0 to 0.05. Within each colorcoded group, an outlier data point may be observed at a worse IQM score. These outliers corresponded to reconstructions made with a regularization weight of 0. Increasing the number of packets used in image reconstruction was only beneficial to a point, beyond which the error in each parameter again grew in magnitude. This outcome shows the expected tradeoff between temporal resolution and image quality.
Figure 5.
Relationship between RMSE scores with the recovery accuracy of quantitative parameters in simulated data of a monoexponentially decaying feature are shown in (A) and (C), whereas the relationship with the SSIM are shown in (B) and (D).
For both the amplitude and the decay parameter, a correlation between the error of the parameter recovery and each IQM score was observed. These correlations were quantified using the Pearson linear correlation coefficient (PLCC), which measures linear correlations, and the Spearman rank correlation coefficient (SRCC), which measures the extent of a more general monotonic relationship. The values of the PLCC and SRCC for the data in Figure 5 are shown in Table 1.
Table 1.
Correlations between IQM Scores and Error in the Recovered Monoexponential Decay Parameters for the Data in Figure 5
i] Abbreviations: PLCC, Pearson linear correlation coefficient; SRCC, Spearman rank correlation coefficient; RMSE, root mean squared error; SSIM, structural similarity index
Simulations of Dynamic Digital Phantom With Tofts Model Signal Evolution
Figure 6 shows the effect of changing the regularization weight and number of packets in a CSMRI reconstruction on the recovered timecourse fidelity of features whose intensity evolves according to equation [2]. The same observations seen in the monoexponentially evolving phantom simulation were observed in this case. In addition, the choice of temporal resolution, made through the choice of the number of packets to include in an image reconstruction, most affected the ability to resolve the period of rapid contrast enhancement.
Figure 6.
Examples of the fidelity of the average reconstructed time course for a voxel evolving according to the Tofts model described by equation [2], where the number of packets was held constant in (A) and the regularization weight was held constant in (B).
Figure 7 shows the relationship between the recovery accuracy ofK t r a n and the image quality as assessed via RMSE and SSIM. Color coding follows the same convention as described for Figure 5. Both RMSE and SSIM show the ability to define a range of IQM scores within which
K t r a n error is minimized, similar to the results for the simple simulation with a monoexponentially evolving phantom. The correlations between the recovery accuracy of
K trans and RMSE and SSIM scores are shown in Table 2.
Figure 7.
Relationship between RMSE scores with the recovery error in quantitative parameters in simulated data of a feature evolving according to the Tofts model described by equation [2] are shown in (A), whereas the relationship with the SSIM is shown in (B).
Table 2.
Correlations between IQM Scores and the Error in the Recovered Tofts Model Parameter for the Data in Figure 7
i] Abbreviations: PLCC, Pearson linear correlation coefficient; SRCC, Spearman rank correlation coefficient; RMSE, root mean squared error; SSIM, structural similarity index
Figure 8 shows the relationship between IQM scores and the recovery accuracy ofK t r a n for the 4 different reference image scenarios described earlier, reconstructed with a regularization weight of 0.05. The choice of reference image was found to change the absolute values of the IQM score corresponding to a given parameter recovery accuracy and, in particular, the score corresponding to minimum parameter error. However, a clear relationship between IQM score and parameter accuracy was observed for all choices, and the observed correlations are shown in Table 3.
Figure 8.
Effect of choosing each of the different reference images, and the relationship that results between each image quality metric and the recovery accuracy ofK t r a n . Results for the RMSE are shown in (A), whereas results for the SSIM are shown in (B). Images were reconstructed with a regularization weight of 0.05.
Table 3.
Correlations between IQM Scores and Error in the RecoveredK t r a n for the Data Shown in Figure 8
i] Abbreviations: PLCC, Pearson linear correlation coefficient; SRCC, Spearman rank correlation coefficient; RMSE, root mean squared error; SSIM, structural similarity index
Discussion
Retrospective Static Pelvic CSMRI
We assessed the response of RMSE and SSIM to wavelet regularization weight and kspace undersampling factor. Both IQMs show that there is a general decline in image quality as the undersampling factor is increased. Furthermore, at each undersampling factor, the effects of under or overregularizing the image are captured by the IQMs as a variation in their reported scores. There was also a tendency for the preferential regularization weight at each undersampling factor to differ slightly between patients and to not always change monotonically. This was true for both RMSE and SSIM, showing their sensitivity to CSMRI reconstruction process for individual patients. The patient results potentially suggest the presence of 2 subpopulations of preferred CS regularization weights, evident in Figure 3; future research can be conducted to assess which image properties caused differences between the responses of the IQMs to each patient on an individual basis that may lead to such subpopulations in CSMRI reconstruction parameters. The present results serve to confirm that each IQM is sensitive to the CSMRI reconstruction process on a patientbypatient basis when reconstructing goldenanglesampled data.
These IQMs could thus be used on a patientbypatient basis to guide the selection of image reconstruction parameters. These results motivate the assessment of IQMs in selecting image acquisition parameters for CSDCEMRI, where the appropriate choice of temporal resolution (ie, undersampling factor) and CS regularization weight is essential.
We note that the magnitude of each IQM, as well as the preferential regularization weight at each undersampling factor, tended to differ in magnitude between the synthetic phantom and the patient data. This is in part owing to the highly nonanthropomorphic nature of our simple synthetic phantom. Although our present work has established a basic understanding of the expected trends in RMSE and SSIM scores as a function of CS regularization weight and undersampling factor, we stress that our current work does not offer a quantitative rule to select these parameters. Future work may use a more anthropomorphically correct phantom that is expected to more closely mirror the reconstruction parameters desirable for real patient data.
We acknowledge that in practice, RMSE and SSIM may not be suitable for the assessment of prospectively accelerated static images. This is owing to the lack of a suitable reference image and the redundancy in generating one. However, the response of RMSE and SSIM IQMs in our simulated study indicates that alternative referencefree IQMs might be a worthwhile area of investigation for this purpose.
Simulations of Dynamic Digital Phantom With Monoexponential Signal Evolution
To validate our simulation framework in the context of CSDCEMRI, we first simulated a phantom containing features that evolved according to a simple monoexponential function. This simple model allowed us to determine the effects of CSDCEMRI reconstructions on the accuracy of recovery of parameters describing contrast enhancement.
Image reconstruction parameters affected the recovered timecourses as we expected they should have. The reconstructed timecourses were most affected by the number of packets included in the reconstruction. While the choice of regularization weight was also found to influence the accuracy of the reconstructed timecourse, in practice, any reasonable nonzero regularization weight was found to be beneficial.
Choosing too high of a temporal resolution (ie, a lower number of packets per reconstruction) degraded the image quality owing to the increased presence of incoherent aliasing artifacts, whereas choosing too low of a temporal resolution (ie, a larger number of packets per reconstruction) reduced the ability to capture the dynamic signal evolution. This outcome shows the expected compromise between temporal resolution and image quality.
We observed that the scalerelated (ie, amplitude) and raterelated (ie, decay constant) parameters in our monoexponential model responded differently to the adjustment of CS regularization weight. In general, the use of a nonzero regularization weight was found to improve the IQM scores without affecting the recovery accuracy of the rate parameter. When recovering the scale parameter with a nonzero regularization weight, the recovery accuracy was found to have a positive offset compared with recovery with a regularization weight of 0. This shows the sensitivity of various model parameters to aspects of the image reconstruction process. For example, the scale parameter should be sensitive to losses in overall image intensity, whereas the rate parameter should be relatively robust to such an effect.
We also note that in the limit of a perfect IQM score (SSIM = 1 or RMSE = 0), the results would not yield perfect recovery accuracy of quantitative parameters. To achieve a perfect IQM score, assuming a noiseless simulated acquisition, a complete sampling of kspace would be required for a given image reconstruction. The temporal span of all the interleaved data for each image reconstruction would then span a period long enough to introduce blurring of the signal evolution. This tradeoff in image reconstruction quality for the appropriate temporal resolution is captured by each IQM in our results, where continuously increasing the number of packets included in each image reconstruction was found to improve the IQM scores while degrading parameter recovery accuracy past a certain point.
The results from the monoexponentially evolving features thus agree with predictable outcomes, verifying the utility of our simulation framework for CSDCEMRI investigations. In particular, these results allow us to confidently determine a more complex contrastenhancement model.
Simulations of Dynamic Digital Phantom with Tofts Model Signal Evolution
The Tofts model exhibits an underlying signal evolution with more complex and physiologically relevant features. This increased temporal signal complexity further motivates the requirement to properly select the image reconstruction parameters in a physiologically realistic CSDCEMRI application.
Similar to the monoexponentially evolving signal model, the reconstructed timecourse for the Tofts model was most sensitive to the number of packets included in the reconstruction. The effect of choosing too low or too high of a temporal resolution was especially pronounced on the recovered timecourses owing to their more complex structure; the period of rapid signal enhancement between 0 and 10 seconds was especially prone to losing its structure.
If RMSE or SSIM were to be implemented in a clinical setting, there would be no “instantaneous reference image,” unlike in our simulated data, and some appropriate choice of a reference image would need to be made. In addition to IQM scores calculated by using the instantaneous reference image, we determined the possible use of 3 realistic methods of obtaining a datagenerated reference image. Despite the more complex model for signal evolution under the Tofts model, the recovery accuracy ofK t r a n was found to correlate very well with both IQMs for all 4 reference images tested. This suggests that RMSE and SSIM may be used to select the appropriate reconstruction parameters in CSDCEMRI to reach a minimum tolerance of parameter recovery accuracy. These results motivate further assessment of the use of these IQMs to guide CSDCEMRI reconstructions for patients in clinical settings. For example, algorithms incorporating RMSE or SSIM could be developed for the automated selection of image reconstruction parameters in goldenanglesampled CSDCEMRI.
We note that we have made direct use of a populationaveraged AIF in our present study. Patientspecific measurement of pharmacokinetic parameters should ideally make use of a patientspecific AIF, which affects theC p ( t ) term in equation [2] through measurement uncertainties introduced by inadequate temporal resolution. Our direct use of a populationaveraged AIF ignores the extra error that would be introduced into parameter reconstructions from a patientspecific AIF measurement; however, the goal of our present study was to establish a baseline from which such factors could be studied. Future work may begin to determine the effect of simulated “in vivo” measurements of AIF on the accuracy of the reconstructed parameters and the resulting correlations with IQMs such as RMSE or SSIM.
Although the extended Tofts model (equation [2]) is parametrized by 3 parameters, namely,v e ,
v p , and
K t r a n , we chose to determine the relationship between RMSE and SSIM scores and the accuracy of recovery of
K t r a n only. We chose to focus exclusively on
K t r a n because previous work has found it to be of primary clinical performance. In the diagnosis of prostate cancer, for example, it was found that the inclusion of pharmacokinetic parameters aside from
K t r a n does not improve performance of datadriven cancer detection using DCEMRI techniques (24), and several successful computeraided diagnosis tools exist that make use of only
K t r a n (25, 26). A subject for future assessment would be the characterization of the relationship of different pharmacokinetic parameters to objective measures of image quality, such as RMSE or SSIM.
In general, our current work has shown preliminary results showing that objective measures of image quality can serve as predictors of dynamic MRI data quality, which opens new research avenues. Future work may determine the use of other IQMs in CSDCEMRI applications or in alternative dynamic MRI applications. This includes both IQMs requiring reference images and those that do not require a reference image, such as the perceptionbased image quality evaluator IQM (27) or the blind/referenceless image spatial quality evaluator IQM (28).
Conclusion
We have determined the relationship between 2 IQMs, RMSE and SSIM, and their correlation with the recovery accuracy ofK t r a n in CSDCEMRI. To accomplish this, we began with a study of the properties of each IQM in retrospectively undersampled static pelvic CSMR images. Each IQM was found to be sensitive to changes in the CSMRI reconstruction process on a patientbypatient basis.
We then studied RMSE and SSIM in CSDCEMRI by designing a simulation framework in Matlab. To validate our simulation framework, we assessed a numerical phantom containing features that evolved according to a simple monoexponential decay. The simulation results agreed with expected outcomes, validating the use of our framework for more complex contrast evolution models such as the extended Tofts model.
Using our simulation framework, we determined the ability of RMSE and of SSIM to predict the recovery accuracy ofK t r a n . Both IQMs require the use of a defined reference image. For 4 different methods of selecting a reference image (one based on the instantaneous image of the simulated phantom, one based on precontrast imaging, one based on latecontrast imaging, and one based on the peak of the AIF used in the Tofts model), we found that both RMSE and SSIM consistently correlate well with the recovery accuracy of
K t r a n . This motivates further assessment into the use of IQMs to guide CSDCEMRI in clinical settings, such as in the design of automated algorithms for tuning the reconstruction parameters of CSDCEMRI acquisitions. More generally, these results suggest that objective measures of image quality have utility in the guidance of dynamic MRI techniques.
Notes
[4] Abbreviations:
IQM
image quality metric
RMSE
root mean square error
SSIM
structural similarity index
ROI
region of interest
PLCC
Pearson linear correlation coefficient
SRCC
Spearman rank correlation coefficient
AIF
arterial input function
MRI
magnetic resonance imaging
DCE
dynamic contrastenhanced
CS
compressed sensing
Acknowledgments
We gratefully acknowledge funding support from the Natural Sciences and Engineering Research Council Discovery Program, Brain Canada Platform Support Grant, the Atlantic Innovation Fund under the Atlantic Canada Opportunities Agency, and an Investigator Sponsored Research Award from GE Healthcare.
Disclosure: No disclosures to report.
Conflict of Interest: The authors have no conflict of interest to declare.
References
Journal Information
Journal ID (nlmta): tom
Journal ID (publisherid): TOMOG
Title: Tomography
Subtitle: A Journal for Imaging Research
Abbreviated Title: Tomog.
ISSN (print): 23791381
ISSN (electronic): 2379139X
Publisher: Grapho Publications, LLC (Ann Abor, Michigan)
Article Information
Self URI: media/vol6/issue4/pdf/tomo06362.pdf
Copyright statement: © 2020 The Authors. Published by Grapho Publications, LLC
Copyright: 2020, Grapho Publications, LLC
License (openaccess, http://creativecommons.org/licenses/byncnd/4.0/):
This is an open access article under the CC BYNCND license (http://creativecommons.org/licenses/byncnd/4.0/).
Publication date (print): December 2020
Volume: 6
Issue: 4
Pages: 362372
Publisher ID: TOMO.2020.00045
DOI: 10.18383/j.tom.2020.00045
PDF
Download the article PDF (2.92 MB)
Download the full issue PDF (5.19 MB)
Mobileready Flipbook
View the full issue as a flipbook (Desktop and Mobileready)