Skip to content

Commit

Permalink
black linting + actions fix
Browse files Browse the repository at this point in the history
  • Loading branch information
@dkazanc
dkazanc committed Dec 21, 2023
1 parent 4920493 commit 99b559e
Show file tree
Hide file tree
Showing 17 changed files with 962 additions and 652 deletions.
134 changes: 70 additions & 64 deletions Demos/4D/Model3D_t.py
View file
Expand Up @@ -14,64 +14,67 @@
import tomophantom
from tomophantom import TomoP3D

print ("Building 4D phantom using TomoPhantom software")
tic=timeit.default_timer()
model = 100 # note that the selected model is temporal (3D + time)
print("Building 4D phantom using TomoPhantom software")
tic = timeit.default_timer()
model = 100 # note that the selected model is temporal (3D + time)
# Define phantom dimensions using a scalar (cubic) or a tuple [N1, N2, N3]
N_size = 256 # or as a tuple of a custom size (256,256,256)
N_size = 256 # or as a tuple of a custom size (256,256,256)
path = os.path.dirname(tomophantom.__file__)
path_library3D = os.path.join(path, "phantomlib", "Phantom3DLibrary.dat")
#This will generate a Time x N_size x N_size x N_size phantom (4D)
# This will generate a Time x N_size x N_size x N_size phantom (4D)
phantom_tm = TomoP3D.ModelTemporal(model, N_size, path_library3D)
toc=timeit.default_timer()
toc = timeit.default_timer()
Run_time = toc - tic
print("Phantom has been built in {} seconds".format(Run_time))

for i in range(0,np.size(phantom_tm,0)):
sliceSel = int(0.5*N_size)
#plt.gray()
plt.figure(1)
for i in range(0, np.size(phantom_tm, 0)):
sliceSel = int(0.5 * N_size)
# plt.gray()
plt.figure(1)
plt.subplot(131)
plt.imshow(phantom_tm[i,sliceSel,:,:],vmin=0, vmax=1)
plt.title('3D Phantom, axial view')
plt.imshow(phantom_tm[i, sliceSel, :, :], vmin=0, vmax=1)
plt.title("3D Phantom, axial view")

plt.subplot(132)
plt.imshow(phantom_tm[i,:,sliceSel,:],vmin=0, vmax=1)
plt.title('3D Phantom, coronal view')
plt.imshow(phantom_tm[i, :, sliceSel, :], vmin=0, vmax=1)
plt.title("3D Phantom, coronal view")

plt.subplot(133)
plt.imshow(phantom_tm[i,:,:,sliceSel],vmin=0, vmax=1)
plt.title('3D Phantom, sagittal view')
plt.imshow(phantom_tm[i, :, :, sliceSel], vmin=0, vmax=1)
plt.title("3D Phantom, sagittal view")
plt.show()
plt.pause(0.3)
#%%
print ("Getting 4D projection data using TomoPhantom software")
# %%
print("Getting 4D projection data using TomoPhantom software")
# Projection geometry related parameters:
Horiz_det = int(np.sqrt(2)*N_size) # detector column count (horizontal)
Vert_det = N_size # detector row count (vertical) (no reason for it to be > N)
angles_num = int(0.5*np.pi*N_size); # angles number
angles = np.linspace(0.0,179.9,angles_num,dtype='float32') # in degrees
angles_rad = angles*(np.pi/180.0)
Horiz_det = int(np.sqrt(2) * N_size) # detector column count (horizontal)
Vert_det = N_size # detector row count (vertical) (no reason for it to be > N)
angles_num = int(0.5 * np.pi * N_size)
# angles number
angles = np.linspace(0.0, 179.9, angles_num, dtype="float32") # in degrees
angles_rad = angles * (np.pi / 180.0)

projData4D_analyt= TomoP3D.ModelSinoTemporal(model, N_size, Horiz_det, Vert_det, angles, path_library3D)
projData4D_analyt = TomoP3D.ModelSinoTemporal(
model, N_size, Horiz_det, Vert_det, angles, path_library3D
)

time_frames = projData4D_analyt.shape[0]
intens_max = 60
sliceSel = 150
for i in range(0,time_frames):
plt.figure(2)
for i in range(0, time_frames):
plt.figure(2)
plt.subplot(131)
plt.imshow(projData4D_analyt[i,:,sliceSel,:],vmin=0, vmax=intens_max)
plt.title('2D Projection (analytical)')
plt.imshow(projData4D_analyt[i, :, sliceSel, :], vmin=0, vmax=intens_max)
plt.title("2D Projection (analytical)")
plt.subplot(132)
plt.imshow(projData4D_analyt[i,sliceSel,:,:],vmin=0, vmax=intens_max)
plt.title('Sinogram view')
plt.imshow(projData4D_analyt[i, sliceSel, :, :], vmin=0, vmax=intens_max)
plt.title("Sinogram view")
plt.subplot(133)
plt.imshow(projData4D_analyt[i,:,:,sliceSel],vmin=0, vmax=intens_max)
plt.title('Tangentogram view')
plt.imshow(projData4D_analyt[i, :, :, sliceSel], vmin=0, vmax=intens_max)
plt.title("Tangentogram view")
plt.show()
plt.pause(0.3)
#%%
# %%
# A capability of building a subset of vertical slices out of 4D phantom (faster)
import timeit
import os
Expand All @@ -80,61 +83,64 @@
import matplotlib.pyplot as plt
import numpy as np

