Development and evaluation of convergent and accelerated penalized SPECT image reconstruction methods for improved dose-volume histogram estimation in radiopharmaceutical therapy.

PURPOSE
Three-dimensional (3D) dosimetry has the potential to provide better prediction of response of normal tissues and tumors and is based on 3D estimates of the activity distribution in the patient obtained from emission tomography. Dose-volume histograms (DVHs) are an important summary measure of 3D dosimetry and a widely used tool for treatment planning in radiation therapy. Accurate estimates of the radioactivity distribution in space and time are desirable for accurate 3D dosimetry. The purpose of this work was to develop and demonstrate the potential of penalized SPECT image reconstruction methods to improve DVHs estimates obtained from 3D dosimetry methods.


METHODS
The authors developed penalized image reconstruction methods, using maximum a posteriori (MAP) formalism, which intrinsically incorporate regularization in order to control noise and, unlike linear filters, are designed to retain sharp edges. Two priors were studied: one is a 3D hyperbolic prior, termed single-time MAP (STMAP), and the second is a 4D hyperbolic prior, termed cross-time MAP (CTMAP), using both the spatial and temporal information to control noise. The CTMAP method assumed perfect registration between the estimated activity distributions and projection datasets from the different time points. Accelerated and convergent algorithms were derived and implemented. A modified NURBS-based cardiac-torso phantom with a multicompartment kidney model and organ activities and parameters derived from clinical studies were used in a Monte Carlo simulation study to evaluate the methods. Cumulative dose-rate volume histograms (CDRVHs) and cumulative DVHs (CDVHs) obtained from the phantom and from SPECT images reconstructed with both the penalized algorithms and OS-EM were calculated and compared both qualitatively and quantitatively. The STMAP method was applied to patient data and CDRVHs obtained with STMAP and OS-EM were compared qualitatively.


RESULTS
The results showed that the penalized algorithms substantially improved the CDRVH and CDVH estimates for large organs such as the liver compared to optimally postfiltered OS-EM. For example, the mean squared errors (MSEs) of the CDRVHs for the liver at 5 h postinjection obtained with CTMAP and STMAP were about 15% and 17%, respectively, of the MSEs obtained with optimally filtered OS-EM. For the CDVH estimates, the MSEs obtained with CTMAP and STMAP were about 16% and 19%, respectively, of the MSEs from OS-EM. For the kidneys and renal cortices, larger residual errors were observed for all algorithms, likely due to partial volume effects. The STMAP method showed promising qualitative results when applied to patient data.


CONCLUSIONS
Penalized image reconstruction methods were developed and evaluated through a simulation study. The study showed that the MAP algorithms substantially improved CDVH estimates for large organs such as the liver compared to optimally postfiltered OS-EM reconstructions. For small organs with fine structural detail such as the kidneys, a large residual error was observed for both MAP algorithms and OS-EM. While CTMAP provided marginally better MSEs than STMAP, given the extra effort needed to handle misregistration of images at different time points in the algorithm and the potential impact of residual misregistration, 3D regularization methods, such as that used in STMAP, appear to be a more practical choice.


