Evaluation of an automatic MR-based gold fiducial marker localisation method for MR-only prostate radiotherapy

MRI-based radiotherapy (RT) planning—also called 'MR-only'—has been proposed for several reasons, including the reduction of systematic spatial uncertainties introduced when registering the CT and MRI, reduction of the costs of treatment and a reduction of patient exposure to ionising radiation (Fraass et al 1987, Lee et al 2003, Karlsson et al 2009, Edmund and Nyholm 2017). To fully exploit the possible advantages offered by MR-only RT, several challenges need to be addressed: (1) performing MRI in the RT treatment position, (2) enabling MR-based dose planning, and (3) generating MR-based reference images for the position verification of image-guided radiotherapy (IGRT) treatment (Kapanen et al 2013).

The first point is addressed by the 'MR-simulator', which utilises dedicated devices such as flat table tops (Mcjury et al 2011), coil supports (Kapanen et al 2013, Sun et al 2014), skin markers and lasers (Schmidt and Payne 2015) to facilitate reproducible patient positioning.

The second point was investigated through the estimation of MR-based electron density maps, which are called pseudo-CT (pCT)⁵ images since they are generated to substitute CT images in the current workflow allowing for the calculation of the attenuation of radiation in the body (Edmund and Nyholm 2017). For prostate cancer radiotherapy, several ways of generating pCTs have been proposed (Lee et al 2003, Chen et al 2004, Dowling et al 2012, Johansson et al 2012, Korhonen et al 2014, Kim et al 2015, Siversson et al 2015, Prior et al 2016), resulting in acceptable treatment plans.

To address the third challenge, we first identify which position verification techniques have been clinically adopted for IGRT (Das et al 2014, Nabavizadeh et al 2016). IGRT for prostate cancer can be performed based on the location of the prostate itself (usually using a cone beam CT), anatomical landmarks, or prostate surrogates, e.g. intraprostatic gold fiducial markers (FMs) and implanted electromagnetic transponders (Bujold et al 2012, Ng et al 2014, Zaorsky et al 2017). Most MR-only oriented investigations have so far focused on the MR-based generation of digitally reconstructed radiographs (DRRs) for bone alignment (Chen et al 2004, Dowling et al 2012, Korhonen et al 2014, Kim et al 2015, Tyagi et al 2016). Nevertheless, it has been shown that bone-based landmarks represent an unreliable surrogate of the true prostate position due to the motion of the prostate with respect to the pelvic bones (Schallenkamp et al 2005), demanding larger PTV margin recipes with respect to FM-based alignment (van Herk et al 1995, Beard et al 1996, Greer et al 2008). Previous research has established that implanted FMs permit significantly improved localisation of the target (Fuller and Scarbrough 2006, van der Heide et al 2007, Kupelian et al 2008). Therefore, the MR-based localisation of intraprostatic gold FMs is crucial for accurate prostate irradiation.

Gold FMs are usually cylindrically shaped, with a diameter ranging between 0.5 and 1.5 mm and a length between 2 and 10 mm. Typically, three markers are implanted in the prostate prior to treatment and permanently remain in the target after completion of the treatment (Habermehl et al 2013, Ng et al 2014, Tanaka et al 2016). During treatment planning, the FMs are easily localised on the CT images thanks to their high signal and specific streaking artefacts (Nederveen et al 2000, Habermehl et al 2013, Chan et al 2015). In MRI, FMs are depicted as signal voids in magnitude images, and the appearance of these voids varies according to the imaging parameters (Vonken et al 2013, Lim et al 2016) and the FM orientation with respect to the magnetic field (Jonsson et al 2012). In current clinical practice, MR-based manual FM localisation is performed by radiation therapy technicians (RTTs) (Parker et al 2003, Ung and Wee 2011, Kapanen et al 2013). As the signal void in the magnitude images is not highly specific, manual FM localisation can be a challenging task (Chan et al 2015), as calcifications or blood clots may resemble FMs (Hong et al 2012, Ng et al 2014). To ease MR-based FM localisation, RTTs usually estimate the position of gold FMs after their initial localisation on the CT, which easily facilitates their distinction from calcifications and bleeding.

However, in an MR-only workflow, CT images are not available, therefore methods that enhance FM visibility, or even enable localisation, are demanded to identify the reference patient position that is used during radiotherapy treatment.

Up to now, several groups have proposed supporting manual FM localisation by enhancing metal object visibility on the MRI through positive contrast (Seevinck et al 2011, Varma et al 2011, Zhu et al 2011, Vonken et al 2013, Dong et al 2015), or sequence parameter optimisation (Lim et al 2016). Attempts to pursue automatic MR-based FM localisation have been scarce. A semi-automatic localisation algorithm was proposed by van Dalen et al (2003), while fully automatic methods for MR-based FM localisation were only recently presented (Ghose et al 2016, Dinis Fernandes et al 2017). A key feature of Ghose et al's method is the use of template matching (Brunelli 2009) between in vivo data and a manifold template on manually segmented FMs. Dinis Fernandes et al (2017) investigated the use of a logistic regression model on individual and combined sequences, obtaining the most promising results when a combination of sequences was employed.

