Optimizing External Surface Sensor Locations for Respiratory Tumor Motion Prediction

Abstract

Real-time tracking of tumor motion due to the patient’s respiratory cycle is a crucial task in radiotherapy treatments. In this work a proof-of-concept setup is presented where real-time tracked external skin attached sensors are used to predict the internal tumor locations. The spatiotemporal relationships between external sensors and targets during the respiratory cycle are modeled using Gaussian Process regression and trained on a preoperative 4D-CT image sequence of the respiratory cycle. A large set ($N \approx 25$) of computer-tomography markers are attached on the patient’s skin before CT acquisition to serve as candidate sensor locations from which a smaller subset ($N \approx 6$) is selected based on their combined predictive power using a genetic algorithm based optimization technique. A custom 3D printed sensor-holder design is used to allow accurate positioning of optical or electromagnetic sensors at the best predictive CT marker locations preoperatively, which are then used for real-time prediction of the internal tumor locations. The method is validated on an artificial respiratory phantom model. The model represents the candidate external locations (fiducials) and internal targets (tumors) with CT markers. A 4D-CT image sequence with 11 time-steps at different phases of the respiratory cycles was acquired. Within this test setup, the CT markers for both internal and external structures are automatically determined by a morphology-based algorithm in the CT images. The method’s in-sample cross validation accuracy in the training set as given by the average root mean-squared error (RMSE) is between 0.00024 and 0.072 mm.

1 Problem

The localization of the internal tumors or structures and the detection of respiratory organ or tumor movement under certain therapies (e.g. in radiation therapy or radiofrequency ablation) and real-time visualization of these movements are an important concern in the safe and effective provision of precision radiotherapy, computer-assisted tumor surgery and biopsy. The information on the movements and tracking of the current positions of targeted tumors, e.g. in prostate, liver, lung or soft tissue tumors are traditionally performed by using abdominal compression, breath hold, respiratory gating, implanted radiation-impermeable markers and real-time motion tracking of these markers with optical (OP) or electromagnetic (EM) tracker [1, 2] or over real-time image processing of interventional 4D-CT, -MRI, and 3D/4D-ultrasound [3, 4]. Because of the difficulties with respect to the intraoperative complexity of these methods, there is currently no clinically established solution for the reliable determination of respiratory movements of internal organs, tumors or soft tissues.

Specifically, in the real-time motion tracking of external sensors with OP or EM tracking systems, it is difficult to determine the optimal amount and location of sensors preoperatively that provides sufficient predictions of the internal tumor locations. The used surface markers are randomly distributed [5] and placed on the patient in the near region of a surgical area and all of them are used for respiratory motion prediction, which can increase the error rate in the real-time prediction [6].

Fig. 1.

Left: Respiratory model with the fixed sensor-holders on it. Right: Interior view of the model during an inhalation. The inside of the balloon (max. inflation ${\oslash }$ 120 cm) is brought out and a rubber band from both ends of the balloon is glued, so that it stays in the middle, if the outside brought in and inflated. To simulate tumor motion; the X-Spot skin markers (${\oslash }$ 1.5 mm) and retro-reflective markers (${\oslash }$ 12 mm) are placed inside of the balloon. A flexible silicone tube (L 200 cm, ${\oslash }$ 20 mm) is used for inflating the balloon (inhalation/exhalation) with the water blaster (82$\,\times \,$5$\,\times \,$15 cm).

2 Materials and Methods

2.1 Phantom Respiratory System Model

The training and evaluation of internal target motion prediction is performed with a custom built phantom model (see Fig. 1). It simulates the most important aspects (with respect to motion prediction) of the human respiratory cycle. The main components consist of a standard rubber hot-water bottle (modeling the abdomen), an internally located spherical rubber balloon (modeling a moving organ). A flexible silicone tube and a water blaster was used to control the amount of air within the model. Sensor-holders within CT skin markers (external input) and skin markers or retro-reflective balls (internal target) are to use in prediction.

2.2 Sensor Holders

The external surface sensors are used in several real-time tumor movement prediction methods [7,8,9,10,11]. In most cases - however - the external sensors are fixed at arbitrary locations that may be sub-optimal for prediction accuracy. Optimizing the spatial distribution and quantity of those surface markers with respect to their prediction power in the preoperative phase therefore can improve the tracking accuracy in the intraoperative phase.

As the two phases typically require different types of external markers, custom made 3D printed sensor-holders were developed (see Fig. 2) to enable switching the sensors while maintaining the same sensor origin. This enables an offline prediction preoperatively using CT markers, and using the pre-trained predictors with a real-time tracking system during the intervention after the known relative transformation between the X-Spot marker and the inserted real-time tracker sensor is applied.

