WaterPreCorrection

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 processing of the simulation is written using SimpleRTK.

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 and the output projections are available in the archive beam_hardening.tgz.

#!/usr/bin/env python
from __future__ import print_function
import SimpleRTK as srtk
import os
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np
import scipy.optimize
import glob

# Script parameters
dir = "output/"
spacingVol = [ 0.5, 0.5, 0.5 ]
sizeOutputVol = [ 512, 1, 512 ]
order = 6;
refMu = 0.02;
originVol = [ -(sizeOutputVol[0]/2+0.5)*spacingVol[0], 0., -(sizeOutputVol[2]/2+0.5)*spacingVol[2]]

# List of filenames
fileNames = srtk.VectorString(0)
for file in os.listdir(dir):
    if file.startswith("attenuation") and file.endswith(".mha"):
        fileNames.push_back(dir + file)

# Read in full geometry
geometryReader = srtk.ThreeDCircularProjectionGeometryXMLFileReader()
geometryReader.SetFileName ( "output/geometry.xml" )
geometry = geometryReader.Execute ( );

projReader = srtk.ProjectionsReader()
projReader.SetFileNames(fileNames)

constantImageSource = srtk.ConstantImageSource()
constantImageSource.SetOrigin( originVol )
constantImageSource.SetSpacing( spacingVol )
constantImageSource.SetSize( sizeOutputVol )
constantImageSource.SetConstant(0.)

writer = srtk.ImageFileWriter()

# Create template
template = constantImageSource.Execute()
drawCyl = srtk.DrawCylinderImageFilter()
drawCyl.SetDensity(refMu)
drawCyl.SetAngle(0)
drawCyl.SetAxis([100, 0, 100])
drawCyl.SetCenter([0, 0, 0])
template = drawCyl.Execute(template)
writer.SetFileName ( 'template.mha')
writer.Execute(template)
template = srtk.GetArrayFromImage(template[:,:,:]).flatten(1)

# Create weights (where the template should match)
weights = constantImageSource.Execute()

drawCyl.SetDensity(1)
drawCyl.SetAngle(0)
drawCyl.SetAxis([125, 0, 125])
weights = drawCyl.Execute(weights)

drawCyl.SetDensity(-1)
drawCyl.SetAngle(0)
drawCyl.SetAxis([102, 0, 102])
weights = drawCyl.Execute(weights)

drawCyl.SetDensity(1)
drawCyl.SetAngle(0)
drawCyl.SetAxis([98, 0, 98])
weights = drawCyl.Execute(weights)

writer.SetFileName ( 'weights.mha')
writer.Execute(weights)
weights = srtk.GetArrayFromImage(weights[:,:,:]).flatten(1)

# Create reconstructed images
waterPre = srtk.WaterPrecorrectionImageFilter()
wpcoeffs = srtk.VectorDouble(order+1, 0)

fdkFilter = srtk.FDKConeBeamReconstructionFilter()
fdkFilter.SetGeometry( geometry )

fdks = [None] * (order+1)
for o in range(0,order+1):
    proj = projReader.Execute();
    wpcoeffs[o-1] = 0.
    wpcoeffs[o] = 1.
    waterPre.SetCoefficients(wpcoeffs)
    proj = waterPre.Execute(proj)  
    sourceVol = constantImageSource.Execute()
    fdk = fdkFilter.Execute( sourceVol, proj )
    writer.SetFileName ( 'fdk%d.mha' % o )
    writer.Execute ( fdk )
    fdks[o] = srtk.GetArrayFromImage(fdk[:,:,:]).flatten(1)

# 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)
proj = projReader.Execute();
waterPre.SetCoefficients(c)
proj = waterPre.Execute(proj)  
sourceVol = constantImageSource.Execute()
fdk = fdkFilter.Execute( sourceVol, proj )
writer.SetFileName ( 'fdk.mha' )
writer.Execute ( fdk )

reader = srtk.ImageFileReader()
reader.SetFileName("fdk.mha")
fdk = reader.Execute()
reader.SetFileName("fdk1.mha")
fdk1 = reader.Execute()

plt.plot(srtk.GetArrayFromImage(fdk)[:,0,256], label='Corrected')
plt.plot(srtk.GetArrayFromImage(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