print ("Building a subset of 4D phantom using TomoPhantom software")
tic=timeit.default_timer()
print("Building a subset of 4D phantom using TomoPhantom software")
tic = timeit.default_timer()
model = 101
N_size = 256 # Define phantom dimensions using a scalar value
DIM_z = (94, 158) # selected vertical subset (a slab) of the phantom
N_size = 256 # Define phantom dimensions using a scalar value
DIM_z = (94, 158) # selected vertical subset (a slab) of the phantom
path = os.path.dirname(tomophantom.__file__)
path_library3D = os.path.join(path, "phantomlib", "Phantom3DLibrary.dat")
phantom_tm = TomoP3D.ModelTemporalSub(model, N_size, DIM_z, path_library3D)
toc=timeit.default_timer()
toc = timeit.default_timer()
Run_time = toc - tic
print("Phantom has been built in {} seconds".format(Run_time))

for i in range(0,np.size(phantom_tm,0)):
for i in range(0, np.size(phantom_tm, 0)):
sliceSel = 32
#plt.gray()
plt.figure(1)
# plt.gray()
plt.figure(1)
plt.subplot(131)
plt.imshow(phantom_tm[i,sliceSel,:,:],vmin=0, vmax=1)
plt.title('4D Phantom, axial view')
plt.imshow(phantom_tm[i, sliceSel, :, :], vmin=0, vmax=1)
plt.title("4D Phantom, axial view")

plt.subplot(132)
plt.imshow(phantom_tm[i,:,70,:],vmin=0, vmax=1)
plt.title('4D Phantom, coronal view')
plt.imshow(phantom_tm[i, :, 70, :], vmin=0, vmax=1)
plt.title("4D Phantom, coronal view")

plt.subplot(133)
plt.imshow(phantom_tm[i,:,:,70],vmin=0, vmax=1)
plt.title('4D Phantom, sagittal view')
plt.imshow(phantom_tm[i, :, :, 70], vmin=0, vmax=1)
plt.title("4D Phantom, sagittal view")
plt.show()
plt.pause(0.5)


print ("Building a subset of 4D projection data using TomoPhantom software")
Horiz_det = int(np.sqrt(2)*N_size) # detector column count (horizontal)
Vert_det = N_size # detector row count (vertical) (no reason for it to be > N)
angles_num = int(0.5*np.pi*N_size); # angles number
angles = np.linspace(0.0,179.9,angles_num,dtype='float32') # in degrees

projData4D_cut = TomoP3D.ModelSinoTemporalSub(model, N_size, Horiz_det, Vert_det, DIM_z, angles, path_library3D)
print("Building a subset of 4D projection data using TomoPhantom software")
Horiz_det = int(np.sqrt(2) * N_size) # detector column count (horizontal)
Vert_det = N_size # detector row count (vertical) (no reason for it to be > N)
angles_num = int(0.5 * np.pi * N_size)
# angles number
angles = np.linspace(0.0, 179.9, angles_num, dtype="float32") # in degrees

projData4D_cut = TomoP3D.ModelSinoTemporalSub(
model, N_size, Horiz_det, Vert_det, DIM_z, angles, path_library3D
)

intens_max = 45
for i in range(0,np.size(projData4D_cut,0)):
for i in range(0, np.size(projData4D_cut, 0)):
sliceSel = 32
#plt.gray()
plt.figure(1)
# plt.gray()
plt.figure(1)
plt.subplot(131)
plt.imshow(projData4D_cut[i,sliceSel,:,:],vmin=0, vmax=intens_max)
plt.title('Sinogram View')
plt.imshow(projData4D_cut[i, sliceSel, :, :], vmin=0, vmax=intens_max)
plt.title("Sinogram View")

plt.subplot(132)
plt.imshow(projData4D_cut[i,:,sliceSel,:],vmin=0, vmax=intens_max)
plt.title('Projection view')
plt.imshow(projData4D_cut[i, :, sliceSel, :], vmin=0, vmax=intens_max)
plt.title("Projection view")

plt.subplot(133)
plt.imshow(projData4D_cut[i,:,:,sliceSel],vmin=0, vmax=intens_max)
plt.title('Tangentogram view')
plt.imshow(projData4D_cut[i, :, :, sliceSel], vmin=0, vmax=intens_max)
plt.title("Tangentogram view")
plt.show()
plt.pause(0.5)
#%%
# %%
98 changes: 63 additions & 35 deletions Demos/RandomPhantoms/RandPhantGen.py
View file
Expand Up @@ -13,12 +13,12 @@
import numpy as np
import matplotlib.pyplot as plt