INTRODUCTION
Radiopharmaceutical therapy (RPT) is a promising treatment method for certain types of cancer such as radioiodine therapy for thyroid cancer, radioimmunotherapy (RIT) for non-Hodgkin's lymphoma (NHL), 1 and peptide receptor radiation therapy (PRRT) for neuroendocrine tumors (NETs). 2 In RPT, the goal is to deliver lethal radiation absorbed doses to tumor cells without producing normal organ toxicities. Since it is hard to estimate the absorbed dose to all tumor cells, treatment to normal organ maximum tolerated doses (MTD) is likely the dosing method of choice. Advanced and accurate dosimetry methods to estimate the radiation absorbed dose to both tumors and normal tissues are thus important in RPT for improved treatment planning.
In RPT, estimating absorbed dose requires knowledge of the radioactivity distribution in the body as a function of time. Emission computed tomography (ECT) imaging modalities, including SPECT and PET, can provide estimates of the three dimensional (3D) activity distribution in the body at the time of imaging; repeating the imaging at multiple time points can provide information about the change in activity distribution over time. Since RPT agents have tended to have relatively long uptake phases in tumors, therapeutic radionuclides with physical decay half-lives of several days have been most appropriate. Most common positron-emitting radionuclides have short half-lives. As a result, radionuclides appropriate for SPECT imaging have received the majority of attention. However, conventional wisdom has been that SPECT imaging requires greater patient imaging time and planar imaging with shorter imaging times has been deemed an adequate alternative. As a result, planar imaging and quantification methods have often been used.
More recently, quantitative SPECT (QSPECT) imaging methods based on statistical iterative reconstruction combined with compensation for the effects that degrade quantitative accuracy, i.e., attenuation, scatter, and collimator-detector response (CDR) have been developed and evaluated for RPT dosimetry by several groups. [3][4][5][6] Most of these studies have focused on organ or tumor based dosimetry, where QSPECT has been demonstrated to improve quantification for total organ activity. QSPECT methods are thus appropriate for organ-based dosimetry, e.g., using MIRD methodology.
On the other hand, 3D RPT dosimetry using 3D activity distributions as inputs has received increased interest in recent years as it can account for the effects of nonuniform dose distributions observed in some tumors and organs. Threedimensional dosimetry also enables patient-specific radiobiological modeling, 7 which has the potential to provide improved predictions of organ and tumor response. Dose-volume histograms (DVHs) and isodose curves, as well as maximum, minimum, and dose percentile, which can be derived from DVHs, are some quantitative measures that summarize the 3D dose distribution information. 8 In addition, the DVH combined with appropriate radiobiological modeling and organ or tumor specific parameters can be used to calculate the tumor control probability (TCP) and normal tissue complication probability (NTCP), which are radiobiological metrics describing the response of tumors and likelihood of late normal organ toxicity. [9][10][11] Even with QSPECT methods, the activity distribution estimates obtained from SPECT imaging are corrupted by noise and partial volume effects (PVEs). These errors in the activity distributions propagate to the DVH estimates and any subsequent summary measures. In Ref. 12, DVHs derived from phantom and OS-EM reconstructed images were calculated for typical organs such as liver, spleen, lungs, and bone marrow. Compared to DVHs obtained from the true absorbed dose distributions in the phantom, large deviations were observed in the DVHs obtained with SPECT images. In a previous simulation study, 13 we found that DVHs estimated from QSPECT images for normal organs such as the liver and kidneys were considerably degraded, mainly as a result of noise, PVEs, and ringing artifacts due to CDR compensation.
SPECT image reconstruction is an ill-posed problem and thus regularization is needed and has often been used to improve the quality of SPECT images, especially when used for visual interpretation tasks. In a previous study, 13 we investigated the ability of simple regularization methods, i.e., selecting an optimal number of iterations and postreconstruction filtering, to improve DVH estimates obtained from OS-EM reconstructed images. However, image reconstruction methods based on penalized maximum-likelihood (PML) or maximum a posteriori (MAP) methods have been derived and provide more sophisticated tools for regularization. [14][15][16] Recent work by Dewaraja et al. 17 showed the potential for using these regularized reconstruction methods to improve activity estimates for RPT applications. In that work, they investigated unfiltered OS-EM and PML reconstruction methods using quadratic and edge-preserving (Huber) penalties and incorporated CT side information into the prior weights. They evaluated the methods in the context of tumor activity quantification in terms of the bias and standard deviation of the total target activity estimate and the root mean square errors (MSEs) of the activity distribution within the target. They found that the anatomic prior produced the best results for targets with uniform activity.
As discussed above, 3D dosimetry allows incorporation of advanced radiobiological modeling into predictions of response. DVH estimates, instead of organ or tumor activity, are more directly relevant to metrics such as the NTCP and TCP. It is not clear whether the optimization criterion for activity inside volumes-of-interest (VOIs) is adequate for optimizing reconstruction in terms of the DVH. For example, in previous work, we found that optimizing for total activity inside organs favored large numbers of iterations and no regularization, 5 while producing the best DVH required a much smaller number of iterations and often postreconstruction filtering. 13 In this work, we developed and evaluated advanced regularized image reconstruction methods based on PML and investigated the efficacy of penalized image reconstruction in obtaining improved DVH estimates. In particular, the algorithms investigated were based on a multitracer hyperbolic prior we developed previously. 18 In this work, we extended that prior for 4D reconstruction of SPECT images from multiple time points to investigate whether simultaneous use of information from multiple time points could improve DVH estimates. Surrogate functions were derived for the prior, and complete data OS-EM (COSEM) and optimization transfer techniques were used to implement a subset-based accelerated and convergent algorithm. 19,20 This 4D reconstruction method was evaluated in comparison to postfiltered OS-EM and PML reconstruction using a 3D hyperbolic prior. The evaluation was performed in terms of reliability (as measured by the mean square error) of organ DVHs using optimized values of the various hyperparameters. MAP methods have previously been studied for lesion detection and quantitation. [21][22][23] However, to the best of our knowledge, this is the first attempt to apply and evaluate the potential of PML based image reconstruction methods in terms of their effect on dose-volume histogram estimates in the context of radiopharmaceutical dosimetry.
This evaluation was performed using a simulation study. A modified NURBS-based cardiac-torso (NCAT) phantom and previously validated Monte Carlo (MC) simulation methods were used to generate projection data. A Monte Carlo generated voxel S kernel was used to compute dose-rate distributions from activity distribution estimates. The reconstruction methods were optimized in terms of lowest mean square error in cumulative dose-rate volume histogram (CDRVH) estimates. Cumulative DVHs (CDVHs) obtained from the phantom and reconstructed images with optimal parameters were calculated, plotted, and compared both qualitatively and quantitatively.
Implementing the cross-time MAP (CTMAP) method clinically would require registering the SPECT/CT images obtained at the different time points and using these registration parameters during evaluation of the 4D prior performed in the reconstruction. The regularization in the CTMAP method is thus sensitive to registration errors, but this sensitivity has not been evaluated in the present work.
Both the single-time MAP (STMAP) and OS-EM methods were applied to patient data. Reconstructed images and CDRVHs qualities are similar to those observed in the simulation study.

