Contour propagation using non-uniform cubic B-splines for lung tumor delineation in 4D-CT

Abstract

Purpose

Accurate target delineation is a critical step in radiotherapy. In this study, a robust contour propagation method is proposed to help physicians delineate lung tumors in four-dimensional computer tomography (4D-CT) images efficiently and accurately.

Methods

The proposed method starts with manually delineated contours on the reference phase. Each contour is fitted by a non-uniform cubic B-spline curve, and its deformation on the target phase is achieved by moving its control vertexes such that the intensity similarity between the two contours is maximized. Since contour is usually the boundary of lesion or tissue which may deform quite differently from the tissues outside the boundary, the proposed method treats each contour as a deformable entity, a non-uniform cubic B-spline curve, and focuses on the registration of contour entity instead of the entire image to avoid the deformation of contour to be smoothed by its surrounding tissues, meanwhile to greatly reduce the time consumption while keeping the accuracy of the contour propagation. Eighteen 4D-CT cases with 444 gross tumor volume (GTV) contours manually delineated slice by slice on the maximal inhale and exhale phases are used to verify the proposed method.

Results

The Jaccard similarity coefficient (JSC) between the propagated GTV and the manually delineated GTV is 0.885 ± 0.026, and the Hausdorff distance (HD) is $2.93\,\pm \,0.93$ mm. In addition, the time for propagating GTV to all the phases is 3.67 ± 3.41 minutes. The results are better than fast adaptive stochastic gradient descent (FASGD) B-spline method, 3D+t B-spline method and diffeomorphic Demons method.

Conclusions

The proposed method is useful to help physicians delineate target volumes efficiently and accurately.

Appendices

Appendix 1: Determination of knot vector

Let ${{\varvec{U}}}=\{u_i , i=0,1,\ldots ,n+6\}$ denote the knot vector, where $u_3 ,u_4 ,\ldots u_{n+3}$ are internal knots since they are inside the domain of definition, and $u_0 ,u_1 ,u_2 ,u_{n+4} ,u_{n+5} ,u_{n+6}$ are external knots since they are outside the domain of definition. As shown in Eq.  1, ${{\varvec{q}}}^{(0)}=\{{{\varvec{q}}}_i^{(0)} ,i=0,1,\ldots ,n\}$ is fitted by ${{\varvec{p}}}^{(0)}$, such that ${{\varvec{q}}}_i^{(0)} ={{\varvec{p}}}^{(0)}(u_{i+3})$, where ${{\varvec{q}}}_0^{(0)}$ and ${{\varvec{q}}}_n^{(0)}$ are the first and last endpoints, respectively. In order to make ${{\varvec{p}}}^{(0)}$ a closed curve and $C^{2}$continuous, the following conditions are satisfied:

