Helical CT issues

[willytown]

Hi all, sorry if this is a easy fix but I cannot figure it out. I am trying to run reconstructions of helical cone beam CT data, my code runs, however the result of the central slice I get is:

when its supposed to be similar to:

my code is as follows:

---

# Python base libraries
import numpy as np
import matplotlib.pyplot as plt
import csv

# TIGRE libraries
import tigre
import tigre.algorithms as algs

# Library to track time
import time

import tifffile as tiff

print("script starting")

#=================================================
# Defining parameters from meta data
#=================================================

# Distance source to detector in mm
D = 390

# Distance source to sample in mm
S = 30

# The number of detector elements in the square detector
N = 3040

#=================================================
# Loading the data of central projections
#=================================================  

# Using a restricted number of projections
start_projection_number = 0000
end_projection_number = 1799
num_projections = end_projection_number - start_projection_number + 1

# Initialize 2D array: rows = projections, cols = detector elements
projections = np.zeros((num_projections, N, N), dtype=np.uint16)
        
for idx, i in enumerate(range(start_projection_number, end_projection_number + 1)):
    file_number = str(i).zfill(4) 
    file_path = f"/mnt/iusers01/fse-ugpgt01/phy01/u63544jw/scratch/tigre_recon_data/full_rotation_data/expt_H00{file_number}.raw"
    
    with open(file_path, 'rb') as file:
        data = np.frombuffer(file.read(), dtype=np.uint16)
        projections[idx] = data.reshape((N, N))

projections = projections.astype(np.float32)

print("projections loaded")

#=================================================
# Reading the meta data
#=================================================

# Open CSV file
with open('/mnt/iusers01/fse-ugpgt01/phy01/u63544jw/scratch/tigre_recon_data/restricted_meta_data.csv', mode='r', newline='', encoding='utf-8') as file:
    reader = csv.DictReader(file)  # Reads as a dictionary where headers are keys
    columns = {}  # Dictionary to store column data

    # Initialize columns based on headers
    # Format is: 'file_ID', 'filename', 'creation_time', 'sample_theta[deg]', 'camera_x_pos[mm]', 
    #'camera_y_pos[mm]', 'sample_y_pos[mm]', 'sample_z_pos[mm]', 'roi_x_pos[mm]', 'roi_y_pos[mm]', 
    #'stage was moving', 'movement duration[ms]', 'exposure time[ms]', 'total time[ms]', 
    #'movement fraction', 'frames averaged', 'binning', 'source values changed'
    for header in reader.fieldnames:
        columns[header] = []

    # Read each row and store data by column
    for row in reader:
        for header in reader.fieldnames:
            columns[header].append(row[header])
            
print("meta data loaded")
            
# Converting meta data to arrays and converting degrees to radians
theta = np.array(columns['sample_theta[deg]'], dtype=np.float32)*(np.pi/(180))
x_cam = np.array(columns['camera_x_pos[mm]'], dtype=np.float32)
z_sample = np.array(columns['sample_z_pos[mm]'], dtype=np.float32)

#=================================================
# Transforming geometry
#=================================================

# Meta data geometry needs to be transformed to a helix relative to the sample
# The real geometry is a static source and detector with a moving sample (for a single helix)
# The detector actually moves but to simplify only data with a static detector will be taken

# Find indices where camera_x_pos[mm] == 0 for a single helix
valid_indices = np.where(x_cam == 0)[0]

# Masking data
theta = theta[valid_indices]
z_sample = z_sample[valid_indices]
projections = projections[valid_indices]

# Checking projections shape
print(projections.shape)

# Number of masked projections
num_projections = len(projections)

# Reordering the projections to (N, N, num_projections) - Maybe needed?
projections_test = np.transpose(projections, (1, 2, 0))

# Making the inital z = 0
z_sample = z_sample - z_sample[0]

# The static positions of the source and detector (mm)
# This is done with the positive y direction from the origin (center of rotation/sample) towards the source
source_static = np.array([0,S,0])
detector_static = np.array([0,-D,0])

# Initialising arrays for transformed source and detector positions
source_positions = np.zeros((num_projections, 3))
detector_positions = np.zeros((num_projections, 3))

