Inconsistencies in CT Reconstruction using SIRT algorithm

[yunsuper123]

Expected Behavior

After implementing SIRT, I expect the CT reconstruction volume to be the same as the ground truth CT.

Actual Behavior

I slightly adjusted demo 1,2,3 to create the geometry and data of a CT .mha file, then implemented SIRT to reconstruct CT volume with demo 7. Although I tried to make sure all the geometry variables were right, I still didn't get the right result shape. The first image is the ground truth CT shown in 3D slicer, and the second is my reconstruction. It's quite possible that I might be overlooking something, so I was wondering if you could shed some light on the potential issues of this inconsistency. Thank you in advance for your assistance and for the amazing contributions you've made to our community.

Code to reproduce the problem (If applicable)

import tigre
import numpy as np
from tigre.utilities import sample_loader
from tigre.utilities import CTnoise
import tigre.algorithms as algs
from matplotlib import pyplot as plt
import imageio
import os
import SimpleITK as sitk

def save_mha(image, filename):
    itkimage = sitk.GetImageFromArray(image)
    sitk.WriteImage(itkimage, filename)

# Lets create a geomtry object
geo = tigre.geometry()

# Distances
geo.DSD = 946.746  
geo.DSO = 538.52  
# Detector parameters
geo.nDetector = np.array([128, 128])  
#eo.dDetector = np.array([500/128, 500/128])  
geo.dDetector = np.array([1000/128, 1000/128])  
geo.sDetector = geo.nDetector * geo.dDetector  
# Image parameters
geo.nVoxel = np.array([225, 512, 512]) 
geo.sVoxel = np.array([281.25, 360, 360]) 
geo.dVoxel = geo.sVoxel / geo.nVoxel 
# Offsets
geo.offOrigin = np.array([0, 0, 0])  
geo.offDetector = np.array([0, 0])  

# Auxiliary
geo.accuracy = 0.5  # Variable to define accuracy of
geo.COR = 0  # y direction displacement for
geo.rotDetector = np.array([0, 0, 0])  # Rotation of the detector, by
geo.mode = "cone"  # Or 'parallel'. Geometry type.

#%% Load data and generate projections
# define angles
angles = np.linspace(0, 2 * np.pi, 72)

groundtruth = sitk.ReadImage("data.mha")
groundtruth = sitk.GetArrayFromImage(groundtruth)

# Simulate forward projection.
# To match with mathematical notation, the projection operation is called Ax
projections = tigre.Ax(groundtruth, geo, angles)

# Optional arguments for all of them
# ==========================================================================
lmbda = 1
lambdared = 0.999
initmode = None
verbose = True
qualmeas = ["RMSE", "SSD"]

# SIRT and SART both have no extra input parameters.
# =========================================================================
imgSIRT, qualitySIRT = algs.sirt(
    projections,
    geo,
    angles,
    20,
    lmbda=lmbda,
    lmbda_red=lambdared,
    verbose=verbose,
    Quameasopts=qualmeas,
    computel2=True,
)

print(imgSIRT.shape, imgSART.shape)
print(qualitySIRT.shape, qualitySART.shape)


# Save the SIRT reconstruction result
save_mha(imgSIRT, 'SIRT.mha')

Specifications

• Python version: 3.9.18
• OS: Windows 10
• CUDA version: nvcc: Build on Tue_Aug_15_22:09:35_Pacific_Daylight_Time_2023 Cuda compilation tools, release 12.2, V12.2.140 Build cuda_12.2.r12.2/compiler.33191640_0