Keywords

1 Introduction

The breast tumor surgery and the following radiotherapy are usually planned using the computed tomography (CT). This fact can be exploited to use the image registration techniques to improve the estimation of the radiation dose margins. However, an alignment of the pre-operative to the post-operative scans is a challenging task. The tumor resection introduces the problem of missing data which leads to several difficulties which are not present in the typical image registration procedures.

Related Work: The majority of recent research about the image registration was dedicated to ensuring that the calculated deformation field is a diffeomorphism. Since the diffeomorphic deformation field is inherently invertible and smooth, it is desirable for majority of the image registration tasks. However, for the tumor resection problem, the deformation field nearby the structure of interest is obviously not invertible and not smooth. There were much less research about the image registration with missing data compared to ensuring the diffeomorphic properties. One of the first works used the thin-plate splines to estimate the brain tumor propagation [1]. However, the author assumed that there is still a partial correspondence between the structures which could be used as control points. For the full tumor resection, it is not true. An interesting work about algorithm based on local affine transformations was presented in [2]. Nonetheless, the main assumption was that the structure of interest is not missing, only its surroundings. The TV-L1 was a promising algorithm [3], however it turned out that its global optimization procedure makes it fail during the tumor bed propagation in the CT. The reason for this was a low influence of the similarity cost gradient nearby the resected tumor. Moreover, the deformation can be very large so an attempt to use local realization of this technique failed. There were even attempts to localize the tumor bed using the rigid registration based on surgical clips which is inherently wrong because the rigid transformation is an isometry preserving the tumor volume [4, 5]. Interesting works were introduced in the context of missing data as a result of the brain tumor resection. A fully automatic method based on the level-set segmentation of intensity disagreements and anisotropic diffusion filter modeling the resection area as a diffusion sink was proposed in [6]. An improved Demons algorithm based on a fourth dimension in order to separate removed tissues from others was proposed in [7]. There was also a work related to an atlas-based segmentation of brain images which used modified bijective Demons to model the tumor growth process with an assumption about a radial growth of the tumor [8]. The model is similar to our method with the difference about ability to reconstruct larger deformations. Our method still uses the Demons force inside the tumor which makes is possible to reconstruct, usually large, breast deformations.

Contribution: In this work, we propose a simple but very effective regularization term for the greedy version of the Demons algorithm. The proposed method is based on the a priori knowledge about the cancer resection. Since we know that the tumor is resected, its volume in the target image significantly decreases. The proposed regularization uses this knowledge to promote the tumor volume reduction. We evaluate the proposed method using artificially generated data, artificial deformations resembling the tumor resection and real, medical data representing women with the breast cancer before and after the tumor surgery. We show that the proposed method greatly decreases the breast tumor volume and improves the tumor bed localization.

2 Methods

The proposed method is based on the greedy, implicit version of the Demons algorithm. For the further reference and the comparison purposes, we present a shortened algorithmic summary of the compositive, symmetric Demons in Algorithm 1.

figure a

The proposed method consists of an empirical, greedy regularization term which enforces the tumor contraction and a scheme to automatically determine an optimal regularization step. Its algorithmic summary is presented in Algorithm 2. The presented algorithm requires a binary mask of the segmented tumor as an additional parameter. The binary mask is usually obtained during the surgery planning and therefore is easily available. The tumor can be segmented manually or automatically. The method does not require any additional parameter tunning (only the smoothing sigmas \(\sigma _d, \sigma _f\) must be chosen). The regularization is added as an additional greedy step during the Demons algorithm, between the deformation field update and the diffusion regularization.

The direction structure \(\mathbf {d_s}\) can be defined as a binary mask with a single point or a structure representing the expected tumor bed shape. Its localization is defined relative to the initial tumor position and transformed during each iteration. In the experiments performed, the direction structure is defined as a center of mass of the tumor because both the source and the target are acquired for the same patient position. Conceptually, it can be compared to the mean shift algorithm. During each iteration, both the tumor and the structure are being transformed, until convergence. In practice, the direction structure is defined using the a priori knowledge, usually with help of the surgical protocol.

The regularization term \(\mathbf {c}\) (a vector field) is based on the direction of a difference between coordinates of the direction structure and the image domain (a grid) with the minimum distance to a given coordinate. It is simply a coordinates subtraction, which can be defined as:

