WaterPreCorrection: Difference between revisions

From Openrtk
Jump to navigation Jump to search
(Created page with "This example illustrates how to apply empirical cupping correction using the [http://onlinelibrary.wiley.com/doi/10.1118/1.2188076/abstract algorithm of Kachelriess et al.] na...")
 
No edit summary
 
(5 intermediate revisions by the same user not shown)
Line 1: Line 1:
This example illustrates how to apply empirical cupping correction using the [http://onlinelibrary.wiley.com/doi/10.1118/1.2188076/abstract algorithm of Kachelriess et al.] named [http://www.openrtk.org/Doxygen/classrtk_1_1WaterPrecorrectionImageFilter.html WaterPrecorrection] in RTK. The example uses a Gate simulation using the [http://wiki.opengatecollaboration.org/index.php/Users_Guide_V7.2:Readout_parameters_for_Radiotherapy_applications:_Actors#Fixed_Forced_Detection_CT Fixed Forced Detection Actor]. It is written using [http://wiki.openrtk.org/index.php?title=SimpleRTK SimpleRTK].
This example illustrates how to apply empirical cupping correction using the [http://onlinelibrary.wiley.com/doi/10.1118/1.2188076/abstract algorithm of Kachelriess et al.] named [http://www.openrtk.org/Doxygen/classrtk_1_1WaterPrecorrectionImageFilter.html WaterPrecorrection] in RTK. The example uses a Gate simulation using the [http://wiki.opengatecollaboration.org/index.php/Users_Guide:Tools_to_Interact_with_the_Simulation_:_Actors#Fixed_Forced_Detection_CT fixed forced detection actor].


The simulation implements a 120 kV beam, a detector with 512x3 with an response curve. Only the primary beam is simulated. The simulation files and the output projections are available here:
The simulation implements a 120 kV beam, a detector with 512x3 pixels and an energy response curve. Only the primary beam is simulated.
[https://midas3.kitware.com/midas/download/item/321827/beam_hardening.tgz]
 
= Python wrapping =
 
This version uses the [https://pypi.org/project/itk-rtk/ Python wrapping]. The simulation files, the output projections and the processing script are available in the archive [https://data.kitware.com/api/v1/file/5d394cea877dfcc9022c922b/download].
 
<source lang="python">
#!/usr/bin/env python
 
import itk
from itk import RTK as rtk
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 = list()
for file in os.listdir(dir):
    if file.startswith("attenuation") and file.endswith(".mha"):
        fileNames.append(dir + file)
 
# Read in full geometry
geometryReader = rtk.ThreeDCircularProjectionGeometryXMLFileReader.New()
geometryReader.SetFilename ( "output/geometry.xml" )
geometryReader.GenerateOutputInformation()
geometry = geometryReader.GetGeometry();
 
ImageType = itk.Image[itk.F,3]
projReader = rtk.ProjectionsReader[ImageType].New()
projReader.SetFileNames(fileNames)
 
constantImageSource = rtk.ConstantImageSource[ImageType].New()
constantImageSource.SetOrigin( originVol )
constantImageSource.SetSpacing( spacingVol )
constantImageSource.SetSize( sizeOutputVol )
constantImageSource.SetConstant(0.)
 
# Create template
drawCyl = rtk.DrawEllipsoidImageFilter[ImageType, ImageType].New()
drawCyl.SetInput(constantImageSource.GetOutput())
drawCyl.SetDensity(refMu)
drawCyl.SetAngle(0)
drawCyl.SetAxis([100, 0, 100])
drawCyl.SetCenter([0, 0, 0])
drawCyl.Update()
template = drawCyl.GetOutput()
itk.imwrite(template, 'template.mha')
template = itk.GetArrayFromImage(template).flatten()
 
# Create weights (where the template should match)
drawCyl.SetDensity(1)
drawCyl.SetAngle(0)
drawCyl.SetAxis([125, 0, 125])
drawCyl.Update()
weights = drawCyl.GetOutput()
weights.DisconnectPipeline()
 
drawCyl.SetInput(weights)
drawCyl.SetDensity(-1)
drawCyl.SetAngle(0)
drawCyl.SetAxis([102, 0, 102])
drawCyl.Update()
weights = drawCyl.GetOutput()
weights.DisconnectPipeline()
 
drawCyl.SetInput(weights)
drawCyl.SetDensity(1)
drawCyl.SetAngle(0)
drawCyl.SetAxis([98, 0, 98])
weights = drawCyl.GetOutput()
 
itk.imwrite(weights, 'weights.mha')
weights = itk.GetArrayFromImage(weights).flatten()
 
# Create reconstructed images
waterPre = rtk.WaterPrecorrectionImageFilter[ImageType,ImageType].New()
wpcoeffs = np.zeros((order+1))
 
fdkFilter = rtk.FDKConeBeamReconstructionFilter[ImageType].New()
fdkFilter.SetGeometry( geometry )
fdkFilter.SetInput(0,constantImageSource.GetOutput())
fdkFilter.SetInput(1,waterPre.GetOutput())
 
fdks = [None] * (order+1)
for o in range(0,order+1):
    wpcoeffs[o-1] = 0.
    wpcoeffs[o] = 1.
    waterPre.SetCoefficients(wpcoeffs)
    waterPre.SetInput(projReader.GetOutput())
    fdkFilter.Update()
    fdk = fdkFilter.GetOutput()
    itk.imwrite(fdk, 'fdk%d.mha' % o )
    fdks[o] = itk.GetArrayFromImage(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)
 
waterPre.SetCoefficients(c)
fdkFilter.Update()
itk.imwrite(fdkFilter.GetOutput(), 'fdk.mha' )
 
fdk = itk.imread('fdk.mha')
fdk1 = itk.imread('fdk1.mha')
 
plt.plot(itk.GetArrayFromImage(fdk)[:,0,256], label='Corrected')
plt.plot(itk.GetArrayFromImage(fdk1)[:,0,256], label='Uncorrected')
plt.legend()
plt.xlabel('Pixel number')
plt.ylabel('Attenuation')
plt.xlim(0,512)
plt.savefig('profile.png')
</source>
 
The resulting central profiles are
 
[[File:profile.png]]
 
= SimpleRTK wrapping =
 
This version is outdated and uses [http://wiki.openrtk.org/index.php?title=SimpleRTK SimpleRTK]. The simulation files, the output projections and the processing script are available in the archive [https://data.kitware.com/api/v1/item/5be99c008d777f2179a2e4ae/download beam_hardening.tgz].


<source lang="python">
<source lang="python">
Line 130: Line 266:
plt.savefig('profile.png')  
plt.savefig('profile.png')  
</source>
</source>
The result is
[[File:profile.png]]

Latest revision as of 02:40, 25 July 2019

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.

Python wrapping

This version uses the Python wrapping. The simulation files, the output projections and the processing script are available in the archive [1].

#!/usr/bin/env python

import itk
from itk import RTK as rtk
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 = list()
for file in os.listdir(dir):
    if file.startswith("attenuation") and file.endswith(".mha"):
        fileNames.append(dir + file)

# Read in full geometry
geometryReader = rtk.ThreeDCircularProjectionGeometryXMLFileReader.New()
geometryReader.SetFilename ( "output/geometry.xml" )
geometryReader.GenerateOutputInformation()
geometry = geometryReader.GetGeometry();

ImageType = itk.Image[itk.F,3]
projReader = rtk.ProjectionsReader[ImageType].New()
projReader.SetFileNames(fileNames)

constantImageSource = rtk.ConstantImageSource[ImageType].New()
constantImageSource.SetOrigin( originVol )
constantImageSource.SetSpacing( spacingVol )
constantImageSource.SetSize( sizeOutputVol )
constantImageSource.SetConstant(0.)

# Create template
drawCyl = rtk.DrawEllipsoidImageFilter[ImageType, ImageType].New()
drawCyl.SetInput(constantImageSource.GetOutput())
drawCyl.SetDensity(refMu)
drawCyl.SetAngle(0)
drawCyl.SetAxis([100, 0, 100])
drawCyl.SetCenter([0, 0, 0])
drawCyl.Update()
template = drawCyl.GetOutput()
itk.imwrite(template, 'template.mha')
template = itk.GetArrayFromImage(template).flatten()

# Create weights (where the template should match)
drawCyl.SetDensity(1)
drawCyl.SetAngle(0)
drawCyl.SetAxis([125, 0, 125])
drawCyl.Update()
weights = drawCyl.GetOutput()
weights.DisconnectPipeline()

drawCyl.SetInput(weights)
drawCyl.SetDensity(-1)
drawCyl.SetAngle(0)
drawCyl.SetAxis([102, 0, 102])
drawCyl.Update()
weights = drawCyl.GetOutput()
weights.DisconnectPipeline()

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

itk.imwrite(weights, 'weights.mha')
weights = itk.GetArrayFromImage(weights).flatten()

# Create reconstructed images
waterPre = rtk.WaterPrecorrectionImageFilter[ImageType,ImageType].New()
wpcoeffs = np.zeros((order+1))

fdkFilter = rtk.FDKConeBeamReconstructionFilter[ImageType].New()
fdkFilter.SetGeometry( geometry )
fdkFilter.SetInput(0,constantImageSource.GetOutput())
fdkFilter.SetInput(1,waterPre.GetOutput())

fdks = [None] * (order+1)
for o in range(0,order+1):
    wpcoeffs[o-1] = 0.
    wpcoeffs[o] = 1.
    waterPre.SetCoefficients(wpcoeffs)
    waterPre.SetInput(projReader.GetOutput())
    fdkFilter.Update()
    fdk = fdkFilter.GetOutput()
    itk.imwrite(fdk, 'fdk%d.mha' % o )
    fdks[o] = itk.GetArrayFromImage(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)

waterPre.SetCoefficients(c)
fdkFilter.Update()
itk.imwrite(fdkFilter.GetOutput(), 'fdk.mha' )

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

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

Profile.png

SimpleRTK wrapping

This version is outdated and uses SimpleRTK. The simulation files, the output projections and the processing script 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')