Skip to content

Commit

Permalink
added modular-beam FBP
Browse files Browse the repository at this point in the history
  • Loading branch information
@kylechampley
kylechampley committed Dec 29, 2023
1 parent c57d219 commit f6540a9
Show file tree
Hide file tree
Showing 15 changed files with 1,097 additions and 89 deletions.
23 changes: 17 additions & 6 deletions demo_leapctype/test_modularbeam.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,13 +7,13 @@

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

# Set the number of detector rows
numRows = numCols
numRows = 1
#numRows = 1

# Make this modular-beam geometry just like a cone-beam dataset
# so let's define sod and sdd when defining our geometry
Expand Down Expand Up @@ -49,6 +49,7 @@

# Set the volume parameters
leapct.set_default_volume()
leapct.set_diameterFOV(leapct.get_numX()*leapct.get_voxelWidth())

# Trouble-Shooting Functions
leapct.print_parameters()
Expand All @@ -65,7 +66,7 @@
leapct.set_FORBILD(f,True)
#leapct.display(f)

leapct.set_gpu(0)
#leapct.set_gpu(0)

# "Simulate" projection data
startTime = time.time()
Expand All @@ -85,8 +86,9 @@
# Reconstruct the data
startTime = time.time()
#leapct.backproject(g,f)
leapct.FBP(g,f)
#leapct.ASDPOCS(g,f,50,10,1,1.0/20.0)
leapct.SART(g,f,20,10)
#leapct.SART(g,f,20,10)
#leapct.OSEM(g,f,10,10)
#leapct.LS(g,f,10,True)
print('Reconstruction Elapsed Time: ' + str(time.time()-startTime))
Expand All @@ -103,4 +105,13 @@
leapct.display(f)



