Back to Research highlights

# Three-dimensional medical X-ray imaging using sparse data

## Research highlight by Finnish Centre of Excellence in Inverse Problems Research

## Introduction

This is a review of a large body of research work on medical X-ray imaging done in the period 2001-2010, resulting in more than a dozen scientific articles, five patents, and one commercial product (www.vt-cube.com). A large number of people have been involved in the research effort; their names appear below in the references to publications.

The research started around the year 2000 from the initiative of a company called Instrumentarium Imaging, specializing in designing and manufacturing medical X-ray imaging devices. Here is a picture of some of the relevant products:

From left to right: dental panoramic imaging device, mammography device, and C-arm surgical imaging device.

The starting point of the research was the desire of using the above X-ray machines to collect projection data from several directions and to compute a three-dimensional (3D) reconstruction of tissue. This is different from traditional computerized tomography (CT), where one collects data from more than a hundred directions from full angle of view, enabling very accurate 3D reconstructions of tissue. The data collected with the above devices has typically a limited angle of view and consists only of at most a dozen images.

The philosophy of the proposed 3D imaging is the following: *Let the doctor collect whatever data is practical to measure with the equipment at hand. Our task is to create mathematical methods that compute 3D reconstructions with high enough quality for the clinical task in question.*

A research project was started in 2001 with funding from Finnish Technology Agency (TEKES), involving researchers from both industry and the CoE in Inverse Problems, with the goal of developing new low-dose, sparse-data 3D reconstruction algorithms.

Traditional filtered back-projection (FBP) methods are not optimal for this kind of imaging, since the incomplete data violates the assumptions of FBP. A better choice is **Bayesian inversion**, a flexible and noise-robust family of algorithms. It turned out that Bayesian inversion is very well-suited for this application. See the following article as well as the examples below.

**Siltanen S, Kolehmainen V, Jarvenpaa S, Kaipio J P, Koistinen P, Lassas M, Pirttila J and Somersalo E**, *Statistical inversion for medical X-ray tomo- graphy with few radiographs I: General theory*, Physics in Medicine and Biology **48** (2003), pp. 1437-1463.

## Single tooth and total variation prior

The study of Bayesian 3D X-ray imaging started with a very simple situation. A single wisdom tooth, was attached to a rotating platform and imaged from full angle of view. Here is a picture of the laboratory setup:

Here is a closer look of the tooth, the rotating platform, and the intraoral digital X-ray detector attached to the wall:

A total of 23 projection images was recorded as shown in the drawing above. Each image has 664 x 872 pixels; below are some examples of how the projection images look like. Note that the data set covers full angle of view, allowing quantitative comparison between the limited angle reconstruction (done using only the angles marked by white discs above) and full-angle reconstruction.

Here are two reconstructed horizontal slices through the tooth. Left slice is computed using full-angle data, and the right slice is computed using the 9 projections from roughly 65 degree angle of view. Note the typical elongation of the limited angle reconstruction towards the direction of the measured rays. A stack of 400 such two-dimensional reconstructions was computed, forming a three-dimensional reconstruction of the whole tooth.

According to Todd Quinto's seminal 1993 article on stable singularities in tomography, certain crisp boundaries can be seen in the horizontal slices. Consequently, stacking all the 400 horizontal reconstructed slices on top of each other and slicing the stack *vertically* leads to the following results:

The limited-angle result is not so different from the full-angle result when sliced this way. This is an example of a Bayesian limited angle reconstruction not being perfect but still showing some aspects of tissue reliably.

The above computations are based on the use of (approximate) total variation prior and the computation of the Maximum a Posteriori (MAP) estimate. For details of computation and more images of reconstructions, see the following article:

**Kolehmainen V, Siltanen S, Jarvenpaa S, Kaipio J P, Koistinen P, Lassas M, Pirttila J and Somersalo E**, *Statistical inversion for medical X-ray tomography with few radiographs II: Application to dental radiology*, Physics in Medicine and Biology **48** (2003), pp. 1465-1490.