$$\begin{aligned} \mathbf {c} = \arg \min (D(\mathbf {d_s}, \mathbf {T_m} \circ \mathbf {u})) - \text {Id}_{\mathrm{u}}, \end{aligned}$$
(1)

where \(\hbox {Id}_{\mathrm{u}}\) denotes the image domain, the \(\mathbf {T_m} \circ \mathbf {u}\) is the transformed tumor and D denotes an Euclidean distance calculation. What is important, the calculation needs to be performed only for the segmented tumor. Then, the vector field \(\mathbf {c}\) is normalized to unit length, which is later compensated by the tumor volume and the average velocity field vector length. The distance calculation must include the knowledge about physical voxel size which is, for clarity, omitted in the algorithm description.

figure b

The optimal regularization step \(\kappa \) controls the tumor contraction size during each iteration. It consists of two steps. Firstly, an average length of velocity vectors s is being calculated. If the velocity field magnitude is significant, the contraction term should have lower influence. It is necessary because enforcing strong contraction during the initial alignment phase leads to medically not reliable deformations. Secondly, the remaining tumor volume is being calculated to make the algorithm more stable. The influence of the regularization term linearly depends on the current tumor volume. The final regularization step is defined as:

$$\begin{aligned} \kappa = \frac{1}{s} (\text {volume}(\mathbf {T_m} \circ \mathbf {u}) - \text {expected}\_\text {volume}), \end{aligned}$$
(2)

where s is the average length of the velocity vectors. Usually, the expected volume after the resection is zero. Adding this term makes to possible to model less idealized scenario, e.g. when the cavity area it not completely filled with the surrounding tissues. However, it is not a common event because the time between the breast tumor surgery and the radiotherapy planning usually is long enough.

The proposed method can be compared to an explicit volume regularization where the gradient of tumor volume with respect to the deformation field is being calculated. However, this method is extremely computationally inefficient. We have compared the results obtained by both approaches and they were very similar, both for the tumor and the tumor bed. However, the proposed method is computationally efficient. It increases the overall Demons algorithm computation time by about 4% using the CPU implementation which is absolutely acceptable. For the explicit volume regularization the computation time depends strongly on the tumor volume and for an average tumor size increases the computation time by about an order of magnitude.

3 Experiments and Results

Three experiments were performed to show the influence of the proposed algorithm on the tumor bed localization. The proposed method was compared to the compositive, symmetric Demons because our previous study has shown that, among the state-of-the-art methods (B-Splines FFD, Demons, TV-L1, and others), the Demons method provided the best results [9]. The data spatial resolution was equal to 0.97 mm x 0.97 mm x 3 mm. The \(\sigma _d, \sigma _f\) were set to 0.5 mm and 3.0 mm respectively. The proposed method is abbreviated as RB Demons (Resection-based Demons).

3.1 Artificial Data

The first experiment used an artificially generated data. This experiment was necessary because the ground-truth about tumor bed localization is unknown for the real, medical data. The artificial data made it possible to state an exact tumor bed localization and shape in both the source and the target. We created 10 artificial cases with breast tumors in different localizations and with different shapes. The artificial cases were synthesized based on the real CT data. Small, big, symmetric, asymmetric convex and concave tumors were introduced. Then, the target image was created with the explicitly stated tumor bed localization and deformed using a known deformation based on the B-Splines transformation. We assumed that the volume of surrounding soft tissues in the source data and the volume of tumor bed in the target data is equal. The artificial data were synthesized using the ASTRA Toolbox [10]. The registration results were evaluated using the relative tumor volume and the Hausdorff distance between the ground-truth and transformed tumor beds. The relative tumor volume (defined as the ratio of the transformed tumor volume to the initial tumor volume) was used because we know that in the ideal case the tumor volume in the target volume is equal to zero. The Hausdorff distance was applied because it directly shows the potential improvement for the radiotherapy planning where the maximum margins are a crucial factor. The results are presented in Fig. 1. The proposed method models the tumor resection and significantly decreases the maximum margin. This confirms that ensuring the tumor resection improves the tumor bed localization. The margin is still above 8 mm because the artificial data uses the same intensity value for the tumor bed and other surrounding soft tissues. This assumption seems valid because in the real CT data the tumor bed in the target is indistinguishable from other soft tissues. As a consequence, the alignment of the tumor bed is not driven by the intensity difference but the proposed regularization term.

Fig. 1.
figure 1

The relative tumor volume and the Hausdorff distance for the artificially generated breast tumor CT data. Initial - the original source and target. Demons - the compositive, symmetric Demons. RB Demons - the proposed method.

