A kernel-based method for markerless tumor tracking in kV fluoroscopic images

Generally, in a kV fluoroscopic image, soft tissues such as blood vessels and tumors have a lower contrast compared to those of hard tissues, such as bones, and their intensity range probably varies with the x-ray energy and irradiation position. Figure 2(a) shows a 16-bit fluoroscopic image in which a tumor target is located at the center of the image. Due to the low contrast, we see that the tumor is relatively indistinguishable and can hardly be tracked by using conventional visual tracking approaches. The normalized histogram of the image is shown in figure 2(b). Note that the dynamical rang of the image intensity is extremely narrow for a 16-bit image whose intensity ranges from 0 to 65535.

Zoom In Zoom Out Reset image size

A desirable high-contrast image should cover a wide range of the intensity scale and the distribution of the image intensities should not be too far from a uniform distribution. Histogram equalization is one of the intensity transformations that can automatically improve the image contrast based only on information available in the histogram of the input image (Gonzalez and Woods 2008).

Given an input image with intensity levels in the range [0, L − 1] (e.g. L = 2¹⁶ for a 16-bit image), the normalized histogram is a discrete function, given by p(r_k) = n_k/N, where r_k is the kth intensity value, n_k is the number of pixels in the image with intensity r_k and N is the total number of pixels in the images. The histogram equalization is used to map each pixel in the input image with intensity r_k into a corresponding pixel with intensity s_k in a new range $\left[0,{{L}_{\text{new}}}-1\right]$ using the transformation T(r_k), given by

$$\begin{eqnarray}&&{{s}_{k}}=T\left({{r}_{k}}\right)=\left({{L}_{\text{new}}}-1\right)\underset{j=0}{\overset{k}{\sum}} p\left({{r}_{k}}\right)=\frac{\left({{L}_{\text{new}}}-1\right)}{N}\underset{j=0}{\overset{k}{\sum}} {{n}_{j}},k=0,1,\ldots ,L-1\end{eqnarray} \tag{ 1 }$$

For reducing the computational cost in the subsequent steps, we convert the input image into an 8-bit image by setting ${{L}_{\text{new}}}=256$. In practice, s_k probably have fractions because they were generated by summing probability values, so we round them to the nearest integer. Figure 3 shows the intensity transformation function T(r_k) for the image in figure 2(a).

Zoom In Zoom Out Reset image size

Figure 4(a) shows the output image obtained from the histogram equalization. Compared with the original image shown in figure 2(a), the output image has higher dynamic range and the tumor target becomes more distinguishable. Figure 4(b) shows the normalized histogram of the output image that exhibits a large variety of gray tones.

Zoom In Zoom Out Reset image size

In practice, given a kV fluoroscopic image sequence, we can use the first frame to generate a lookup table according to equation (1) for the histogram equalization and apply this lookup table to subsequent frames for improving the image contrast.

In order to track a tumor target in a fluoroscopic image sequence, the tumor target should be chosen in advance. Since the contrast of the observed image has been enhanced in the previous step, the tumor target can more easily be delineated in a rectangle by clinicians. In figure 4, a red rectangle covering the tumor area is delineated by a clinician. The tumor target is then represented by a feature vector $\mathbf{\hat{q}}$ which consists of m elements ${{\left\{{{\hat{q}}_{u}}\right\}}_{u=1,\ldots ,m}}$, where ${{\hat{q}}_{u}}$ denotes the quantity of the u-th element. The feature vector is generated based on a weighted histogram described as follows.

Suppose that the tumor target is delineated within a rectangle $T\left(\mathbf{x}_{i}^{*}\right)$ of size n pixels, where T is the intensity of a pixel located at ${{\left\{\mathbf{x}_{i}^{*}=\left(x_{i}^{*},y_{i}^{*}\right)\right\}}_{i=1,2,\ldots ,n}}$, and that the center of the rectangle is 0. An m-bin histogram index function is defined by

$$\begin{eqnarray}&&b\left(\mathbf{x}_{i}^{*}\right)=\frac{mT\left(\mathbf{x}_{i}^{*}\right)}{255}+1\end{eqnarray} \tag{ 2 }$$

