Volume-of-change cone-beam CT for image-guided surgery

Since for many surgical scenarios, the changes imparted to the patient in any given step of the procedure may be local and relatively small, there is a significant amount of redundant information, i.e. unchanged patient anatomy, between repeated scans. This implies that, if we know how to relate the already-known information to new information in the current scan, we can highlight the regional changes and restore them from a much smaller amount of data than is needed to reconstruct an entirely new image volume, which leads us to extensively utilize available patient-specific prior knowledge for intraoperative 3D reconstructions. In this paper, we use a preoperative CT or CBCT scan as patient-specific prior knowledge to augment sparse-scan data obtained during surgery. In particular, the prior anatomical model is used to generate difference projection images, ideally, showing surgical changes only, thus connecting patient-specific prior knowledge to high-quality VOC reconstructions. Specifically, instead of performing a complete CBCT scan that involves a full radiation dose, the VOC approach attempts to reconstruct only the aspects of the region that have changed since the previous scan from a sparse set of projections. Only new information, i.e. the VOC, is reconstructed, and the result may be fused with the prior volume, to present a complete and up-to-date view of the surgical field.

To fully utilize the information available from the prior volume, an accurate 3D–2D registration between the prior volume and the observed projections is necessary. Once a sparse set of intraoperative projections are obtained during the surgery, the 3D prior volume is registered to them. Since patient anatomy contains both structures that deform such as soft tissue, and ones that remain rigid such as bone, between scans, we first obtain a bone-only prior volume by thresholding. During the registration process, digitally reconstructed radiographs (DRRs) of the prior volume are repeatedly generated and compared to the selected sparse subset of 2D intraoperative projections. To measure data similarity between DRRs and the observed projections, we used the gradient difference similarity metric because of its known robustness to the slowly-varying soft tissue background in the projections and local data inconsistency introduced by surgical tools and implants (Penney et al 1998). The gradient difference of two images is computed as

$$\begin{equation} G ( s ) = \mathop \sum \limits_{i,j} \frac{{A_v }}{{A_v + ( {I^{V{\rm diff}} ( {i,j} )} )^2 }} + \mathop \sum \limits_{i,j} \frac{{A_h }}{{A_h + ( {I^{H{\rm diff}} ( {i,j} )} )^2 }},\end{equation} \tag{ 1 }$$

$$\begin{eqnarray} &&I^{V{\rm diff}} ( {i,j} ) = \frac{{{\rm d}I^{{\rm proj}} }}{{{\rm d}i}} - s\frac{{{\rm d}I^{{\rm DRR}} }}{{{\rm d}i}}, \\&&\ms I^{H{\rm diff}} \left( {i,j} \right) = \frac{{{\rm d}I^{{\rm proj}} }}{{{\rm d}j}} - s\frac{{{\rm d}I^{{\rm DRR}} }}{{{\rm d}j}}, \end{eqnarray} \tag{ 2 }$$

where I^proj is the observed x-ray projection, I^DRR is the generated DRR, I^Vdiff and I^Hdiff are vertical and horizontal gradient differences, respectively, s is a scaling factor that considers the intensity difference between the x-ray image and the DRR, and A_v and A_h are regularizing constants, which in our experiments were the variance of the respective gradient x-ray image. Since a new set of DRRs has to be created at every iteration during the optimization, Siddon's forward projection method (Siddon 1985) was implemented on GPU for generating DRRs. The registration was optimized for pose by using the downhill simplex method (Nelder and Mead 1965).

Once the prior volume is registered to the observed 2D projections, we compute a series of difference images that highlight surgical changes. The idea of computing difference images to highlight the structure of interest is not new, and has been used in digital subtraction angiography (Heautot et al 1998, Bidaut et al 1998, Anxionnat et al 2001) where contrast-enhanced vascular structures are highlighted by subtracting mask images obtained prior to the contrast agent injection at the same image poses. Since such mask images are not generally available in most surgical procedures, we instead use the patient-specific prior volume to create DRRs that are used as mask images, thus computing difference images by subtracting DRRs from the intraoperative x-rays.

