LLNL / LEAP

comprehensive library of 3D transmission Computed Tomography (CT) algorithms with Python and C++ APIs, a PyQt GUI, and fully integrated with PyTorch
https://leapct.readthedocs.io
MIT License
104 stars 10 forks source link

Half-Fan Weighting Artifacts #92

Open lc82111 opened 1 month ago

lc82111 commented 1 month ago

I observed a ring artifact in the overlapping region when using half-fan weighting. Are there any solutions?

image

image

cfg = get_config(sdd=990.805, sod=537.5, nAngles=720,
                 nBins=2, nPad=0, detPhyOffsetX=-138.23, detPhyOffsetY=-96,
                 volZDim=128, volXDim=512, volYDim=512, volXSpacing=1, volYSpacing=1, volZSpacing=1.28)
projs1, angles1 = readDcmProj(cfg, './20240717_504/AVG_504_120_30_air_cw.png', './20240717_504/504_120_30_cw/')
angles = np.linspace(-180, 180, cfg.nAngles).astype(np.float32)

leapct = tomographicModels()
leapct.set_conebeam(cfg.nAngles, cfg.detRows, cfg.detCols, cfg.detSpacing, cfg.detSpacing, cfg.detRowCenter, cfg.detColCenter, angles, cfg.SOD, cfg.SDD)
leapct.set_offsetScan(True)
leapct.set_truncatedScan(True)
# leapct.set_default_volume()
leapct.set_volume(cfg.volXDim, cfg.volYDim, cfg.volZDim, cfg.volXSpacing, cfg.volZSpacing, cfg.volOffsetX, cfg.volOffsetY, cfg.volOffsetZ)
leapct.print_parameters()

vol = leapct.allocateVolume(0)
leapct.FBP(projs1,vol)

print(vol.max(), vol.min(), vol.shape) # 0.030418077 -0.008564696 (128, 512, 512)
print(vol[40, 120:400, 120:400].max(), vol[40, 120:400, 120:400].min())  # 0.019195706 -0.00047323745
print(projs1[4].max(), projs1[4].min(), projs1[4].mean())  # 3.6975892 0.0 1.7670326
vol[vol<0] = 0

fig, axe = plt.subplots(2, 4, figsize=(20, 10))
# vmax = 0.018
vmax = 0.019
vmin = 0.015
_ = axe[0,0].imshow(vol[10, 100:412, 100:412], cmap='gray', vmin=vmin, vmax=vmax) #; fig.colorbar(_, ax=axe[0,0])
_ = axe[0,1].imshow(vol[20, 100:412, 100:412], cmap='gray', vmin=vmin, vmax=vmax)# ; fig.colorbar(_, ax=axe[0,1])
_ = axe[0,2].imshow(vol[40, 100:412, 100:412], cmap='gray', vmin=vmin, vmax=vmax)# ; fig.colorbar(_, ax=axe[0,2])
_ = axe[0,3].imshow(vol[50, 100:412, 100:412], cmap='gray', vmin=vmin, vmax=vmax)# ; fig.colorbar(_, ax=axe[0,3])
_ = axe[1,0].imshow(vol[60, 100:412, 100:412],  cmap='gray', vmax=vmax)# ; fig.colorbar(_, ax=axe[1,3])
_ = axe[1,1].imshow(vol[70, 100:412, 100:412],  cmap='gray', vmax=vmax)# ; fig.colorbar(_, ax=axe[1,3])
_ = axe[1,2].imshow(vol[80, 100:412, 100:412],  cmap='gray', vmax=vmax)# ; fig.colorbar(_, ax=axe[1,3])
_ = axe[1,3].imshow(vol[110, 100:412, 100:412],  cmap='gray', vmax=vmax)# ; fig.colorbar(_, ax=axe[1,3])

======== CT Cone-Beam Geometry ======== number of angles: 720 number of detector elements (rows, cols): 704 x 704 angular range: 360.500702 degrees detector pixel size: 0.615057 mm x 0.615057 mm center detector pixel: 507.583130, 576.243469 sod = 537.500000 mm sdd = 990.804993 mm

======== CT Volume ======== number of voxels (x, y, z): 512 x 512 x 128 voxel size: 1.000000 mm x 1.000000 mm x 1.280000 mm FOV: [-256.000000, 256.000000] x [-256.000000, 256.000000] x [-81.919999, 81.919999]