where ⌊·⌋ is a ceiling function that returns the smallest integer greater than or equal to $mT\left(\mathbf{x}_{i}^{*}\right)/255$. The function $b\left(\mathbf{x}_{i}^{*}\right)$ indicates the index of the histogram bin corresponding to the intensity of the pixel at location $\mathbf{x}_{i}^{*}$. We also define a profile function of a Gaussian kernel, denoted by k(x), given by

$$\begin{eqnarray}&&k\left(x\right)={{\text{e}}^{-\frac{x}{2}}}\end{eqnarray} \tag{ 3 }$$

The details of the profile function of a kernel can be found in Comaniciu et al (2003).

The feature vector representing the tumor target is defined as

$$\begin{eqnarray}&&{\hat q_u} = C\mathop \sum \limits_{i = 1}^n k\left( {{{\left\| {\frac{{{\bf{x}}_i^*}}{{\bf{h}}}} \right\|}^2}} \right)\delta \left[ {b({\bf{x}}_i^*) - u} \right]\,\,\,\,u = 1,2, \ldots ,m.\end{eqnarray} \tag{ 4 }$$

where $\left\| {{\bf{x}}_i^*/{\bf{h}}} \right\|$ is the normalized distance between $\mathbf{x}_{i}^{*}$ and the center of $T\left(\mathbf{x}_{i}^{*}\right)$, $\mathbf{h}=\left({{h}_{x}},{{h}_{y}}\right)$ is the kernel bandwidth that normalizes the coordinate $\mathbf{x}_{i}^{*}$ within a unit circle, δ is the Kronecker delta function defined by