Fig. 2.

(a) Two M2 screws that are placed from both directions to establish a rigid setup once an OP or EM sensor is placed in the main part. (b) The main part within a X-Spot marker. It used in offline prediction step. (c) The EM-Sensor-Holder to fix the EM sensor in it. (d) View of sensor-holder within an OP marker. (e) View of sensor-holder within an EM sensor. The parts a, c, d and e are intended for the real-time prediction.

The main part of the sensor-holder consists of an X-Spot CT marker centered in a sensor attachment point. During the preoperative phase, 10–25 of these empty sensor-holders are fixed to the phantom patient. During the intraoperative phase, the sensor-holder can optionally hold an optical- or magnetic-tracking sensor. When used with optical tracking, the sensor-holder can hold an active IRED tracker sensor (11$\,\times \,$7 $\times$ 5 mm, NDI Optotrak Certus) (Fig. 2(d)). When used with magnetic (EM) tracking, it can hold a 5-DOF NDI Aurora sensor (Length: 8 mm, ${\oslash }$ 0.6 mm) (Fig. 2(e)) concentrically.

2.3 Data Acquisition

In order to validate sensor optimization, a 4D-CT scan of phantom patient is acquired. The phantom model with 27 external (candidate locations) and 5 internal (target locations) markers in the balloon is placed into the CT device. During the imaging the respiration cycle is simulated by manually adjusting different air amount within the balloon using the water blaster connected to the model with a flexible tube. For the 4D-CT a scanner in Univ. Clinic for Radiology (Siemens healthcare Austria) in Medical Univ. of Innsbruck is used. The scan consists of 11 discrete time steps of a breathing cycle (see Fig. 3). Each axial CT slice (512$\,\times \,$512 px) has a thickness of 1.0 mm and the 11 discrete CT phases consist of 261 images with 0.488$\,\times \,$0.488$\,\times \,$0.488 mm pixel spacing.

Fig. 3.

Detected external markers in the first phase of the visualized 4D-CT images. Left frame is used to set parameters for the automatic marker detection algorithm and contains a marker list for the marker management. The detected markers are inserted into the marker list and visualized in the standard DICOM views (axial, sagittal, coronal) in the tiled right window. The green blobs are accepted automatically as external markers (based on given geometrical properties in the detection algorithm) and blue blobs are possible candidates to be accepted manually. The geometry view (right bottom) represents the distances of all detected markers. (Color figure online)

2.4 Automatic Marker Detection

In order to learn the respiratory cycle of the patient from the observed CT images, where the optimal sensor locations for prediction are determined, a regressor is trained to predict internal motion given this data. The precise locations of external and internal fiducials in the CT image space are detected by a GPU accelerated volumetric detection method [12].

For this purpose, each of 11 CT phases are thresholded and binarized to determine the centroid of the fiducials. On the resulting image the 3D fiducials are filtered for spherical structures using morphological opening with a spherical structuring ball element of the appropriate scale given the voxel size of the dataset and the physical dimensions of the markers. Using a geometry filter on the resulting spherical blobs best candidates are selected based on shape and size. From the best candidates the blob centroids are calculated and stored to be used in the prediction step (see Fig. 4).

Fig. 4.

Top Left: Discrete chronological 3D movement positions of all external surface fiducials during in-/exhalation over 11 timesteps (Timestep 1: Fully inhaled, Timestep 11: Fully exhaled). The location coordinates are obtained by automatic marker detection. Bottom left: The positions of one external fiducial in 4D-CT image space. Movement positions of internal fiducials in top right and bottom right.

2.5 Determining Surface Sensor Locations with GA

In order to train an accurate prediction of tumor motion from a few optimally positioned fiducials, a multi-objective genetic algorithm (GA) based feature selection method (similar to [13, 14]) is proposed.

Before the imaging step, a larger set of CT markers $\lbrace c_k \rbrace$, $k \approx 25$ are fixed on the surface at randomized candidate locations.

After imaging, all external surface marker locations are detected. This results in k time-series, each with T timesteps and 3 output dimensions (the spatial coordinates of the sensor in the CT reference frame). The target marker locations over the T timesteps yield the time-series $y \in R^{T \times 3}$.

An individual during the GA search is represented by an element of a k-dimensional binary vector $I = \lbrace 0, 1 \rbrace ^{k}$, where the nth bit represents whether the nth external sensor is used for prediction (1) or not (0).