SDD: 990.805 SOD: 537.5 detColCenter: 576.2434642032332 detCols: 704 detPhyOffsetX: -138.23 detPhyOffsetY: -96 detPhySize: 433 detPixOffsetX: -224.74346420323323 detPixOffsetY: -156.08314087759814 detRowCenter: 507.5831408775981 detRows: 704 detSpacing: 0.6150568181818182 device: cuda:0 nAngles: 720 nBins: 2 nPad: 0 volOffsetX: 0 volOffsetY: 0 volOffsetZ: 0 volPhyXDim: 512 volPhyYDim: 512 volPhyZDim: 163.84 volXDim: 512 volXSpacing: 1 volYDim: 512 volYSpacing: 1 volZDim: 128 volZSpacing: 1.28

kylechampley commented 1 month ago

That's weird. I simulated data with your geometry and I did not get an artifact. If you are able to send your data to me, I can take a look it at.

hws203 commented 1 month ago

@lc82111 I guess that your thinking ring artifact is a real ring of your target sample. As you can see the projection image, there is a circle tap on top of cylinder. Those border area of this tap may be a ring.

lc82111 commented 1 month ago

@kylechampley Thanks for your kind help. The raw projections can be found at this link:

https://drive.google.com/file/d/1xb7X9mQo2kcIhWD_UGi40JXHCmd-6zQq/view?usp=sharing

lc82111 commented 1 month ago

@lc82111 I guess that your thinking ring artifact is a real ring of your target sample. As you can see the projection image, there is a circle tap on top of cylinder. Those border area of this tap may be a ring.

@hws203 Thanks for the suggestion. However, I believe the ring artifact I'm observing is not a physical characteristic of the sample. The phantom is flawless, and the ring and central dark region are absent in the full-scan data.

kylechampley commented 1 month ago

@lc82111, the data in the link has no air scan projection. Could you send this and any other calibration files you have?

kylechampley commented 1 month ago

Well, I tried to reconstruct without the air scan data and I think I got something good. It is different from what you showed, but there is still that inner circle. I do think this inner circle is real. Here is my result.

image

And you can see that the inner cylinder is indeed in the raw data, so it is indeed real. See this projection:

image

kylechampley commented 1 month ago

Here is the script I used. I turned your dicom files into tiff files before running the reconstruction.

import numpy as np
from leap_preprocessing_algorithms import *
from leapctype import tomographicModels
leapct = tomographicModels()

angles = np.linspace(-180, 180, 720).astype(np.float32)
leapct.set_conebeam(720, 704*2, 704*2, 0.615057*0.5, 0.615057*0.5, 507.583130*2.0, 576.243469*2.0, angles, 537.5, 990.804993)
leapct.set_offsetScan(True)
leapct.set_truncatedScan(True)
leapct.set_default_volume(4.0) # just show a low resolution example
leapct.print_parameters()

g = leapct.load_projections(r'D:\tomography\issue92\raw.tif')
ROI = [30, 80, 950, 1000]
makeAttenuationRadiographs(leapct, g, None, None, ROI)

vol = leapct.allocateVolume(0)
leapct.FBP(g,vol)

leapct.display(vol)
lc82111 commented 1 month ago

AVG_504_120_30_air_cw

Above is the averaged air scan projection in png format.

I will try your script and report back. Thanks for your help.

lc82111 commented 1 month ago

image

image

To avoid ambiguity, I have added a yellow circle to highlight the artifact. Apologies for any confusion.

lc82111 commented 1 month ago

This ring artifact appears to be more prominent in RTK reconstruction as below:

image

Furthermore, I believe the ring artifact is likely attributed to the RTK DDF (Displaced Detector Filter) weighting, as evidenced by the abrupt changes in the weighted image:

image

Therefore, I utilized LEAP CT to validate my hypothesis.

kylechampley commented 1 month ago

I believe this is a cone-beam artifact caused by the void cylinder at the bottom of that projection. Notice that the diameter of your dark circle is about the same size as the air cylinder at the bottom of the projection.

I think you also have some beam hardening artifacts here. If you know something about the spectra, you could use LEAP's BHC algorithms to help mitigate this.

hws203 commented 1 month ago

I think so too.