$$\begin{eqnarray}\delta (x) = \left\{ {\begin{array}{*{20}{c}} {1,}&{{\rm{if}}\,x = 0}\\ {0,}&{{\rm{otherwise}}} \end{array}} \right.\end{eqnarray} \tag{ 5 }$$

and C is the normalization constant given by

$$\begin{eqnarray}&&C = \frac{1}{{\sum\nolimits_{i = 1}^n {k\left( {{{\left\| {\frac{{{\bf{x}}_i^*}}{{\bf{h}}}} \right\|}^2}} \right)} }}\end{eqnarray} \tag{ 6 }$$

Figure 5(a) shows the tumor target of size 80 × 90 pixels delineated by the red rectangle shown in figure 4(a). Figure 5(b) is the corresponding Gaussian kernel obtained from equation (3) where h_x = 40 and h_y = 45. Figure 5(c) shows the feature vector ${{\left\{{{\hat{q}}_{u}}\right\}}_{u=1,\ldots ,m}}$ representing the tumor target where m is 25.

Zoom In Zoom Out Reset image size

Here, the feature vector ${{\left\{{{\hat{q}}_{u}}\right\}}_{u=1,\ldots ,m}}$ can be considered as a weighted intensity probability density function (PDF) of the tumor target. The Gaussian kernel function gives a larger weight for the pixel whose location is closed to the center of the target, while gives a smaller weight for the pixel whose location is far away from the center. Therefore, the feature vector not only contains intensity information of tumor target, but also contains spatial information.

Tracking the tumor target in the subsequent frame is to find the location of a target candidate whose feature vector denoted by $\mathbf{\hat{p}}$, is closest to the feature vector of the tumor target $\mathbf{\hat{q}}$.

Let $T\left({{\mathbf{x}}_{i}}\right)$ be a target candidate, centered at $\mathbf{y}$ in the current frame, where i = 1, 2, ..., n. Using the Gaussian kernel defined by equation (3), the feature vector of this target candidate, denoted by $\mathbf{\hat{p}}\left(\mathbf{y}\right)={{\left\{{{\hat{p}}_{u}}\left(\mathbf{y}\right)\right\}}_{u=1,\ldots ,m}}$, is given by

$$\begin{eqnarray}&&{\hat p_u}({\bf{y}}) = C\mathop \sum \limits_{i = 1}^n k(\frac{{{\bf{y}} - {{\bf{x}}_i}}}{{\bf{h}}}{^2})\delta [b({{\bf{x}}_i}) - u],\,\,\,\,u = 1,2, \ldots ,m.\end{eqnarray} \tag{ 7 }$$

Tracking the target can be formulated by finding the coordinate $\mathbf{y}$ where the similarity between two feature vectors $\mathbf{\hat{q}}$ and $\mathbf{\hat{p}}\left(\mathbf{y}\right)$ reaches the maximum. In this paper, we employ the Bhattacharyya coefficient as the similarity metric to measure the similarity between the target model $\mathbf{\hat{q}}$ and the target candidate $\mathbf{\hat{p}}\left(\mathbf{y}\right)$ (Comaniciu et al 2003). The Bhattacharyya coefficient, denoted by $\rho \left[\mathbf{\hat{p}}\left(\mathbf{y}\right),\mathbf{\hat{q}}\right]$, is defined by

$$\begin{eqnarray}&&\rho \left[\mathbf{\hat{p}}\left(\mathbf{y}\right),\mathbf{\hat{q}}\right]=\underset{u=1}{\overset{m}{\sum}} \sqrt{{{{\hat{p}}}_{u}}\left(\mathbf{y}\right){{{\hat{q}}}_{u}}}\end{eqnarray} \tag{ 8 }$$

The Bhattacharyya coefficient can be interpreted as the correlation between the vectors $\left(\sqrt{{{{\hat{p}}}_{1}}},\ldots ,\sqrt{{{{\hat{p}}}_{m}}}\right)$ and $\left(\sqrt{{{{\hat{q}}}_{1}}},\ldots ,\sqrt{{{{\hat{q}}}_{m}}}\right)$. Since the feature vector of the target candidate $\mathbf{\hat{p}}\left(\mathbf{y}\right)$ carries continuous spatial information, the Bhattacharyya coefficient between $\mathbf{\hat{p}}\left(\mathbf{y}\right)$ and $\mathbf{\hat{q}}$ is also a smooth function of $\mathbf{y}$. Suppose that the center of the tumor target in the previous frame is ${{\mathbf{\hat{y}}}_{0}}$. The Bhattacharyya coefficient in equation (8) can be approximated by the following Taylor expansion around the location ${{\mathbf{\hat{y}}}_{0}}$ :

$$\begin{eqnarray}&&\rho [{\bf{\hat p}}({\bf{y}}),{\bf{\hat q}}] \approx \frac{1}{2}\mathop \sum \limits_{u = 1}^m \sqrt {{{\hat p}_u}({{{\bf{\hat y}}}_0}){{\hat q}_u}} + \frac{C}{2}\mathop \sum \limits_{u = 1}^n {w_i}k\left( {{{\left\| {\frac{{{\bf{y}} - {{\bf{x}}_i}}}{{\bf{h}}}} \right\|}^2}} \right)\end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray}&&{{w}_{i}}=\underset{u=1}{\overset{m}{\sum}} \sqrt{\frac{{{{\hat{q}}}_{u}}}{{{{\hat{p}}}_{u}}\left({{{\mathbf{\hat{y}}}}_{0}}\right)}}\delta \left[b\left({{\mathbf{x}}_{i}}\right)-u\right]\end{eqnarray} \tag{ 10 }$$

In equation (9), the first term is independent of $\mathbf{y}$ and the second term represents a density estimation computed with the kernel profile k (·) at $\mathbf{y}$ in the current image.

According to the theorem of density gradient estimation (Comaniciu et al 2003), the maximization of the Bhattachayya coefficient $\rho \left[\mathbf{\hat{p}}\left(\mathbf{y}\right),\mathbf{\hat{q}}\right]$ can be solved by the following mean-shift algorithm.

Given a tumor target whose center is located at ${{\mathbf{\hat{y}}}_{0}}$ in the previous frame, its feature vector $\mathbf{\hat{q}}={{\left\{{{\hat{q}}_{u}}\right\}}_{u=1,\ldots ,m}}$ is computed according to equation (4). In the current frame, the tracking process consists of the following three steps:

Step 1: Compute the feature vector $\mathbf{\hat{p}}\left({{\mathbf{\hat{y}}}_{0}}\right)={{\left\{{{\hat{p}}_{u}}\left({{\mathbf{\hat{y}}}_{0}}\right)\right\}}_{u=1,\ldots ,m}}$ of a target candidate located at ${{\mathbf{\hat{y}}}_{\mathbf{0}}}$ in the current frame and the Bhattachayya coefficient between $\mathbf{\hat{p}}\left({{\mathbf{\hat{y}}}_{0}}\right)$ and $\mathbf{\hat{q}}$