''' Compare to cone-beam
ct2 = tomographicModels()
ct2.set_conebeam(numAngles, numRows, numCols, pixelSize, pixelSize, 0.5*(numRows-1), 0.5*(numCols-1), leapct.setAngleArray(numAngles, 360.0), 1100, 1400)
ct2.set_default_volume()
f_cone = ct2.allocateVolume()
ct2.FBP(g,f_cone)
leapct.display(f_cone)
leapct.display(f_cone-f)
quit()
#'''
11 changes: 7 additions & 4 deletions src/filtered_backprojection.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,9 +106,9 @@ bool filteredBackprojection::convolve1D(float* g, parameters* params, bool data_

bool filteredBackprojection::filterProjections(float* g, parameters* params, bool data_on_cpu)
{
if (params->geometry == parameters::MODULAR)
if (params->geometry == parameters::MODULAR && params->modularbeamIsAxiallyAligned() == false)
{
printf("Error: not implemented for modular geometries\n");
printf("Error: projection filtering only implemented for modular geometries whose rowVectors are aligned with the z-axis\n");
return false;
}

Expand Down Expand Up @@ -165,8 +165,11 @@ bool filteredBackprojection::execute(float* g, float* f, parameters* params, boo
{
if (params->geometry == parameters::MODULAR)
{
printf("Error: FBP not implemented for modular geometries\n");
return false;
if (params->modularbeamIsAxiallyAligned() == false || params->min_T_phi() < 1.0e-16)
{
printf("Error: FBP only implemented for modular geometries whose rowVectors are aligned with the z-axis\n");
return false;
}
}
if (params->geometry == parameters::CONE && params->helicalPitch != 0.0 && params->whichGPU < 0)
{
Expand Down
85 changes: 71 additions & 14 deletions src/leapctype.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,17 +75,61 @@ def __init__(self, param_id=None, lib_dir=""):
if _platform == "linux" or _platform == "linux2":
import readline
from ctypes import cdll
self.libprojectors = cdll.LoadLibrary(os.path.join(current_dir, "../build/lib/libleap.so"))

fullPath = os.path.join(current_dir, 'libleap.so')
fullPath_backup = os.path.join(current_dir, '../build/lib/libleap.so')

if os.path.isfile(fullPath):
self.libprojectors = cdll.LoadLibrary(fullPath)
elif os.path.isfile(fullPath_backup):
self.libprojectors = cdll.LoadLibrary(fullPath_backup)
else:
print('Error: could not find LEAP dynamic library at')
print(fullPath)
print('or')
print(fullPath_backup)
self.libprojectors = None

elif _platform == "win32":
from ctypes import windll
try:
self.libprojectors = windll.LoadLibrary(os.path.join(current_dir, r'..\win_build\bin\Release\libleap.dll'))
except:
self.libprojectors = ctypes.CDLL(os.path.join(current_dir, r'..\win_build\bin\Release\libleap.dll'), winmode=0)

fullPath = os.path.join(current_dir, 'libleap.dll')
fullPath_backup = os.path.join(current_dir, r'..\win_build\bin\Release\libleap.dll')

if os.path.isfile(fullPath):
try:
self.libprojectors = windll.LoadLibrary(fullPath)
except:
self.libprojectors = ctypes.CDLL(fullPath, winmode=0)
elif os.path.isfile(fullPath_backup):
try:
self.libprojectors = windll.LoadLibrary(fullPath_backup)
except:
self.libprojectors = ctypes.CDLL(fullPath_backup, winmode=0)
else:
print('Error: could not find LEAP dynamic library at')
print(fullPath)
print('or')
print(fullPath_backup)
self.libprojectors = None

elif _platform == "darwin": # Darwin is the name for MacOS in Python's platform module
from ctypes import cdll
# there is current no support for LEAP on Mac, but maybe someone can figure this out
self.libprojectors = cdll.LoadLibrary(os.path.join(current_dir, "../build/lib/libleap.dylib"))
from ctypes import cdll

fullPath = os.path.join(current_dir, 'libleap.dylib')
fullPath_backup = os.path.join(current_dir, '../build/lib/libleap.dylib')

if os.path.isfile(fullPath):
self.libprojectors = cdll.LoadLibrary(fullPath)
elif os.path.isfile(fullPath_backup):
self.libprojectors = cdll.LoadLibrary(fullPath_backup)
else:
print('Error: could not find LEAP dynamic library at')
print(fullPath)
print('or')
print(fullPath_backup)
self.libprojectors = None

if param_id is not None:
self.param_id = param_id
Expand Down Expand Up @@ -918,7 +962,7 @@ def get_FBPscalar(self):
self.set_model()
return self.libprojectors.get_FBPscalar()

def FBP(self, g, f, inplace=False):
def FBP(self, g, f=None, inplace=False):
"""Performs a Filtered Backprojection (FBP) reconstruction of the projection data, g, and stores the result in f
The CT geometry parameters and the CT volume parameters must be set prior to running this function.
Expand All @@ -942,9 +986,13 @@ def FBP(self, g, f, inplace=False):
self.libprojectors.FBP.restype = ctypes.c_bool
self.set_model()
if has_torch == True and type(q) is torch.Tensor:
if f is None:
f = self.allocateVolume(0.0,True)
self.libprojectors.FBP.argtypes = [ctypes.c_void_p, ctypes.c_void_p, ctypes.c_bool]
self.libprojectors.FBP(q.data_ptr(), f.data_ptr(), q.is_cuda == False)
else:
if f is None:
f = self.allocateVolume()
self.libprojectors.FBP.argtypes = [ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ctypes.c_bool]
self.libprojectors.FBP(q, f, True)
return f
Expand Down Expand Up @@ -2367,7 +2415,11 @@ def sketch_system(self,whichView=None):

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
self.drawCT(ax,whichView)
if whichView is None or len(whichView) == 1:
self.drawCT(ax,whichView)
else:
for i in range(len(whichView)):
self.drawCT(ax,whichView[i])
self.drawVolume(ax)

ax.set_xlabel('X (mm)')
Expand Down Expand Up @@ -2877,6 +2929,12 @@ def loadData(self, fileName):
try:
from PIL import Image
import glob
hasPIL = True
except:
print('Error: Failed to load PIL or glob library!')
print('To install PIL do: pip install Pillow')
return None
if hasPIL == True:
currentWorkingDirectory = os.getcwd()
dataFolder, baseFileName = os.path.split(fileName)
if len(dataFolder) > 0:
Expand All @@ -2894,18 +2952,17 @@ def loadData(self, fileName):
justDigits.append(int(digitStr))
ind = np.argsort(justDigits)

#print('found ' + str(len(fileList)) + ' images')
#print('reading first image: ' + str(fileList[0]))
firstImg = np.array(Image.open(fileList[0]))
x = np.zeros((len(fileList), firstImg.shape[0], firstImg.shape[1]), dtype=np.float32)
print('found ' + str(x.shape[0]) + ' images of size ' + str(x.shape[1]) + ' x ' + str(x.shape[2]))
for i in range(len(fileList)):
anImage = np.array(Image.open(fileList[ind[i]]))
#anImage = np.array(Image.open(fileList[ind[i]]).rotate(-0.5))
x[i,:,:] = anImage[:,:]
os.chdir(currentWorkingDirectory)
return x

except:
print('Error: Failed to load PIL or glob library!')
print('To install PIL do: pip install Pillow')
return None
'''
try:
from PIL import Image
Expand Down
Loading

0 comments on commit f6540a9

Please sign in to comment.