FDK based reconstruction

[karthik162]

Friends,

I am trying to make a cone beam based reconstruction.

Below is the code, however when executed, the slice outline is distorted like projection is taken as parallel beam, instead of cone beam.

import tigre
import numpy as np
import imageio
import os

# Configuration
# -------------------------------------------------------------------------------------
# Geometry parameters
distance_source_detector = 860                        # Distance from source to detector (DSD) in mm
distance_source_origin = 740                          # Distance from source to origin (DSO) in mm

# Detector parameters
n_detector = np.array((470, 1669))                     # Number of pixels (nDetector)
d_detector = np.array((0.429, 0.429))                 # Size of each pixel (dDetector) in mm
s_detector = n_detector * d_detector                  # Total size of the detector (sDetector) in mm

# Image parameters
n_voxel = np.array((470, 1669, 1669))                   # Number of voxels (nVoxel)
s_voxel = np.array((107, 382, 382))           # Size of each voxel (sVoxel) in mm
d_voxel = s_voxel / n_voxel                           # Size of each voxel (dVoxel) in mm

# Offset correction parameters
off_origin = np.array((0, 0, 0))                      # Offset of image from origin (offOrigin) in mm
off_detector = np.array((0, 0))                       # Offset of detector (offDetector) in mm

# Geometry setup
tigre_geo = tigre.geometry()
tigre_geo.DSD = distance_source_detector                # Distance from source to detector (mm)
tigre_geo.DSO = distance_source_origin                  # Distance from source to object (mm)
tigre_geo.nDetector = n_detector                        # Number of pixels (detector)
tigre_geo.dDetector = d_detector                        # Size of each pixel (mm)
tigre_geo.sDetector = s_detector                        # Total size of the detector (mm)
tigre_geo.nVoxel = n_voxel                              # Number of voxels
tigre_geo.sVoxel = s_voxel                              # Total size of the image (mm)
tigre_geo.dVoxel = d_voxel                              # Size of each voxel (mm)
tigre_geo.offOrigin = off_origin                        # Offset of image from origin (mm)
tigre_geo.offDetector = off_detector                    # Offset of detector (mm)
tigre_geo.mode = 'cone'

# Print all parameters to verify before reconstruction
print("Verification of geometry and reconstruction parameters:")
print(f"Distance Source to Detector (DSD): {tigre_geo.DSD} mm")
print(f"Distance Source to Origin (DSO): {tigre_geo.DSO} mm")
print(f"Detector - Number of pixels: {tigre_geo.nDetector}")
print(f"Detector - Pixel size: {tigre_geo.dDetector} mm")
print(f"Detector - Total size: {tigre_geo.sDetector} mm")
print(f"Image - Number of voxels: {tigre_geo.nVoxel}")
print(f"Image - Voxel size: {tigre_geo.dVoxel} mm")
print(f"Image - Total size: {tigre_geo.sVoxel} mm")
print(f"Offset - Image offset from origin: {tigre_geo.offOrigin}")
print(f"Offset - Detector offset: {tigre_geo.offDetector}")
print(f"Geometry Mode: {tigre_geo.mode}")

# Load projections
input_dir = 'C:/ALL_PROJECTIONS/2ND-SET/PUMPKIN-2/SCALED'  # Input directory for projections
num_of_projections = 1800                                 # Number of projections
angles = np.linspace(0, 2 * np.pi, num=num_of_projections, endpoint=False)  # Angles

projections = []
for i in range(num_of_projections):
    filename = os.path.join(input_dir, f'{i:04d}.tif')  # Files named 0000.tif to 1799.tif
    img = imageio.v2.imread(filename)                     # Read each projection image
    projections.append(img)
projections = np.array(projections)

# Reconstruction using FDK
print("Starting FDK reconstruction...")
img_sirt = tigre.algorithms.fdk(projections, tigre_geo, angles)  # Set desired number of iterations
print("FDK reconstruction completed.")

# Saving the reconstructed image slices in 32-bit format
output_dir = 'C:/RECONSTRUCTED_SLICES/'                    # Output directory
os.makedirs(output_dir, exist_ok=True)                     # Create directory if it doesn't exist

# Save each slice of the reconstructed image as 32-bit TIFF
for i in range(img_sirt.shape[0]):
    output_filename = os.path.join(output_dir, f'reconstructed_slice_{i:04d}.tif')
    imageio.v2.imwrite(output_filename, img_sirt[i, :, :].astype(np.float32))  # Save each slice as 32-bit
    print(f"Saved {output_filename}")