2.A. NCAT phantom and Monte Carlo simulation
A modified version of the 3D NCAT phantom that incorporated a more detailed kidney model was used in this work to provide a realistic and flexible model of human anatomy. 24 Each kidney in the new model was comprised of a renal cortex, medulla, and pelvis. This enabled realistic modeling of the nonuniform activity distribution in the kidneys, which is especially relevant in PRRT. 25 The anatomical parameters used in the NCAT phantom for major organs were adjusted to match those of a patient in a high-dose 90 Y Zevalin trial. 26 The activities and kinetics for all organs except the kidneys were from the same patient and measured after an 111 In Zevalin injection, which served as the imaging surrogate for 90 Y Zevalin therapy. The activity concentrations for the kidney, including the cortex, medulla, and pelvis, were from a different patient imaged after administration of 111 In-DTPAoctreotide in a PRRT dosimetry study. 27 The renal medulla and pelvis were assumed to have the same activity concentrations. In RIT with Zevalin, the liver is the major doselimiting organ. By contrast, in PRRT, a substantial fraction of the radiolabeled peptides are excreted by and retained in the kidneys, especially in the renal cortices. We combined organ uptakes from two different patients undergoing different therapies to enable us to investigate dosimetry of both the liver and kidneys in a single simulation study. There is some cross talk between the organs because of partial volume effects, so the results reported are not identical to what would have been observed independently. Nevertheless, this approach was computationally expedient and enabled study of the effects of the parameters on the DVH of each organ.
We simulated time activity curves at five time points, i.e., 1, 5, 24, 72, and 144 h postadministration, assuming monoexponential decay that modeled only the biological clearance phase. The uptake phase was ignored because the biological clearance phase tends to dominate the dose. For each organ, the activity at the 1 h time point, volume, and the half-life are listed in Table I.
The projection data were simulated using a modified version of the SimSET code 28 in combination with an angular response function (ARF) based simulation of the interactions in the collimator and detector. 5,29 The ARF models the distance and rotational variation of the collimator point response function. Parameters appropriate for a GE Millennium VG/Hawkeye SPECT/CT system were used to perform the Monte Carlo simulation. We simulated a 15.875 mm thick crystal, an energy resolution of 9.5% at 140.5 keV with inverse square root dependence on energy, and a medium energy general purpose (MEGP) collimator. Both 111 In photopeaks (171 and 245 keV) were separately simulated for all organs and regions listed in Table I. Low noise projection data were generated in 128 transaxial and 170 axial projection bins at 120 views over 360 • using a 0.442 cm projection bin size. The counts in two 20% wide energy windows centered at 171 and 245 keV were summed with the appropriate abundance (0.90 T I. Organ volumes, activities, and half-lives for the simulated NCAT phantom.

Heart
Lung Liver  for 171 keV and 0.94 for 245 keV). In order to reduce the effects of discretization, the NCAT phantom was discretized with a voxel size of 0.221 cm, half of the projection bin size, for the simulation. The projection data simulated from all organs were scaled and summed based on the activity distribution in Table I. The low noise projection data were scaled to a clinical count level appropriate for SPECT with a 30 min acquisition time. The average number of counts per transaxial slice at the 24 h time point was 1.355 × 10 5 . Finally, the projection data were filtered with a 2D Gaussian filter having a FWHM of 0.40 cm to simulate the detector intrinsic spatial resolution. Poisson noise was simulated to produce projections with realistic noise levels. Figure 1 shows one slice of the phantom, attenuation map, and the simulated noisy projection data. Figure 2 shows a slice of the detailed kidney model.