$$\begin{eqnarray}&&\rho \left[\mathbf{\hat{p}}\left({{\mathbf{\hat{y}}}_{\mathbf{0}}}\right),\mathbf{\hat{q}}\right]=\underset{u=1}{\overset{m}{\sum}} \sqrt{{{{\hat{p}}}_{u}}\left({{{\mathbf{\hat{y}}}}_{0}}\right){{{\hat{q}}}_{u}}}\end{eqnarray} \tag{ 11 }$$

Step 2: Compute the weight {w_i}_{i = 1,...,n} according to equation (10).

Step 3: Compute the new location of the target

$$\begin{eqnarray}&&{{\bf{\hat y}}_1} = \frac{{\sum\nolimits_{i = 1}^n {{{\bf{x}}_i}{w_i}g} \left( {{{\left\| {\frac{{{{{\bf{\hat y}}}_0} - {{\bf{x}}_i}}}{{\bf{h}}}} \right\|}^2}} \right)}}{{\sum\nolimits_{i = 1}^n {{w_i}g} \left( {{{\left\| {\frac{{{{{\bf{\hat y}}}_0} - {{\bf{x}}_i}}}{{\bf{h}}}} \right\|}^2}} \right)}}\end{eqnarray} \tag{ 12 }$$

where g(x) = − k'(x). In this paper, since k(x) is the profile function of a Gaussian kernel, g(x) is then given byg(x) = 0.5e^−x/2.

Step 4: If $\left\| {{{{\bf{\hat y}}}_1} - {{{\bf{\hat y}}}_0}} \right\| <$, stop the iterations, otherwise set ${{\mathbf{\hat{y}}}_{0}}~\leftarrow ~{{\mathbf{\hat{y}}}_{1}}$ and return to Step 1.

Once the iteration is stopped in current frame, we replace the tumor target by using the tracking result and continue the tracking process to the next frame. Therefore, we only need to decide the tumor target in the first frame during whole tracking process. In practice, it is difficult to delineate a tumor target during the treatment. Therefore, we consider that this step can be conducted in patient positioning.

To demonstrate the efficiency of the mean-shift algorithm, figure 6 gives an example to show the tracking process. In figure 6(a), the red rectangle whose center is ${{\mathbf{\hat{y}}}_{0}}$ indicates the tumor's location in the previous frame and the green rectangle shows an estimated location of the tumor target obtained by the mean-shift iteration processing. Figure 6(b) shows the Bhattacharyys similarity surface and the mean-shift iteration that converges at the maximum of the similarity surface within about 20 iterations. In this figure, the similarity surface reaches its maximum at the coordinates (−43, −21) which indicates the values of tumor motion along the vertical and horizontal directions.

Zoom In Zoom Out Reset image size

Figure 9 shows the contrast-enhanced images obtained from the histogram equalization described in section 2.1. Compared with the original images in figure 7, the tumor target delineated by the red rectangles are more distinguishable. Specially, for the third (Patient 3) and fourth (Patient 4) images, the visibilities of tumor targets have been significantly improved.

Zoom In Zoom Out Reset image size

To assess the contrast enhancement numerically, we employ the image entropy to measure the degree of randomness of image. The image entropy is defined by

$$\begin{eqnarray}&&H=-\underset{k=0}{\overset{L-1}{\sum}} p\left({{s}_{k}}\right) {{\log}_{2}}~p\left({{s}_{k}}\right)\end{eqnarray} \tag{ 13 }$$

where L and p(s_k) are defined in section 2.1. A high entropy implies better contrast in image, while low entropy implies low contrast in image. Figure 10 shows the comparison of the entropies between the original images and the enhanced images. For the first and second patients, the image entropies are increased by about 20% . For the third and fourth patients, the image entropies are increased by about 100% . These results demonstrate that the histogram equalization is capable of increasing the image entropies effectively.

Zoom In Zoom Out Reset image size

In figure 9, four tumor targets within the red rectangles are delineated manually in the first frame before the target tracking. These tumor targets are represented by 25-element feature vectors. We then apply the tracking algorithm to localize the tumor target in the subsequent frames.