print(f"Reconstructed slices saved to: {output_dir}")

[AnderBiguri]

there are 2 things that I would say are going wrong here:

1- You have definitely not applied the log transform to the projection data. The background should be black, not white.
2-its plausible (given the ghosting artifacts in the boundary of the object) that your angles should be reversed, or some other similar mismatch of axis direction w.r.t TIGRE axis.

[karthik162]

Thanks very much for the reply. It is very useful.
With regard to the ghosting artifacts in the boundary, I tried reconstruction using OsCAR also.
I found that, there is no ghosting artifact. (I think)

From Oscar

From TIGRE using the above program, as you can see. The Ghosting artifact is not present in Oscar, while it is present in TIGRE.
Is there any solution?
Compressed dataset - link

I have uploaded the dataset in Google drive and given the link also. Please review the same.
https://drive.google.com/drive/folders/1EXlriObYIQYXfIAgFz9x5BGfbSutgQ0f?usp=sharing

[karthik162]

(tigre_env) C:\Users\admin>python C:\ASTRA-TOOLBOX\ALL-OLD\TIGRE\33TIGRE-PUMPKIN-BPF.py
Verification of geometry and reconstruction parameters:
Distance Source to Detector (DSD): 860 mm
Distance Source to Origin (DSO): 740 mm
Detector - Number of pixels: [ 470 1669]
Detector - Pixel size: [0.4 0.4] mm
Detector - Total size: [188. 667.6] mm
Image - Number of voxels: [ 470 1669 1669]
Image - Voxel size: [1. 1. 1.] mm
Image - Total size: [ 470 1669 1669] mm
Offset - Image offset from origin: [0 0 0]
Offset - Detector offset: [0 0]
Geometry Mode: cone
Starting FDK reconstruction...
FDK reconstruction completed.
Saving reconstructed slices as 32-bit TIFF files...
Saved slice 1/470 to C:/RECONSTRUCTED_SLICES/reconstructed_slice_0000.tif
Saved slice 2/470 to C:/RECONSTRUCTED_SLICES/reconstructed_slice_0001.tif

[karthik162]

(tigre_env) C:\Users\admin>python --version
Python 3.10.15
(tigre_env) C:\Users\admin>conda list tigre

packages in environment at C:\Users\admin.conda\envs\tigre_env:

Name Version Build Channel

pytigre 2.4.0 pypi_0 pypi

(tigre_env) C:\Users\admin>

[AnderBiguri]

As I said before, you certainly are not putting everything correctly into TIGRE, but its something you will need to find out, as I don't know how your scanner is, so I cant translate that. Almost surely you just need angles=-angles; as I said before. Similarly I am not sure if you fixed it already, but the log trasnform tends to be needed proj=-log(proj/max(proj(:)) when using real data.

I don't know how OSCAR works, but maybe you just happen to have the correct axis direction by chance compared to Oscar, or maybe it does the log transform for you, I don't know. Just to clarify that I myself (and the community) have reconstructed hundreds of CT scans with TIGRE, so there is no issue with TIGRE itself, but in how you are inputting the data.

If I have some time free next week I can try to test your code with your data, but I am almost sure its just a minor change, like the ones I already suggested.

Is there any solution?

Yes, to fix the way you define the geometry and pass the data to TIGRE.

[karthik162]

Thanks a lot for quick suggestions.
I made the modifications, the slice improved, in terms of reductions in ghosting and better contrast within the slice.

However, is it possible to reduce further?
import tigre
import numpy as np
import imageio
import os

Configuration

-------------------------------------------------------------------------------------

Geometry parameters

distance_source_detector = 860 # Distance from source to detector (DSD) in mm
distance_source_origin = 740 # Distance from source to origin (DSO) in mm

Detector parameters

n_detector = np.array((470, 1669)) # Number of pixels (nDetector)
d_detector = np.array((0.4, 0.4)) # Size of each pixel (dDetector) in mm
s_detector = n_detector * d_detector # Total size of the detector (sDetector) in mm

Image parameters

n_voxel = np.array((470, 1669, 1669)) # Number of voxels (nVoxel)
s_voxel = np.array((188, 667, 667)) # Size of each voxel (sVoxel) in mm
d_voxel = s_voxel / n_voxel # Size of each voxel (dVoxel) in mm

Offset correction parameters

off_origin = np.array((0, 0, 0)) # Offset of image from origin (offOrigin) in mm
off_detector = np.array((0, 0)) # Offset of detector (offDetector) in mm

Geometry setup