If a marker is used, its x, y, z coordinates within the CT reference frame are added to the input coordinate set used for prediction. This yields a $3 \times p$ dimensional input feature for each time-step, where p is the number of enabled markers within the individual.

For each individual I, the fitness function is defined by a multi-objective function $F(I) = \left( F_1(I), S(I) \right)$.

The primary component is given by the weighted sum

$$ F_1(I) = E(I) + \alpha * \min (0, S(I) - K) $$

where E(I) is the average RMS error between the predicted and target locations using X as the input feature set over a 3-fold cross-validation on the T timesteps and S(I) is the number of features enabled, K is the maximum preferred number of enabled fiducials and $\alpha$ is a scaling parameter, which balances the trade-off between additional prediction error and the number of enabled fiducials. This setup leads to an optimization goal of finding the minimum achievable prediction error with as few sensors as possible, but softly punishing configurations that have more than K enabled sensors.

For each individual, the predictions are evaluated using 3 Gaussian Process Regressors (GPR) ($G_i: X \rightarrow t_i$, $i = {1,2,3}$) for each coordinate of the target, with C*SE + W where C is constant kernel, SE is squared exponential and W is white noise kernel [15]. (see Fig. 5).

The C kernel is configured with the constant value: 1.0, constant value bounds: 1e−3, 1e3, SE with length scale: 10.0, length scale bounds: 1e−2, 1e2 and W with noise level: 0.1 and noise level bounds: 1e−10, 1e$+$0.5. The Gaussian Process Regressor is configured with normalized target value without an optimizer. In GA, the parameters for population: 600, cv-proba: 0.5, mu-proba: 0.2, generation: 50, cv-independent-proba: 0.5 and mu-independent-proba: 0.05 are used.

Fig. 5.

Prediction result (green dotted lines) for target 1 (red) based on surface marker 8 (blue) from the recommended group (8 and 14) in 11 different time-steps. (Color figure online)

Table 1. Overview of the offline prediction results. The RMS is the average mean squared error of the deviations in each timestep when using the recommended surface sensors.

3 Results

Table 1 represents the prediction results for 5 internal target markers in the balloon using automatically recommended surface marker list. Each input marker in the recommended surface sensor group is processed with the listed individual target respectively. The best result is obtained from the C*SE + W of the GPR algorithm. The prediction is performed with a raw data for each surface and target marker respectively.

4 Discussion

In many treatments, where respiratory motion prediction and tracking is a necessary approach to apply, the success of a treatment strongly depends on the accuracy of prediction, which is related to the detection and placement of the external surface marker locations on the patient.

Therefore, the quantity, location and distribution of the sensors is important and challenging during real-time tumor motion prediction. In the most clinical approaches, the location of external markers is chosen empirically; that is, operator-dependent [5]. Within this work, we examined that the distribution and selection of the best placement locations automatically and intelligently, which gives better prediction results with using less numbers of external sensors to use e.g. in the thorax or abdominal regions. The proposed method is tested through a custom built respiratory phantom model to provide the required dataset, that surrogates same breathing circle of a real patient [16] (Table 2).

Table 2. The 3D positional movement variations of the internal target markers (marker centroids) from fully exhaled to fully inhaled time-steps after 4D-CT.

Fig. 6.

Current (preoperative phase) and future (intraoperative phase) workflow of the presented work.

The primary aim of this work was to separate the preoperative and intraoperative phase based on offline prediction method using sensor-holders. (see Fig. 6). Having already optimal sensor amount and locations with the GPR and GA before intervention also optimizes the whole complex workflow of the real-time respiratory motion prediction and serves a reliable solution to this complexity. All processes in intraoperative phase workflow will be evaluated in the future on the phantom model and real patients to validate the presented approach. The OP or EM tracker will serve sensor data to be used as input for the real-time GPR.

5 Conclusion

In this work, a method is presented to automatically determine the optimal locations of the external OP/EM sensors on the patient’s surface to predict internal tumor motions with high accuracy. The preoperatively determined optimal sensor locations can be used in real-time tumor motion tracking using the sensors of NDI Optotrak Certus or NDI Aurora tracking systems.

The experiments and evaluations on the built realistic respiratory phantom model have shown that by using our method, best possible locations of surface sensors can be determined with high accuracy and serves a reliable prediction, based on the selected tumor inside the phantom. The experiments give also information the distributing of recommended sensors locations have a high correlation between the surface motion and the internal tumor motion. With this procedure, EM or OP 3D measurement technologies can be used for real-time prediction and it is suitable for use in the medical environment.