Although the results of these approaches have been promising, we see room for improvement, hypothesising that a more accurate fiducial localisation may be achieved by increasing the information available, e.g. using complex MR images. In this sense, recently, a template matching method based on phase correlation for metal object localisation was proposed (Zijlstra et al 2017), inspired by Wachowicz et al (2006) and Kuglin and Hines (1975), which was applied in several brachytherapy applications (Zijlstra et al 2016, Beld et al 2016). Instead of merely using the spin dephasing information presented in magnitude images, this method exploits the phase dispersion induced by susceptibility markers, which is hopefully detectable when employing complex-valued MR images. Furthermore, detailed knowledge of the objects to be localised (geometry, orientation, magnetic susceptibility and B₀ distortion), as well as of the MR imaging parameters, are incorporated in this approach. In this way, we aim to exploit all available prior knowledge.

In this study, we will apply this approach to FM localisation after performing the adaptations necessary to enable gold FM localisation. In particular, we generated a new template specific to gold FMs, and designed an image processing pipeline that allows the selection of candidates and eventually their identification. We will present a clinical evaluation of the method against MR-based manual localisation and the gold standard in the clinical practice: CT-based localisation.

We conducted a study to assess the performance of the MR-based automatic FM localisation approach proposed here. The study was performed on prostate cancer patients undergoing external beam radiotherapy treatment to evaluate the performance in a clinical setting. This study is divided into four parts. First, we acquired the CT and MR images (section 2.1), including a dedicated MR sequence for FM localisation purposes. Second, we obtained the locations of the FMs as they were estimated from the CT images in clinical practice (section 2.2). Third, we performed both manual and automatic MR-based localisation (section 2.3). Finally, we evaluated the detection rate, accuracy and precision of automatic FM localisation and compared it to the outcome of the clinical CT-based and MR-based manual localisation method (section 2.4).

2.1. Patient data collection

2.2. Clinical fiducial marker localisation

2.3. MR-based fiducial marker localisation

MR-based manual and automatic FM localisation was performed on the second echo of the GRE sequence.

2.3.1. Manual localisation.

An experienced RTT manually localised the FMs by clicking the apparent top and bottom of the marker on the magnitude images. Subsequently, the centre of the FM was calculated as the mean position between the top and bottom, and the orientation of the vector connecting the top and bottom. Note that the apparent location of the top and bottom of each FM coincided with the centre of a voxel in the GRE sequence; nevertheless, the calculated centre of the FM may have been positioned with subvoxel resolution. The actual position of the top and bottom was finally assigned knowing the position of the centre, the orientation and the length of the FM. The term 'MR_man' was used to refer to the FM locations performed on the MR images by the manual observer.

2.3.2. Automatic localisation.

Here, we describe a method that localises FMs using complex images (real and imaginary) of the GRE sequence. Zijlstra et al (2017) recently developed and presented this method for brachytherapy seed localisation. Here, we present an adaptation of the approach for gold FM localisation purposes. For clarity, all the steps of the method will be introduced in this section and the details that may be useful to reproduce the method have been provided in the appendix. During the algorithm design, we assumed that the number of implanted FMs (three) and the target delineation (CTV) were available prior to FM localisation. To assess the specificity of the approach, the algorithm was allowed to search for a maximum of four FMs, even if only three had been implanted. The approach was subdivided into four main steps: (1) library generation, (2) template matching, (3) reduction of candidates and (4) selection of FMs.

With respect to the previous work Zijlstra et al (2017), a new template library specific to the gold fiducials and scan parameters was generated. Furthermore, the image processing pipeline described at points (3) and (4) is novel. Note that in these steps, some parameters involved in candidate reduction and selection have been determined empirically after heuristic optimisation, which may have influenced the results in a positive way due to overfitting. Figure 1 shows these steps organised according to their chronological order within the patient study. Here follows a description of each of the four steps:

• (1)  
  Library generation. First of all, the artefact generated by the FM in the GRE sequence was simulated in two steps: (i) calculation of the B₀ perturbation induced by an FM and (ii) simulating the resulting MR signal in proximity of the FM within the frameworks presented by Bouwman and Bakker (2012) and Zijlstra et al (2016). In particular, during the calculations, a 3D analytical model of the gold FM in a uniform background was used. The term 'template' was used to refer to the output of the simulation. Note that the output of the simulation was complex-valued data. To allow localisation of the FM along with its orientation, several templates were generated orienting the FM along different orientations with respect to the main magnetic field direction. This collection of templates was referred to as the 'library'. Note that the library generation step was performed only once during the whole study and is specific for the sequence and the FMs employed.