# Applying transformation
for i in range(num_projections):
    
    # Rotation matrix
    Rz = np.array([
        [np.cos(theta[i]), -np.sin(theta[i]), 0],
        [np.sin(theta[i]),  np.cos(theta[i]), 0],
        [0,                 0,                1]
    ])
    
    # Rotate and translate
    source_positions[i] = (Rz @ source_static) + np.array([0, 0, -z_sample[i]])
    detector_positions[i] = (Rz @ detector_static) + np.array([0, 0, -z_sample[i]])


#=================================================
# Reconstructions using Tigre
#=================================================

# Set up geometry
geo = tigre.geometry()

# Set detector parameters
geo.nDetector = np.array([N, N])      # Detector elements (width, height)
geo.dDetector = np.array([0.139, 0.139])      # mm per detector element 
geo.sDetector = geo.nDetector * geo.dDetector

# Set volume parameters
geo.nVoxel = np.array([1024, 1024, 512]) # Number of pixels
geo.sVoxel = np.array([200, 200, 100])  # mm size of the volume
geo.dVoxel = geo.sVoxel / geo.nVoxel # mm size per voxel

# Offset parameters
geo.offOrigin = np.array([0, 0, 0])
geo.offDetector = np.array([0, 0])

# Distance source to origin and source to detector
geo.DSO = S
geo.DSD = D

# Assign source and detector positions
geo.SourcePositions = source_positions.astype(np.float32)
geo.DetectorPositions = detector_positions.astype(np.float32)

# FDK reconstruction

# Check
print("projections_test.shape:", projections_test.shape)
print("theta.shape:", theta.shape)
print("geo.SourcePositions.shape:", geo.SourcePositions.shape)
print("geo.DetectorPositions.shape:", geo.DetectorPositions.shape)

# Tracking the time taken to reconstruct by FDK
# Start time
start_fdk = time.perf_counter()

# Run FDK (helical mode will be inferred from position arrays)
recon_fdk = algs.fdk(projections, geo, theta)

# End time
end_fdk = time.perf_counter()

# Printing time taken
print(f"FDK compute time: {end_fdk - start_fdk:.6f} seconds")

# Saving to a .tiff stack
# First need to rescale data 
recon_fdk_norm = (recon_fdk - recon_fdk.min()) / (recon_fdk.max() - recon_fdk.min())
recon_fdk_uint16 = (recon_fdk_norm * 65535).astype('uint16')

# Writing to file
tiff.imwrite('/mnt/iusers01/fse-ugpgt01/phy01/u63544jw/scratch/tigre_recon/out_test/recon_fdk_stack_highres.tiff', np.transpose(recon_fdk_uint16, (2, 1, 0)))

print("data written to .tiff stack")

---

And I know that the projections load properly as I can visualise them, e.g.

The meta data is in:
restricted_meta_data.csv
and my code outputs:

---

script starting
projections loaded
meta data loaded
(100, 3040, 3040)
projections_test.shape: (3040, 3040, 100)
theta.shape: (100,)
geo.SourcePositions.shape: (100, 3)
geo.DetectorPositions.shape: (100, 3)
FDK compute time: 125.237112 seconds
data written to .tiff stack

I have tried everything I can think of, and cannot figure it out, I have even committed the sin of asking chatGPT, however nothing has worked. Any help or pointers would be greatly appreciated.

[AnderBiguri]

Certainly chatGPT won't know about TIGRE, of course! hehe.

It clearly is something wrong with the geometry. Is there any chance you can reduce this largie peice of code as much as possible? That will make my job of reviewing your code a bit easier, and you more likely to get help :)

Ander

[willytown]

DSO = 30mm, DSD = 390mm, in the simulated data, OS-SART produced

after one iteration, but after 10:

where the phantom is 5 nested spheres. But when I changed the DSO to 60mm it produced:

after 10 iterations

[AnderBiguri]

Thats quite unexpected actually, can you produce code for the simulation, so I can have a try?

[willytown]

Yeah sure! The only other things I changed was the image size, as a larger image would fit when I increased the DSO

# Python base libraries
import numpy as np
import matplotlib.pyplot as plt

# TIGRE libraries
import tigre
import tigre.algorithms as algs

# Library to track time
import time

import tifffile as tiff

print("script starting")

