DICOM-CT-PD (3rd Generation CT)#
[1]:
from __future__ import annotations
from pytomography.metadata import ObjectMeta
from pytomography.algorithms import SART
from pytomography.io.CT import dicom_ct_pd
from pytomography.projectors.CT import CTGen3SystemMatrix
import matplotlib.pyplot as plt
import os
First, lets specify the paths of all the projection data:
[2]:
save_path = '/disk1/ct_proj/DICOM-CT-PD_FD' # Path where you downloaded the CT-PD data
paths = [os.path.join(save_path, file) for file in os.listdir(save_path) if '.dcm' in file]
We can load the projections and corresponding metadata as follows
This may take a couple minutes as the dataset is huge
[3]:
proj, proj_meta = dicom_ct_pd.get_projections_and_metadata_gen3(paths)
We can view some of the projections:
[4]:
plt.subplots(5,1,figsize=(5,5))
plt.subplot(511)
plt.pcolormesh(proj[0].cpu().T, cmap='Greys_r')
plt.axis('off')
plt.subplot(512)
plt.pcolormesh(proj[1400].cpu().T, cmap='Greys_r')
plt.axis('off')
plt.subplot(513)
plt.pcolormesh(proj[4800].cpu().T, cmap='Greys_r')
plt.axis('off')
plt.subplot(514)
plt.pcolormesh(proj[6400].cpu().T, cmap='Greys_r')
plt.axis('off')
plt.subplot(515)
plt.pcolormesh(proj[8800].cpu().T, cmap='Greys_r')
plt.axis('off')
plt.show()
Now lets set up our object space metadata
[5]:
dx = dy = dz = 1.0 # mm
Nx = Ny = Nz = 512
object_meta = ObjectMeta(dr=(dx,dx,dx), shape=(Nx,Ny,Nz))
We can create the system matrix as follows. The N_splits argument can save on memory if you are using a small number of subsets or have little GPU RAM. The device='cpu' means that the projections are expected on CPU (the internal projector code will still use pytomography.device and transfer the projections to GPU if needed)
[6]:
system_matrix = CTGen3SystemMatrix(object_meta, proj_meta, N_splits=1, device='cpu')
Here we reconstruct using SART with 3 iterations and 40 subsets (this takes around 10 mins)
[7]:
recon_algorithm = SART(system_matrix, proj)
recon = recon_algorithm(n_iters=3, n_subsets=40)
Now we can look at an axial slice:
[8]:
plt.figure(figsize=(5.5,5))
plt.imshow(recon[:,255].cpu().numpy().T * 10, cmap='Greys_r', vmin=0.15, vmax=0.26, interpolation='gaussian', origin='lower')
plt.ylim(10,500)
plt.axis('off')
plt.colorbar(label='cm$^{-1}$')
plt.title('PyTomography: OS-SART(3it40ss)')
If anyone wants to implement FBP for the geometry of this scanner, please do and submit a pull request! (see the
CTConeBeamFlatPanelSystemMatrixbackward method for how filtered back projection should be implemented)If anyone wants to write a DICOM exporter for these CT images, please do, and submit a pull request!