from tomophantom import TomoP2D
from tomophantom import TomoP3D
from tomophantom.generator import foam2D,foam3D
from tomophantom import TomoP2D
from tomophantom import TomoP3D
from tomophantom.generator import foam2D, foam3D

N_size = 256 # define the grid size
tot_objects = 300 # the total number of objects to generate
N_size = 256 # define the grid size
tot_objects = 300 # the total number of objects to generate

# define ranges for parameters
x0min = -0.9
Expand All @@ -33,57 +33,85 @@
ab_max = 0.25

# 2D example
(Objfoam2D,myObjects) = foam2D(x0min, x0max, y0min, y0max, c0min, c0max, ab_min, ab_max, N_size, tot_objects, object_type = 'mix')
(Objfoam2D, myObjects) = foam2D(
x0min,
x0max,
y0min,
y0max,
c0min,
c0max,
ab_min,
ab_max,
N_size,
tot_objects,
object_type="mix",
)
plt.figure()
plt.imshow(Objfoam2D, vmin=0, vmax=0.3, cmap="gray")
plt.imshow(Objfoam2D, vmin=0, vmax=0.3, cmap="gray")

# Generate a sinogram
angles_num = int(0.5*np.pi*N_size); # angles number
angles = np.linspace(0,180,angles_num,dtype='float32')
angles_rad = angles*(np.pi/180)
P = int(np.sqrt(2)*N_size) #detectors
angles_num = int(0.5 * np.pi * N_size)
# angles number
angles = np.linspace(0, 180, angles_num, dtype="float32")
angles_rad = angles * (np.pi / 180)
P = int(np.sqrt(2) * N_size) # detectors

sino_Objfoam2D = TomoP2D.ObjectSino(N_size, P, angles, myObjects)
plt.figure()
plt.imshow(sino_Objfoam2D, cmap="gray")
plt.imshow(sino_Objfoam2D, cmap="gray")

#%%
# %%
# 3D example
print ("Generating 3D random phantom")
(Objfoam3D, myObjects) = foam3D(x0min, x0max, y0min, y0max, z0min, z0max, c0min, c0max, ab_min, ab_max, N_size, tot_objects, object_type = 'mix')
print("Generating 3D random phantom")
(Objfoam3D, myObjects) = foam3D(
x0min,
x0max,
y0min,
y0max,
z0min,
z0max,
c0min,
c0max,
ab_min,
ab_max,
N_size,
tot_objects,
object_type="mix",
)
plt.figure()
sliceSel = int(0.5*N_size)
sliceSel = int(0.5 * N_size)
plt.subplot(131)
plt.imshow(Objfoam3D[sliceSel,:,:],vmin=0, vmax=0.5)
plt.title('3D Object, axial view')
plt.imshow(Objfoam3D[sliceSel, :, :], vmin=0, vmax=0.5)
plt.title("3D Object, axial view")

plt.subplot(132)
plt.imshow(Objfoam3D[:,sliceSel,:],vmin=0, vmax=0.5)
plt.title('3D Object, coronal view')
plt.imshow(Objfoam3D[:, sliceSel, :], vmin=0, vmax=0.5)
plt.title("3D Object, coronal view")

plt.subplot(133)
plt.imshow(Objfoam3D[:,:,sliceSel],vmin=0, vmax=0.5)
plt.title('3D Object, sagittal view')
plt.imshow(Objfoam3D[:, :, sliceSel], vmin=0, vmax=0.5)
plt.title("3D Object, sagittal view")
plt.show()

Horiz_det = int(np.sqrt(2)*N_size) # detector column count (horizontal)
Vert_det = N_size # detector row count (vertical) (no reason for it to be > N)
angles_num = int(0.5*np.pi*N_size); # angles number
angles = np.linspace(0.0,179.9,angles_num,dtype='float32') # in degrees
Horiz_det = int(np.sqrt(2) * N_size) # detector column count (horizontal)
Vert_det = N_size # detector row count (vertical) (no reason for it to be > N)
angles_num = int(0.5 * np.pi * N_size)
# angles number
angles = np.linspace(0.0, 179.9, angles_num, dtype="float32") # in degrees

print ("Building 3D analytical projection data with TomoPhantom")
print("Building 3D analytical projection data with TomoPhantom")
ProjData3D = TomoP3D.ObjectSino(N_size, Horiz_det, Vert_det, angles, myObjects)

sliceSel = 150
plt.figure()
plt.figure()
plt.subplot(131)
plt.imshow(ProjData3D[:,sliceSel,:])
plt.title('2D Projection (analytical)')
plt.imshow(ProjData3D[:, sliceSel, :])
plt.title("2D Projection (analytical)")
plt.subplot(132)
plt.imshow(ProjData3D[sliceSel,:,:])
plt.title('Sinogram view')
plt.imshow(ProjData3D[sliceSel, :, :])
plt.title("Sinogram view")
plt.subplot(133)
plt.imshow(ProjData3D[:,:,sliceSel])
plt.title('Tangentogram view')
plt.imshow(ProjData3D[:, :, sliceSel])
plt.title("Tangentogram view")
plt.show()
#%%
# %%

0 comments on commit 99b559e

Please sign in to comment.