tigre_geo = tigre.geometry()
tigre_geo.DSD = distance_source_detector # Distance from source to detector (mm)
tigre_geo.DSO = distance_source_origin # Distance from source to object (mm)
tigre_geo.nDetector = n_detector # Number of pixels (detector)
tigre_geo.dDetector = d_detector # Size of each pixel (mm)
tigre_geo.sDetector = s_detector # Total size of the detector (mm)
tigre_geo.nVoxel = n_voxel # Number of voxels
tigre_geo.sVoxel = s_voxel # Total size of the image (mm)
tigre_geo.dVoxel = d_voxel # Size of each voxel (mm)
tigre_geo.offOrigin = off_origin # Offset of image from origin (mm)
tigre_geo.offDetector = off_detector # Offset of detector (mm)
tigre_geo.mode = 'cone'
tigre_geo.accuracy = 0.5
tigre_geo.COR = 0
tigre_geo.rotDetector = np.array([0, 0, 0])

Print all parameters to verify before reconstruction

print("Verification of geometry and reconstruction parameters:")
print(f"Distance Source to Detector (DSD): {tigre_geo.DSD} mm")
print(f"Distance Source to Origin (DSO): {tigre_geo.DSO} mm")
print(f"Detector - Number of pixels: {tigre_geo.nDetector}")
print(f"Detector - Pixel size: {tigre_geo.dDetector} mm")
print(f"Detector - Total size: {tigre_geo.sDetector} mm")
print(f"Image - Number of voxels: {tigre_geo.nVoxel}")
print(f"Image - Voxel size: {tigre_geo.dVoxel} mm")
print(f"Image - Total size: {tigre_geo.sVoxel} mm")
print(f"Offset - Image offset from origin: {tigre_geo.offOrigin}")
print(f"Offset - Detector offset: {tigre_geo.offDetector}")
print(f"Geometry Mode: {tigre_geo.mode}")

Load projections

input_dir = 'C:/ALL_PROJECTIONS/2ND-SET/PUMPKIN-2/SCALED' # Input directory for projections
num_of_projections = 1800 # Number of projections

Angles - without tilt correction

angles = np.linspace(0, 2 * np.pi, num=num_of_projections, endpoint=False)

Define tilt correction in radians (e.g., 10 degrees to radians)

tilt_correction = 10 * np.pi / 180 # Example: 10 degrees tilt correction (in radians)

Apply tilt correction to each angle

corrected_angles = angles + tilt_correction

Ensure angles are within the range [0, 2*pi]

corrected_angles = np.mod(corrected_angles, 2 * np.pi)

angles = -angles

Load projections with error handling for missing files

projections = []
for i in range(num_of_projections):
filename = os.path.join(input_dir, f'{i:04d}.tif') # Files named 0000.tif to 1799.tif
if os.path.isfile(filename):
img = imageio.v2.imread(filename) # Read each projection image
# Apply log transfer: proj = -log(proj/max(proj))
img_log = -np.log(img / np.max(img)) # Logarithmic transformation
projections.append(img_log)
else:
print(f"Warning: {filename} not found. Skipping this projection.")
projections = np.array(projections)

Check if projections were loaded successfully

if projections.size == 0:
raise ValueError("No projections were loaded. Please check the input directory and file naming.")

Reconstruction using FDK

print("Starting FDK reconstruction...")
img_fdk = tigre.algorithms.fdk(projections, tigre_geo, angles)
print("FDK reconstruction completed.")

Saving the reconstructed image slices in 32-bit format

output_dir = 'C:/RECONSTRUCTED_SLICES/' # Output directory
os.makedirs(output_dir, exist_ok=True) # Create directory if it doesn't exist

Save each slice of the reconstructed image as 32-bit TIFF

print("Saving reconstructed slices as 32-bit TIFF files...")
for i in range(img_fdk.shape[0]):
output_filename = os.path.join(output_dir, f'reconstructed_slice_{i:04d}.tif')
imageio.v2.imwrite(output_filename, img_fdk[i, :, :].astype(np.float32)) # Save each slice as 32-bit
print(f"Saved slice {i + 1}/{img_fdk.shape[0]} to {output_filename}")

print(f"All reconstructed slices have been saved in the directory: {output_dir}")

[AnderBiguri]

I'll check the code whan I have time.

However, to be clear: yes, you can get exactly the same reconstructions that you get from the machine, from Astra or from Oscar using Tigre, and if you are not is because there are still mistakes in the code. Keep playing with the parameters, something is still off.