## Intraoral imaging, head phantom and total variation prior

This measurement models chair-side intraoral X-ray imaging modality. The idea is to insert a small digital X-ray detector inside the mouth of a patient, let the patient hold the detector against her teeth, and take for example five images from different directions. Using a small steel ball attached to the crown of one of the teeth enables the calibration of the projection directions.

The yellow device on the left bottom corner in the image below is the X-ray source, which was moved into different positions for recording projection data.

Seven radiographs were recorded from a 60 degree angle of view according to the drawing below on the left. On the right is one of the 664x872 pixel images. The data consists of roughly 4 million pixel values. The region of interest was divided into some 42 million voxels, so the data matrix has size 4000000x42000000. It is not advisable to construct such a large matrix explicitly, so matrix-free iterative reconstruction methods were designed.

Here are three vertical slices through the teeth showing that the method can separate features at different depths inside the tissue. This is crucial for example in the detection of alveolar disease, or bone loss between teeth. The loss may not be visible in a single X-ray image, but a slice at the correct depth will reveal the condition.

The above computations are based on the use of (approximate) total variation prior and the computation of the Maximum a Posteriori (MAP) estimate. For details of computation and more images of reconstructions, see the following articles:

**Kolehmainen V, Siltanen S, Jarvenpaa S, Kaipio J P, Koistinen P, Lassas M, Pirttila J and Somersalo E**, *Statistical inversion for medical X-ray tomography with few radiographs II: Application to dental radiology*, Physics in Medicine and Biology **48** (2003), pp. 1465-1490.

**Kolehmainen V, Vanne A, Siltanen S, Jarvenpaa S, Kaipio J P, Lassas M and Kalke M**, Bayesian Inversion Method for 3-D Dental X-ray Imaging. Elektrotechnik & Informationstechnik **124** (2007), pp. 248-253.

## Highlight: 3D imaging using a panoramic device

The 3D X-ray imaging project has so far resulted in one commercial product, the VT tomograph (www.vt-cube.com ). The idea is to use a panoramic dental imaging device in a novel way to produce 3D information helpful for dental implant planning.

Often the information contained in a 2D X-ray image, panoramic or otherwise, is not enough for dental implant planning. Namely, the aim is to replace a missing tooth by an artificial one which is screwed into bone.

The important thing is to measure the distance to the nerve located inside the mandible. The hole should be drilled deep enough for the implant to hold, but not so deep that the nerve is damaged as the result can be facial paralysis. This measurement cannot be reliably made from a projection image or a panoramic image because of geometric distortions.

The panoramic device was reprogrammed so that it collects limited-angle data as shown in the drawing below. Using a Bayesian 3D reconstruction method, the three-dimensional structure of the region of interest can be estimated. Then, a two dimensional slice can be extracted, enabling a reliable measurement of the desired distance.

The crucial thing about the above 3D imaging method is that a dental clinic can use the existing equipment: a software update turns a 2D panoramic device into a 3D tomographic instrument. See www.vt-cube.com for more information.

Furthermore, the projection data can be complemented by an additional panoramic image and the reconstruction quality improved that way:

See **Hyvonen N, Kalke M, Lassas M, Setala H and Siltanen S**, *Three-dimensional X-ray imaging using hybrid data collected with a digital panoramic device.* To appear in Inverse Problems and Imaging.

## Modified level set method for X-ray tomography

Reconstructions from limited angle data are often stretched along the directions of the X-rays, almost regardless of the inversion algorithm. Removing those stretching artifacts has been one of the central goals in the research project. The most successful approach so far is the modified level set method.

The core object in every level set method is the so-called level set function, whose zero level set defines the boundary between two different components. In inverse problems, the unknown object to be reconstructed is represented differently (often as constants) inside and outside the domain defined by the level set. The level set function deforms according to a certain time evolution scheme, resulting in flexible and topologically variable domains determined by the zero level set.