• (2)  
  Template matching. For each patient, a region of interest (ROI) was defined as the smallest volume containing the clinical delineation of the prostate after an isotropic expansion of 1.5 cm. A phase correlation (PC) (Kuglin and Hines 1975) was performed between the image in the ROI and each template in the library. The template with the highest PC within the library was associated with each voxel in the ROI, such that each voxel was associated with an FM orientation with respect to the main magnetic field.
• (3)  
  Reduction of candidates. In principle, all the voxels in the ROI may be selected as FM candidates. To reduce the number of candidates two sub-steps were performed. (i) We hypothesised that the voxels with the highest PC have the highest probability of containing an FM. The first 250 voxels with the highest PC were selected as FM candidates. (ii) To further reduce the number of candidates, we discarded the voxels that were located in hypointense regions, with a volume much larger than the hypointense region induced by the FM.
• (4)  
  Selection of FMs. The remaining candidates were sorted, minimising the error calculated by a local linear regression between the in vivo data and the template. Afterwards, for each FM candidate, a smoothness criterion was applied to the phase image starting with the candidate with the lowest linear regression error until the criterion was satisfied for a maximum of four voxels.

Zoom In Zoom Out Reset image size

The output of the algorithm was a centre, top and bottom for each FM. The total time to simulate a library of templates was below 7 s, and the maximum time needed to locate the FMs for a single patient was about 20 s in Matlab (R2015a, MathWorks, Natick, MA, USA) and using a quad core 3.4 GHz commercial CPU. The term 'MR_auto' was used to refer to the FM locations performed on the MR images by the method proposed here.

2.4. Evaluation

The performance of the proposed MR-based automatic FM localisation approach (MR_auto) was compared with the performance of clinical CT-based (CT), and manual MR-based (MR_man) localisation. The detection rate and inter-marker distances (IDs) between the FM positions were determined.

Figure 2 shows an example of the FM location for one patient, as obtained by the three approaches in their respective frame of reference. Since the frames of reference differ between the localisation methods, the FMs were labelled (from 1 to 3) according to the position of their centre along the craniocaudal direction prior to comparison. In cases in which the FMs were in the same transverse plane, the left–right direction was used for labelling.

Zoom In Zoom Out Reset image size

For the three FM localisation approaches, the apparent detection rate was calculated as the number of located FMs divided by the number of implanted FMs (51, since three FMs were implanted for each patient). We define such a detection rate as 'apparent' since, at this stage, no check has been performed in terms of false or positive FM detection.

2.4.1. Geometrical evaluation and detection rate.

The inter-marker distance (ID) between two FMs was calculated, referring to ${\rm ID}_{ij}$ as the distance between FMs with the label i and j (with $i, j=1, 2, 3$ and $i\neq j$). To compare the FM location of the three localisation approaches, the differences between the IDs with the same label between all the FM localisation modalities were calculated as ID${\hspace{0pt}}_{MR_{\rm man}}$-ID_CT, ID${\hspace{0pt}}_{MR_{\rm auto}}$-ID_CT, and ID${\hspace{0pt}}_{MR_{\rm man}}$-ID${\hspace{0pt}}_{MR_{\rm auto}}$. As these differences can be both positive and negative, a value of zero is to be expected when calculating the mean, assuming no non-rigid deformations. The absolute values for the ID differences were also calculated. To test the precision of the FM detection with respect to CT, each FM was considered as correctly localised when the associated absolute difference IDs⁶ were <5 mm (=the length of an FM). We considered the mean and STD of ID${\hspace{0pt}}_{MR_{\rm man}}$-ID_CT and ID${\hspace{0pt}}_{MR_{\rm auto}}$-ID_CT, excluding the IDs associated with incorrectly localised FMs, to verify whether MR-based localisation was biased with respect to the clinical gold standard (CT). Furthermore, we considered the STD of the absolute value of ID${\hspace{0pt}}_{MR_{\rm man}}$-ID_CT and ID${\hspace{0pt}}_{MR_{\rm auto}}$-ID_CT as representative of the precision of the MR-based localisation with respect to the clinical gold standard (CT). To allow for a complete comparison of the localisation approaches, the mean and STD were also calculated between ID${\hspace{0pt}}_{MR_{\rm auto}}$-ID${\hspace{0pt}}_{MR_{\rm man}}$. Note that the precision of the IDs between MR-based and CT-based FM localisation are expected to be larger than the IDs between MR-based FM localisation approaches due to the possible presence of a non-rigid inter-scan motion.

The benefit of these metrics is that IDs are independent of the absolute position of the FMs in the coordinates used during the localisation; this therefore allows a comparison of the imaging modality with different frames of reference without recurring to image registration. However, there are certain drawbacks associated with the use of these metrics, since the mislocalisation of one FM may bias two inter-marker distances. To provide additional metrics that are independent of the inter-marker distances, the euclidean distance between FMs with the same label was calculated between MR_auto and MR_man, as well as the difference of the FM positions in each direction ($X, Y, Z$). The advantage of these metrics is that MR-based localisation is intrinsically performed in an identical frame of reference. This enables the accuracy of the FM localisation of MR_auto to be estimated against MR_man. The mean, maximum difference and STD among all the patients and corresponding FMs were calculated. We considered the mean and STD as representative of the accuracy and precision, respectively, of localisation between the MR-based localisation approaches. Finally, a t-test at a 1% significance level was conducted on the difference between each of the components of MR_auto-MR_man to verify whether any bias was observed in terms of FM localisation.

