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.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.0, axis=[125, 0, 125]
)
weights = rtk.draw_ellipsoid_image_filter(weights, density=-1.0, axis=[102, 0, 102])
weights = rtk.draw_ellipsoid_image_filter(weights, density=1.0, 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.0
wpcoeffs[o] = 1.0
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