In transmission tomography, the measurement of x-ray photons after being attenuated by an object under a monochromatic x-ray assumption can be approximated as

$$\begin{equation} \bar y_i = I_{i0} {\rm exp}( - [ {{\bf A}\mu } ]_i ),\end{equation} \tag{ 3 }$$

where I_i0 is the number of unattenuated photons, µ is the digitized linear attenuation map of the object, A is the system matrix, and [Aμ]_i is the line integral of the linear attenuation map µ measured at the detector element i. Assuming an ideal detector, each measurement y_i at the detector element i is independent of measurement in neighbouring pixels and follows a Poisson distribution with the expectation value of $\bar y_i$. Under this assumption, the joint probability of this measurement is derived as

$$\begin{equation} P(y|\mu ) = \mathop \prod \limits_{i = 1}^N \frac{{\bar y_i ^{y_i } }}{{y_i !}}{\rm exp}( { - \bar y_i } )\end{equation} \tag{ 4 }$$

and the corresponding negative log-likelihood is given by

$$\begin{equation} {-}L( {y;\mu } ) = \mathop \sum \limits_{i = 1}^N \left\{ {I_{i0} {\rm exp}\left( { - \left[ {{\bf A}\mu } \right]_i } \right) + y_i \left( {\left[ {{\bf A}\mu } \right]_i } \right)} \right\}\end{equation} \tag{ 5 }$$

ignoring the constant term. Direct minimization of equation (5) yields the maximum likelihood (ML) estimate $\hat{\mu}$ which is very noisy. Alternatively, a penalized likelihood (PL) or a maximum a posteriori method regularizes the problem by applying roughness penalty to the objective function, thus reducing the noise. In this approach, the following PL objective function modified from equation (5) is minimized:

$$\begin{equation} \hat \mu = {\rm arg}\mathop {\min }\limits_\mu \{ { - L( {y;\mu } ) + \beta R(\mu )} \},\end{equation} \tag{ 6 }$$

where R(μ) is a smoothness prior and β controls the relative contribution of the likelihood and the prior terms. In order to penalize difference in the neighbouring pixels and preserve edges, we typically use the following penalty:

