Half-fan cone-beam CT

[jbgeb]

Hi!

Thanks for maintaining this cool and helpful package!

I have a question concerning FBP for sparse-view offset-scan (half-fan) cbct (axial scan). In particular, I am looking at the case where I want to group my projections into a given amount of bins (these could e.g. be breathing phases in 4D phase-correlated FDK).

In the simplest case, I have two bins. The first n projections are assigned to bin A, the following n to bin B and then the next projection to A again, etc. I then want to reconstruct an image from all projections in bin A and a second image using all projections in bin B.

What I have now observed in leapct is that offset-scan FBP only works if a bin contains projections with corresponding angles around phi=0, so that angle[-1]-angle[0] > some_theshold (e.g. 359). If a bin does not happen to have projections around phi=0, the FBP fails with the error Error: offsetScan requires at least 360 degrees of projections!. I would expect the algorithm to be independent of the offset of the phis, as the number and distribution of projections and their angles are the same for both case mentioned above.

I can only solve this by shifting the angles in a bin to be around phi=0, reconstructing, and then rotating the final image back. This, however, does not result in the exact same image.

Is there another way to go about this with leapct?

Thanks in advance!

Here is also a code example based on your offset-scan demo. FBP works in the current form (of course with streaking artifacts) but fails when setting e.g. start_idx=0 in line 87.

import sys
import os
import time
import numpy as np
import matplotlib.pyplot as plt
from leapctype import *

leapct = tomographicModels()


# Specify the number of detector columns which is used below
# Scale the number of angles and the detector pixel size with N
numCols = 512
numAngles = 4 * int(360 * numCols / 1024)
pixelSize = 0.65 * 512 / numCols

# Set the number of detector rows
numRows = numCols

# Reduce the number of columns so that the phantom extends past either the
# right or left side of the projections (but not both)
# Note that we translate the centerCol parameter so that only one side is truncated
numCols = 300
centerCol = 0.5 * (numCols - 1) + 100

angles = leapct.setAngleArray(numAngles, 360.0)

# Set the scanner geometry
leapct.set_conebeam(
    numAngles,
    numRows,
    numCols,
    pixelSize,
    pixelSize,
    0.5 * (numRows - 1),
    centerCol,
    angles,
    1100,
    1400,
)

# Set the offsetScan flag to True
# This needs to be done after the geometry is set
# If you plan on setting the volume with the set_default_volume command, please set this flag first
# Setting this flag will automatically enlarge the circular field of view window
leapct.set_offsetScan(True)

# Set the volume parameters.
# It is best to do this after the CT geometry is set
leapct.set_default_volume()
# leapct.set_diameterFOV(leapct.get_voxelWidth()*leapct.get_numX())

# Trouble-Shooting Functions
leapct.print_parameters()
# leapct.sketch_system()

# Allocate space for the projections and the volume
g = leapct.allocateProjections()
f = leapct.allocateVolume()

# Specify simplified FORBILD head phantom
# One could easily do this in Python, but Python is soooooo slow for these types of operations,
# so we implemented this feature with multi-threaded C++
leapct.set_FORBILD(f, True)

# "Simulate" projection data
startTime = time.time()
leapct.project(g, f)
print("Forward Projection Elapsed Time: " + str(time.time() - startTime))

# Reconstruct the data
startTime = time.time()
leapct.FBP(g, f)
print("Reconstruction Elapsed Time: " + str(time.time() - startTime))


# Reconstruct the data for limited angles
num_phases = 2
num_cycles = 10
num_angles_per_cycle = numAngles // num_cycles
print(f"Number of angles per cycle: {num_angles_per_cycle}")


angles_phases = []
g_phases = []

start_idx = -int(num_angles_per_cycle / 4)
# start_idx = 0
for i in range(num_cycles):
    end_idx = start_idx + num_angles_per_cycle // num_phases
    print(f"Start index: {start_idx}, End index: {end_idx}")

    if start_idx < 0:
        angles_cycle = list(angles[start_idx:]) + list(angles[:end_idx])
        g_cycle = list(g[start_idx:]) + list(g[:end_idx])
    else:
        angles_cycle = list(angles[start_idx:end_idx])
        g_cycle = list(g[start_idx:end_idx])

    angles_phases += angles_cycle
    g_phases += g_cycle

    start_idx = start_idx + num_angles_per_cycle

# Sort angles_phases and g_phases based on angles
sorted_indices = np.argsort(angles_phases)
angles_phases = np.array(angles_phases)[sorted_indices]
g_phases = np.array(g_phases)[sorted_indices]

angles_phases = np.ascontiguousarray(angles_phases, dtype=np.float32)
g_phases = np.ascontiguousarray(g_phases, dtype=np.float32)
numAngles = g_phases.shape[0]


leapct.set_conebeam(
    numAngles,
    numRows,
    numCols,
    pixelSize,
    pixelSize,
    0.5 * (numRows - 1),
    centerCol,
    angles_phases,
    1100,
    1400,
)

f_phases = leapct.allocateVolume()

startTime = time.time()
leapct.FBP(g_phases, f_phases)
print("Reconstruction Elapsed Time: " + str(time.time() - startTime))


# Display the result with napari
leapct.display(f)
leapct.display(f_phases)

[kylechampley]

Thanks for including the code. I was able to reproduce your issue. I'll offer you a quick fix and a reason for this issue.

The quick fix is this: Take your first projection, add 360 to it and then put this projection at the end of the projection data. If you do this you will not have to rotate your reconstruction.

Now here is an explanation of why you are seeing this. First of all, half-fan FBP really does require projections over 360 degrees. Your projections have large gaps between sets of projections, about 18 degrees, while consecutive projections are only 0.5 degrees apart. Strickly speaking, FBP is not designed to handle such a large deviation of projection angle spacing. Such situations really should be handled by iterative reconstruction. Iterative reconstruction algorithms don't really care if the angular spacing is grossly unequal.

For your situation, LEAP sees that most of your angles are 0.5 degrees apart. So, when it looks at the difference of the last projection (which is at 341.5) and the first projection (at 0.0), it feels that this difference is much further from 360
than the average angular sampling distance, so it claims that you do not have a full revolution of angular spacing.

And thus when you map the projection at 0.0 degrees to 360, then the difference of the last and first projection is 360 - 0.5 = 359.5 which is close enough to 360 (within one mean angular sampling distance).

Hopefully this explanation makes sense. If not, feel free to respond.