Full size image

3.2 Artificial Deformations

The second experiment was based on an artificial but real-like deformation fields resembling the tumor resection process. The applied deformation fields were proposed in [9]. The applied deformations significantly decreased the tumor volume and smoothly deformed the surrounding soft tissues. However, please notice that a real tumor resection cannot be modeled by the artificial deformation field because then it should point outside the image. The applied vector lengths and the calculated RMSE for both the original Demons and the proposed RB Demons are shown in Fig. 2a. The RMSE for the proposed method is lower and the relative improvement on average is equal to 15.25%.

3.3 Real Data

The final experiment used a real, medical data. We acquired 20 CT scans acquired before the breast cancer surgery, and after, during the radiotherapy planning. The tumors were manually segmented by a medical expert with more than 20 years experience in the breast cancer radiotherapy. The images were firstly rigidly registered based on bones segmentation, SIFT and RANSAC algorithms. For the real data, the ground-truth about the tumor bed is unknown. Therefore, a quantitative assessment of the tumor bed propagation is impossible. However, we can still use the relative tumor volume reduction to evaluate the calculated deformation fields. If the method is unable to model the complete tumor resection, it is certainly incorrect. The calculated relative tumor volumes are shown in Fig. 2b. The original Demons algorithm is incorrect (for one case it even doubled the tumor volume) while the proposed method resembled all the resections well. Moreover, a visual assessment is crucial for the correctness evaluation. An example of the tumor propagation using different methods is shown in Fig. 3. An exemplary 3D visualization of the tumor propagation is shown in Fig. 4.

Fig. 2.
figure 2

The RMSE [mm] for artificially applied deformation fields and the relative tumor volume [%] for real data using different registration methods. The applied deformation fields resemble the tumor resection process. The rigid registration is shown to present the influence of warping error which is negligible. Please note that (a) and (b) are from different experiments.

Full size image
Fig. 3.
figure 3

An exemplary visualization of the propagated tumor in the CT target using three image registration algorithms. The data is shown on the sagittal, coronal and transversal planes respectively. The colors indicate following methods: - Rigid Registration - Original Demons - Proposed RB Demons. Please note that for the rigid registration the tumor is not even inside the body. (Color figure online)

Full size image
Fig. 4.
figure 4

An exemplary 3D visualization of the propagated tumor in the CT target for real, medical data. The colors indicate following methods: - Rigid Registration - Original Demons - Proposed RB Demons. Please note the significant volume decrease of the tumor using the proposed method. (Color figure online)

Full size image

4 Discussion and Conclusion

The experiments presented a significant volume reduction for the propagated tumors with different sizes and shapes. The evaluation based on the tumor volume reduction is able to show just the evaluated algorithm wrongness. A volume greater than zero clearly indicates the incorrectness of the registration method. The proposed method modeled a complete resection of tumors for both the artificial and real data while the original Demons failed. The proposed algorithm not only improved the tumor propagation but also the tumor bed localization which is the structure of interest for the radiotherapy planning.

The Hausdorff distance shows the potential improvement on the radiotherapy. The 1.5 mm decrease of Hausdorff distance (compared to the original Demons algorithm), for an average tumor size, can lead to the reduction of the irradiated volume by 32.5%. This value is significant and can decrease the risk of both the secondary carcinogenesis and tumor recurrence.

Moreover, the proposed method is independent of the registration forces calculation which makes it useful in the multi-modal registration. Therefore, this technique can be easily extended to multi-modal problems, like e.g. a real-time imaging system for the tumor surgery. This method, without any modifications, can be used in the MIND-based [11] or the NMI-based Demons [12].

The proposed method can be extended to incorporate the mechanical properties of the tissues being deformed. A biomechanical model of the breast will be used in the further study to better model the propagation destination for large deformations. Moreover, the calculated vector directions can be changed to more complex shape to make it useful for complex, concave and anisotropic shapes.

To conclude, we proposed a greedy, empirical regularization term for the Demons algorithm which ensures an appropriate tumor resection. The proposed method not only improved the tumor propagation but also the tumor bed localization. In the further research, we will incorporate this technique into multi-modal registration e.g. for the real-time surgery based on MRI-USG scans. Moreover, we will extend this technique to a structure-based contraction which should improve the tumor bed localization even more. We will perform a research about reliable, quantitative evaluation methods for image registration algorithms dedicated to the missing data problem.