$$\begin{equation} R( \mu ) = \left\{ \begin{array}{@{}ll@{}} {\dsty\frac{1}{{2\delta }}( {{\rm \Psi }\mu } )^2 } & {\rm for}\quad {| {{\rm \Psi }\mu } | \le \delta } \\\ms {| {{\rm \Psi }\mu } | - \dsty\frac{\delta }{2}} & {\rm for}\quad {| {{\rm \Psi }\mu }| > \delta }, \end{array} \right. \end{equation} \tag{ 7 }$$

where Ψ is the spatial gradient operator and δ is a small constant. This penalty is similar to the 1-norm to the spatial gradient of the image estimate (μ) and also Geman prior (Geman and Reynolds 1992, Geman et al 1992) which is closely related to the total variation (TV) norm (Tang et al 2009), but is differentiable everywhere.

It has been reported that this PL reconstruction with a proper choice of smoothness prior can provide an improved reconstruction compared to TV-based CS reconstruction (Tang et al 2009). In particular when the dose level is very low, TV-based CS reconstruction results in patchy images while mitigating streak artefacts, which could be harmful. Since our goal is to achieve low-dose and high-contrast VOC reconstruction, the PL approach is an appropriate choice.

In the VOC reconstruction method, we have two measurements; one is the projections ($\bar y^{{\rm intraop}}$) obtained during the surgical procedure and the other is the DRRs ($\bar y^{{\rm prior}} )$ generated from the registered prior volume. The measured x-ray image $\bar y^{{\rm intraop}}$ is created from the intraoperative patient anatomy (μ^intraop) that contains the regional change made by the surgical procedure, but the DRR is generated from the prior volume (μ^prior) that does not involve the surgical change. However, when we compute the difference between these two datasets (after taking the logarithm) assuming no registration error, the resulting difference can be considered as a projection of the difference between the intraoperative patient volume and the prior volume, i.e. μ^VOC = μ^intraop − μ^prior. Therefore, the difference projection can be described from equation (3) as

$$\begin{equation} \bar y_i^{{\rm diff}} = I_{i0} {\rm exp}( { - [ {{\bf A}(\mu ^{{\rm intraop}} - \mu ^{{\rm prior}} )} ]_i } ) = I_{i0} {\rm exp}( { - [ {{\bf A}\mu ^{{\rm VOC}} }]_i } ),\end{equation} \tag{ 8 }$$

which is the projection of the change. Although the difference images can be described as above, the likelihood has changed because this difference image computation involves a nonlinear subtraction of Poisson random variables. However, due to the complexity of the nonlinear subtraction of Poisson random variables, we assume that it is approximately Poisson in this paper, and a more accurate noise model is currently under investigation. Under the Poisson approximation, μ^VOC can be estimated by minimizing the following objective function using the likelihood in equation (5) and the penalty in equation (7):

$$\begin{equation} \hat \mu ^{{\rm VOC}} = {\rm arg}\mathop {\min }\limits_\mu \{ { - L( {y^{{\rm diff}} ;\mu }) + \beta R(\mu )} \}.\end{equation} \tag{ 9 }$$

Once the VOC is reconstructed, we can recover the entire intraoperative volume by fusing the VOC reconstruction with the prior volume:

$$\begin{equation} \hat \mu ^{{\rm intraop}} = \hat \mu ^{{\rm VOC}} + \mu ^{{\rm prior}} .\end{equation} \tag{ 10 }$$

The PL estimation problem in equation (9) can be solved in various ways, and we used the ordered subset transmission tomography approach proposed by Erdogan and Fessler (1999).

An experimental x-ray bench that provides a precise and flexible platform for investigation of CBCT imaging performance was used for a phantom study. The bench consists of an x-ray tube (DU694 in an EA10 housing; Dunlee, Aurora, IL), a FPD (2048 × 1536 format 4030 CB detector with 0.194 mm pixel pitch and CsI:Tl converter; Varian Medical Systems, Salt Lake City, UT) and a motion control system (Compumotor 6k8 and Dynaserv G3 servo drive with a DR 1060B motor; Parker Hannifin, OH). The bench offers a wide range of source–detector motions and x-ray techniques, emulating a broad range of x-ray systems including C-arms with precise control and reproducibility of system geometry and a high degree of user control over the acquisition process. The phantom images were acquired at 100 kVp and 125 mA (500 mAs). Each image has the size of 1024 × 768 pixels with a pixel pitch of 0.388 mm after a 2 × 2 binning. The phantom imaging setup on the CBCT bench is shown in figure 2.

Zoom In Zoom Out Reset image size

A sacroplasty was simulated on a dry cadaveric sacrum phantom without soft tissue. A preoperative full sweep CBCT scan (360 cone-beam projections over a 360° arc) was first acquired from which a CBCT volume (400 × 400 × 400 voxels with a voxel size of 0.5 × 0.5 × 0.5 mm³) was reconstructed as a prior volume. Bone cement (Spineplex, Stryker Spine, Allendale, NJ) was injected into the upper portion of the sacrum to create a surgical change. After the bone cement injection, another full sweep CBCT scan was performed, and different subsets of 15, 30, 45, 60 projections that were uniformly sampled over 180° among the full scan data were extracted to generate sparse post-injection scans. We placed six stainless steel beads on the sacral surface to ensure the registration and validate the VOC reconstruction method on an ideal case, i.e. with an accurate registration and using the same imaging system at the same x-ray energy. Once the prior volume was registered to the post-injection x-ray images using these six beads, DRRs were generated from the registered prior volume and subtracted from the post-injection x-ray images to create difference images. A series of VOCs (300 × 300 × 90 voxels with a voxel size of 0.5 × 0.5 × 0.5 mm³) were reconstructed from the different subsets of 15, 30, 45, 60 difference images by the PL reconstruction method with different β values.

To assess the VOC reconstruction quality, we also reconstructed a post-injection CBCT volume (400 × 400 × 400 voxels with a voxel size of 0.5 × 0.5 × 0.5 mm³) by using all 360 post-injection projections and computed a ground truth VOC by subtracting the prior CBCT volume from the post-injection CBCT volume. Since pure subtraction of two CBCT volumes yields noisy images of the change and the bone cement has uniform density, we applied a median filter to this subtracted volume to remove quantum noise while maintaining the overall shape and edge. We compared the sparse VOC reconstructions to this ground truth using the similarity metrics described in section 4.

A prototype FPD-based mobile C-arm for high-performance CBCT (C14, Siemens Medical Solutions, Erlangen, Germany) was used for imaging cadaver specimens during simulated vertebroplasty procedures. The C-arm is based on a Siemens Powermobil and employs a 30 × 30 cm² CsI:TI FPD with an active matrix of 1536 × 1536 at 194 µm pixel size (PaxScan 3030+, Varian Imaging Products, Palo Alto, CA). The C-arm operated at the tube voltage of 120 kVp and current of 2.3 mA (230 mAs), and a total of 200 CBCT images were acquired over an angular range of 178°. Figure 6 shows the cadaver imaging setup.

Zoom In Zoom Out Reset image size

Vertebroplasty was performed on a thoracic vertebra in a cadaver torso. A pre-injection C-arm CBCT scan was first performed to get 200 pre-injection x-ray projections over a 178° arc, from which a prior volume of 512 × 512 × 512 voxels with a voxel size of 0.2930 × 0.2930 × 0.2930 mm³ was reconstructed. A surgeon carefully injected bone cement (Confidence Spinal Cement, DePuySpine, Raynham, MA) into a thoracic vertebra without moving the cadaver torso so that the CBCT scan could be repeated on the same C-arm orbit at the same subject position. After the bone cement injection, a post-injection CBCT scan was performed to get 200 post-injection x-ray projections at the same poses as the pre-injection scan and a post-injection volume with the same volume and voxel sizes to the pre-injection volume was reconstructed. Since both the cadaver and the C-arm did not move between these two scans and the C-arm image acquisition was reproducible, we assumed that these two scans were registered. Therefore, difference images were computed by directly subtracting the pre-injection x-ray images from the corresponding post-injection x-ray images.

Among 200 difference images, we uniformly sampled subsets of 10, 20, 40, 66 projections to simulate sparse intraoperative scans. VOCs of 250 × 250 × 250 voxels with a voxel size of 0.2930 × 0.2930 × 0.2930 mm³ were reconstructed by the PL method from these difference images. The reconstructed VOCs were compared to a ground truth VOC that was computed by subtracting the pre-injection CBCT volume from the post-injection CBCT volume followed by the median filtering to remove quantum noise as we did for the phantom experiment.

We have also performed vertebroplasty on a lumbar spine of the torso cadaver. In this case, the cadaver and the C-arm moved between the pre- and post-injection scans. A preoperative CT volume (512 × 512 × 232 voxels with a voxel size of 0.9062 × 0.9062 × 3.0 mm³) was first obtained before injecting the bone cement, and both pre- and post-injection C-arm CBCT scans were performed before and after the cement injection. In the post-injection scan, we obtained 200 projections over a 178° arc, from which we extracted 10, 20, 40, 66 projections to create sparse projections. Since the shape of the cadaver torso was maintained, rigid 3D–2D registrations between the prior CT and the pre- and post-injection x-ray projections were first computed. Difference images were computed by subtracting the DRRs of the registered CT from the post-injection x-ray images as described in section 2.2.

VOCs of 250 × 250 × 400 voxels were reconstructed from different image sets by the PL method, and merged with the prior volume to show a full surgical field. A ground truth VOC was computed by subtracting pre-injection CBCT volume from the post-injection CBCT volume (both have 512 × 512 × 512 voxels with a voxel size of 0.2930 × 0.2930 × 0.2930 mm³) followed by the median filtering, and the VOC reconstructions were compared to the ground truth.

Figure 7 shows example images of the pre- and post-injection x-ray projections and the computed difference images. Since there was only negligible motion of the cadaver and the C-arm between these two scans, the difference image clearly shows the bone cement while other structures are eliminated. Note that this subtraction process is commonly used in digital subtraction angiography (Heautot et al 1998, Bidaut et al 1998, Anxionnat et al 2001). Similarity measures between the reconstructed VOCs and the ground truth for different β values are plotted in figure 8. Although there was noticeable reconstruction quality degradation when only 10 images were used, it was minor when 20, 40 and 66 images were used. Figure 9 shows example slice images of the VOC reconstruction computed from 20 projections and merged with the prior volume for β = 10⁴ that was chosen based on the similarity measures shown in figure 8. Unlike the FBP reconstruction computed from the same 20 projections, both the surgical change including extravasation and the unchanged patient anatomy are clearly visible in the VOC reconstruction.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the lumbar vertebroplasty case, a prior CT volume was used as a prior volume and registered to the post-injection CBCT x-ray projections using the method in section 2.2. Figure 10 shows example images of the post-injection CBCT projections, DRR computed from the registered prior CT and the resulting difference image. Since the prior CT was rigidly registered to the CBCT projections, there are intensity mismatches especially in the soft tissue and air regions. Also DRR looks blurred compared to the post-injection x-ray image due to the low spatial resolution of the CT data. However, the bone cement is clearly highlighted in the difference image. Similarity measures between the reconstructed VOCs and the ground truth as a function of different β values are computed and plotted in figure 11. Similar to the thoracic vertebroplasty case, there were only minor reconstruction quality variations between 20 and 66 projections although image quality degradation was relatively large when only 10 projections were used. Figure 12 shows example slice images of the VOC reconstruction computed from 20 projections with β = 10⁷ (chosen from the similarity measures shown in figure 11) and merged with the pre-injection CBCT prior volume. Both the surgical change and the background patient anatomy are clearly visible and the extravasation can be detected in the final VOC reconstruction although it is very hard to see the injected bone cement in relation to the patient anatomy in the FBP reconstruction from the same 20 post-injection x-ray images.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The key advantage of the proposed VOC reconstruction approach is that one can achieve low-dose, high-quality imaging as well as compensate for the truncation that is common in CBCT due to the small lateral FOV. Therefore, we assume that the prior volume has all the information of the patient anatomy except for the surgical change-–i.e. the prior volume is not truncated. This assumption is important because we subtract all redundant information in the prior volume from the intraoperative x-ray scan data to highlight the surgical change only. This difference image computation can be considered as a form of truncation compensation because the ray sum from all the non-changed patient anatomical structures will be subtracted in the projection domain and the resulting difference images therefore show only the changes without truncation. This is a reasonable assumption, since a non-truncated, patient-specific prior volume is often available in image-guided surgical procedures from the preoperative, diagnostic and/or planning images. CT is a very common preoperative imaging modality for image-guided orthopaedic/spine surgery and radiation therapy. For anatomical sites for which CBCT can be used to provide a non-truncated prior volume (e.g. head and neck imaging and/or CBCT data acquired from an offset-detector geometry), that image may also serve as a prior volume. As one can easily imagine, there are intensity inconsistencies between the DRRs created from the prior volume and the intraoperative x-ray projections at views for which the truncated anatomy is involved in DRR generation. Theoretically, the VOC reconstruction approach will reconstruct not only the changes made by the surgical procedure but also any difference between the current target volume and the prior volume. However, accurate reconstruction can be computed only when we consider both positive and negative intensity differences during the reconstruction, as discussed below in section 9.3. In this case, the full CBCT FOV must be reconstructed because the change will no longer be local but spread over the entire FOV.

In this paper, we demonstrated the sparse-scan scenario to achieve low-dose imaging. This is a reasonable choice because we can simply reduce the radiation dose by reducing the number of projections and a sparse scan is known to be effectively handled by an iterative reconstruction technique with CS regularization. Another possible way of dose reduction is a reduced-dose scan where dose per projection is reduced. The reduced-dose scan can be better handled by the statistical reconstruction approach than conventional FBP or other iterative reconstruction techniques such as POCS or TV-based CS approach because the statistical approach takes into account the noise property in the observed image while others do not.

One may also achieve dose reduction by limiting the exposure to a small FOV (using, for example, dynamic collimation) when the exact target location is known (Heuscher and Noo 2011). This approach may require hardware changes such as collimators, beam filters and controller to dynamically change the FOV during the scan. Additionally, a smaller FOV may cause significant errors in 3D–2D registration of the prior volume to intraoperative x-ray projections due to the lack of common information between them. (Note that the prior volume does not reflect the surgical change.) Therefore, we focused the current work on the sparse-scan scenario in which the number of projections collected over a semicircular arc is significantly reduced, thus reducing the overall dose. The other dose reduction scenarios are areas of on-going research.

We assumed that the noise in the difference image obeys a Poisson distribution, and enforced positivity constraints because the surgical change we considered here, i.e. bone cement, has larger linear attenuation coefficient than the patient anatomy. Although this is partly true in an ideal case where the 3D–2D registration between the prior volume and the observed x-rays is perfect, there exist both positive and negative changes in the difference images in reality due to registration errors. Even when the 3D–2D registration is perfect, the noise in the difference images does not obey a Poisson distribution although Poisson approximation shows reasonably good reconstruction. Therefore, further modification of the noise model and the constraints may be needed to improve reconstruction and allow not only positive changes but also negative changes, which is our future research. For example, one may use the Gaussian noise model without positivity constraints or experimentally estimate the noise distribution rather than assuming Poisson to handle this issue.

In our experiments, we only performed 3D–2D rigid registration because the target was not moved or significantly deformed between the prior scan and the intraoperative scans. However, in reality, there could be a considerable amount of soft-tissue deformation between the prior scans (usually preoperative CT scan or initial intraoperative CBCT scan) and the intraoperative scans. Therefore, an accurate 3D–2D deformable registration might be necessary to compute accurate difference images showing surgical change only. However, the 3D–2D deformable registration problem is very challenging and is an active area of research that deserves separate papers. One viable approach we can consider for image-guided surgery applications is to perform 3D–3D deformable registration as soon as we get an initial intraoperative CBCT volume for visualizing the surgical field and registering surgical tools within a surgical navigation system at the beginning of the surgery (Nithiananthan et al 2011, Uneri et al 2011). Since the patient does not move significantly during the surgical procedure, the successive CBCT scans will already be roughly registered to the initial CBCT volume as well as the prior images. Therefore, the successive 3D–2D registration problem will be relatively easier than direct registration of the 3D prior volume to 2D intraoperative scans.

Our registration and reconstruction algorithms were implemented on a combination of Matlab and C++ with the forward and backprojection implemented on a GPU. Since they were not optimized for the best performance, the overall process still requires significant amount of computation time (∼10 min for 10–66 projections and volume sizes shown in our experiments). We believe that this can be significantly improved by setting the VOI more tighter using the surgical plan (note that the VOIs in our experiments are large) and designing the algorithm for GPU and optimizing the code using C++, which is under development. For example, Jia et al (2010) reported that a CBCT volume can be reconstructed using a TV-based compressed sensing reconstruction from 20–40 projections in a few minutes using their GPU-friendly version of forward–backward splitting algorithm.