In addition, the detection rate was calculated for MR_man and MR_auto as the number of correctly localised FMs divided by the total number of implanted ones (51). A case study based on the observation of CT and MR images for all the patients was chosen to obtain further in-depth information on the causes of incorrect detection and to verify whether a false positive FM had been detected. Cases with further positive FMs were reported; furthermore, patient $\#14$ also highlighted the performance of FM localisation in the presence of hip implantation.

2.4.2. Impact on patient alignment.

To compare the two MR-based localisation methods (MR_man and MR_auto) and investigate the possible influence of their difference on the patient alignment, for each patient we performed registration between the CT and both the MR-based localisation methods. This was done by matching the centre of the FMs between the CT and MR_man and between the CT and MR_auto. The matching of the FMs was performed by only applying translations to realistically mimic patient alignment at the linear accelerator. More specifically, for each patient, the centre of mass (the average among the three FMs) of the FMs localised by MR_man and MR_auto was aligned (translation only) with the centre of mass of the FMs localised on the CT. To compare the differences between the registrations, we assess their quality in terms of the root mean square of the residual mismatch between the FMs with the same label (Ung and Wee 2011), referring to the root mean square error (RMSE). The RMSE was calculated for each patient after the registration. In cases in which the RMSE was larger than 2.5 mm, the quality of registration was considered poor, and the registered FMs were visually checked to assess the origin of the poor registration. In addition, the difference between the translation vectors obtained aligning MR_man to CT and MR_auto to CT in each direction ($X, Y, Z$) was calculated for each patient using the mean, STD and range over the patient population. A t-test at 1% significance level was calculated for each component of the difference between the translation vectors over the patient population to verify whether any systematic shift had been introduced between registration based on the two MR-based localisation modalities.

3.1. Detection rates

