from cil.io import TIFFStackReader
from cil.utilities.jupyter import islicer
from cil.processors import TransmissionAbsorptionConverter
from cil.framework import AcquisitionGeometry, AcquisitionData
from cil.utilities.display import show_geometry, show2D
from cil.recon import FDK
import numpy as np
"""
The gvxr json file
{
"WindowSize": [800, 450],
"Source": {
"Position": [0, 200, 0, "mm"],
"Shape": "Point",
"Beam": [
{
"Energy": 500,
"Unit": "keV",
"PhotonCount": 1000
}
]
},
"Detector": {
"Position": [0, -150, 0, "mm"],
"UpVector": [0, 0, 1],
"NumberOfPixels": [900, 900],
"Spacing": [0.2, 0.2, "mm"]
},
"Samples": [
{
"Label": "internals",
"Path": "input_data/TurboPump/internals.stl",
"Unit": "mm",
"Material": ["Element", "Ti"],
"Density": 4.506,
"Transform": [
["Rotation", 90, 1, 0, 0],
["Scaling", 0.2, 0.2, 0.2]
]
},
{
"Label": "front_flange",
"Path": "input_data/TurboPump/front_flange.stl",
"Unit": "mm",
"Material": ["Element", "Al"],
"Density": 2.7,
"Transform": [
["Rotation", 90, 1, 0, 0],
["Scaling", 0.2, 0.2, 0.2]
]
},
{
"Label": "rear_flage",
"Path": "input_data/TurboPump/rear_flange.stl",
"Unit": "mm",
"Material": ["Element", "Al"],
"Density": 2.7,
"Transform": [
["Rotation", 90, 1, 0, 0],
["Scaling", 0.2, 0.2, 0.2]
]
},
{
"Label": "housing",
"Path": "input_data/TurboPump/housing.stl",
"Unit": "mm",
"Material": ["Element", "Fe"],
"Density": 7.874,
"Transform": [
["Rotation", 90, 1, 0, 0],
["Scaling", 0.2, 0.2, 0.2]
]
},
{
"Label": "roller_bearing",
"Path": "input_data/TurboPump/ThrustRollerBearing.stl",
"Unit": "mm",
"Material": ["Element", "Ti"],
"Density": 4.506,
"Transform": [
["Rotation", 90, 1, 0, 0],
["Scaling", 0.2, 0.2, 0.2]
]
}
],
"Scan": {
"NumberOfProjections":721,
"FinalAngle": 360,
"IncludeFinalAngle": false,
"CenterOfRotation": [0,0,0],
"OutFolder": "./input_data/TurboPump/scan/scan"
}
}
"""
reader = TIFFStackReader(file_name='/mnt/materials/IBSim/TurboPump-scan/projections/')
data_in = reader.read()
islicer(data_in)
It can be hard to create the right geometry. Each system will use it's own standard definitions so you can't just plug the same numbers in. CIL defines a right hand cordinate system in which you can place your objects in space. The good news is CIL will convert the geometry to ASTRA and TIGRE for you, so you can compare reconstruction backends without worrying about redefining the geometry.
show_geometry will help you visualise the geometry you've created with CIL's definitions. Is this what you expect?
ag_3d = AcquisitionGeometry.create_Cone3D([0, -200, 0],[0, 150, 0] )
ag_3d.set_angles(np.linspace(0, 360,721, False))
ag_3d.set_panel([900, 900],[0.2,0.2])
ag_3d.set_labels(['angle','vertical','horizontal'])
print(ag_3d)
show_geometry(ag_3d)
AcquisitionData from the geometry and the raw data array¶data = AcquisitionData(data_in, deep_copy=False, geometry=ag_3d)
islicer(data)
data_corr = TransmissionAbsorptionConverter(white_level=data_in.max())(data)
islicer(data_corr)
data_ss = data_corr.get_slice(vertical='centre')
show2D(data_ss)
show_geometry(data_ss.geometry)
Tune the reconstruction roi by hand. Set the number of voxels to reconstruct and the window position
ig_2D = data_ss.geometry.get_ImageGeometry()
print(ig_2D)
ig_2D.voxel_num_x = 600
ig_2D.voxel_num_y = 600
ig_2D.center_x = -40 * ig_2D.voxel_size_x
ig_2D.center_y = 0 * ig_2D.voxel_size_y
print(ig_2D)
reco1 = FDK(data_ss, ig_2D).run()
show2D(reco1)
Start with the default image geometry. This has the voxel sizes set to pixel_sixe/magnification. From the 2D dataset we can use the same roi for x and y, but will have to tune the z direction manually. Let's reconstruct it at a quatar resolution first. We need to scale everything by our reconstruction volume resolution.
#resolution = 1
resolution = 1/4
ig = ag_3d.get_ImageGeometry(resolution)
print(ig)
ig.voxel_num_x = 600*resolution
ig.voxel_num_y = 600*resolution
ig.voxel_num_z = 400*resolution
ig.center_x = -40*resolution * ig.voxel_size_x
ig.center_y = 0*resolution * ig.voxel_size_y
ig.center_z = -100*resolution * ig.voxel_size_z
print(ig)
recon_full = FDK(data_corr,ig).run()
# look at the maximum value in each direction to help us constrain the reconstruction window
show2D([recon_full.array.max(axis=0),recon_full.array.max(axis=1)])
Once we are happy with the speicied roi let's reconstruct at full resolution, go back and set the resolution to 1 (the default value).
islicer(recon_full, size=25)