demo.ct_mace

This demo illustrates a 2D tomography example using a high-dynamic range phantom. First the radon transform is applied using relatively sparse views. Then proportional noise is added and filtered backprojection applied to produce an initial reconstruction.

The forward model is the forward radon transform, the backprojection can be chosen as either filtered or unfiltered backprojection, and the prior agent can be bm3d, a version of TV, wavelet-based denoising, or bilateral.

Functions

add_noise(sinogram)	Add realistic noise to the sinogram.
baseline(sinogram, noisy_sinogram, ...)	Use filtered backprojection for baseline recon.
ct_demo()	Illustrate a MACE reconstruction on a CT example with high dynamic range.
display_images(image_list, image_titles, ...)	Display images to capture high dynamic range.
get_image_and_mask(image_path[, imsize])	Load high dynamic range image and specify recon region.
get_scaled_sinogram(ground_truth, theta[, ...])	Apply physically realistic scaling to sinogram.

demo.ct_mace.ct_demo()

Illustrate a MACE reconstruction on a CT example with high dynamic range.

demo.ct_mace.get_image_and_mask(image_path, imsize=None)

Load high dynamic range image and specify recon region.

The images in this demo are designed to mimic single-energy (~100 KeV) CT images with high dynamic range. The pixel values are in Hounsfield units, with air as 0 and water as 1000. Hounsfield units are closely related to the atomic weight of the associated material. We then scale so that water is 1 and air is 0. In these scaled units, steel and other dense metals are about 12 to 15, but may go up to 20.

The demo image has high dynamic range with some characteristics like those seen in CT scans of baggage.

With parallel beam CT, only the center circular region can be properly reconstructed. Here we mask to that region.

Note that if the image is resized here, then the pixel pitch (physical distance between pixels as measured on the object being imaged) must be adjusted later as part of scaling.

Parameters

• image_path – path to image file

• imsize – Leave imsize = [] to use the natural image size or set to an ordered pair to resize to that shape

Returns

Ground truth image in scaled Hounsfield units mask to circular region of interest image_scale (1 if the image size is unchanged)

demo.ct_mace.get_scaled_sinogram(ground_truth, theta, image_scale=1)

Apply physically realistic scaling to sinogram.

Scaling: All CT scans are relative to a baseline scan with no objects - i.e., a scan of air, which makes air 0.

The raw projection operator (radon function in python or matlab) sums along pixels, with assumed distance 1 between pixels. To get the correct units, we need to scale. The first step is raw projection:

\(\text{Raw projection} = Ax\) (\(x\) is the image as scaled above, \(Ax\) is the output of the radon method)

This result needs to be scaled according to physical units.

We scale by the pixel pitch (distance between pixels) times the x-ray density of water.

We assume pixel pitch of 0.93 mm = 0.093 cm and water density at 100KeV of 0.17 \(\text{cm}^{-1}\). With an image size of 512x512 pixels, this corresponds to a width of 512*0.93/10 = 47.16 cm or 18.75 in. Scaling of the projection is then

\[y \text{ (or scaled projection) } = Ax * \text{pixel_pitch} * \text{water_xray_density}\]

Note that if the original image was resized, then we have to adjust the pixel pitch accordingly so that the image has the same physical dimensions. For example, if the image was resized from 512x512 to 256x256, then the pixel pitch is increased by a factor of two.

The sinogram values should be on the order of 0-5 for most images. Here are approximate X-ray densities for common materials you might have in an image

water ~0.17 \(\text{cm}^{-1}\) steel ~3.0 \(\text{cm}^{-1}\) aluminum ~1.0 \(\text{cm}^{-1}\) plastic ~0.1 \(\text{cm}^{-1}\)

Parameters

• ground_truth – Ground truth image from get_image_and_mask

• theta – Vector of angles of views to use

• image_scale – image scale from get_image_and_mask

Returns

sinogram of the ground truth, scaled as described above. sino_scale so that the sinogram equals radon transform * sino_scale

demo.ct_mace.add_noise(sinogram)

Add realistic noise to the sinogram.

Noise modeling: The noise is modeled as additive noise but with variance that depends on the signal.

\[y_{noise} = y + w\]

For an entry \(y(i)\), the noise \(w(i)\) is \(N(0, \sigma^2(i))\) with

\[\sigma^2(i) = \frac{1}{\lambda_0} \exp(y(i)),\]

where \(\lambda_0\) is the photon count at the detector for an empty scan. For convenience, we define \(\alpha = 1/\lambda_0\).

Parameters

sinogram – sinogram from get_scaled_sinogram

Returns

noisy sinogram

demo.ct_mace.baseline(sinogram, noisy_sinogram, sino_scale, theta)

Use filtered backprojection for baseline recon.

Invert the radon transform using filtered backprojection for both the noise-free and noisy sinograms.

Since the sinogram was scaled using sino_scale as described above, we need to apply the inverse scaling to the reconstructions.

Parameters

• sinogram – sinogram from get_scaled_sinogram

• noisy_sinogram – noisy sinogram from add_noise

• sino_scale – scale from get_scaled_sinogram

• theta – angles of the views in the sinogram

Returns

Filtered backprojection for sinogram and noisy_sinogram

demo.ct_mace.display_images(image_list, image_titles, ground_truth)

Display images to capture high dynamic range.

Note that these images have a scale bar to indicate the intensity and that all the images are scaled to have the same intensity range. We show several intensity bands to highlight the high dynamic range. The reconstructions have values outside the given range, so the intensities are clipped, but the full range is shown in the titles.

Parameters

• image_list – list of images to display

• image_titles – title for the images

• ground_truth – ground truth image for calculating NRMSE

Returns

None