Table 1 provides the summary statistics for detection rates obtained by the three approaches introduced in sections 2.2 and 2.3. From the first column of the table, it appears that the localisation clinically performed on the CT, as well as that manually performed on the MR (MR_man), detected all the implanted FMs. Interestingly, the MR-based automatic approach proposed here (MR_auto) never located $> 3$ FMs for any of the patients, even if the detection of a maximum of four FMs was allowed. For one patient (#1), two FMs were localised resulting in an apparent detection rate of 0.98. Considering the incorrectly detected FMs, the result was that one of them had not been correctly localised for patient $\#4$ by either the MR_man or MR_auto, lowering the detection rate for both MR-based localisation approaches with respect to their apparent detection rate.

Table 1. Apparent detection rate (first column), and detection rate (second column) for the CT-based clinical localisation (CT), MR-based manual (MR_man) and MR-based automatic (MR_auto).

	Apparent detection rate	Detection rate
CT	51/51 $=$ 1	Ground truth
MRman	51/51 $=$ 1	50/51 $=$ 0.98
MRauto	50/51 $=$ 0.98	49/51 $=$ 0.96

3.2. Geometrical accuracy and precision

3.2.1. MR_man versus MR_auto.

The euclidean distance between FMs with the same label between MR_auto and MR_man was 0.91, 0.53 and 3.14 mm in terms of the mean, STD and maximum among all the patients and FMs, respectively. Both the spatial accuracy (mean) and precision (STD) between the MR-based localisation approaches appears to be below the image resolution ($1.1 \times 1.1 \times 1.2$ mm³). The maximum difference value occurred for patient $\#4$ and FM 2, but as outlined in the previous section, incorrect detection occurred for this FM. Excluding this FM from the statistics, the maximum value obtained was reduced to 1.73 mm, which is a distance <2.5 mm (2 voxels). Figure 3 shows a histogram of the difference MR_auto-MR_man of the position of the corresponding FMs (with the same label) for all three directions. From the figure, it is obvious that the maximum difference observed between the manual and automatic MR-based FM localisation approaches was within 2 voxels. The null hypothesis was rejected ($p <0.001$) in the anterior–posterior ($X=$ readout direction) direction according to the result of the t-test. This result suggests that FM localisation between the manual observer and the automatic approach along the X direction significantly differed. In particular, in terms of the average among all the FMs, the absolute difference in X was 0.26  ±  0.49 mm (1 STD).

Zoom In Zoom Out Reset image size

3.2.2. Inter-marker distances.

Figure 4 shows the difference between the IDs of the corresponding (having the same label) FMs between localisation approaches. From this figure, we can see that patient #4 has ID differences >5 mm for FM $\#2$ for both the MR-based localisation approaches. This suggests that such an FM was incorrectly located by MR_man and MR_auto, resulting in a false positive detection. Furthermore, ID differences relative to FM $\#3$ localised by MR_auto are not plotted for patient $\#1$ in figure 4, since such an FM was not located (false negative detection).

Zoom In Zoom Out Reset image size

In summary, figure 4 highlights that: for patient $\#1$ MR_auto, FM $\#3$ resulted in false negative detection for the MR_auto approach, while for patient $\#4$, FM $\#3$ resulted in false positive detection for the MR_man and the MR_auto approaches. These FMs are considered as outliers and a more detailed analysis is given in the following section. As reported in table 2, averaging the ID differences among the true positive detected FMs in all the patients, the mean and STD were about 0 mm and <1.2 mm for all the compared FM localisation approaches. We conclude that the ID metric was not biased over the patient population. We identified the STD over the absolute values of the ID differences as indicative of the precision of the FM localisation. It is apparent from table 2 that the precision of the localisation was comparable between approaches. Therewith, the precision resulting from a comparison of all three methods was lower than the voxel size ($1.1 \times 1.1 \times 1.2$ mm²) of the GRE sequence. From the table, it can be seen that the maximum ID difference obtained was <3 mm when the localisation was performed between different imaging modalities, and <2 mm between the MR-based FM localisation approaches.

Table 2. An overview of the statistics of the correctly detected FM in terms of inter-marker distance (ID) differences and the absolute value between the localisation approaches (CT, MR_man and MR_man). The results are expressed in mm.

	ID difference		|ID difference|
	Mean	STD	Mean	STD	Maximum
ID${\hspace{0pt}}_{MR_{\rm man}}$-IDCT	0.01	1.13	0.90	0.66	2.80
ID${\hspace{0pt}}_{MR_{\rm auto}}$-IDCT	0.00	1.14	0.89	0.70	2.66
ID${\hspace{0pt}}_{MR_{\rm man}}$-ID${\hspace{0pt}}_{MR_{\rm auto}}$	0.03	0.77	0.62	0.45	1.71

3.3. Outliers of MR-based FM localisation

To understand the possible origins of false negative and false positive detections, we present the CT and MR images of patients $\#1$ and $\#4$. As shown in figure 5 (left), in patient $\#1$ false negative detection occurred in the proximity of a large hypointense signal around the FM on the GRE (circle on the bottom left image). To further understand the origin of such a large hypointense signal, we checked all the MR sequences acquired during the imaging session. From a 3D T1 weighted sequence, we concluded that bleeding results in an enlargement of the dimension of the hypointense region in proximity to the FM. This may explain the false positive detection, since the artefact generated by the FM may not, in this case, match the simulation due to the presence of blood. Figure 5 (right) shows the scenario for patient $\#4$: a hyperintense area on the CT (within the square), which is not characterised by the CT streaking artefacts, was localised as an FM on the MR images, but it was not detected as an FM on the CT. Such a hyperintense area on the CT occurred in the proximity of a calcification. We hypothesise that the presence of calcifications contributed to the false positive and negative detection for patient $\#4$. During the inspection of all the patients, we did not notice any further false positive detection.

Zoom In Zoom Out Reset image size

3.4. Impact on patient alignment

As shown in figure 6 (left), the quality of registrations, only taking into account translation, was comparable regarding the MR-based methods for all the patients. For patient $\#4$, poor registration quality was found. As already observed in the previous section 3.3, for such patients false positive and negative detection occurred. Interestingly, even for this patient ($\#4$), when considering the centre of mass, patient alignment seems to be acceptable (figure 6, right).

Zoom In Zoom Out Reset image size

As can be seen in the table 3, the range of differences between the translation vector is within the voxel resolution, and therefore we can conclude that patient alignment based on the FM location obtained with the method proposed here was comparable (within 1 voxel) to the patient alignment based on FM location obtained by a manual observer. Nevertheless, the p-value highlights that a systematic shift occurred in the anterior–posterior direction (X), which is the frequency encoding direction. This is in accordance with the results reported in section 3.2, and it reveals that the FM localisation difference also impacts patient alignment along the readout direction. In the following section, we will elaborate on the possible causes of such a reported difference.

Table 3. The statistics in terms of mean, STD, range and t-test over the patient population of the difference between the translation vectors obtained with the two MR-based localisation approaches. The statistics, except for the p-values, are expressed in mm.

	Translation vector
	MRauto-MRman
	X	$Y$	Z
Mean	0.28	$-0.08$	0.17
STD	0.33	0.39	0.42
Minimum	$-0.37$	$-0.75$	$-0.20$
Maximum	0.75	0.56	1.20
p-value	0.003a	0.41	0.12

^ap-value $<0.01$ $=$ the null hypothesis is rejected.

The state of the art IGRT of prostate cancer within an MR-only workflow requires the localisation of gold FMs (Fuller and Scarbrough 2006, van der Heide et al 2007, Kupelian et al 2008).

To facilitate MR-based FM localisation, we proposed an MR-based automatic gold FM localisation approach, which relies on a template matching method applied to the complex MRI data, aiming to fully exploit signal behaviour in the vicinity of an FM. The method proposed here is performed in a similar way to an operator resulting in 50/51 FMs automatically localised within 2 voxels (2.4 mm) of the manually localised FMs. In addition, when compared to the clinical standard (CT), automatic MRI-based FM localisation results in all three FMs being correctly detected for 15 out of 17 patients with a precision $<$$1$ mm. When comparing our approach to the clinical standard, inter-scan motion may have occurred between the two imaging sessions (van Herk et al 1995, Beard et al 1996, McPartlin et al 2016, Kupelian et al 2008). To minimise the effect of inter-scan motion, we used inter-marker distances (IDs) as a metric, which are insensitive to rigid-body transformation but sensitive to prostate deformation (Kupelian et al 2005, van der Heide et al 2007). The difference between the inter-marker distances (IDs) of the automatic MR-based localisation and clinical standard resulted in a mean of 0 mm and an STD of 1.14 mm. This result may indicate that prostate deformation was not observed in our study over the patient population. This is in line with previous research when considering ID variation between irradiation sessions (Kupelian et al 2005).

To understand the performance of FM localisation in the presence of a hip implant, we manually inspected the FM location for the patient $\#14$. No evidence was found for the systematic difference between MR-based FM localisation and between MR- and CT-based localisation. The result is surprising and it encourages future research for these cases, nevertheless, these findings cannot be extrapolated to all hip implant patients.

When considering false detection, bleeding and calcifications were among the reported causes. This is in line with previous research (Hong et al 2012, Ng et al 2014). More specifically, calcifications are expected to occur in about 40% of the cases (Hong et al 2012). Given the high rate, finding the correct FM location for such cases is highly necessary for facilitating clinical implementation. Previous studies showed that susceptibility-weighted imaging may help to identify calcifications and bleeding (Wu et al 2009, Bai et al 2013, Chen et al 2014, Berberat et al 2014), and thus an approach that foresees the integration of susceptibility-weighted imaging with the approach presented here may be of interest for a future study.

To date, two studies have investigated the performance of MR-based approaches for gold FM localisation. Ghose et al (2016) proposed to exploit manifold learning and spectral clustering to localise gold markers from manually segmented volumes. Their method showed promising results, however, it required manual interaction to generate the references for data training. Manual interaction may need to be repeated when varying the MR image parameters or type of FMs. Dinis Fernandes et al (2017) investigated the use of logistic regression, obtaining the most promising results when several MR sequences were employed. In terms of detection rate, our study resulted in a higher number of correctly detected FMs with respect to the studies just presented.

An aspect that is important for MR-only radiotherapy is the accuracy of FM localisation. MR images may be subject to geometric distortion caused by the system and/or the patient (Fransson et al 2001, Walker et al 2014, Weygand et al 2016). When considering MR-simulators, strict geometrical fidelity tests are mandatory to ensure that system-dependent geometric distortion is minimised (Kapanen et al 2013). Note that since system-related distortion increases proportionally with the distance from the isocentre of the MR scanner, we expect the geometric fidelity to be higher than 1 mm within a distance of 10 cm from the isocentre (Wang et al 2004, Sun et al 2014, Xing et al 2016). If the patient is positioned with the prostate in the centre of the FOV, we might expect accuracy of the same order. Nevertheless, patient-related distortion may also contribute, especially when metallic objects are present (Hargreaves et al 2011), or when the rectum is filled with air, as both metal and air induce magnetic susceptibility-related field changes. Previous studies have demonstrated that gold FMs may cause small deviations in the position ($<1$ mm) of the imaged marker with respect to the actual marker (Jonsson et al 2012). As Jonsson et al (2012) showed, deviations can be mitigated by the MR imaging parameters. We believe that geometrical accuracy is crucial for MR-only RT, therefore we utilised a method that provides the original location of the FM by taking into account the FM-induced B₀ distortion. By doing so, potential geometrical distortion related to the interference between the local field inhomogeneities caused by the fiducials and the readout gradients are mitigated. In addition, the dedicated GRE sequence used in our study had a high BW (1142 Hz/voxel), minimising possible distortion in the readout direction to a sub-voxel level. Surprisingly, a statistically significant difference in the readout direction (anterior–posterior, or X) between MR-based localisation approaches was found. A possible explanation for this might be that the observer perceived that the signal void was shifted with respect to the real position of the FM, whereas the automatic method incorporated these effects. As Jonsson et al (2012) underlined, a deviation of about 0.5 mm was reported for the susceptibility induced by a gold marker and a bandwidth $=$ 1200 Hz/voxel. This is compatible with the observed deviation along X, and it is in line with the previously accepted distortion for MR-based FM localisation (Tyagi et al 2016, Dinis Fernandes et al 2017). In our study, we considered the MR-based manual FM localisation as a standard without further questioning its precision or geometrical accuracy. Unfortunately, the availability of studies to quantify inter-observer variability in gold FM localisation is limited when excluding x-ray-based studies. Huisman et al (2005) reported that inter-observer variability resulted in uncertainties within 0.4–0.6 mm. The uncertainties were calculated as STD at the centre of the prostate between the CT-MR registration. In our study, uncertainties (STD) impacting the translation were in the range of 0.3–0.4 mm, which is in line with the results of the previous study (Huisman et al 2005).

The FM localisation has only been demonstrated with a single centre dataset; hence, a generalisation of the method across scanners and different shapes of gold FMs requires further tests. In any case, the method proposed here is, in principle, applicable to different FMs as long as the shape, dimensions and composition of the FMs are known. In addition, our method is flexible since it may be employed for different MR sequences. In principle, a different sequence may possibly impact the appearance of gold FMs and the performance of marker localisation. In this sense, optimisation of the image parameters might result, which has the benefit of reducing false positive or false negative detection, and might be the object of future studies. Note that if a different combination of MR sequence and FM is employed, a new template library will need to be generated.

We believe that the proposed approach will facilitate MR-only RT workflows by reducing manual interaction. Moreover, such a method could also be applied to a standard MR-CT workflow to reduce the workload for registration purposes. To ensure patient alignment, at least three FMs need to be localised (Ung and Wee 2011). The results obtained in terms of false positive detection would not allow patient alignment for 2/15 patients. This rate is considered too high to accept a fully automatic procedure. In this sense, we believe that a semi-automatic scenario would still be favourable to reduce workload and limit FM localisation variability (van Dalen et al 2003). In particular, the proposed sequence has already been shown to facilitate the manual depiction of FM with respect to single echo GRE and fast spin echo sequences (Schieda et al 2015). Further speculating, introducing a fourth intraprostatic FM is a scenario that could facilitate automatic FM localisation. This can be carried out without creating additional patient discomfort, since for each needle insertion, two FMs can be released, as this has already currently been done for two of the FMs. In this way, the automatic algorithm can still ensure target localisation on the reference images in case one out of four FMs is falsely detected. This may enable the safer introduction of automatic FM localisation. Another approach that it is worth investigation is automatic quantification of the correctness of FM localisation. Interestingly, Ghose et al designed an automatic warning to call attention to the possibly incorrectly detected FMs. We believe that such a feature facilitates clinical implementation of MR-based automatic FM localisation. In this sense, our experience suggests that phase correlation (PC) could help to quantify the confidence of detection since a high PC corresponds to a high match, but further investigation is required.

As a final note, we would like to stress the fact that, in this article, FM localisation was performed in a separate sequence with respect to the sequences used for pCT generation or target delineation. Ideally, to minimise intra-scan motion and to ensure geometric fidelity, the sequence used for target delineation and FM localisation should be performed in a single scan with an acquisition time on the order of 2/3 min, since intra-scan motion increases with the treatment time (McPartlin et al 2016). In this sense, we believe that MR-only would not benefit from a multi-sequence approach (Dinis Fernandes et al 2017). The proposed method can be, in principle, applied to any GRE sequence, and can therefore be integrated with the sequence used for pCT generation or for target delineation. Actually, the sequence presented here is similar to what has already been used by other groups (Tyagi et al 2016) and by our group (Maspero et al 2017) for pCT generation. In any case, when it is not possible to incorporate FM localisation with the sequence used for pCT generation or target delineation, our suggestion is to acquire three sequences one after another, such that the intra-scan motion is as low as possible and the relative position of the prostate in the sequence is consistent.

This project was undertaken to design an automatic MR-based method to localise intraprostatic gold fiducial markers and evaluate its performance. The results of this study indicate that gold FMs were automatically localised with a detection rate, accuracy and precision comparable to a human expert. We envision that the method may facilitate the introduction of MR-only RT workflow, reducing manual workload and ensuring geometrical fidelity.

The research is funded by the ZonMw IMDI Programme, project number: 1040030. The project is co-funded by Philips Healthcare. The study received approval of the medical ethical commission (Medisch Ethische Toetsingscommissie) and was classified under protocol number 15-444/C. Cornelis A T van den Berg declares to be a minority shareholder of MRCode BV. Max A Viergever discloses research contracts with Philips Healthcare, and with Pie Medical Imaging BV. Peter R Seevinck declares to be a majority shareholder of MRIGuidance BV.

In this section, details that may be useful for reproducing the work are presented along with an evaluation of the discriminative power of the fiducial marker (FM) selection. In particular, steps (1) and (2) have already been presented elsewhere, Zijlstra et al (2017), and represent an adaptation of the previous method for matching the gold FMs. The details follow, maintaining the division into steps presented already in section 2.3.2:

• (1)  
  Library generation. (i) The B₀ perturbation ($\Delta B_0$) induced by one FM was obtained by fast-forward field-shift calculations (Bouwman and Bakker 2012). During the calculations, the 3D analytical model of the gold FM consisted of a cylinder with diameter $=$ 1 mm and length $=$ 5 mm. The uniform simulated background was water. To perform the $\Delta B_0$ calculation, the magnetic susceptibility (χ) was set to $-9$ ppm in the background (water) and to $-34$ ppm in the FM (gold). (ii) To simulate the signal in the proximity of an FM for the GRE sequence employed during the in vivo data acquisition, we used the FORECAST framework (Zijlstra et al 2016). This framework simulates the off-resonance artefacts for steady-state GRE sequences with a lower computational complexity than full Bloch simulations. During the simulation, the proton density relative to water was set to 1 in the background and to 0 in the FM. In addition, the T₂ value of the background was set to 50 ms. Fast forward $\Delta B_0$ simulations were performed on a grid with a resolution $0.07 \times 0.07 \times 0.08$ mm³ and an FOV of $23.6 \times 23.7 \times 25.2$ mm³. The resolution of the $\Delta B_0$ calculation was chosen to be higher than the typical MR scan resolution, to avoid—at least at this stage—the partial volume effect. When employing the FORECAST framework, the resolution of the simulation was matched to the GRE sequence. During the generation of the library, a total of 321 templates were generated, resulting in an angular resolution of approximately $8^{\circ}$. Note that the FM is cylindrically symmetric, therefore during library generation any orientations with respect to B₀ with symmetry degeneracy were not considered.
• (2)  
  Template matching. Phase correlation (PC) (Kuglin and Hines 1975) was performed between the image (I) in the ROI and each template (T) in the library. Mathematically, PC was implemented as follows:

  $$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle {\rm PC} = \mathfrak{F}^{-1} \frac{\mathfrak{F} (I)^{\ast} \cdot \mathfrak{F}(T)}{\vert \mathfrak{F} (I)^{\ast} \cdot \mathfrak{F}(T) \vert + \epsilon} \nonumber \end{align*}$$

  where $\mathfrak{F}$ refers to the fast Fourier transform, ^* refers to complex conjugation and $\newcommand{\e}{{\rm e}} \epsilon=2.2\cdot 10^{-16}$ is the floating point accuracy during the calculation that was added to the denominator to avoid singularities in case both T and I were zero.
• (3)  
  Reduction of candidates. After the selection of the first 250 candidates (i) we performed further reduction of the image using a method inspired by Vonken et al (2013). In particular, a heuristic approach was used to define the volume of the hypointense region induced by an FM. First of all, the magnitude images were thresholded (25% of the maximum intensity in the ROI) and the regions were defined in the binary masks according to their connectivity (26 neighbourhood connected components). Second, the volume of one hypointense region was calculated from the simulated library since the size of the region may change according to the orientation (Jonsson et al 2012): for each template the volume of the hypointense region was calculated as the number of voxels with an intensity lower than the average value in the template, minus five times its standard deviation (STD). Averaging the volumes over the whole library, we estimated that the hypointense region was about 46 voxels (∼70 mm³). Finally, the candidates in regions with a volume larger than 184 voxels, which correspond to four times (=the maximum amount of FMs that we searched for) the volume of the hypointense region in the proximity of an FM, were excluded.
• (4)  
  Selection of FMs. The linear regression was locally performed on 3D patches with a template that was the same size. The smoothness criterion was applied to patches with equally sized templates and the candidate as the centre. In particular, the criterion was defined as follows: FMs were selected only if the local smoothness (S) of the image phase ($\phi_{I(x, y, z)}$) was increased after subtracting the template phase voxel-wise ($\phi_{T(x, y, z)}$). More specifically, S over the image (I) was calculated as:

  $$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle {\rm S}( \phi_{I(x,y,z)}) = \vert \frac{\partial {\rm e}^{{\rm i} \phi_{I(x,y,z)}}}{{\partial x} } \vert \nonumber \end{align*}$$

  where here (x) represents the phase encoding direction (transverse plane, left–right). The criterion can be mathematically written as:

  $$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \sum [ {\rm S}( \phi_{I(x,y,z)}) - {\rm S}( \phi_{I(x,y,z)} - \phi_{T(x,y,z)}) ] > 0 \nonumber \end{align*}$$

  where $\sum$ indicates the voxel-wise sum over the patch. The smoothness criterion was only calculated over the phase encoding direction after heuristically observing that this direction maximised the specificity of the selection.

The discriminative power of the candidate selection during steps three and four is reported in table A1 in terms of the number of voxels considered to be an FM candidate. As shown in the table, step four (i.e. the smoothness criterion) is revealed to have the highest discriminative power. Furthermore, the STD of step three (ii) is large since such a step is only applied in cases where a patient presents large signal voids that are usually generated by calcifications. Since not all the patients present intraprostatic calcifications, the discriminative power of such a step is variable.