Hello，everone！I built a simple XCIST phantom made of water and titanium cylindrical rods. During simulation, discontinuous dark edges occur along the metal boundaries on the output sinogram. This dark fringe artifact is extremely obvious after computing the ratio of low-energy sinogram to high-energy sinogram.
I provide five figures:
（1）Reconstructed CT slice
（2）Full original sinogram
（3）Windowed local zoom of raw sinogram showing metal dark edges
（4）Noise-free dual-energy division sinogram
（5）Dual-energy division sinogram with simulated noise
My simulation code is attached. May I ask what leads to such discontinuous dark shadows at metal edges in XCIST simulation, and how to mitigate this artifact?
# Copyright 2024, GE Precision HealthCare. All rights reserved. See https://github.com/xcist/main/tree/master/license
###------------ import XCIST-CatSim
import gecatsim as xc
import gecatsim.reconstruction.pyfiles.recon as recon
##--------- Initialize
ct = xc.CatSim() # initialization
##--------- Make changes to parameters (optional)
Phantom
ct.phantom.callback = "Phantom_Voxelized" # name of function that reads and models phantom
ct.phantom.projectorCallback = "C_Projector_Voxelized" # name of function that performs projection through phantom
#ct.phantom.filename = 'CatSim_logo_1024.json' # phantom filename, not actually used in AAPM Recon
ct.phantom.filename = 'C:/Users\Administrator\Desktop\XCIST\CatSim\gecatsim\phantom/xjx_phantom/phantom7_single_1024.json' # phantom filename
ct.phantom.centerOffset = [0.0, 0.0, 0.0] # offset of phantom center relative to origin (in mm)
ct.phantom.scale = 0.8 # re-scale the size of phantom
ct.phantom.projectorNumThreads = 1
physics
ct.physics.energyCount = 12 # number of energy bins
ct.physics.monochromatic = -1 # -1 for polychromatic (see protocol.ct);
ct.physics.colSampleCount = 1 # number of samples of detector cells in lateral direction
ct.physics.rowSampleCount = 1 # number of samples of detector cells in longitudinal direction
ct.physics.srcXSampleCount = 2 # number of samples of focal spot in lateral direction
ct.physics.srcYSampleCount = 2 # number of samples of focal spot cells in longitudinal direction
ct.physics.viewSampleCount = 2 # number of samples of each view angle range in rotational direction
ct.physics.recalcDet = 0 # recalculate detector geometry
ct.physics.recalcSrc = 0 # recalculate source geometry and relative intensity
ct.physics.recalcRayAngle = 0 # recalculate source-to-detector-cell ray angles
ct.physics.recalcSpec = 0 # recalculate spectrum
ct.physics.recalcFilt = 0 # recalculate filters
ct.physics.recalcFlux = 0 # recalculate flux
ct.physics.recalcPht = 0 # recalculate phantom
ct.physics.enableQuantumNoise = 0 # enable quantum noise
ct.physics.enableElectronicNoise = 0 # enable electronic noise
ct.physics.disable_low_signal_correction = 0 # Disable the low signal correction
ct.physics.rayAngleCallback = "Detector_RayAngles_2D" # name of function to calculate source-to-detector-cell ray angles
ct.physics.fluxCallback = "Detection_Flux" # name of function to calculate flux
ct.physics.scatterCallback = "Scatter_ConvolutionModel" # name of function to calculate scatter
ct.physics.scatterKernelCallback = "" # name of function to calculate scatter kernel ("" for default kernel)
ct.physics.scatterScaleFactor = 1 # scale factor, 1 appropriate for 64-mm detector and 20-cm water
ct.physics.callback_pre_log = "Scatter_Correction"
ct.physics.prefilterCallback = "Detection_prefilter" # name of function to calculate detection pre-filter
ct.physics.crosstalkCallback = "CalcCrossTalk" # name of function to calculate X-ray crosstalk in the detector
ct.physics.col_crosstalk = 0.025
ct.physics.row_crosstalk = 0.02
ct.physics.opticalCrosstalkCallback = "CalcOptCrossTalk" # name of function to calculate X-ray crosstalk in the detector
ct.physics.col_crosstalk_opt = 0.04
ct.physics.row_crosstalk_opt = 0.045
ct.physics.lagCallback = "" # name of function to calculate detector lag
ct.physics.opticalCrosstalkCallback = "" # name of function to calculate optical crosstalk in the detector
ct.physics.DASCallback = "Detection_DAS" # name of function to calculate the detection process
ct.physics.outputCallback = "WriteRawView" # name of function to produce the simulation output
ct.physics.callback_post_log = ''
ct.physics.EffectiveMu = 0.2
ct.physics.BHC_poly_order = 5
ct.physics.BHC_max_length_mm = 300
ct.physics.BHC_length_step_mm = 10
protocol
ct.protocol.scanTypes = [1, 1, 1] # flags for airscan, offset scan, phantom scan
ct.protocol.scanTrajectory = "Gantry_Helical" # name of the function that defines the scanning trajectory and model
ct.protocol.viewsPerRotation = 1000 # total numbers of view per rotation
ct.protocol.viewCount = 1000 # total number of views in scan
ct.protocol.startViewId = 0 # index of the first view in the scan
ct.protocol.stopViewId = ct.protocol.startViewId + ct.protocol.viewCount - 1 # index of the last view in the scan
ct.protocol.airViewCount = 1 # number of views averaged for air scan
ct.protocol.offsetViewCount = 1 # number of views averaged for offset scan
ct.protocol.rotationTime = 1.0 # gantry rotation period (in seconds)
ct.protocol.rotationDirection = 1 # gantry rotation direction (1=CW, -1 CCW, seen from table foot-end)
ct.protocol.startAngle = 0 # relative to vertical y-axis (n degrees)
ct.protocol.tableSpeed = 0 # speed of table translation along positive z-axis (in mm/sec)
ct.protocol.startZ = 0 # start z-position of table
ct.protocol.tiltAngle = 0 # gantry tilt angle towards negative z-axis (in degrees)
ct.protocol.wobbleDistance = 0.0 # focalspot wobble distance
ct.protocol.focalspotOffset = [0, 0, 0] # focalspot position offset
ct.protocol.mA = 500 # tube current (in mA)
ct.protocol.spectrumCallback = "Spectrum" # name of function that reads and models the X-ray spectrum
ct.protocol.spectrumFilename = "xcist_kVp100_tar7_bin1.dat" # name of the spectrum file
ct.protocol.spectrumUnit_mm = 1 # Is the spectrum file in units of photons/sec/mm^2/?
ct.protocol.spectrumUnit_mA = 1 # Is the spectrum file in units of photons/sec//mA?
ct.protocol.spectrumScaling = 1 # scaling factor, works for both mono- and poly-chromatic spectra
ct.protocol.bowtie = "large.txt" # name of the bowtie file (or [] for no bowtie)
ct.protocol.filterCallback = "Xray_Filter" # name of function to compute additional filtration
ct.protocol.flatFilter = ['air', 0.001] # additional filtration - materials and thicknesses (in mm)
ct.protocol.dutyRatio = 1.0 # tube ON time fraction (for pulsed tubes)
ct.protocol.maxPrep = -1 # set the upper limit of prep, non-positive will disable this feature
Scanner
ct.scanner.detectorCallback = "Detector_ThirdgenCurved" # name of function that defines the detector shape and model
ct.scanner.sid = 550.0 # source-to-iso distance (in mm)
ct.scanner.sdd = 950.0 # source-to-detector distance (in mm)
ct.scanner.detectorColsPerMod = 1 # number of detector columns per module
ct.scanner.detectorRowsPerMod = 1 # number of detector rows per module
ct.scanner.detectorColOffset = -1.25 # detector column offset relative to centered position (in detector columns)
ct.scanner.detectorRowOffset = 0.0 # detector row offset relative to centered position (in detector rows)
ct.scanner.detectorColSize = 1.0 # detector column pitch or size (in mm)
ct.scanner.detectorRowSize = 1.0 # detector row pitch or size (in mm)
ct.scanner.detectorColCount = 900 # total number of detector columns
ct.scanner.detectorRowCount = ct.scanner.detectorRowsPerMod # total number of detector rows
ct.scanner.detectorPrefilter = [] # detector filter
ct.scanner.focalspotCallback = "SetFocalspot" # name of function that defines the focal spot shape and model
ct.scanner.focalspotData = "vct_large_fs.npz" # Parameterize the model
ct.scanner.targetAngle = 7.0 # target angle relative to scanner XY-plane (in degrees)
ct.scanner.focalspotWidth = 1.0
ct.scanner.focalspotLength = 1.0
ct.scanner.focalspotWidthThreshold = 0.5
ct.scanner.focalspotLengthThreshold = 0.5
Detector
ct.scanner.detectorMaterial = "GOS" # detector sensor material
ct.scanner.detectorDepth = 3.0 # detector sensor depth (in mm)
ct.scanner.detectionCallback = "Detection_EI" # name of function that defines the detection process (conversion from X-rays to detector signal)
ct.scanner.detectionGain = 0.1 # factor to convert energy to electrons (electrons / keV)
ct.scanner.detectorColFillFraction = 0.9 # active fraction of each detector cell in the column direction
ct.scanner.detectorRowFillFraction = 0.9 # active fraction of each detector cell in the row direction
ct.scanner.eNoise = 25 # standard deviation of Gaussian electronic noise (in electrons)
recon
FOV = 400
ct.recon.fov = FOV # diameter of the reconstruction field-of-view (in mm)
ct.recon.imageSize = 512 # number of columns and rows to be reconstructed (square)
ct.recon.sliceCount = 1 # number of slices to reconstruct
ct.recon.sliceThickness = 0.579 # reconstruction slice thickness AND inter-slice interval (in mm)
ct.recon.centerOffset = [0.0, 0.0, 0.0] # reconstruction offset relative to center of rotation (in mm)
ct.recon.reconType = 'fdk_equiAngle' # Name of the recon function to call
ct.recon.kernelType = 'standard' # 'R-L' for the Ramachandran-Lakshminarayanan (R-L) filter, rectangular window function
ct.recon.startAngle = 0 # in degrees; 0 is with the X-ray source at the top
ct.recon.unit = 'HU' # '/mm', '/cm', or 'HU'
ct.recon.mu = 0.02 # in /mm; typically around 0.02/mm
ct.recon.huOffset = -1000 # unit is HU, -1000 HU by definition but sometimes something else is preferable
ct.recon.printReconParameters = False # Flag to print the recon parameters
ct.recon.saveImageVolume = True # Flag to save recon results as one big file
ct.recon.saveSingleImages = False # Flag to save recon results as individual imagesrecon.printReconParameters = False # Flag to print the recon parameters
ct.recon.displayImagePictures = False # Flag to display the recon results as .png images
ct.recon.saveImagePictureFiles = False # Flag to save the recon results as .png images
ct.recon.displayImagePictureAxes = False # Flag to display the axes on the .png images
ct.recon.displayImagePictureTitles = False # Flag to display the titles on the .png images
name of result
ct.resultsName = "xjx_phantom7_single_noise_500mA_2D_kVp100_tar7_bin1"
##--------- Run simulation
ct.run_all() # run the scans defined by protocol.scanTypes
##--------- Reconstruction
ct.do_Recon = 1
recon.recon(ct)
##--------- Show results
import matplotlib.pyplot as plt
imgFname = "%s_%dx%dx%d.raw" %(ct.resultsName, ct.recon.imageSize, ct.recon.imageSize, ct.recon.sliceCount)
img = xc.rawread(imgFname, [ct.recon.sliceCount, ct.recon.imageSize, ct.recon.imageSize], 'float')
plt.imshow(img[0,:,:], cmap='gray', vmin=-400, vmax=400)
plt.show()
Hello，everone！I built a simple XCIST phantom made of water and titanium cylindrical rods. During simulation, discontinuous dark edges occur along the metal boundaries on the output sinogram. This dark fringe artifact is extremely obvious after computing the ratio of low-energy sinogram to high-energy sinogram.
I provide five figures:
（1）Reconstructed CT slice
（2）Full original sinogram
（3）Windowed local zoom of raw sinogram showing metal dark edges
（4）Noise-free dual-energy division sinogram
（5）Dual-energy division sinogram with simulated noise
My simulation code is attached. May I ask what leads to such discontinuous dark shadows at metal edges in XCIST simulation, and how to mitigate this artifact?
###------------ import XCIST-CatSim
import gecatsim as xc
import gecatsim.reconstruction.pyfiles.recon as recon
##--------- Initialize
ct = xc.CatSim() # initialization
##--------- Make changes to parameters (optional)
Phantom
ct.phantom.callback = "Phantom_Voxelized" # name of function that reads and models phantom
ct.phantom.projectorCallback = "C_Projector_Voxelized" # name of function that performs projection through phantom
#ct.phantom.filename = 'CatSim_logo_1024.json' # phantom filename, not actually used in AAPM Recon
ct.phantom.filename = 'C:/Users\Administrator\Desktop\XCIST\CatSim\gecatsim\phantom/xjx_phantom/phantom7_single_1024.json' # phantom filename
ct.phantom.centerOffset = [0.0, 0.0, 0.0] # offset of phantom center relative to origin (in mm)
ct.phantom.scale = 0.8 # re-scale the size of phantom
ct.phantom.projectorNumThreads = 1
physics
ct.physics.energyCount = 12 # number of energy bins
ct.physics.monochromatic = -1 # -1 for polychromatic (see protocol.ct);
ct.physics.colSampleCount = 1 # number of samples of detector cells in lateral direction
ct.physics.rowSampleCount = 1 # number of samples of detector cells in longitudinal direction
ct.physics.srcXSampleCount = 2 # number of samples of focal spot in lateral direction
ct.physics.srcYSampleCount = 2 # number of samples of focal spot cells in longitudinal direction
ct.physics.viewSampleCount = 2 # number of samples of each view angle range in rotational direction
ct.physics.recalcDet = 0 # recalculate detector geometry
ct.physics.recalcSrc = 0 # recalculate source geometry and relative intensity
ct.physics.recalcRayAngle = 0 # recalculate source-to-detector-cell ray angles
ct.physics.recalcSpec = 0 # recalculate spectrum
ct.physics.recalcFilt = 0 # recalculate filters
ct.physics.recalcFlux = 0 # recalculate flux
ct.physics.recalcPht = 0 # recalculate phantom
ct.physics.enableQuantumNoise = 0 # enable quantum noise
ct.physics.enableElectronicNoise = 0 # enable electronic noise
ct.physics.disable_low_signal_correction = 0 # Disable the low signal correction
ct.physics.rayAngleCallback = "Detector_RayAngles_2D" # name of function to calculate source-to-detector-cell ray angles
ct.physics.fluxCallback = "Detection_Flux" # name of function to calculate flux
ct.physics.scatterCallback = "Scatter_ConvolutionModel" # name of function to calculate scatter
ct.physics.scatterKernelCallback = "" # name of function to calculate scatter kernel ("" for default kernel)
ct.physics.scatterScaleFactor = 1 # scale factor, 1 appropriate for 64-mm detector and 20-cm water
ct.physics.callback_pre_log = "Scatter_Correction"
ct.physics.prefilterCallback = "Detection_prefilter" # name of function to calculate detection pre-filter
ct.physics.crosstalkCallback = "CalcCrossTalk" # name of function to calculate X-ray crosstalk in the detector
ct.physics.col_crosstalk = 0.025
ct.physics.row_crosstalk = 0.02
ct.physics.opticalCrosstalkCallback = "CalcOptCrossTalk" # name of function to calculate X-ray crosstalk in the detector
ct.physics.col_crosstalk_opt = 0.04
ct.physics.row_crosstalk_opt = 0.045
ct.physics.lagCallback = "" # name of function to calculate detector lag
ct.physics.opticalCrosstalkCallback = "" # name of function to calculate optical crosstalk in the detector
ct.physics.DASCallback = "Detection_DAS" # name of function to calculate the detection process
ct.physics.outputCallback = "WriteRawView" # name of function to produce the simulation output
ct.physics.callback_post_log = ''
ct.physics.EffectiveMu = 0.2
ct.physics.BHC_poly_order = 5
ct.physics.BHC_max_length_mm = 300
ct.physics.BHC_length_step_mm = 10
protocol
ct.protocol.scanTypes = [1, 1, 1] # flags for airscan, offset scan, phantom scan
ct.protocol.scanTrajectory = "Gantry_Helical" # name of the function that defines the scanning trajectory and model
ct.protocol.viewsPerRotation = 1000 # total numbers of view per rotation
ct.protocol.viewCount = 1000 # total number of views in scan
ct.protocol.startViewId = 0 # index of the first view in the scan
ct.protocol.stopViewId = ct.protocol.startViewId + ct.protocol.viewCount - 1 # index of the last view in the scan
ct.protocol.airViewCount = 1 # number of views averaged for air scan
ct.protocol.offsetViewCount = 1 # number of views averaged for offset scan
ct.protocol.rotationTime = 1.0 # gantry rotation period (in seconds)
ct.protocol.rotationDirection = 1 # gantry rotation direction (1=CW, -1 CCW, seen from table foot-end)
ct.protocol.startAngle = 0 # relative to vertical y-axis (n degrees)
ct.protocol.tableSpeed = 0 # speed of table translation along positive z-axis (in mm/sec)
ct.protocol.startZ = 0 # start z-position of table
ct.protocol.tiltAngle = 0 # gantry tilt angle towards negative z-axis (in degrees)
ct.protocol.wobbleDistance = 0.0 # focalspot wobble distance
ct.protocol.focalspotOffset = [0, 0, 0] # focalspot position offset
ct.protocol.mA = 500 # tube current (in mA)
ct.protocol.spectrumCallback = "Spectrum" # name of function that reads and models the X-ray spectrum
ct.protocol.spectrumFilename = "xcist_kVp100_tar7_bin1.dat" # name of the spectrum file
ct.protocol.spectrumUnit_mm = 1 # Is the spectrum file in units of photons/sec/mm^2/?
ct.protocol.spectrumUnit_mA = 1 # Is the spectrum file in units of photons/sec//mA?
ct.protocol.spectrumScaling = 1 # scaling factor, works for both mono- and poly-chromatic spectra
ct.protocol.bowtie = "large.txt" # name of the bowtie file (or [] for no bowtie)
ct.protocol.filterCallback = "Xray_Filter" # name of function to compute additional filtration
ct.protocol.flatFilter = ['air', 0.001] # additional filtration - materials and thicknesses (in mm)
ct.protocol.dutyRatio = 1.0 # tube ON time fraction (for pulsed tubes)
ct.protocol.maxPrep = -1 # set the upper limit of prep, non-positive will disable this feature
Scanner
ct.scanner.detectorCallback = "Detector_ThirdgenCurved" # name of function that defines the detector shape and model
ct.scanner.sid = 550.0 # source-to-iso distance (in mm)
ct.scanner.sdd = 950.0 # source-to-detector distance (in mm)
ct.scanner.detectorColsPerMod = 1 # number of detector columns per module
ct.scanner.detectorRowsPerMod = 1 # number of detector rows per module
ct.scanner.detectorColOffset = -1.25 # detector column offset relative to centered position (in detector columns)
ct.scanner.detectorRowOffset = 0.0 # detector row offset relative to centered position (in detector rows)
ct.scanner.detectorColSize = 1.0 # detector column pitch or size (in mm)
ct.scanner.detectorRowSize = 1.0 # detector row pitch or size (in mm)
ct.scanner.detectorColCount = 900 # total number of detector columns
ct.scanner.detectorRowCount = ct.scanner.detectorRowsPerMod # total number of detector rows
ct.scanner.detectorPrefilter = [] # detector filter
ct.scanner.focalspotCallback = "SetFocalspot" # name of function that defines the focal spot shape and model
ct.scanner.focalspotData = "vct_large_fs.npz" # Parameterize the model
ct.scanner.targetAngle = 7.0 # target angle relative to scanner XY-plane (in degrees)
ct.scanner.focalspotWidth = 1.0
ct.scanner.focalspotLength = 1.0
ct.scanner.focalspotWidthThreshold = 0.5
ct.scanner.focalspotLengthThreshold = 0.5
Detector
ct.scanner.detectorMaterial = "GOS" # detector sensor material
ct.scanner.detectorDepth = 3.0 # detector sensor depth (in mm)
ct.scanner.detectionCallback = "Detection_EI" # name of function that defines the detection process (conversion from X-rays to detector signal)
ct.scanner.detectionGain = 0.1 # factor to convert energy to electrons (electrons / keV)
ct.scanner.detectorColFillFraction = 0.9 # active fraction of each detector cell in the column direction
ct.scanner.detectorRowFillFraction = 0.9 # active fraction of each detector cell in the row direction
ct.scanner.eNoise = 25 # standard deviation of Gaussian electronic noise (in electrons)
recon
FOV = 400
ct.recon.fov = FOV # diameter of the reconstruction field-of-view (in mm)
ct.recon.imageSize = 512 # number of columns and rows to be reconstructed (square)
ct.recon.sliceCount = 1 # number of slices to reconstruct
ct.recon.sliceThickness = 0.579 # reconstruction slice thickness AND inter-slice interval (in mm)
ct.recon.centerOffset = [0.0, 0.0, 0.0] # reconstruction offset relative to center of rotation (in mm)
ct.recon.reconType = 'fdk_equiAngle' # Name of the recon function to call
ct.recon.kernelType = 'standard' # 'R-L' for the Ramachandran-Lakshminarayanan (R-L) filter, rectangular window function
ct.recon.startAngle = 0 # in degrees; 0 is with the X-ray source at the top
ct.recon.unit = 'HU' # '/mm', '/cm', or 'HU'
ct.recon.mu = 0.02 # in /mm; typically around 0.02/mm
ct.recon.huOffset = -1000 # unit is HU, -1000 HU by definition but sometimes something else is preferable
ct.recon.printReconParameters = False # Flag to print the recon parameters
ct.recon.saveImageVolume = True # Flag to save recon results as one big file
ct.recon.saveSingleImages = False # Flag to save recon results as individual imagesrecon.printReconParameters = False # Flag to print the recon parameters
ct.recon.displayImagePictures = False # Flag to display the recon results as .png images
ct.recon.saveImagePictureFiles = False # Flag to save the recon results as .png images
ct.recon.displayImagePictureAxes = False # Flag to display the axes on the .png images
ct.recon.displayImagePictureTitles = False # Flag to display the titles on the .png images
name of result
ct.resultsName = "xjx_phantom7_single_noise_500mA_2D_kVp100_tar7_bin1"
##--------- Run simulation
ct.run_all() # run the scans defined by protocol.scanTypes
##--------- Reconstruction
ct.do_Recon = 1
recon.recon(ct)
##--------- Show results
import matplotlib.pyplot as plt
imgFname = "%s_%dx%dx%d.raw" %(ct.resultsName, ct.recon.imageSize, ct.recon.imageSize, ct.recon.sliceCount)
img = xc.rawread(imgFname, [ct.recon.sliceCount, ct.recon.imageSize, ct.recon.imageSize], 'float')
plt.imshow(img[0,:,:], cmap='gray', vmin=-400, vmax=400)
plt.show()