While implementing the classical piecewise constant level set method for our tomographic application we noticed that the level set function itself looks like a reconstruction of the X-ray attenuation coefficient. We decided to represent the attenuation coefficient by zero outside the level set (providing natural enforcement of non-negativity) and by *the level set function itself* inside the level set. The result is an inversion method with very efficient truncation of the reconstruction to zero.

Here is an example using the synthetic Shepp-Logan phantom. We collected 10 projection images with full angle of view and reconstructed the phatntom (left) using filtered back-projection (middle) and our level set method (right).

We repeated the above experiment with real data measured using a surgical C-arm device from a knee specimen. Left: full-angle reconstruction from 60 projections using the level set method, middle: filtered back-projection recon- struction from 10 projections, right: level set reconstruction from 10 projections.

In the case of limited angle tomography, our novel level set method turns out to be quite effective in avoiding stretching artifacts. Left: level set reconstruction from full-angle data, middle: back-projection (without filtering) from limited- angle data, right: level set reconstruction from limited-angle data.

Using the level set function itself to represent the attenuation coefficient inside the level set has an additional advantage. Namely, we can give a convergence proof for our method. See the article **Kolehmainen V, Lassas M and Siltanen S**, Limited data X-ray tomography using nonlinear evolution equations, SIAM Journal of Scientific Computation **30** (2008), pp. 1413-1429.

## Wavelet-based reconstructions using discretization- invariant Besov space priors

As demonstrated above, using the total variation prior and MAP estimation leads to useful 3D reconstructions from sparse projection data. However, the conditional mean (CM) estimates corresponding to the total variation prior do not behave well, as shown in the article

**Lassas M and Siltanen S**, Can one use total variation prior for edge-preserving Bayesian inversion? Inverse Problems **20** (2004), pp. 1537-1563.

The problem with Bayesian interpretation of the total variation prior is two-fold: when discretization is refined, the CM estimates either diverge or approach a Gaussian smoothness prior, depending on scaling. Also, there is no scaling that would ensure the convergence of both MAP and CM estimates to meaningful limits as the discretization is refined. Our conclusion is that

*Total variation prior should not be used in Bayesian inversion*.

The above negative result was our motivation for designing a prior that would (a) have the desirable edge-preserving properties of total variation prior, and (b) would be discretization-invariant so that both MAP and CM estimates make sense in the fine discretization limit. It turns out that wavelet-based Besov space priors fulfill both (a) and (b).

Concerning property (b), see the article **Lassas M, Saksman E and Siltanen S**, Discretization invariant Bayesian inversion and Besov space priors. Inverse Problems and Imaging **3** (2009), pp. 87-122.

Concerning property (a), let us look at some reconstructions. Here the target was a mammographic specimen. Remarkably, 90% of the wavelet coefficienst were left out in the optimization problem of computing the MAP estimate, redusing the degrees of freedom significantly and providing automatic compression of the result.

**Rantala M, Vanska S, Jarvenpaa S, Kalke M, Lassas M, Moberg J and Siltanen S**, Wavelet-based reconstruction for limited angle X-ray tomography, IEEE Transactions on Medical Imaging **25** (2006), pp. 210-217.

Encouraged by the promising results reported in the above paper, we proceeded to test the wavelet-based reconstruction approach with local tomography. Here the target is a dry mandible specimen, and the whole objet was not visible in any of the projection images. We tried using all wavelet scales in all of the reconstruction area (below left) and only the coarse wavelets outside the region of interest (below right).

We compared the region of interest to Lambda-tomographic reconstruction from the same data set:

These results were published in **Niinimaki K, Siltanen S and Kolehmainen V**, Bayesian multiresolution method for local tomography in dental X-ray imaging. Physics in Medicine and Biology **52** (2007), pp. 6663-6678.

The Besov space prior has three parameters (p,q,s), whose optimal choice is not well-understood. One possible solution, namely empirical determination of those parameters from measure data, is described in the article **Vanska S, Lassas M and Siltanen S**, Statistical X-ray tomography using empirical Besov priors, International Journal of Tomography & Statistics **11** (2009), pp. 3-32. Here is an example on choosing parameter s empirically when p=1.5=q: