Effect of magnetic resonance imaging pre-processing on the performance of model-based prostate tumor probability mapping

This study is based on 31 PCa patients from a single institution enrolled in a prospective, research ethics board approved clinical trial (NCT03378856 and NCT03367702). This research was approved by Comité d'éthique du CHUM, No 17.032 and 17.244. The research was conducted in accordance with the principles embodied in the Declaration of Helsinki, as well as in accordance with local statutory requirements. All participants gave written informed consent to participate in the study. Written informed consent was also given for publication by all participants.

Multi-parametric imaging consisting of a T2w (TR = 3377 ms, TE = 80 ms, FOV = 0.185 × 0.185 × 4.5 mm³), diffusion-weighted imaging (DWI) (TR = 3280 ms, TE = 79 ms, FOV = 1.774 × 1.774 × 4.5 mm³, b-values = 0–2000 s mm⁻²), and dynamic contrast-enhanced (DCE) (TR = 5.76 ms, TE = 0 ms, FOV = 0.917 × 0.917 × 3.0 mm³, dynamic scan time = 15.55 s, scan duration = 311 s, contrast agent = ProHance (gadoteridol, Bracco Imaging), dose = 0.2 cc kg = up to maximum of 20 cc, injection rate = 2 cc s⁻¹) compliant with PI-RADS v2.1 (Weinreb et al 2016, Turkbey et al 2019) was performed on a 3 T MRI scanner (Philips Ingenia) prior to radiotherapy.

A genitourinary radiation oncologist manually delineated contours of the prostate, peripheral zone (PZ), and tumors in both the PZ and transition zone (TZ) with a PI-RADS score of 4 or 5 via the aid of institutional PI-RADS reports. ADC maps were computed from three b-value DW images (b = 50, 600, 800 s mm⁻²) using the scanner proprietary software. K^trans maps were derived from the DCE-MRI images according to the generalized kinetic model (Tofts et al 1999). A population-based arterial input function (Parker et al 2006) was used to generate K^trans maps after adjusting the functional form to fit the temporal features of this dataset. In the absence of a pre-contrast T1 acquisition, a constant value of 1597 ms was used (Fennessy et al 2012). K^trans generation was performed using the PkModeling module in 3D Slicer (v 4.10.1) (http://slicer.org) (Fedorov et al 2012). K^trans images were normalized according to the median value in the PZ. The ADC images and K^trans images were resampled to a voxel size of 0.19 × 0.19 × 4.5 mm³ to match the resolution of the T2w images. ADC images were not resampled when performing multi-modal deformable registration to the T2w images. No motion correction was performed prior to ADC or K^trans image generation.

Zoom In Zoom Out Reset image size

Prior to the application of each pre-processing step, optimal parameter options were determined for each potential method. Each parameter option was initially set to the default or recommended value for a given method. A range of values were then tested for each parameter option individually while keeping all other parameters at the default or recommended settings. The performance of a particular parameter option was then evaluated by one or more of the relevant individual metrics for that pre-processing method. Once all parameter options were determined, a final implementation of each method was carried out using these optimal parameter values. The performance of this implementation was assessed by all relevant individual metrics, as well as the application-specific task of tumor localization, since the performance of a pre-processing tool has been shown to vary according to the individual metric used for evaluation (Viswanath et al 2011, Rohlfing 2012). Where applicable, results for the individual metrics are represented by violin plots, which are a combination of a box plot and a kernel density plot. As such, they display both summary statistics and the density of the variable. Like a box plot, the white dot depicts the median, and the thick gray bar shows the interquartile range. In addition, the thin gray line represents the distribution of the variable. Either side of the line is a kernel density estimation in which wider sections correspond to a higher probability of a given value, and narrow sections correspond to a lower probability of that value.

MR images are often susceptible to intensity non-uniformity artifacts. These artifacts may occur as a result of inhomogeneities in the radio frequency field, eddy currents, or attributes that are inherent to the particular biological sample being imaged. The effect of the bias field manifests as a nonlinear, smooth variation of signal intensity values across the image. The result is an image for which intensity values of a given tissue type may vary with location in the image. This spatial invariance has been shown to adversely affect image processing algorithms for performing tasks such as segmentation (Li et al 2014), registration (Xu et al 2013), and classification (Li et al 2009), among others.

The aim of bias field correction is thus to remove this intensity inhomogeneity, and numerous methods have been proposed. We chose to explore the more commonly used retrospective models, so intensity inhomogeneity in the T2w images was therefore corrected using two well-established, publicly available methods: N4 (v 2.2.0) (Tustison et al 2010) and PABIC (v 1.0) (Styner 2000). Optimal parameter settings were first determined for each method by testing a range of values for each available parameter option, including number of iterations, shrink factors, and order. The number of iterations refers to the maximum number of iterations that the algorithm is run and is employed as at least one of the stopping conditions for the algorithm. In general, a higher number of iterations can produce a more reliable result, though this will increase the runtime. A shrink factor refers to a parameter that is associated with a multi-resolution implementation of the algorithm; in this approach, the amount of data to process is reduced by initially running the algorithm on a coarse resolution of the data. In addition to data reduction, the number of local minimums/maximums that could potentially confound the performance of the algorithm is also reduced. Fine detail is then gradually added to the image by gradually increasing the image resolution at subsequent levels of the multi-resolution approach. The shrink factor is the degree to which the image is downsampled at each resolution level. The order parameter refers to the order of the b-spline in the case of N4 and to the order of the Legendre polynomials in the case of PABIC. Additional parameters that are specific to N4 consisted of the mesh grid resolution, number of histogram bins, Wiener noise, and full width at half maximum (FWHM). As for PABIC, the initial step size, as well as the mutation factors c_grow and c_shrink, were also optimized.

The performance of each parameter setting was evaluated individually using three commonly employed measures of the effects of bias field correction: the percent coefficient of variation (CoV), entropy, and the coefficient of joint variation (CJV). The CoV is expressed as a quotient between the standard deviation and mean of a given tissue class, t:

$$\begin{eqnarray*}&&{\rm{C}}{\rm{o}}{\rm{V}}=\displaystyle \frac{\sigma \left(t\right)}{\mu \left(t\right)},\end{eqnarray*}$$

where t represents the whole prostate gland. The CJV takes two tissue classes into account, which in our case consisted of the central gland (CG) and PZ:

$$\begin{eqnarray*}&&{\rm{C}}{\rm{J}}{\rm{V}}=\displaystyle \frac{\sigma \left({\rm{C}}G\right)+\sigma \left({\rm{P}}{\rm{Z}}\right)}{\left|\mu \left({\rm{C}}G\right)-\mu \left({\rm{P}}{\rm{Z}}\right)\right|}.\end{eqnarray*}$$

Each of these metrics should exhibit decreased values for bias field-corrected images. Moreover, the intensity inhomogeneity that occurs as a result of the bias field is thought to raise the entropy of the image data distribution. Entropy can be described as:

$$\begin{eqnarray*}&&{\rm{E}}{\rm{n}}{\rm{t}}{\rm{r}}{\rm{o}}{\rm{p}}{\rm{y}}=\displaystyle \sum _{i=1}^{G}\left[P\left({g}_{i}\right)\times \,\mathrm{log}\left[P\left({g}_{i}\right)\right]\right],\end{eqnarray*}$$

where G represents the total number of unique gray level intensities ${g}_{i}$ in the region of interest (ROI) and $P\left({g}_{i}\right)$ is the probability of occurrence associated with each ${g}_{i}.$ Thus, the entropy of bias field-corrected images should also decrease. The final implementation of each bias field correction method was then performed using a combination of all optimal parameter settings.

One well-recognized limitation of MR images is the lack of tissue-specific meaning in intensity values. This can lead to a drift in intensity values across both inter- and intra-patient acquisitions. Normalization aims to correct this drift by standardizing the intensity values across a dataset and is typically performed by one of three methods: centering the image intensity range according to its mean value, normalization of gray level values according to a specific ROI, and transformation of the histogram of the image to a reference intensity distribution. We refer to these methods as basic normalization, ROI normalization, and histogram matching normalization, respectively. Basic normalization refers here to z-score normalization. ROI normalization was performed by dividing all intensities by the median value in the PZ. Histogram matching normalization was implemented according to the method introduced by Nyú and Udupa (1999). This approach consists of two steps: (1) a training stage in which the landmarks of a standard histogram are determined from a training set of images and (2) transformation of each image histogram to the standard histogram. All bias field-corrected T2w images were included in the training step. A value of 4095 was chosen for s₂, and quartiles were chosen as histogram landmarks. The bias field-corrected T2w images were thus normalized according to each of these techniques. The performance of each normalization method was assessed with the CoV, which is routinely used to evaluate this type of pre-processing task.

Image distortion on DWI of the prostate due to chemical shift and susceptibility artifacts is a well-recognized phenomenon (Donato et al 2014, Hao et al 2017). Such distortions manifest as either local stretching or compression in the acquired image (Gallichan et al 2010). As a result, rigid registration is insufficient for the correction of local distortions on prostate DWI (Giannini et al 2015). Non-rigid registration, then, is often employed for more accurate alignment of T2w- and DW-images of the prostate (Buerger et al 2015). Thus, multi-modal deformable registration of the ADC image to the T2w image, representing the last step before a tumor map generation, was carried out using three popular, open source toolkits: elastix (Klein et al 2010), NiftyReg (Modat et al 2014), and ANTs (Avants et al 2011). In order to compare the performance of each toolkit, the same registration approach was followed. The registration of DW images is typically performed using non-rigid registration due to local deformations in these images that occur as a result of motion and susceptibility artifacts. In view of this, an intensity-based, non-rigid registration method that employs cubic b-splines to parameterize the deformation field was adopted as this has been shown to be a relatively simple, robust approach for both mono- and multi-modal implementations. The images were first aligned using an initial translation and then underwent a global transformation (rigid + affine). Non-rigid registration was then carried out by means of free-form deformation that utilizes a standard gradient-based optimization of cubic b-splines with mutual information as the similarity metric. Common parameter manipulations, including number of resolution levels, maximum number of iterations, interpolation order, number of samples, degree of smoothing, and regularization weight, were initially assessed using the Dice similarity coefficient (DSC) to measure the overlap of prostate gland delineations before performing the final registrations using optimal parameter values. The DSC is defined as follows:

$$\begin{eqnarray*}&&{\rm{D}}{\rm{S}}{\rm{C}}=\displaystyle \frac{2\left|{\rm{A}}\cap {\rm{B}}\right|}{\left|{\rm{A}}\right|+\left|{\rm{B}}\right|},\end{eqnarray*}$$

where A and B represent the segmentations of the T2w image and ADC image, respectively. The registrations were first estimated between T2w and b = 0 (b0) images, and the resulting transformation was then applied to the ADC images. While the fundamental structure of this approach is the same for each of the three toolkits, the algorithmic implementations differ, potentially resulting in disparities in registration outcomes. The quality of registration for each toolkit was evaluated by three different metrics: (1) a tissue overlap measure represented by the DSC, (2) a deformation field assessment represented by the determinant of the Jacobian, and (3) the mean target registration error (TRE). The calculation of the TRE was based on natural landmarks identified in both the T2w and ADC images, such as BPH nodules or calcifications. The TRE was computed as the distance between the geometric centers of corresponding landmarks in the T2w and ADC images after registration.

Once the pre-processing steps were evaluated, tumor probability maps were generated using a previously developed tumor localization model that was trained on mpMRI data acquired on both Siemens and Philips scanners (Dinh et al 2017). The aim of this model is to reduce the number of false positive results by combining mpMRI features with prior knowledge features. A total of 31 features are extracted per voxel, which include 29 mpMRI features (3 intensity, 24 texture, and 2 blobness), one from a population-based tumor probability atlas (Ou et al 2009), and one from a map constructed from the biopsy report for that patient. Processing of all mpMRI and prior knowledge data before feature extraction is summarized in figure 2. The mpMRI features were extracted from the T2w, ADC, and K^trans images according to the methods presented in Dinh et al (2017). Multi-parametric MRI-based features included both intensity and texture features. Texture features followed those proposed by Litjens et al (2014) and Vos et al (2012). Gaussian derivatives up to second order were extracted from the T2w images at exponentially increasing scales (σ = 1.5, 2.4, 3.8 and 6.0 mm). As for the ADC and K^trans images, normalized multiscale blobness features (Lindeberg 1998) were computed over the same four scales. The first element of prior information is a population-based tumor probability atlas presented by Ou et al (2009). The tumor probability atlas consists of voxels that correspond to the number of specimens that contained a tumor at that voxel. The tumor probability atlas was registered to the T2w images using the non-rigid registration algorithm implemented with elastix. The second element of prior knowledge is a biopsy map that was constructed using a detailed biopsy scheme in which the prostate was divided into six regions: apex, middle, base, and left/right. Each voxel value in the biopsy map was assigned a value of zero or one, depending on whether or not a positive biopsy was reported in that location (1 = positive). Once constructed, the biopsy map was smoothed by convolving it with a Gaussian kernel (σ = 4.5 mm).

Zoom In Zoom Out Reset image size

The tumor localization model was previously created and trained by Dinh et al (2017). Prostatectomy specimens obtained from each patient after mpMRI acquisition were used to establish ground truth tumor locations in a training cohort of patients. Slices of the prostatectomy specimens were stained with hematoxylin and eosin (H&E) and then registered to the T2w images. Each voxel of the training cohort was labeled as healthy or tumor tissue according to delineations on the H&E-stained slices by a pathologist. The aforementioned 31 features were then extracted for each voxel, and a logistic regression model was fitted to the histologic data obtained from the training cohort of patients according to Dinh et al (2017). The logistic regression function can be expressed as follows:

$$\begin{eqnarray*}&&F\left(x\right)=\displaystyle \frac{1}{1+{\,{\rm{e}}}^{-\left({\beta }_{0}+\,\displaystyle \sum {\beta }_{i}{x}_{i}\right)}},\end{eqnarray*}$$

where $F\left(x\right)$ is the probability that a voxel $x$ belongs to a tumor, ${x}_{i}$ is the value of feature $i,$ ${\beta }_{i}$ is the regression coefficient of feature $i,$ and ${\beta }_{0}$ is the offset. The model was validated in both a single-center and cross-center setting. Training and testing datasets for the single-center setting were established using a leave-one-out cross-validation approach. All patients were used to fit the model for the cross-center validation.

The classification model presented by Dinh et al (2017) was used here as presented in that study without any modifications to the model. The only difference is that instead of using non pre-processed imaging data, we are using pre-processed data as an input for testing the model in order to assess the effect of pre-processing on tumor classification accuracy. The pre-trained model was then used to generate tumor probability maps for four pre-processing variations: no pre-processing, bias field correction alone, bias field correction followed by normalization, and bias field correction and normalization followed by registration. The resultant tumor probability maps were evaluated by comparison with manually contoured tumors using the area under the curve (AUC). A Wilcoxon signed-rank test was performed with Bonferroni multiple comparison correction to determine whether tumor localization results differed significantly between various combinations of pre-processing steps (table 5). As four tests were performed in this study, the new significant p-value after Bonferroni correction is 1.25 × 10⁻².

To test the generalizability of this approach, the 31 patients were randomly divided into a training cohort of 20 patients and a testing cohort of 11 patients. We then repeated the entire analysis by first estimating the optimal parameter values for each pre-processing step on the training cohort and applying those pre-processing steps to the testing cohort. We then generated tumor probability maps on the testing cohort for each of the four pre-processing variations and evaluated the results using the AUC.

Clinically significant PCa has a characteristic appearance on mpMRI. In general, it is hypointense to surrounding tissue on T2w images and has an 'erased charcoal' appearance in the TZ. Likewise, since PCa is known to restrict diffusion, it also appears as hypointense on ADC images. The intensity values of PCa on DCE images, on the other hand, are expected to exhibit early enhancement, which manifests itself as hyperintense signal compared to healthy tissue. Image artifacts present on mpMRI can thus impact the expected appearance of PCa and ultimately affect any features that are extracted. Intensity inhomogeneities and intensity drift may lead to irregularities in intensity values that are not representative of the underlying tissue. This can affect both healthy and tumor tissue alike, resulting, for example, in healthy tissue that is more hypointense or tumor tissue that is more hyperintense than it should be. A tumor classification model that includes intensity value features may then be more prone to false positives and false negatives as a result if these artifacts are not mitigated.

Texture features that aim to characterize spatial patterns of image intensity can also be affected by the presence of unwanted artifacts. Materka and Strzelecki (2013) note that these features are susceptible to intensity inhomogeneities that require bias field correction. Moreover, it has also been shown that texture features depend not only on texture, but also on the mean and variance within an ROI. Thus, these features are also susceptible to intensity drift, and normalization is necessary to ensure that accurate texture features are calculated.

In order to assess the impact of each pre-processing step on the features extracted from mpMRI, we compared the value of each of the features, both within healthy tissue and in the suspected tumor, for each pre-processing step to the corresponding values in the non pre-processed images using the concordance correlation coefficient (CCC) (Isaksson et al 2020). The CCC is a measure of reproducibility between two sets x and y and is expressed as

$$\begin{eqnarray*}&&{\rm{C}}{\rm{C}}{\rm{C}}=\displaystyle \frac{2{s}_{xy}}{{\sigma }_{x}^{2}+{\sigma }_{y}^{2}+{\left({\mu }_{x}-{\mu }_{y}\right)}^{2}}\end{eqnarray*}$$

where $\sigma$ is the standard deviation, $\mu$ is the mean, and ${s}_{xy}$ is the covariance. The features evaluated with the CCC included only those extracted from the T2w and ADC images because only those are affected by the pre-processing steps. They included the T2w and ADC intensity values, the 24 texture features extracted from the T2w images, and the blobness feature calculated for the ADC images.

The optimal parameter values for each bias field correction method are shown in table 1. Of all the parameters available for N4, only the value of the spline order is equal to the default value. All of the optimized values for the PABIC parameters differ from the default values. The bias field correction applied to the T2w images improved image quality when compared to unprocessed images, but the results varied according to each method and the individual metrics used for evaluation (figure 3). N4 achieved the lowest value for CoV and entropy while PABIC marginally outperformed those images corrected with N4 when considering the CJV. In order to assess the benefit of bias field correction within the context of tumor localization, tumor probability maps were generated using T2w images corrected with either N4 or PABIC. Tumor classification performance in terms of the AUC was highest using PABIC bias field correction, whereas classification using N4 was equivalent to that of the T2w images before any bias field correction (table 2). As a result, PABIC was the bias field correction method selected for this study. A closer inspection of the effect of bias field correction on tumor localization showed that removal of the bias field led to intensity values that more accurately characterized PCa, i.e. an increase in true positive values in the tumor probability maps (figure 4).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

After bias field correction with PABIC, each of the three normalization methods was applied to the T2w images. Neither basic nor ROI normalization achieved a CoV value less than that achieved by PABIC alone (figure 5). Histogram matching normalization, however, did minimize the CoV, by comparison. While all three methods led to an improvement in tumor localization compared to bias field correction alone, basic and histogram matching normalization slightly outperformed ROI normalization in terms of the AUC (table 2). Histogram matching was thus chosen as the method for normalization as it was found to score best for both the CoV and tumor localization. In general, bias field-corrected images that were normalized according to the histogram matching technique demonstrated a slight reduction in the likelihood of false positives present in tumor probability maps generated on T2w images that had only been corrected using PABIC, though the difference was not substantial (figure 6).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Identical parameter options were assessed for each registration method where available. The resulting optimal values are given in table 3. Of the three methods, NiftyReg adhered most closely to the default or recommended values. All methods performed best at a single resolution level with no smoothing while both elastix and NiftyReg achieved the highest DSC value for a regularization weight of 0.05. All other parameter values differed according to the method used for registration. DSC values demonstrated an improvement in registration quality for elastix when compared to NiftyReg and ANTs (figure 7). The mean range of Jacobian determinant values over all subjects showed that elastix is the only method within the acceptable range of zero to slightly more than one (table 4). The mean TRE values, however, showed that NiftyReg attained a higher registration quality over elastix and ANTs. Tumor localization performance, though, was greatest for elastix compared to NiftyReg and ANTs (table 2). Thus, elastix was ultimately chosen as the preferred method for registration of the ADC images to the T2w images. When successful, this registration improved alignment of the tumor location across modalities and thereby improved the accuracy of the tumor localization model in that region (figure 8). Registration of the ADC using elastix resulted in poor alignment with the T2w image for a single patient (figure 9). In this case, the bladder of the ADC image was improperly aligned with the prostate gland of the T2w image, which resulted in misidentification of the tumor in the resultant probability map. Removal of this patient led to an increase in the AUC to a value of 0.87, though the result was still not significant (table 5). Tumor probability results for all three cases are shown in (see figures S1–S3 in supplementary materials).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When considering the cumulative effect of multiple pre-processing steps, the application of bias field correction along with normalization achieved the highest level of tumor localization performance, followed closely by the full pipeline of bias field correction, normalization, and registration (figure 10). In general, the application of these pre-processing steps improves tumor localization by increasing the number of true positive results, as well as reducing the number of false positives. This improvement is greater for tumors located within the PZ compared to those in the TZ (figure 11), which are traditionally more difficult to identify. Similar results were obtained for the analysis of the generalizability of using an optimized pre-processing pipeline to improve tumor localization accuracy (see S5–S13 in supplementary materials).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The effect of each pre-processing step on the features extracted from the T2w and ADC images is shown in figure 12. The 24 texture features extracted from the T2w images are designated by ${G}_{i}\left(s\right)$ where G is the Gaussian derivative in direction i at scale s (mm). Gaussian derivatives up to second order were calculated at four exponentially increasing scales ($\sigma$ = 1.5, 2.4, 3.8 and 6.0 mm). Any feature with a CCC value above an absolute threshold of 0.8 was considered to be unchanged (Akoglu 2018). In general, the features extracted from the T2w images that have undergone bias field correction are more highly correlated with non pre-processed images than those images that have been both bias field corrected and normalized. In fact, there is a large decrease in the mean CCC values of the T2w features after both bias field correction and intensity normalization have been performed (table 6). This is consistent with the change in tumor localization accuracy that is found with intensity normalization in addition to bias field correction. Additionally, both features extracted from the ADC images were found to differ from non pre-processed images in healthy tissue as well as tumor tissue. Moreover, features extracted from mpMRI of tumor tissue were less correlated with non pre-processed images than those extracted from healthy tissue.

Zoom In Zoom Out Reset image size