experiment_04_analyticArtifactPhantom_SimAndRecon

[xie-pixel]

Hi everyone, I'm trying to reproduce the sources of metal artifacts using XCIST. I ran the experiment_04_analyticArtifactPhantom_SimAndRecon.py script in the evaluation directory without any modifications, but the results I got differ from those in the paper. Here's what I obtained:

[MingyeWu]

Happy new year!

The Root Cause:
We identified that the new streak artifact was caused by the low-signal correction, which was introduced after the evaluation study and is enabled by default.
In this simulation, there are metal rods that cause extremely high attenuation. However, since this is an ideal case with no noise and sufficient precision, the low-signal correction is unnecessary. In fact, because the correction replaces low-signal pixels with the average of neighboring pixels, it introduces error and causes streak artifacts.

The Solution:
We added a cfg flag (cfg.physics.disable_low_signal_correction) to disable the low-signal correction when needed. We updated the following code:
pyfiles/PrepView.py
The evaluation sample script (experiment_analyticArtifactPhantom.py)

[xie-pixel]

Happy New Year to you too! 🙏
Thank you so much for the quick fix! The low-signal correction was exactly what was causing the issue, and your solution resolved it immediately. I really appreciate your help.

[xie-pixel]

Dear MingyeWu,
First of all, thank you very much for your detailed reply to my previous question about metal artifact simulation. Your guidance was extremely helpful for my research.
Now I encounter a new issue that I cannot resolve on my own, and I hope to get your insights:In your published work (the ideal artifact-free image in Figure 1 of your paper), I noticed faint dark stripe artifacts that should not exist under ideal conditions. To verify this, I constructed my own Fe+water phantom and ran simulations under strictly ideal, artifact-free conditions (no noise, no beam hardening, no scatter, no View aliasing, etc.), but the resulting image still shows the same type of dark stripe artifacts.

My core goal is to generate a completely artifact-free ideal image of the Fe+water phantom under simulation, which will serve as the gold standard for subsequent metal artifact correction research. Below is the key code snippet I used for the ideal simulation. I suspect the issue may lie in the phantom modeling logic, numerical calculation precision, or unconsidered physical assumptions in the simulation pipeline, but I cannot pinpoint the exact cause.
Any suggestions or guidance on this issue would be greatly appreciated. Thank you very much for your time and help!

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)

modify Phantom

#re-scale the size of phantom

ct.phantom.callback = "Phantom_Analytic" # name of function that reads and models phantom

ct.phantom.projectorCallback = "C_Projector_Analytic" # name of function that performs projection through phantom

ct.phantom.filename = "analyticArtifactPhantom.ppm" # phantom filename

ct.phantom.scale = 1 # re-scale the size of phantom

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

ct.phantom.filename = 'C:/Users\Administrator\Desktop\XCIST\CatSim\gecatsim\phantom/xjx_phantom/phantom1_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.5 # re-scale the size of phantom

modify Scanner

ct.scanner.detectorRowsPerMod = 1 # Will be revised for certain experiments.
ct.scanner.detectorRowCount = ct.scanner.detectorRowsPerMod # Will be revised for certain experiments.
ct.scanner.detectorColOffset = 0.25 # Will be revised for certain experiments.

ct.protocol.viewsPerRotation = 2000 # Will be revised for certain experiments.
ct.protocol.viewCount = ct.protocol.viewsPerRotation # Will be revised for certain experiments.
ct.protocol.stopViewId = ct.protocol.viewCount - 1 # Will be revised for certain experiments.
ct.protocol.spectrumScaling = 0.931
ct.protocol.mA = 500
ct.protocol.maxPrep = -1

ct.physics.energyCount = 12 # Using 120 kVp, so 10 kV/bin.
ct.physics.monochromatic = 88 # Will be revised for certain experiments.
ct.physics.colSampleCount = 4
ct.physics.rowSampleCount = 1
ct.physics.srcXSampleCount = 4
ct.physics.srcYSampleCount = 1
ct.physics.viewSampleCount = 1
ct.physics.enableQuantumNoise = 0 # Will be revised for certain experiments.
ct.physics.enableElectronicNoise = 0 # Will be revised for certain experiments.
ct.physics.disable_low_signal_correction = 1 # Disable the low signal correction

ct.recon.fov = 300.0
ct.recon.sliceThickness = ct.scanner.detectorRowSize*ct.scanner.sid/ct.scanner.sdd
ct.recon.sliceCount = 1 # Will be revised for certain experiments.
ct.recon.kernelType = 'standard'
ct.recon.mu = 0.0178

change

ct.protocol.viewsPerRotation = 360

ct.protocol.viewCount = ct.protocol.viewsPerRotation

ct.protocol.stopViewId = ct.protocol.viewCount - 1

ct.physics.monochromatic = -1

ct.recon.mu = 0.02

ct.scanner.detectorRowsPerMod = 64

ct.scanner.detectorRowCount = ct.scanner.detectorRowsPerMod

ct.physics.scatterCallback = "Scatter_ConvolutionModel"

ct.physics.scatterKernelCallback = ""

ct.physics.scatterScaleFactor = 1

ct.physics.enableQuantumNoise = 1

ct.physics.enableElectronicNoise = 1

ct.physics.disable_low_signal_correction = 0

name of result

ct.resultsName = "test_ideal"

##--------- 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()