Note
Go to the end to download the full example code.
EasyCT device creation#
This is an example how to create a new device. In this case a new system for EasyCT The device should be run only one time to create a new device. A folder with a unique identifier will be created Afterwars the device can be read from the folder and added to the new TOR files created
# # sphinx_gallery_thumbnail_path = ‘examples/EasyPETCT/easyPETCT.png’
Imports to create the device
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
from toor.Geometry.easyPETBased import EasyCTGeometry
from toor.DetectionLayout.RadiationProducer import GenericRadiativeSource
from toor.DetectionLayout.Modules import easyPETModule
from toor.Designer import DeviceDesignerStandalone
from toor.Device import StoreDeviceInFo, EnergyResolutionFunction
# from toor.TORFilesReader import ToRFile
SYSTEM ENERGY RESPONSE FUNCTION (Not mandatory)
def systemEnergyResponseFunction(E, Er, p1, p2):
"""
Energy response function of the system
:param energy: energy of the photon
:param Er: energy resolution
"""
fwhm = np.sqrt((p1 / E) ** 2 + (p2) ** 2)
return fwhm / E
energies = np.array([30, 59.6, 511])
energy_resolution = np.array([0.63, 0.33, 0.14])
fit = curve_fit(systemEnergyResponseFunction, energies, energy_resolution)
plt.figure()
plt.plot(energies, energy_resolution * 100, 'ro', label='Data')
plt.plot(np.arange(25, 600, 10), systemEnergyResponseFunction(np.arange(25, 600, 10), *fit[0]) * 100, 'b-', label='Fit')
plt.xlabel('Energy (keV)')
plt.ylabel('Energy Resolution (%)')
plt.legend()
plt.savefig("../../images/system_energy_response_function.png")
# .. image:: ../../images/system_energy_response_function.png
# :alt: EasyCT Diagram
# :width: 400px
# :align: center
energy_window = [energies - energies * systemEnergyResponseFunction(energies, *fit[0]),
energies + energies * systemEnergyResponseFunction(energies, *fit[0])]
systemEnergyResolution = EnergyResolutionFunction(p1=fit[0][1], p2=fit[0][2])
C:\Users\pedro\OneDrive\Documentos\GitHub\Infinity-Tomographic-Reconstruction\docs\source\examples\EasyPETCT\1_createNewDevice.py:51: OptimizeWarning: Covariance of the parameters could not be estimated
fit = curve_fit(systemEnergyResponseFunction, energies, energy_resolution)
Setup the type of the detector module. You should not call the PETModule class directly. This object should entry as argument in the geometry class type for proper setting. This allows to set multiple cells. Number of modules, rotations and translations are set after the geometry class is created.
_module = easyPETModule
# _module = PETModule
Setup the x-ray source
Now we define the characteristics of the x-ray source using the GenericRadiativeSource class. The source is set to be an Am-241 source with a focal spot diameter of 1 mm, and the shielding is set to be a cylinder made of lead with a density of 11.34 g/cm³ and a thickness of 0.5 mm. Set x-ray producer object
xrayproducer = GenericRadiativeSource()
xrayproducer.setSourceName("Am-241")
xrayproducer.setSourceActivity(1.0 * 37000)
xrayproducer.setFocalSpotDiameter(1)
xrayproducer.setShieldingShape("Cylinder")
xrayproducer.setShieldingMaterial("Lead")
xrayproducer.setShieldingDensity(11.34)
xrayproducer.setShieldingThickness(0.5)
xrayproducer.setShieldingHeight(4)
xrayproducer.setShieldingRadius(12.5)
xrayproducer.setMainEmissions({1: {"energy": 59.54, "intensity": 0.36},
2: {"energy": 26.34, "intensity": 0.024},
})
The next step is to choose the geometry type, which is DualRotationCTGeometry in this case. This function is inherited from the DualRotationGeometry class which is an Device Object. Here we set the distance between the two points of rotation, the distance between the fan motor and the detector modules (closest side) and the distance between the fan motor and the detector modules (far side). as well as the initial position of the x-ray source.
newDevice = EasyCTGeometry(detector_moduleA=_module, detector_moduleB=_module, x_ray_producer=xrayproducer)
newDevice.setDeviceName("EasyCT_simulation_16_2_special")
newDevice.setDeviceType("CT")
newDevice.setEnergyResolutionFunction(systemEnergyResolution) # use to apply energy cuts
newDevice.setDistanceBetweenMotors(30) # Distance between the two points of rotation
newDevice.setDistanceFanMotorToDetectorModulesOnSideA(
0) # Distance between the fan motor and the detector modules (closest side)
newDevice.setDistanceFanMotorToDetectorModulesOnSideB(
60) # Distance between the fan motor and the detector modules (far side)
newDevice.xRayProducer.setFocalSpotInitialPositionWKSystem([12.55, 3, 0]) # simulation 31 _1
newDevice.xRayProducer.setFocalSpotInitialPositionWKSystem([-2, 0, 0]) # simulation 16_2 e real
# newDevice.xRayProducer.setFocalSpotInitialPositionWKSystem([2.45, 7.7, 0]) # simulation 16_2 special
newDevice.evaluateInitialSourcePosition() # evaluate the initial position of the source
Calculating source position for all events detected...
Focal spot initial position set to: [[28. 0. 0.]]
Set modules Side A. For each module, should be in the list the equivalent rotation and translation variables. If for example two modules are set, the variables should be in the list as follows:
moduleSideA_X_translation = np.array([15, 20], dtype=np.float32) moduleSideA_Y_translation = np.array([0, 0], dtype=np.float32)
…
Very important. The translations are regarding the fan motor center. The rotations are regarding the center of the module.
newDevice.setNumberOfDetectorModulesSideA(1)
moduleSideA_X_translation = np.array([15], dtype=np.float32)
moduleSideA_Y_translation = np.array([0], dtype=np.float32)
moduleSideA_Z_translation = np.array([0], dtype=np.float32)
moduleSideA_alpha_rotation = np.array([0], dtype=np.float32)
moduleSideA_beta_rotation = np.array([0], dtype=np.float32)
moduleSideA_sigma_rotation = np.array([0], dtype=np.float32)
for i in range(newDevice.numberOfDetectorModulesSideA):
# newDevice.detectorModulesSideA[i].model32()
newDevice.detectorModulesSideA[i].model16_2()
newDevice.detectorModulesSideA[i].setXTranslation(moduleSideA_X_translation[i])
newDevice.detectorModulesSideA[i].setYTranslation(moduleSideA_Y_translation[i])
newDevice.detectorModulesSideA[i].setZTranslation(moduleSideA_Z_translation[i])
newDevice.detectorModulesSideA[i].setAlphaRotation(moduleSideA_alpha_rotation[i])
newDevice.detectorModulesSideA[i].setBetaRotation(moduleSideA_beta_rotation[i])
newDevice.detectorModulesSideA[i].setSigmaRotation(moduleSideA_sigma_rotation[i])
Set modules Side B.
newDevice.setNumberOfDetectorModulesSideB(1)
moduleSideB_X_translation = np.array([-75], dtype=np.float32)
moduleSideB_Y_translation = np.array([0], dtype=np.float32)
moduleSideB_Z_translation = np.array([0], dtype=np.float32)
moduleSideB_alpha_rotation = np.array([0], dtype=np.float32)
moduleSideB_beta_rotation = np.array([0], dtype=np.float32)
moduleSideB_sigma_rotation = np.array([180], dtype=np.float32)
for i in range(newDevice.numberOfDetectorModulesSideB):
# newDevice.detectorModulesSideB[i].model32()
newDevice.detectorModulesSideB[i].model16_2()
newDevice.detectorModulesSideB[i].setXTranslation(moduleSideB_X_translation[i])
newDevice.detectorModulesSideB[i].setYTranslation(moduleSideB_Y_translation[i])
newDevice.detectorModulesSideB[i].setZTranslation(moduleSideB_Z_translation[i])
newDevice.detectorModulesSideB[i].setAlphaRotation(moduleSideB_alpha_rotation[i])
newDevice.detectorModulesSideB[i].setBetaRotation(moduleSideB_beta_rotation[i])
newDevice.detectorModulesSideB[i].setSigmaRotation(moduleSideB_sigma_rotation[i])
Set the inital coordinates of the system. In both coordinate
newDevice.generateInitialCoordinatesWKSystem()
newDevice.generateInitialCoordinatesXYSystem()
Calculating source position for all events detected...
Save the device in a folder with a unique identifier. The folder will be created in the current directory.
modifyDevice = True
if not modifyDevice:
print("Creating new device")
newDevice.generateDeviceUUID() # one time only
newDevice.createDirectory() # one time only
print("Device created in: ", newDevice.deviceDirectory)
storeDevice = StoreDeviceInFo(device_directory=newDevice.deviceDirectory) # one time only
device_path = newDevice.deviceDirectory
else:
device_path = "C:\\Users\\pedro\\OneDrive\\Documentos\\GitHub\\Infinity-Tomographic-Reconstruction\\configurations\\08d98d7f-a3c1-4cdf-a037-54655c7bdbb7_EasyCT"
device_path = "C:\\Users\\pedro\\OneDrive\\Documentos\\GitHub\\Infinity-Tomographic-Reconstruction\\configurations\\5433bf80-06b5-468f-9692-674f4b007605_EasyCT_simulation_16_2_special"
storeDevice = StoreDeviceInFo(device_directory=device_path) # one time only
storeDevice.createDeviceInDirectory(object=newDevice)
Device created successfully
readDevice = StoreDeviceInFo(device_directory=device_path) # read the device saved
newDevice_Read = readDevice.readDeviceFromDirectory()
designer = DeviceDesignerStandalone(device=newDevice)
designer.addDevice()
designer.addxRayProducerSource()
designer.startRender()
Test some initial positions of the source and the detectors
unique_header = np.repeat(np.arange(0, 32), 13)
axial_motor_angles = (np.zeros(32 * 13))
fan_motor_angles = np.tile(np.arange(-90, 105, 15), 32)
Calculate the coordinates for the previous angles
newDevice.detectorSideBCoordinatesAfterMovement(axial_motor_angles, fan_motor_angles, unique_header)
axial_motor_angles = np.array([0, 0], dtype=np.float32)
fan_motor_angles = np.array([0, 0], dtype=np.float32)
newDevice.sourcePositionAfterMovement(axial_motor_angles, fan_motor_angles)
plt.figure(figsize=(10, 10))
plt.plot(newDevice.originSystemWZ[0], newDevice.originSystemWZ[1], 'ro', label='Origin Fan Motor')
# plot source center
plt.plot(newDevice.sourceCenter[:, 0], newDevice.sourceCenter[:, 1], 'bo', label='Source Center')
plt.plot(newDevice.originSystemXY[0], newDevice.originSystemXY[1], 'ko', label='Origin FOV')
plt.plot(newDevice.centerFace[:, 0], newDevice.centerFace[:, 1], 'go', label='Center Face Detector Module B')
plt.plot(newDevice._verticesB[:, :, 0], newDevice._verticesB[:, :, 1], 'mo', label='Vertices Base Detector Module B')
plt.plot([np.ones(newDevice.centerFace.shape[0]) * newDevice.originSystemWZ[0, 0],
newDevice.centerFace[:, 0]], [np.ones(newDevice.centerFace.shape[0]) * newDevice.originSystemWZ[1, 0],
newDevice.centerFace[:, 1]], '-')
plt.xlabel('X (mm)')
plt.ylabel('Z (mm)')
plt.legend()
plt.figure(figsize=(10, 10))
# x an Z direction
plt.plot(newDevice.originSystemWZ[0], newDevice.originSystemWZ[2], 'ro', label='Origin Fan Motor')
# plot source center
plt.plot(newDevice.sourceCenter[:, 0], newDevice.sourceCenter[:, 2], 'bo', label='Source Center')
plt.plot(newDevice.originSystemXY[0], newDevice.originSystemXY[2], 'ko', label='Origin FOV')
plt.plot(newDevice.centerFace[:, 0], newDevice.centerFace[:, 2], 'go', label='Center Face Detector Module B')
plt.plot(newDevice._verticesB[:, :, 0], newDevice._verticesB[:, :, 2], 'mo', label='Vertices Base Detector Module B')
plt.xlabel('X (mm)')
plt.ylabel('Y (mm)')
plt.legend()
plt.show()
Calculating parametric positions of the center and vertices of the detector for all events...
Centroid calculated for all events...
Vertice 0 calculated for all events...
Vertice 1 calculated for all events...
Vertice 2 calculated for all events...
Vertice 3 calculated for all events...
Vertice 4 calculated for all events...
Vertice 5 calculated for all events...
Vertice 6 calculated for all events...
Vertice 7 calculated for all events...
[[[ 29.825 60. -18.24 ]
[ 29.824999 90. -18.24 ]
[ 27.824152 89.99998 -18.24 ]
...
[ 29.824999 90. -15.96 ]
[ 27.824152 89.99998 -15.96 ]
[ 27.823095 59.99993 -15.96 ]]
[[ 14.301821 57.910255 -18.24 ]
[ 6.5372486 86.88803 -18.24 ]
[ 4.6045837 86.370155 -18.24 ]
...
[ 6.5372486 86.88803 -15.96 ]
[ 4.6045837 86.370155 -15.96 ]
[ 12.368147 57.39206 -15.96 ]]
[[ -0.15155411 51.874023 -18.24 ]
[-15.151554 77.85478 -18.24 ]
[-16.884327 76.85434 -18.24 ]
...
[-15.151554 77.85478 -15.96 ]
[-16.884327 76.85434 -15.96 ]
[ -1.8852196 50.873013 -15.96 ]]
...
[[ -1.8852215 -50.873013 15.96 ]
[-16.884327 -76.85435 15.96 ]
[-15.151554 -77.85478 15.96 ]
...
[-16.884327 -76.85435 18.24 ]
[-15.151554 -77.85478 18.24 ]
[ -0.15155602 -51.874023 18.24 ]]
[[ 12.368145 -57.39206 15.96 ]
[ 4.6045856 -86.370155 15.96 ]
[ 6.5372505 -86.88803 15.96 ]
...
[ 4.6045856 -86.370155 18.24 ]
[ 6.5372505 -86.88803 18.24 ]
[ 14.301818 -57.910255 18.24 ]]
[[ 27.823093 -59.99993 15.96 ]
[ 27.824154 -89.99998 15.96 ]
[ 29.825 -90. 15.96 ]
...
[ 27.824154 -89.99998 18.24 ]
[ 29.825 -90. 18.24 ]
[ 29.824997 -60. 18.24 ]]]
Calculating source position for all events detected...
Total running time of the script: (0 minutes 4.361 seconds)