2.B. Image reconstruction method
Previously, we derived and validated a regularized image reconstruction method for simultaneous dual-isotope myocardial perfusion SPECT imaging using a hyperbolic prior. 18 The prior makes use of the automatic spatial registration of the images of the two radionuclides guaranteed by simultaneous dual-isotope acquisition. The image reconstruction was formulated as a joint estimation problem. In this work, F. 2. One coronal slice through the nonuniform kidney model. we extend the prior to SPECT image reconstruction in RPT dosimetry using the fact that projection data for the same patient are usually acquired at several time points for dosimetry purposes. We assumed that the images from the different time points were spatially registered. One goal of this investigation was to see the potential benefit that could be obtained if this perfect registration was achieved, though achieving this in practice would be quite difficult. This 4D reconstruction method will be referred to as the CTMAP. CTMAP applies a 4D prior that penalizes differences in time activity curves among neighboring voxels; by this means it incorporates temporal information into the reconstruction process.
We also investigated a 3D version of the algorithm, referred to as STMAP, which reconstructs the image for each time point independently without using the cross-time penalty, and can thus be applied even if the images at multiple time-points are not registered. Comparing results between CTMAP and STMAP shows the potential benefit that could be obtained if the 4D data were perfectly registered prior to or during the reconstruction.
The objective function for the two regularized methods is defined as follows: where X t is the image vector for the tth time point, − → X = (X 1 ,X 2 ,...,X T ) is the set of image vectors for the T time points, and β is the prior weight. In this study, T was 5. The first term is the negative sum of the log likelihood at the various time points and the second term is the prior defined as The potential function, ψ, has the general form In Eq. (2), − → X j is the vector of image values at voxel j across all time points, J is the number of voxels in each image, N j is the neighborhood of voxel j, the constants {w j k } are non-negative weights applied to differences between voxels j and k. In Eqs. (3) and (4), x j, t and x k, t are the values of voxels j and k at the tth time point. Eq. (3) defines the CTMAP prior as the square root of the difference in time activity curves between neighboring voxels j and k. The variables c and δ t are adjustable parameters; when c = 0 and δ t = 1 for all time points t, Eq. (3) reduces to the Euclidean distance between the two vectors representing the discrete time activity curves. Eq. (4) defines the STMAP prior as the sum of the square root of the difference between neighboring voxels at each time point. In this work, we fixed c to be 1.
The hyperbolic prior is convex and edge preserving. It is a compromise between quadratic and absolute value (i.e., total variation) priors. Similar forms of hyperbolic priors have been used in CT image reconstruction, 30-32 breast tomosynthesis, 33 and CT image restoration. 34 We redefine the prior below as it only depends on the difference between neighboring voxels. Let where In Ref. 18, we derived a separable surrogate function for the hyperbolic prior and implemented a convergent and accelerated algorithm using the COSEM framework 19 and optimization transfer technique. 20 We followed a similar procedure and derived the update equations for CTMAP; the update equations for STMAP are similar and its derivation is already given in Ref. 18. The first step of the optimization is to update the complete data 19 where the superscript (n,l) denotes the nth iteration and lth subset. In Eq. (7), C i j, t is analogous to the complete data concept used in conventional ML-EM, and represents the number of photons from voxel j detected in bin i at time point t; 19 p i, t is the projection data at bin i and time point t; P i j is the probability that a photon from voxel j is recorded in projection bin i; and S l represents the set of projections in subset l.
The second step is to update the activity estimates for all time points simultaneously. This is done using where In this work, images were reconstructed using both CTMAP and STMAP as well as OS-EM. Both the MAP algorithms and OS-EM were implemented with model-based compensations for attenuation, scatter, and the collimator-detector response function (CDRF). 5,35,36 The abundance-weighted average of the true attenuation maps for the 171 and 245 keV photons emitted by 111 In was used in the reconstruction. Scatter compensation was performed using the effective source scatter estimation (ESSE) method. 37 The set of full CDRFs, including interactions in and penetration through the collimator and detector, was estimated by Monte Carlo simulation of a point source at various distances from the surface of the collimator. 38 In all reconstructions, we used 24 subsets (five views per subset) per iteration. The objective function includes the hyperparameters β and a set of δ values (one for each time point). We used a grid search method to find the optimal β and investigated four sets of δ values. A total of seventy β values ranging from 10 −7 to 2.0 were used. The four sets of δ used are given in Table II. For the first two sets of δ values, the differences at five time points had equal weights, i.e., 0.1 and 1. The third set of δ values was chosen to take into account that the change in voxel values over time due to the radioactive decay of 111 In, i.e., the weights for activity differences were scaled in inverse proportion to the reduction in activity due to decay at those time points. The last set of δ values simply used a larger δ at each time point than that used in the third set to approximately model the effects of biological decay, which would result in a faster decrease in activity than radioactive decay alone. For the convergent MAP algorithms, we used images obtained with ten iterations of OS-EM as initialization and iterated to effective convergence, i.e., when we observed only small changes in the dose-rate volume histogram obtained from reconstructed images. For images reconstructed using OS-EM, a postreconstruction Gaussian filter was applied to the image obtained after 1-150 iterations. We investigated 15 FWHM values ranging from 0.5 to 5.0 pixels in increments of 0.5 pixels and from 5.0 to 10.0 pixels in increments of 1 pixel. Images from each combination of iteration number and FWHM were used to produce dose-rate distribution estimates.

2.C. Dose rate and absorbed dose estimate method
The methods were evaluated in terms of their effects on dose rate and dose for two reasons. First, these are the quantities that are important for the target application of therapy planning. Second, since dose-rate estimation is essentially a smoothing of the activity estimate, in principle, the optimal parameters and performance of the method could be different for activity estimation rather than dose-rate estimation.
We estimated dose rate and then absorbed dose images from the reconstructed activity estimates obtained from both the MAP and OS-EM methods using the following procedure. First, the reconstructed images were interpolated from a voxel size of 0.442 cm to 0.221 cm using trilinear interpolation in order to be consistent with the phantom. The voxel values in the SPECT images were then converted to activities by scaling with a calibration factor that was calculated from the geometric sensitivity of the imaging system. The geometric sensitivity is defined as the number of detected geometric photons per emitted photon and was calculated by Monte Carlo simulation of a point source in air using SimSET and an ARF simulation that included only geometric photons.
The focus of the study was on 90 Y dosimetry. Thus the activity estimates from 111 In imaging were converted to 90 Y activity estimates by compensating for the difference in physical decay half-lives.
The dose-rate distributions were estimated from the 90 Y distribution by convolution with a voxel S kernel. The voxel S kernel was calculated from MC simulations using the 3D-RD software package. 39 A water-filled sphere with a diameter of 33 × 0.221 cm was used. A voxel source was placed at the center of the sphere and 90 Y beta emission from it was simulated. We recorded the energy deposited in the voxels inside the sphere. A total of 10 × 10 6 particle histories were used in the simulation. The voxel size, 0.221 cm, was the same as the NCAT phantom and the interpolated activity distribution images. Thus convolution with the voxel S kernel gave the dose-rate distribution. This voxel S kernel convolution technique assumes a homogeneous medium in the body and is less accurate for non-soft-tissues and at the boundaries of these tissues. The importance of this limitation is reduced in our study because we focused on dose-limiting organs such as the liver and kidneys. This approach might have led to some errors at the upper boundary of the liver where it is close to the lung. However, only a very small proportion of the liver volume was thus affected, and this would have a relatively minor effect on the DVH for the entire liver. In addition, the same absorbed dose estimation method was used for all reconstruction methods compared. Thus differences between reconstruction methods due to this source of error were expected to be small.
The absorbed dose images were estimated from the set of dose-rate images as follows. The trapezoidal rule was applied to integrate the dose-rate data in each voxel for the intervals between the five simulated time points. Exponential fits were adopted to calculate the dose estimates before the first and after the last time points. In particular, exponential fits of the data at the five time points were computed and the resulting half-life was compared to the half-life of physical decay. Since a half-life longer than physical decay indicates a failure of the exponential fitting due to noise, the smaller of the fitted or physical half-life was used. This half-life was used to extrapolate the activity from time zero to the first time point and the last time point to infinity. The resulting exponentials were integrated over these times and added to the results of the trapezoidal integration between the five time points. A similar computational scheme was used in Ref. 40.

