Cupping correction

This example illustrates how to apply empirical cupping correction using the algorithm of Kachelriess et al. named WaterPrecorrection in RTK. The example uses a Gate simulation using the fixed forced detection actor.

The simulation implements a 120 kV beam, a detector with 512x3 pixels and an energy response curve. Only the primary beam is simulated.

The simulation files, the output projections and the processing script are available here.

#!/usr/bin/env python

import itk
from itk import RTK as rtk
import os
import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize
import glob

# Script parameters
directory = "output/"
spacing = [ 0.5 ] * 3
size = [ 512, 1, 512 ]
order = 6
reference_attenuation = 0.02
origin = [ -(size[0]/2+0.5)*spacing[0], 0., -(size[2]/2+0.5)*spacing[2]]

# List of filenames
file_names = list()
for file in os.listdir(directory):
    if file.startswith("attenuation") and file.endswith(".mha"):
        file_names.append(directory + file)

# Read in full geometry
geometry = rtk.read_geometry("output/geometry.xml")

# Crate template image
ImageType = itk.Image[itk.F,3]
projections_reader = rtk.ProjectionsReader[ImageType].New(file_names=file_names)
projections = projections_reader.GetOutput()
constant_image_filter = rtk.ConstantImageSource[ImageType].New(origin=origin, spacing=spacing, size=size)
constant_image = constant_image_filter.GetOutput()
template = rtk.draw_ellipsoid_image_filter(constant_image, density=reference_attenuation, axis=[100, 0, 100])
itk.imwrite(template, 'template.mha')
template = itk.array_from_image(template).flatten()

# Create weights (where the template should match)
weights = rtk.draw_ellipsoid_image_filter(constant_image, density=1., axis=[125, 0, 125])
weights = rtk.draw_ellipsoid_image_filter(weights, density=-1., axis=[102, 0, 102])
weights = rtk.draw_ellipsoid_image_filter(weights, density=1., axis=[98, 0, 98])
itk.imwrite(weights, 'weights.mha')
weights = itk.array_from_image(weights).flatten()

# Create reconstructed images
wpcoeffs = np.zeros((order+1))
fdks = [None] * (order+1)
for o in range(0,order+1):
    wpcoeffs[o-1] = 0.
    wpcoeffs[o] = 1.
    water_precorrection = rtk.water_precorrection_image_filter(projections, coefficients=wpcoeffs, in_place=False)
    fdk = rtk.fdk_cone_beam_reconstruction_filter(constant_image, water_precorrection, geometry=geometry)
    itk.imwrite(fdk, f'fdk{o}.mha')
    fdks[o] = itk.array_from_image(fdk).flatten()

# Create and solve the linear system of equation B.c= a to find the coeffs c
a = np.zeros((order+1))
B = np.zeros((order+1,order+1))
for i in range(0,order+1):
    a[i] = np.sum(weights * fdks[i] * template)
    for j in np.arange(i,order+1):
        B[i,j] = np.sum(weights * fdks[i] * fdks[j])
        B[j,i] = B[i,j]
c = np.linalg.solve(B,a)

water_precorrection = rtk.water_precorrection_image_filter(projections, coefficients=c)
fdk = rtk.fdk_cone_beam_reconstruction_filter(constant_image, water_precorrection, geometry=geometry)
itk.imwrite(fdk, 'fdk.mha' )

fdk = itk.imread('fdk.mha')
fdk1 = itk.imread('fdk1.mha')

plt.plot(itk.array_from_image(fdk)[:,0,256], label='Corrected')
plt.plot(itk.array_from_image(fdk1)[:,0,256], label='Uncorrected')
plt.legend()
plt.xlabel('Pixel number')
plt.ylabel('Attenuation')
plt.xlim(0,512)
plt.savefig('profile.png')

The resulting central profiles are

Profile