$$\begin{aligned} \left\{ {\begin{array}{l} {{\varvec{q}}}_0^{(0)} ={{\varvec{q}}}_n^{(0)} \\ {{\varvec{d}}}_n^{(0)} ={{\varvec{d}}}_0^{(0)} \\ {{\varvec{d}}}_{n+1}^{(0)} ={{\varvec{d}}}_1^{(0)} \\ {{\varvec{d}}}_{n+2}^{(0)} ={{\varvec{d}}}_2^{(0)} \\ \end{array}} \right. . \end{aligned}$$

As shown in Eq. 5, the $n+1$ internal knots $u_3 ,u_4 ,\ldots u_{n+3}$ are defined by the chord length method [35], which is widely used and usually performs well.

$$\begin{aligned} \left\{ {\begin{array}{l} u_3 =0 \\ u_b =u_{b-1} +{\left| {{{\varvec{q}}}_{b-3}^{(0)} -{{\varvec{q}}}_{b-4}^{(0)} } \right| }/L,\quad b=4,5,\ldots ,n+3 \\ \\ L =\sum _{i=1}^n {\left| {{{\varvec{q}}}_i^{(0)} -{{\varvec{q}}}_{i-1}^{(0)} } \right| } \\ \end{array}} \right. \end{aligned}$$

(5)

And the six outside knots are

$$\begin{aligned} \left\{ {\begin{array}{l} u_0 =u_n -1 \\ u_1 =u_{n+1} -1 \\ u_2 =u_{n+2} -1 \\ u_{n+4} =u_4 +1 \\ u_{n+5} =u_5 +1 \\ u_{n+6} =u_6 +1 \\ \end{array}} \right. . \end{aligned}$$

When the knot vector ${{\varvec{U}}}$ is defined, $N_{k,3} (u),k=0,1,\ldots ,n+2$ can be defined by the de Boor–Cox recursion formula [36], namely,

$$\begin{aligned} \left\{ {\begin{array}{l} N_{k,0} (u)=\left\{ {\begin{array}{l} 1\quad \mathrm{if} u\in [u_k ,u_{k+1} ), \\ 0\quad \mathrm{otherwise}. \\ \end{array}} \right. \\ N_{k,h} (u)=\frac{u-u_k }{u_{k+h} -u_k }N_{k,h-1} (u)+\frac{u_{k+h+1} -u}{u_{k+h+1} -u_{k+1} }N_{k+1,h-1} (u)\quad \mathrm{if}h\ge 1 \\ \frac{0}{0}=0 \\ \end{array}} \right. . \end{aligned}$$

Fig. 11

Schematic diagram of how to derive the propagated contour on the z-th slice, where u = 0.3, $\varepsilon$ = 0.25, K = 4

Appendix 2: Slice contours derived from a set of propagated contours

Assume ${{\varvec{p}}}_s^{(0)}$ and ${{\varvec{p}}}_{s+1}^{(0)}$ are the contours on the s-th and (s+1)-th slice of the reference phase, respectively, and their propagated contours are ${{\varvec{p}}}_s$ and ${{\varvec{p}}}_{s+1}$, respectively. Search a point ${{\varvec{p}}}_{s+1} (u)$ on the contour ${{\varvec{p}}}_{s+1}$ such that it is the closest point to ${{\varvec{p}}}_s (0)$. Let K be a sufficiently large positive integer, $\varepsilon =1/K,k=0,1,2,\ldots K-1$, then K line segments between the points ${{\varvec{p}}}_s (\varepsilon k)$ and points ${{\varvec{p}}}_{s+1} (u+\varepsilon k\hbox { }\bmod \hbox { }1)$ are generated, where mod is the modulo operator. Assume that the target tumor is composed of M contours. Then for every two adjacent contours, K line segments are generated in accordance with the above method, and $K(M-1)$ line segments are generated in total. The $K(M-1)$ line segments intersect the z-th slice plane at a series of intersections which are taken as the sampling points, and the contour on the z-th slice is obtained by fitting a non-uniform cubic B-spline curve to these sampling points. The schematic diagram of how to derive the propagated contour on the z-th slice is shown in Fig. 11.

Fig. 12

An example of keeping the derived contour closed. a The derived contour is unclosed, and b the closed contour is obtained by replacing line segments with rays. The green curves are manually delineated contours, and the red curves are delineated by the proposed method

Near the top or bottom of the tumor, sometimes the derived contours are unclosed, as shown in Fig. 12a. This can be explained in Fig. 11: if ${{\varvec{p}}}_s$ did not exist, i.e., ${{\varvec{p}}}_{s+1}$ is the highest contour of the tumor, then some line segments between the points ${{\varvec{p}}}_s (\varepsilon k)$ and points ${{\varvec{p}}}_{s+1} (u+\varepsilon k\hbox { }\bmod \hbox { }1)$ would not exist, which would lead to a unclosed contour on the z-th slice. The solution is to replace the line segments between ${{\varvec{p}}}_{s\hbox {+}2} (u+\varepsilon k\hbox { }\bmod \hbox { }1)$ and ${{\varvec{p}}}_{s\hbox {+}1} (\varepsilon k)$ with rays whenever ${{\varvec{p}}}_{s+1}$ is the highest contour. Then the closed contour is obtained as shown in Fig. 12b. A similar strategy is also applied to the lowest contour.