Support for plotting helical trajectory with curved detector?

[ProjX0]

Can the d09_curved_polygon_detector_array.py script plot a helical trajectory along with the curved detector using matplotlib (plt)?

[ProjX0]

Below is the modified d09_curved_polygon_detector_array.py where I have changed the geometry variables. I would like to visualize the geometry by plotting it along with a helical trajectory in the CT geometry of the "true" system.

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

# First we define the geometry
numCols_per_module = 80
numSticks = 48
numCols = numSticks*numCols_per_module
numRows = 128
numAngles = 720
pixelSize = 0.15 # mm

# Set source-to-object and source-to-detector distances
sod = 330.0
sdd = sod + 280.0

leapct_regular = tomographicModels()

# Define the angular rotation samples
phis = leapct_regular.setAngleArray(numAngles, 360.0)

# Set the regular geometry as a curved cone-beam system
leapct_regular.set_conebeam(numAngles, numRows, numCols, pixelSize, pixelSize, 0.5*(numRows-1), 0.5*(numCols-1), phis, sod, sdd)
leapct_regular.set_curvedDetector()
leapct_regular.set_normalizedHelicalPitch(0.5)
leapct_regular.set_volume(numCols, numCols, numRows, pixelSize*sod/sdd, pixelSize*sod/sdd)
# leapct_regular.set_volume(100, 100, 100, 0.17, 0.17)

g = leapct_regular.allocate_projections()


# Make another tomographicModels object to track the true geometry
leapct_true = tomographicModels()

# Define the angular spacing between detector modules
#angularSpacing = np.arctan(pixelSize*numCols_per_module / sdd)
angularSpacing = np.arctan((pixelSize*numCols_per_module + 1.0) / sdd)

# Define the angles of each module
alphas = (np.array(range(numSticks))-(numSticks-1)/2.0)*angularSpacing

# Now we define the CT geometry of the "true" system
leapct_true.set_volume(numCols, numCols, numRows, pixelSize*sod/sdd, pixelSize*sod/sdd)
f = leapct_true.allocateVolume() # shape is numZ, numY, numX
leapct_true.set_FORBILD(f,True)

for n in range(numAngles):
    phi = phis[n]*np.pi/180.0 - np.pi/2.0
    sourcePositions = np.zeros((numSticks,3),dtype=np.float32)
    sourcePositions[:,0] = sod*np.cos(phi)
    sourcePositions[:,1] = sod*np.sin(phi)
    moduleCenters = np.zeros((numSticks,3),dtype=np.float32)
    colVectors = np.zeros((numSticks,3), dtype=np.float32)
    rowVectors = np.zeros((numSticks,3), dtype=np.float32)
    rowVectors[:,2] = 1.0
    for m in range(numSticks):
        colVectors[m,0] = -np.sin(phi-alphas[m])
        colVectors[m,1] = np.cos(phi-alphas[m])
        moduleCenters[m,0] = -sdd*np.cos(phi-alphas[m])+sourcePositions[m,0]
        moduleCenters[m,1] = -sdd*np.sin(phi-alphas[m])+sourcePositions[m,1]
    
    leapct_true.set_modularbeam(numSticks, numRows, numCols_per_module, pixelSize, pixelSize, sourcePositions, moduleCenters, rowVectors, colVectors)
    leapct_true.set_normalizedHelicalPitch(0.5)
    g_sim = leapct_true.allocateProjections() # shape is numAngles, numRows, numCols
    

    leapct_true.sketch_system()
    
    # "Simulate" the data
    leapct_true.project(g_sim,f)

    # Concatenate the projections onto each modular into one projection
    for m in range(numSticks):
        g[n,:,m*numCols_per_module:(m+1)*numCols_per_module] = g_sim[m,:,:]


# To enable the rebinning, we must specify the angle (in degrees) of each detector pixel along a single row
fanAngles = np.array(range(numCols),dtype=np.float32)
fanAngles_per_module = (np.array(range(numCols_per_module))-(numCols_per_module-1)/2.0)*np.arctan(pixelSize/sdd)
for m in range(numSticks):
    fanAngles[m*numCols_per_module:(m+1)*numCols_per_module] = alphas[m]+fanAngles_per_module
fanAngles *= 180.0/np.pi

# Now perform the rebinning.  The "order" argument specifies the order of the interpolating polynomial.
startTime = time.time()
leapct_regular.rebin_curved(g, fanAngles, order=6)
leapct_regular.sketch_system()
print('Rebinning Elapsed Time: ' + str(time.time()-startTime) + ' seconds')

# Finally, perform FBP reconstruction and display the result
f = leapct_regular.FBP(g)
f[f<0.0] = 0.0
plt.imshow(np.squeeze(f[f.shape[0]//2,:,:]), cmap=plt.get_cmap('gray'))
plt.show()

[kylechampley]

Yes, I agree that the sketch_system function needs improvements. I am interested in improving it, but not sure when I will have time to get to this.

[ProjX0]

I understand. I'll be looking forward to your improvements whenever you have the time to work on them! Thank you

[ProjX0]

I modified the internal source code for sketch_system(). So, it's not perfect like the picture below, but I made a temporary representation of the geometry with helical trajectory and curved detector.

[kylechampley]

Looks good!