reconstruction.py

import os
from PIL import Image
import numpy as np
import astra
import matplotlib.pyplot as plt

## Load Projection Data:
# List of paths to your projection images
directory_path = r"E:\Tiankuo\Dr_Chow_gnptab_mice\1. raw\month 0\870" #["path1.tif", "path2.tif", ...]
image_files = [f for f in os.listdir(directory_path) if f.endswith('.tif')]
image_paths = [os.path.join(directory_path,f) for f in image_files]

projections = [np.array(Image.open(p)) for p in image_paths]
projections = np.transpose(projections, (1, 0, 2))
projections = np.stack(projections, axis=0)
# Transpose the projections
#projections = np.transpose(projections, (1, 0, 2))


## Set Up the Geometries:
# Geometry parameters
det_spacing = [0.010609024, 0.010609024]  # pixel_size_y, pixel_size_x in mm
det_count = projections.shape[1:3]  # detector pixel count
angles = np.linspace(0, 2*np.pi, projections.shape[0], endpoint=False)  # projection angles
source_origin = 93.055 # distance from the X-ray source to the origin (rotation center)  mm
origin_det = 65.793  # distance from the origin to the detector  mm 

# Create projection geometry
proj_geom = astra.create_proj_geom('cone', det_spacing[0], det_spacing[1], det_count[0], det_count[1], angles, source_origin, origin_det)
 
## Use the ASTRA toolbox to perform the reconstruction:
# Create a data object for the projection data
proj_id = astra.data3d.create('-sino', proj_geom, projections)
 
# Create a volume geometry and data object
vol_geom = astra.create_vol_geom(det_count[0], det_count[1], det_count[0])  # Assuming a cubic volume
rec_id = astra.data3d.create('-vol', vol_geom)

# Set up and run the FDK algorithm
cfg = astra.astra_dict('FDK_CUDA_3D')
cfg['ProjectionDataId'] = proj_id
cfg['ReconstructionDataId'] = rec_id
alg_id = astra.algorithm.create(cfg)
astra.algorithm.run(alg_id)

# Get the reconstruction data
reconstruction = astra.data3d.get(rec_id)

# Cleanup
astra.algorithm.delete(alg_id)
astra.data3d.delete(proj_id)
astra.data3d.delete(rec_id)

 
## Visualize the Reconstruction:
plt.imshow(reconstruction[reconstruction.shape[0] // 2, :, :], cmap='gray')

#plt.imshow(reconstruction, cmap='gray')
plt.show()