Figures 11(a)–(d) show the tracking results in four image sequences. For comparison, we also plot the ground truths of the tumor motion. The ground truths are the mean of manually localized positions generated by three observers. We can see that our tracking results are in good agreement with the ground truths.

Zoom In Zoom Out Reset image size

We also implemented four basic template matching (TM)-based tracking methods by using different similarity metrics: mean squared error (MSE), correlation ratio (CR), correlation coefficient (CC) and mutual information (MI) for comparison. All the template matching methods are performed on the same contrast-enhanced image sequences. In addition, in the template matching-based methods, we limited the searching region within a pre-defined region of size 100 × 100 pixels, which includes the tumor positions at all breathing phases. Here, the template matching methods for comparison do not represent other previously published methods, but are instead a basic implementation of template matching. Table 3 summarizes the absolute means of errors (mm) and the standard deviations of errors (mm), receptively denoted by $\bar{e}$ and σ, between the tracking results and the ground truths. For each image sequence, we calculate the errors along two directions corresponding to the image acquisition angle. We also summarize the averages of errors along three directions. This table demonstrates that the accuracy of our tracking method is higher than those of the template matching methods.

As mentioned in section 2.2, the proposed method is to use a histogram-based feature vector to represent a tumor target. If the intensities of a tumor target have a narrow dynamical range, for example, when the tumor overlaps with heart or spine, the accuracy of our method will be degraded inevitably. In future work, we will focus on tracking the tumor when it overlaps with the heart or spine. We plan to utilize an adaptive histogram equalization or motion enhancement technique (Berbeco et al 2006) as a preprocessing technique to enhance the tumor target and extend the validity of the proposed method.

Table 3. Absolute mean of error $\left(\bar{e}\right)$ and standard deviation of errors (σ) (mm) of tumor tracking.

Patient No.		Patient 1		Patient 2		Patient 3		Patient 4		Average
Direction		SI	LR	SI	LR	SI	AP	SI	AP	SI	LR	AP
Moving range		16.9	10.5	2.9	1.3	12.7	3.8	22.4	9.6	1.5	1.0	1.4
Method 1	$\bar{e}$	1.4	1.3	1.0	0.7	2.1	1.3	1.4	1.5	2.6	1.3	3.6
	σ	2.1	1.8	1.0	0.7	2.6	2.1	3.1	3.2	2.2	1.3	2.7
Method 2	$\bar{e}$	1.9	1.4	0.6	0.4	2.0	1.1	2.8	2.1	1.6	0.9	1.6
	σ	2.3	1.9	1.2	0.9	3.4	2.1	2.9	4.5	2.5	1.3	3.3
Method 3	$\bar{e}$	1.0	0.7	0.3	0.3	2.1	1.2	2.9	2.0	1.6	0.5	1.6
	σ	2.0	1.2	0.5	0.9	3.6	2.8	3.0	3.2	2.3	1.1	3.0
Method 4	$\bar{e}$	1.7	1.1	0.7	0.4	2.8	2.1	2.5	2.1	1.9	0.8	2.1
	σ	2.5	1.3	1.0	0.9	3.0	2.9	2.9	2.3	2.3	1.1	2.6
Proposed	$\bar{e}$	1.0	0.6	0.1	0.2	0.4	0.1	1.6	0.7	0.8	0.4	0.4
	σ	1.9	1.1	0.3	0.5	2.2	0.6	2.4	2.0	1.7	0.8	1.3

Note: Method 1: MSE-based TM, Method 2: CR-based TM, Method 3: CC-based TM, Method 4: MI-based TM.

As described in section 3.1, the processing of contrast enhancement is performed by using an intensity look-up table that is obtained from the first frame of an image sequence. The contrast enhancement for the subsequent frame is computationally inexpensive. In fact, the contrast enhancement for such a frame was about 3 ms in our experiments.

Therefore, we mainly focus on evaluating the computational time of the tracking process. Figure 12 shows the average computational times (s/frame) of the tracking process of the proposed method with comparison of four template matching methods. The average tracking time of the proposed method is about 45 ms per frame, while the lowest computational time of template matching method is about 200 ms per frame that is about 4 times of the proposed method. Consider that the fps of the kV imaging system in our experiments is about 15 fps (about 67 ms/frame), the proposed method is capable of providing a real-time performance in practice.

Zoom In Zoom Out Reset image size