2.D. Histogram generation and figure of merit (FOM)
The CDRVHs for the liver, kidney, and the renal cortices were computed and plotted from the dose-rate images at each time point. In order to make the histograms more relevant to the final estimated dose, we used a bin size equivalent to 0.25% of the median value of dose/activity ratio obtained for 90 Y Zevalin for the liver and 0.2% of a literature value of the dose/activity ratio for 90 Y-DOTATOC for the kidneys. Specifically, the dose/activity ratios used were 4.8 mGy/MBq for the liver 41 and 3.3 mGy/MBq for the kidneys. 42 The CDRVHs were computed for each combination of β and δ for the MAP algorithms and each iteration number and FWHM of the Gaussian filter for OS-EM. The CDRVHs from the phantom were also computed and served as the ground truth. The MSE of the CDRVH estimate served as the FOM in this study and was computed for each histogram using the following formula: where h r b is the histogram value in the bth bin obtained from the reconstructed images, h p b is the histogram value in the bth bin for the phantom, and B is the total number of bins in the histogram.
The dose-rate images that gave the minimum MSE of the CDRVH at each of the five time points were integrated to obtain the optimal absorbed dose distribution estimate. The optimal parameters for each time point were usually different. A discussion of the hyperparameters and their effects on MSE is given in Sec. 3.D. Tables III, IV, and V show the optimal parameters for the CDRVH at each time point for each organ with OS-EM, CTMAP, and STMAP, respectively. The CDVH was then computed from the optimal dose distribution estimate analogously to the CDVRH. This was performed for both the phantom and the estimates obtained from SPECT. The MSE of the CDVH for the reconstruction relative to the phantom was computed as described above for the CDRVH.

RESULTS FROM SIMULATION STUDY
In this section, we present the evaluation results from the NCAT simulation study.