#=================================================
# Defining parameters from meta data
#=================================================

# Distance source to detector in mm
D = 390 

# Distance source to sample in mm
S = 30 # Changed to 60 for second run 

# The number of detector elements in the square detector
N = 3040

num_projections = 100

#=================================================
# Tigre geometry
#=================================================

# Set up geometry
geo = tigre.geometry()

# Set detector parameters
geo.nDetector = np.array([N, N])      # Detector elements (width, height)
geo.dDetector = np.array([0.139, 0.139])      # mm per detector element 
geo.sDetector = geo.nDetector * geo.dDetector

# Set volume parameters
geo.nVoxel = np.array([500, 500, 100]) # Number of voxels
geo.sVoxel = np.array([25, 25, 50])  # mm size of the volume - changed to 50 for all on second run
geo.dVoxel = geo.sVoxel / geo.nVoxel # mm size per voxel

# Offset parameter for offset of detector
geo.offDetector = np.array([0, 0])

# Doing it manually
z_test = np.linspace(-17.176, 17.176, num_projections)
angle = np.linspace(0, 2*np.pi, num_projections)

# Image offset - makes it helical
geo.offOrigin = np.zeros((num_projections, 3))
geo.offOrigin[:,2] = z_test

# Distance source to origin and source to detector
geo.DSO = S
geo.DSD = D

#=================================================
# Creating a sphere phantom with rings
#=================================================

# Function to create a spherical phantom with concentric rings inside
def stalagmite_phantom(n_voxels):
    
    # Creating grid of points
    nx, ny, nz = n_voxels
    x = np.linspace(-1, 1, nx)
    y = np.linspace(-1, 1, ny)
    z = np.linspace(-1, 1, nz) 
    X, Y, Z = np.meshgrid(x, y, z, indexing='ij')

    # Spherical radius 
    R = 0.75
    r = np.sqrt(X**2 + Y**2 + Z**2)

    # Base sphere mask
    sphere = r <= R
    
    # Number of concentric rings
    num_rings = 10

    # Create concentric ring pattern
    radial_fraction = r / R
    #ring_pattern = 0.5 * (1 + np.sin(2 * np.pi * num_rings * radial_fraction))
    ring_pattern = ((np.floor(num_rings * radial_fraction) % 2) == 0).astype(float)

    # Apply rings only inside sphere
    phantom = sphere * ring_pattern.astype(np.float32)

    return phantom

phantom = stalagmite_phantom(geo.nVoxel)

# project data
data = tigre.Ax(phantom, geo, angle)

#=================================================
# Reconstruction
#=================================================

# Reconstruct
recon_fdk = algs.fdk(data, geo, angle)

# Saving to a .tiff stack
# First need to rescale data 
recon_fdk = (recon_fdk - recon_fdk.min()) / (recon_fdk.max() - recon_fdk.min())
recon_fdk = (recon_fdk * 65535).astype('uint16')

# Writing to file
tiff.imwrite('/mnt/iusers01/fse-ugpgt01/phy01/u63544jw/scratch/tigre_recon/out_test/isometric_fdk.tiff', np.transpose(recon_fdk, (2, 1, 0)))

# Number of iterations
n_iter = np.array([1,5,10, 30])

# for loop over number of iterations
for n in n_iter:
      
    # Start time
    start = time.perf_counter()

    # Running algorithm
    recon = algs.ossart(data, geo, angle, n)

    # End time
    end = time.perf_counter()
    
    print(f"OS-SART compute time for {n} iterations: {end - start:.6f} seconds")
    
    # Rescaling
    recon = (recon - recon.min()) / (recon.max() - recon.min())
    recon = np.clip((recon * 65535), 0, 65535).astype('uint16')
    
    # Writing to file
    tiff.imwrite(f"/mnt/iusers01/fse-ugpgt01/phy01/u63544jw/scratch/tigre_recon/out/isometric_cgls_{n}.tiff", np.transpose(recon, (2, 1, 0)))
    
    # Deleting large variables for memory management
    del recon

[AnderBiguri]

Sorry, I've been a bit busy. Any progress?

[willytown]

This was for my Masters project at university which has been submitted, and so I have stopped working on it for now. I may revisit it over summer, and if I figure anything out, I will update the thread.