Chunked data shows an error only when offOrigin[0] is positive.

[X-Strahl]

Expected Behavior

With FDK reconstructing subgroups of z-slices and then stacking them should lead to the same result as reconstructing the same volume.

Actual Behavior

When comparing a non chunked and chunked and then combined reconstruction, there is only a difference in the chunked and non-chunked recon, in the bottom chunk. The top chunk shows no error while the bottom chunk shows an error on the order of 0.001-0.003 to the full recon.

Code to reproduce the problem

from tigre.utilities import sample_loader
import tigre
import tigre.algorithms as algs
import numpy as np

####FUNCTIONS####
def split_geometry_in_two_chunks(full_volume_tigre_geometry):
    """Calculate a top and bottom geometry object based on the full size volume"""
  
    full_size_nVox = full_volume_tigre_geometry.nVoxel
    full_size_sVox = full_volume_tigre_geometry.sVoxel
    full_size_dVox = full_volume_tigre_geometry.dVoxel

    n_vox_z_top_chunk = full_size_nVox[0] // 2
    n_vox_z_bottom_chunk = full_size_nVox[0] // 2

    assert n_vox_z_top_chunk + n_vox_z_bottom_chunk == full_size_nVox[0]

    n_vox_z_slice = np.array([n_vox_z_top_chunk, n_vox_z_bottom_chunk])

    # half way point in mm of the full size volume
    # we need this because offOrigin is calculated relative to this position
    half_z_full_volume_mm = full_size_dVox[0] * full_size_nVox[0] / 2.0

    z_offsets_each_chunk = []
    sum_previous_chunk_z_mm = np.double(0.0)
    # this loop calculates the z offset (in mm) for each chunk
    for idx, num_z_vox in enumerate(n_vox_z_slice):
        half_z_this_chunk_mm = np.double(full_size_dVox[0]) * np.double(num_z_vox) / np.double(2.0)
        z_offset = (np.double(sum_previous_chunk_z_mm) + np.double(half_z_this_chunk_mm)) - np.double(half_z_full_volume_mm)
        z_offsets_each_chunk.append(z_offset)
        sum_previous_chunk_z_mm += np.double(full_size_dVox[0]) * np.double(num_z_vox)

    # now use the offsets we calculated to make two new geometry objects
    top_chunk_geometry = tigre.geometry_default(nVoxel=full_size_nVox)
    top_chunk_geometry.nVoxel = np.array([ n_vox_z_top_chunk, full_size_nVox[1], full_size_nVox[2]])
    top_chunk_geometry.sVoxel = top_chunk_geometry.nVoxel * full_size_dVox
    top_chunk_geometry.offOrigin = np.array([ z_offsets_each_chunk[0], 0.0,0.0])
    top_chunk_geometry.dVoxel = full_size_dVox

    bottom_chunk_geometry = tigre.geometry_default(nVoxel=full_size_nVox)
    bottom_chunk_geometry.nVoxel = np.array([ n_vox_z_top_chunk, full_size_nVox[1], full_size_nVox[2]])
    bottom_chunk_geometry.sVoxel = top_chunk_geometry.nVoxel * full_size_dVox
    bottom_chunk_geometry.offOrigin = np.array([z_offsets_each_chunk[1],0.0, 0.0])
    bottom_chunk_geometry.dVoxel = full_size_dVox

    return top_chunk_geometry, bottom_chunk_geometry
####FUNCTIONS_END####    
      

#Load data
head_phantom = sample_loader.load_head_phantom(np.array([512,512,512]))
tigre_geometry = tigre.geometry_default(nVoxel=head_phantom.shape)
angles = np.linspace(0, 2*np.pi, 360) #Create 360 angles

#Define output geometry
# Inverse of magnification - 1 voxel maps to 1 pixel on detector
voxel_scale = tigre_geometry.DSO/tigre_geometry.DSD
tigre_geometry.nDetector=np.array([512,512])
tigre_geometry.sDetector=np.array([409.6,409.6])

tigre_geometry.dVoxel =  voxel_scale * np.array(
                [
                    tigre_geometry.dDetector[0],
                    tigre_geometry.dDetector[1],
                    tigre_geometry.dDetector[1],
                ]
            )

tigre_geometry.sVoxel = tigre_geometry.dVoxel * tigre_geometry.nVoxel
projections_default = tigre.Ax(head_phantom, tigre_geometry, angles) #Create projections

#Reconstruct
reconstruction_raw = algs.fdk(projections_default, tigre_geometry, angles, filter='ram_lak')
nz, ny, nx = reconstruction_raw.shape

#Create chunks
top_chunk_geometry, bottom_chunk_geometry = split_geometry_in_two_chunks(tigre_geometry)
#reconstruct
reconstruction_top_chunk = algs.fdk(projections_default, top_chunk_geometry, angles, filter='ram_lak')
reconstruction_bottom_chunk = algs.fdk(projections_default, bottom_chunk_geometry, angles, filter='ram_lak')
chunked_volume = np.concatenate([reconstruction_top_chunk, reconstruction_bottom_chunk])

#Create difference
diff = (chunked_volume - reconstruction_raw)

Specifications

• MATLAB/python version: Python
• OS: Ubunut
• CUDA version: 12.0

[AnderBiguri]

While this looks like numerical error, its extremely suspicious that it only happens in one side, so I need to investigate.