3.A. Qualitative image quality
The primary goal of this paper was to investigate the effects of the reconstruction methods on the quality of the DVHs. However, the reconstructed images themselves shed some light on the effects of the regularization methods. Figure 3 shows images of the phantom and the reconstructed images obtained with CTMAP and OS-EM. The MAP images shown were obtained with various values of β and single value of δ. We can see that with increasing β the reconstructed images are smoother. Also, for later time points the images show increased noise with the same β due to the decreased activity. We also observed in the experiments that images tended to be smoother with a larger value of δ. Analysis of the effects of the methods and parameters on the DVH estimates are presented in Secs. 3.B and 3.C. More analysis of the effects of β and δ on the DVH estimates is given in Sec. 3.D. Figure 4 shows the CDRVH of the liver at one time point, in this case 5 h, obtained from the phantom and optimal reconstructions; the MSEs with the optimal parameters are shown in Table VI. Both CTMAP and STMAP substantially improved the shape of the estimated histogram both above and below the median dose rate (corresponding to the dose T IV. The optimal parameters with CTMAP for each organ at five time points. i.e., optimal OS-EM. The results for CTMAP and STMAP were close both qualitatively and quantitatively, although the MSE of the CDRVH from CTMAP was about 13% less than that from STMAP. Similar results were observed for the CDVH of the liver, as shown in Fig. 5 Figures 6 and 7 show plots of the CDVH for the kidneys and renal cortices, respectively, obtained from the phantom and reconstructed images obtained using CTMAP, STMAP, and OS-EM. Differences in the shape of CDVHs between MAP and OS-EM are evident, but not large. Table VIII lists the MSEs obtained with the parameters optimized for the kidneys and renal cortices. Note that the MSEs for the kidneys and renal cortices from the three reconstruction methods were very close.

3.D. The selection of hyperparameters
The MAP algorithms use two hyperparameters, β and δ. In this study, we used seventy values of β and four sets of δ T V. The optimal parameters with STMAP for each organ at five time points. values, as described above. Figures 8 and 9 show the MSE of the CDRVH as a function of β for the various sets of δ values for the liver and kidneys at 72 h. These graphs illustrate the effects of these parameters on the quality of the estimated CDRVH. Note that the optimal combination of β and δ was different for each time point and organ. For instance, a larger β value was required for the liver than for the kidneys in order to achieve minimal MSE. Larger β and δ values result in stronger smoothing of noise and also in reduced resolution (and hence larger PVEs). Because of the coupled effect of δ and β, the optimal β was different for each set of δ values. For example, with a larger δ, the optimal β required was often smaller. In addition, with smaller δ values (e.g., set F. 4. The CDRVH of the liver at 5 h from phantom and SPECT images using CTMAP, STMAP and OS-EM with optimal parameters. 1 where all δ values were 0.1) the MSE changed less over a larger range of β. Also notice that CTMAP was less sensitive to changes in β than STMAP, possibly because of the stabilizing effects of the information at other time points. We did not find that any of the four sets of δ values gave consistently better results. However, since the first set of δ gave near-optimal MSE for a wider range of β, it might be preferred for practical reasons. Choosing a different optimal combination of parameters for each time point might be impractical. For example, in our experiments using CTMAP with five time points requires specification of the values of 30 hyperparameters: one β and five δ values for each of the five time points. However, in the F. 5. The CDVH of the liver from phantom and SPECT images using CTMAP, STMAP and OS-EM with optimal parameters. above results, we noticed that the MSE was not as sensitive to β for the first and second sets of δ values investigated. This was especially true for the liver, as can be seen from Fig. 8. For instance, the MSE was relatively flat for β in the range [0.2, 0.5] for both CTMAP and STMAP and for all five time points when using the second δ set, i.e., all ones. In fact, the first and second sets were very simple parameter choices and reduced parameter selection to a single value, i.e., β. Further, there was a range of β values that produced near-optimal results. To demonstrate this, we used β = 0.4 and β = 0.5 for all five time points with all δ values set to 1 for all time points. The resulting histograms and MSE results for the liver are shown in Fig. 10 and Table IX. Note the results from both MAP algorithms were still substantially better than the optimal OS-EM results. In fact, the MSEs for CTMAP and STMAP were about 17% and 23%, respectively, of the MSE for OS-EM. These results indicate near-optimal results can be obtained with a much simpler selection of hyperparameters.

APPLICATION TO PATIENT DATA
One practical advantage of STMAP is that it does not require that the SPECT data from multiple times be spatially registered. Thus, to demonstrate qualitatively the difference in images and dose-volume histograms from the edge preserving prior, we applied the STMAP method to a set of patient data. Results from one patient are presented in this section. The SPECT/CT data were acquired at 24 h postinjection of 111 In Zevalin and were part of a previous clinical trial. Approval was obtained from the Johns Hopkins Institutional Review Board to use these data for the evaluation of reconstruction algorithms. Unlike the simulation where we combined the activity distribution from one patient with Zevalin therapy and the other patient from a PRRT dosimetry study using 111 In-DTPA-octreotide, the patient data represent only the Zevalin distribution. The scan was acquired on a GE Millenium VG/Hawkeye SPECT/CT system with a 15.875 mm thick crystal and a MEGP collimator. The SPECT projection data were acquired into 128 × 128 projection matrices with a 0.442 cm bin size at 120 views over 360 • . The x-ray CT data were acquired with the SPECT pixel size in-plane and a 1 cm axial slice thickness. The registered CT and SPECT images were used to manually define organ VOIs. A detailed T VII. The minimal MSEs (×10 −3 ) of the liver CDVH using optimal parameters. description of the therapeutic and imaging protocol is given in Ref. 43. Both the STMAP and OS-EM methods were used to reconstruct SPECT images with compensations for attenuation, scatter, and the CDRF. Dose-rate images were calculated using the same procedure as in the simulation. Reconstructed SPECT images from the two methods are shown in Fig. 11 for several sets of parameters. The optimal parameters for STMAP for the simulated liver CDRVH at 24 h were β = 0.1 and δ = 10. The corresponding optimal parameters for OS-EM were 3 iterations and a FWHM of 1.0 pixels. The CDRVHs of the liver, the dose-limiting organ in myeloablative Zevalin therapy, are shown with those optimal parameters in Fig. 12. In Sec. 3.D, we showed that near-optimal results with CTMAP and STMAP can be obtained with a much simpler selection of hyperparameters such as 0.5 for β and 1 for all δ. So, the CDRVH of patient liver obtained with STMAP is also shown with β = 0.5 and δ = 1. In addition, a liver CDRVH obtained with OS-EM using 2 iterations and a FWHM of 4.0 pixels is shown. For comparison, Fig. 13 shows the simulation results for CDRVHs of the liver at 24 h using the same sets of parameters. Note that the dose-rate images used to obtain the CDRVHs shown here were computed from image counts and were not converted to activity. This does not affect the shape of the histogram. Also note the count level for the patient at the abdomen region was approximately 20% lower than for the simulated phantom. This results in a shift of the patient histograms on the horizontal axis compared to those in simulation.

CTMAP
From the simulation results in Fig. 13, we observed that STMAP with the two sets of parameters produced a very similar CDRVH for the liver. In the simulation, OS-EM with 2 iterations and an FWHM of 4 pixels showed a reduction of noise at the upper end of CDRVH, but larger deviations from the phantom CDRVH compared to that obtained for OS-EM with 3 iterations and an FWHM of 1.5 pixels. The true activity distribution for the patient data is not available so only a qualitative comparison can be made. However, Fig. 12 illustrates that the CDRVHs from the patient with STMAP and OS-EM had similar trends compared to those in the simulation. This indicates that the advantages of STMAP in terms of CDRVH can potentially be realized in patients.

DISCUSSION
The DVH is a valuable tool for treatment planning and TCP and NTCP computation. Accurate DVH estimation depends on accurate dose distribution estimates, while dose distribution estimates are based on activity distribution estimates obtained from ECT imaging. Dose distribution estimates are degraded by noise and limited spatial resolution. In previous studies, 13 we found that noise, PVEs, and ringing artifacts due to collimator-detector response compensation are the major factors that contribute to the errors in dose histograms obtained from OS-EM reconstructed images. In this work, we developed and evaluated a convergent and accelerated 4D MAP algorithm using a hyperbolic prior, CTMAP, and compared it to results from postfiltered OS-EM F. 6. The CDVH of kidneys from phantom and SPECT images using CTMAP, STMAP, and OS-EM with optimal parameters. and a previously developed 3D MAP method, STMAP, also using a hyperbolic prior. The results from the simulation study showed that both MAP algorithms produced qualitatively and quantitatively superior CDVH estimates compared to OS-EM for the liver, which is the major dose-limiting organ in myeloablative Zevalin RIT. For the kidneys and renal cortices, which are dose-limiting organs in PRRT therapy, both MAP algorithms and OS-EM gave similar results and large residual errors in the CDVHs were still evident.
As pointed out in Ref. 13, for large organs such as the liver, noise, and ringing artifacts were the major factors degrading the histograms; for small organs with detailed features, such as the kidney and renal cortex, PVEs were the dominant factor degrading the histograms. The liver simu-F. 7. The CDVH of renal cortices from phantom and SPECT images using CTMAP, STMAP, and OS-EM with optimal parameters. lated in this study is a large organ and noise and ringing artifacts are the major factors limiting the accuracy of the DVH for the count levels studied. The edge-preserving hyperbolic priors, using either 3D or 4D information, were effective in reducing noise and ringing artifacts without increasing PVEs, resulting in substantial improvements in the CDVH compared to OS-EM with optimal filtering. However, these priors, while potentially producing sharper edges, did not reverse the degradations in the CDVH due to partial volume effects for small organs such as the kidneys. In fact, a larger β value was often optimal for the liver in order to reduce noise while a smaller β was optimal for the kidneys and renal cortices in order to preserve resolution. This indicates that less prior weight is optimal for small structures with fine details. it might be helpful to employ a prior that uses anatomical information to reduce PVEs, e.g., Refs. 17 and 44-47.
The CTMAP algorithm uses temporal information to try to improve the histogram estimates. However, this requires registration of images across multiple time points. In our simulation experiments, CTMAP marginally improved the CDVH estimates compared to STMAP: the MSE in the CDVH with CTMAP was about 16% less for the liver and 9% less for kidneys and renal cortices than that with STMAP. Further, there was little visible difference in the shape of the CDVH. A full implementation of CTMAP with clinical data would require registering the reconstructed images at each time point and using this registration information during calculation of the MAP prior inside the reconstruction process. Since registration is likely to be imperfect, residual misregistration errors would impact the computation of the MAP prior, and thus the improvement due to the 4D prior would be expected to be smaller in practice than observed in the simulation. As a result, for applications such as RPT where image registration may be difficult to achieve, STMAP might be preferred in practice. It should be noted, however, that calculation of dose-volume histograms from SPECT images at multiple time points also requires image registration. Residual misregistration thus degrades dose-volume histograms obtained even from 3D reconstruction methods. For other dynamic imaging applications, where the data at different time points are acquired over a short period of time without moving the patient from the imaging table, the proposed 4D method might prove practical and worthwhile. The selection of hyperparameters is an important issue for MAP algorithms in general, and the situation is no different for the algorithms used in this work. From Figs. 8 and 9, we observe that the values of both β and δ have a substantial effect on the fidelity of the histogram. However, we observed that the MSE was less sensitive to the β value when all δ values for the different time points were the same. Furthermore, we noticed that there was a range of β values that gave near-optimal results for all time points and that near-optimal MSEs could be obtained using a single β and set of δ values for all time points. These results suggest that the hyperparameter selection problem may be no more difficult than selecting an optimal iteration and filter for OS-EM and that the proposed MAP algorithms have the potential for practical use in RPT dosimetry. One issue that has not been investigated is whether the optimal hyperparameters are patient dependent. A comprehensive study of these effects on MAP algorithms using a population of phantoms, e.g., as in Ref. 26, is beyond the scope of this work, but is an important topic for future study. As a preliminary investigation, we used the same phantom and increased the activity in the kidneys by a factor of 3. We found that a smaller β value was needed to obtain the minimal MSE for the kidney CDRVH. For example, the optimal β value for the kidneys at 1 h decreased from 0.03 to 0.007 with the same set of δ values. This is because there was less noise in the kidneys, thus favoring less smoothing to avoid increasing PVEs. This indicates that the optimal hyperparameter values are count-dependent, though the amount of change for smaller changes in organ activities or changes in organ shapes remain unclear.
A practical disadvantage of the current implementation of the proposed MAP reconstruction methods is the large computational cost required. In this work, we used the COSEM framework to implement the MAP methods as convergent and accelerated algorithms. As observed in this work and elsewhere in the literature, the convergence of the MAP algorithms using COSEM is slow. 19 This computational cost may not be acceptable in clinical practice. On the other hand, the improvements in histogram estimates obtained with MAP algorithms were gained not as a result of the iterative algorithms, but due to the prior's property of noise control and edge-preservation. Therefore, faster implementations of MAP algorithms, 48 either convergent such as modified block sequential regularized EM (BSREM) 49 or nonconvergent, such as an ordered-subsets version of De Pierro's modified EM algorithm, 50 might prove more useful in practice.
There are several factors that have not been investigated in this paper and are important for practical applications of the proposed methods to 3D dosimetry. These include the effects of errors in image registration and VOI misdefinition on the accuracy of the dose-volume histograms. Studies with multiple noise realizations would be helpful to further understand the precision of the dose-volume histograms. Demonstration that registration of data acquired at different time points at the voxel level is feasible is also essential. However, the goal of this study was to demonstrate the potential of the proposed 3D and 4D MAP reconstruction methods and these other factors are beyond the focus of the present study. We have studied these factors in the context of organ activity estimation and dosimetry, 6,26,43,51,52 and anticipate doing so for voxel dosimetry in the future.
The major purpose of this paper was to investigate the ability of more sophisticated regularization methods to reduce errors in dose-volume histograms. The purpose was not to provide a comprehensive evaluation of the errors in DVH that would be inherent in those from human studies. While every effort was made to make the simulation as realistic as possible, there are a number of effects that were not modeled. It is, however, of interest to estimate the errors that other effects would have on the DVHs obtained here.
The simulation does explicitly and realistically model the effects of reconstruction and compensation for image degrading factors, partial volume effects, Poisson noise, and integrating the time activity curve in voxels to obtain voxel dose values. Thus the reported results from simulation do include the impact of these effects in the errors observed in the DVHs.
The evaluation study, did use simulation not a real SPECT system. The simulation modeled the acquisition process very realistically using anatomically realistic phantoms and validated Monte Carlo simulations. For 111 In, we have previously reported errors in organ activities from simulation and phantom studies that were similar in magnitude, and expect the same for the DVH. The simulation did not model variations in performance over time, and it is possible that there is some variability due to acquisition. We are currently studying these effects by repeating phantom and sensitivity measurements over time. In that study, we have found that these variations are a relatively small effect with a coefficient of variation of less than 2%, and expect that there would be a similar error for the DVH. The study did not simulate the effects of deadtime. However, for the small pretherapy administered activities simulated in this study, the effects would be small, especially if the system calibration was performed using similar count rates. Finally, this paper did not include the effects of inaccuracy and imprecision of the calibration factor used to convert reconstructed voxel values to activity concentration. The fact that there was previously good agreement between simulation and experiment in terms of accuracy of organ activity estimates indicates that this factor is small. In addition, we are currently investigating these effects and a paper in preparation indicates that, when care is used in source preparation and measurement, errors in the calibration factor are less than 2%. These errors would not affect the shape of the DVH and would affect all of the reconstruction methods in the same way.
While the simulation does model partial volume effects for a single anatomy, we have previously reported that variations in anatomy and uptake result in a variation in errors over a population. We found that these errors were on the order of a few percent. These effects would primarily affect the precision of the DVH near the edges of objects, regions where the DVH accuracy is already the poorest. The effect would thus be the greatest for the kidneys and the magnitude should be bounded by the difference between the DVH for the phantom and the reconstructed image.

CONCLUSIONS
In this work, we developed and evaluated the potential of penalized image reconstruction methods for improved histogram estimates in RPT dosimetry. An edge-preserving hyperbolic prior was used, and both 3D and 4D convergent and accelerated algorithms were implemented. The results showed that MAP algorithms substantially improved the cumulative dose and dose-rate volume histogram estimates for large organs such as the liver compared to optimally postfiltered OS-EM reconstructions, reducing the mean square error in the histograms by a factor of about 5. For small organs with fine structural detail such as the kidneys and renal cortices, a large residual error was observed for both MAP algorithms and OS-EM, likely due to partial volume effects. These results suggest that incorporation of anatomical information into priors may be important for small objects. An important practical problem with MAP algorithms is the selection of hyperparameters. The results of this study suggest that a range of hyperparameters and the use of the same values for each time point can produce near-optimal results. However, methods for hyperparameter selection as a function of patient anatomy and count level remain an important topic for future research before MAP algorithms can be routinely practical for RPT dosimetry applications.
The 4D CTMAP method produced marginally better dosevolume histograms than the 3D STMAP method. However, it requires registered images at each time point for calculation of the 4D prior. Since residual misregistration errors would likely reduce the benefit of the 4D prior, it is unclear whether the relatively small benefit of the 4D prior would remain. Thus the 3D STMAP prior is likely preferred for clinical dosimetry applications. The 3D STMAP method, which does not use temporal information and thus does not require registration across time points, was applied to patient data at a single time point. The image quality and obtained DVHs were similar in character to those obtained in the simulations.

ACKNOWLEDGMENTS
This work was supported by the National Cancer Institute of the National Institutes of Health under award number R01-CA109234. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. The authors thank the editors and reviewers for their comments which have improved the presentation of the paper. Under a licensing agreement between GE Healthcare and the Johns Hopkins University, one of the authors (Eric Frey) is entitled to a share of royalty received by the University on sales of iterative reconstruction software. Portions of that software were used to obtain some results in this paper. The terms of this arrangement are being managed by the Johns Hopkins University in accordance with its conflict of interest policies.