diff --git a/Demos/4D/Model3D_t.py b/Demos/4D/Model3D_t.py index 3b185f5..401d817 100644 --- a/Demos/4D/Model3D_t.py +++ b/Demos/4D/Model3D_t.py @@ -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 @@ -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) -#%% \ No newline at end of file +# %% diff --git a/Demos/RandomPhantoms/RandPhantGen.py b/Demos/RandomPhantoms/RandPhantGen.py index 6f99cb4..0f8ff64 100644 --- a/Demos/RandomPhantoms/RandPhantGen.py +++ b/Demos/RandomPhantoms/RandPhantGen.py @@ -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 @@ -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() -#%% \ No newline at end of file +# %% diff --git a/Demos/RandomPhantoms/RandPhantGen_crystal.py b/Demos/RandomPhantoms/RandPhantGen_crystal.py index b853a06..c174fc8 100755 --- a/Demos/RandomPhantoms/RandPhantGen_crystal.py +++ b/Demos/RandomPhantoms/RandPhantGen_crystal.py @@ -6,7 +6,7 @@ import numpy as np import matplotlib.pyplot as plt -from tomophantom import TomoP2D +from tomophantom import TomoP2D from tomophantom.TomoP2D import Objects2D from scipy.ndimage import gaussian_filter from tomobar.methodsDIR import RecToolsDIR @@ -15,7 +15,7 @@ # create a 2D object explicitly without using parameters file -N_size = 512 # define the grid +N_size = 512 # define the grid # A PHANTOM WITHOUT ARTEFACTS a_el1_min = 0.7 @@ -25,13 +25,15 @@ b_el1_max = 0.75 b_el1 = random.uniform(b_el1_min, b_el1_max) -el1 = {'Obj': Objects2D.ELLIPSE, - 'C0' : 0.7, - 'x0' : 0.0, - 'y0' : 0.0, - 'a' : a_el1, - 'b' : b_el1, - 'phi': 0.0} +el1 = { + "Obj": Objects2D.ELLIPSE, + "C0": 0.7, + "x0": 0.0, + "y0": 0.0, + "a": a_el1, + "b": b_el1, + "phi": 0.0, +} a_el2_min = 0.6 a_el2_max = a_el1 @@ -40,13 +42,15 @@ b_el2_max = b_el1 b_el2 = random.uniform(b_el2_min, b_el2_max) -el2 = {'Obj': Objects2D.ELLIPSE, - 'C0' : -0.4, - 'x0' : 0.0, - 'y0' : 0.0, - 'a' : a_el2, - 'b' : b_el2, - 'phi': 0.0} +el2 = { + "Obj": Objects2D.ELLIPSE, + "C0": -0.4, + "x0": 0.0, + "y0": 0.0, + "a": a_el2, + "b": b_el2, + "phi": 0.0, +} C0_min = 0.01 C0_max = 0.2 @@ -61,56 +65,64 @@ phi_max = 90.0 phi = random.uniform(phi_min, phi_max) -el3 = {'Obj': Objects2D.RECTANGLE, - 'C0' : C_0, - 'x0' : 0.0, - 'y0' : 0.0, - 'a' : a_el3, - 'b' : b_el3, - 'phi': phi} +el3 = { + "Obj": Objects2D.RECTANGLE, + "C0": C_0, + "x0": 0.0, + "y0": 0.0, + "a": a_el3, + "b": b_el3, + "phi": phi, +} -GROUND_TRUTH = TomoP2D.Object(N_size, [el1,el2,el3]) +GROUND_TRUTH = TomoP2D.Object(N_size, [el1, el2, el3]) -plt.close('all') +plt.close("all") plt.figure(1) -plt.rcParams.update({'font.size': 21}) +plt.rcParams.update({"font.size": 21}) plt.imshow(GROUND_TRUTH, vmin=0, vmax=1, cmap="BuPu") -plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical') -plt.title('{}'.format('Ground Truth Object')) -#%% +plt.colorbar(ticks=[0, 0.5, 1], orientation="vertical") +plt.title("{}".format("Ground Truth Object")) +# %% # define all objects bellow: -el1 = {'Obj': Objects2D.ELLIPSE, - 'C0' : 0.7, - 'x0' : 0.0, - 'y0' : 0.0, - 'a' : a_el1, - 'b' : b_el1, - 'phi': 0.0} - -el2 = {'Obj': Objects2D.ELLIPSE, - 'C0' : 0.7, - 'x0' : 0.0, - 'y0' : 0.0, - 'a' : a_el1-0.015, - 'b' : b_el1-0.015, - 'phi': 0.0} +el1 = { + "Obj": Objects2D.ELLIPSE, + "C0": 0.7, + "x0": 0.0, + "y0": 0.0, + "a": a_el1, + "b": b_el1, + "phi": 0.0, +} + +el2 = { + "Obj": Objects2D.ELLIPSE, + "C0": 0.7, + "x0": 0.0, + "y0": 0.0, + "a": a_el1 - 0.015, + "b": b_el1 - 0.015, + "phi": 0.0, +} Object1 = TomoP2D.Object(N_size, [el1]) Object2 = TomoP2D.Object(N_size, [el2]) Object_loop = Object1 - Object2 Object_loop = gaussian_filter(Object_loop, sigma=3) -Object_loop = gaussian_filter(Object1+Object_loop, sigma=3) +Object_loop = gaussian_filter(Object1 + Object_loop, sigma=3) -el3 = {'Obj': Objects2D.ELLIPSE, - 'C0' : -0.4, - 'x0' : 0.0, - 'y0' : 0.0, - 'a' : a_el2, - 'b' : b_el2, - 'phi': 0.0} +el3 = { + "Obj": Objects2D.ELLIPSE, + "C0": -0.4, + "x0": 0.0, + "y0": 0.0, + "a": a_el2, + "b": b_el2, + "phi": 0.0, +} Object3 = TomoP2D.Object(N_size, [el3]) @@ -118,25 +130,29 @@ Object = Object_loop + Object_liquor -el6 = {'Obj': Objects2D.RECTANGLE, - 'C0' : 2.0*C_0, - 'x0' : 0.0, - 'y0' : 0.0, - 'a' : a_el3+0.02, - 'b' : b_el3+0.02, - 'phi': phi} +el6 = { + "Obj": Objects2D.RECTANGLE, + "C0": 2.0 * C_0, + "x0": 0.0, + "y0": 0.0, + "a": a_el3 + 0.02, + "b": b_el3 + 0.02, + "phi": phi, +} Object6 = TomoP2D.Object(N_size, [el6]) Object_crystal_edge = gaussian_filter(Object6, sigma=3) Object += Object_crystal_edge -el5 = {'Obj': Objects2D.RECTANGLE, - 'C0' : C_0, - 'x0' : 0.0, - 'y0' : 0.0, - 'a' : a_el3, - 'b' : b_el3, - 'phi': phi} +el5 = { + "Obj": Objects2D.RECTANGLE, + "C0": C_0, + "x0": 0.0, + "y0": 0.0, + "a": a_el3, + "b": b_el3, + "phi": phi, +} Object5 = TomoP2D.Object(N_size, [el5]) Object_crystal = gaussian_filter(Object5, sigma=3) @@ -144,106 +160,118 @@ Object -= Object_crystal # forming dictionaries with artifact types -_noise_ = {'type' : 'Gaussian', - 'sigma' : 0.1, # noise amplitude - 'seed' : None} +_noise_ = {"type": "Gaussian", "sigma": 0.1, "seed": None} # noise amplitude # adding zingers and stripes -_zingers_ = {'percentage' : 0.5, - 'modulus' : 50} +_zingers_ = {"percentage": 0.5, "modulus": 50} -_stripes_ = {'percentage' : 1.2, - 'maxthickness' : 3.0, - 'intensity' : 0.25, - 'type' : 'partial', - 'variability' : 0.005} +_stripes_ = { + "percentage": 1.2, + "maxthickness": 3.0, + "intensity": 0.25, + "type": "partial", + "variability": 0.005, +} -Object = artefacts_mix(Object, _noise_, _zingers_, _stripes_, _sinoshifts_= {}) +Object = artefacts_mix(Object, _noise_, _zingers_, _stripes_, _sinoshifts_={}) Object = gaussian_filter(Object, sigma=3) plt.figure(2) -plt.rcParams.update({'font.size': 21}) +plt.rcParams.update({"font.size": 21}) plt.imshow(Object, vmin=0, vmax=1, cmap="BuPu") -plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical') -plt.title('{}'.format('Distorted Phantom')) - -#%% -# Generate projection data of distorted phantom - -angles_num = int(np.pi*N_size); # angles number -angles = np.linspace(0.0,179.9,angles_num,dtype='float32') -angles_rad = angles*(np.pi/180.0) -P = N_size #detectors - -Rectools = RecToolsDIR(DetectorsDimH = P, # Horizontal detector dimension - DetectorsDimV = None, # Vertical detector dimension (3D case) - CenterRotOffset = 0.0, # Center of Rotation scalar - AnglesVec = angles_rad, # A vector of projection angles in radians - ObjSize = N_size, # Reconstructed object dimensions (scalar) - device_projector='gpu') +plt.colorbar(ticks=[0, 0.5, 1], orientation="vertical") +plt.title("{}".format("Distorted Phantom")) + +# %% +# Generate projection data of distorted phantom + +angles_num = int(np.pi * N_size) +# angles number +angles = np.linspace(0.0, 179.9, angles_num, dtype="float32") +angles_rad = angles * (np.pi / 180.0) +P = N_size # detectors + +Rectools = RecToolsDIR( + DetectorsDimH=P, # Horizontal detector dimension + DetectorsDimV=None, # Vertical detector dimension (3D case) + CenterRotOffset=0.0, # Center of Rotation scalar + AnglesVec=angles_rad, # A vector of projection angles in radians + ObjSize=N_size, # Reconstructed object dimensions (scalar) + device_projector="gpu", +) sino_num = Rectools.FORWPROJ(Object) -_noise_={} -_zingers_={} -_sinoshifts_={} +_noise_ = {} +_zingers_ = {} +_sinoshifts_ = {} -_stripes_ = {'percentage' : 0.75, - 'maxthickness' : 2.0, - 'intensity' : 0.15, - 'type' : 'mix', - 'variability' : 0.005} +_stripes_ = { + "percentage": 0.75, + "maxthickness": 2.0, + "intensity": 0.15, + "type": "mix", + "variability": 0.005, +} -sino_num_artifacts = artefacts_mix(sino_num, _noise_, _zingers_, _stripes_, _sinoshifts_) +sino_num_artifacts = artefacts_mix( + sino_num, _noise_, _zingers_, _stripes_, _sinoshifts_ +) plt.figure(3) -plt.rcParams.update({'font.size': 21}) +plt.rcParams.update({"font.size": 21}) plt.imshow(sino_num_artifacts, cmap="BuPu") -plt.title('{}'.format('Distorted Phantom')) -#%% -FBPrec = Rectools.FBP(sino_num_artifacts) # perform FBP reconstruction +plt.title("{}".format("Distorted Phantom")) +# %% +FBPrec = Rectools.FBP(sino_num_artifacts) # perform FBP reconstruction plt.figure() -plt.rcParams.update({'font.size': 20}) +plt.rcParams.update({"font.size": 20}) plt.imshow(FBPrec, vmin=0, vmax=1, cmap="gray") -plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical') -plt.title('FBP reconstruction') -#%% +plt.colorbar(ticks=[0, 0.5, 1], orientation="vertical") +plt.title("FBP reconstruction") +# %% from tomobar.methodsIR import RecToolsIR -RectoolsIR = RecToolsIR(DetectorsDimH = P, # DetectorsDimH # detector dimension (horizontal) - DetectorsDimV = None, # DetectorsDimV # detector dimension (vertical) for 3D case only - CenterRotOffset = None, # Center of Rotation (CoR) scalar (for 3D case only) - AnglesVec = angles_rad, # array of angles in radians - ObjSize = N_size, # a scalar to define reconstructed object dimensions - datafidelity='LS',# data fidelity, choose LS, PWLS (wip), GH (wip), Student (wip) - device_projector='gpu') - -_data_ = {'projection_norm_data' : sino_num_artifacts} # data dictionary -lc = RectoolsIR.powermethod(_data_) # calculate Lipschitz constant (run once to initialise) -_algorithm_ = {'iterations' : 350, - 'lipschitz_const' : lc} -# adding regularisation using the CCPi regularisation toolkit -_regularisation_ = {'method' : 'PD_TV', - 'regul_param' : 0.001, - 'iterations' : 150, - 'device_regulariser': 'gpu'} +RectoolsIR = RecToolsIR( + DetectorsDimH=P, # DetectorsDimH # detector dimension (horizontal) + DetectorsDimV=None, # DetectorsDimV # detector dimension (vertical) for 3D case only + CenterRotOffset=None, # Center of Rotation (CoR) scalar (for 3D case only) + AnglesVec=angles_rad, # array of angles in radians + ObjSize=N_size, # a scalar to define reconstructed object dimensions + datafidelity="LS", # data fidelity, choose LS, PWLS (wip), GH (wip), Student (wip) + device_projector="gpu", +) + +_data_ = {"projection_norm_data": sino_num_artifacts} # data dictionary +lc = RectoolsIR.powermethod( + _data_ +) # calculate Lipschitz constant (run once to initialise) +_algorithm_ = {"iterations": 350, "lipschitz_const": lc} -# Run FISTA reconstrucion algorithm with regularisation +# adding regularisation using the CCPi regularisation toolkit +_regularisation_ = { + "method": "PD_TV", + "regul_param": 0.001, + "iterations": 150, + "device_regulariser": "gpu", +} + +# Run FISTA reconstrucion algorithm with regularisation RecFISTA_reg = RectoolsIR.FISTA(_data_, _algorithm_, _regularisation_) plt.figure() plt.imshow(RecFISTA_reg, vmin=0, vmax=1, cmap="gray") -plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical') -plt.title('TV-Regularised FISTA reconstruction') +plt.colorbar(ticks=[0, 0.5, 1], orientation="vertical") +plt.title("TV-Regularised FISTA reconstruction") plt.show() -#%% +# %% # Segment the result of FBPrec while using the GROUND_TRUTH as ideally segmented phantom -#%% +# %% ############################################################################## ####################3D phantom generator###################################### ############################################################################## @@ -270,15 +298,17 @@ c_el1_max = 0.85 c_el1 = random.uniform(c_el1_min, c_el1_max) -el1 = {'Obj': Objects3D.ELLIPSOID, - 'C0' : 0.7, - 'x0' : 0.0, - 'y0' : 0.0, - 'z0' : 0.0, - 'a' : a_el1, - 'b' : b_el1, - 'c' : c_el1, - 'phi1': 0.0} +el1 = { + "Obj": Objects3D.ELLIPSOID, + "C0": 0.7, + "x0": 0.0, + "y0": 0.0, + "z0": 0.0, + "a": a_el1, + "b": b_el1, + "c": c_el1, + "phi1": 0.0, +} a_el2_min = 0.6 @@ -291,15 +321,17 @@ c_el2_max = c_el1 c_el2 = random.uniform(c_el2_min, c_el2_max) -el2 = {'Obj': Objects3D.ELLIPSOID, - 'C0' : -0.4, - 'x0' : 0.0, - 'y0' : 0.0, - 'z0' : 0.0, - 'a' : a_el2, - 'b' : b_el2, - 'c' : c_el2, - 'phi1' : 0.0} +el2 = { + "Obj": Objects3D.ELLIPSOID, + "C0": -0.4, + "x0": 0.0, + "y0": 0.0, + "z0": 0.0, + "a": a_el2, + "b": b_el2, + "c": c_el2, + "phi1": 0.0, +} C0_min = 0.01 C0_max = 0.2 @@ -320,34 +352,36 @@ phi_max = 180.0 phi1 = random.uniform(phi_min, phi_max) -el3 = {'Obj': Objects3D.CUBOID, - 'C0' : C_0, - 'x0' : x0_rand, - 'y0' : y0_rand, - 'z0' : z0_rand, - 'a' : a_el3, - 'b' : b_el3, - 'c' : c_el3, - 'phi1': phi1} +el3 = { + "Obj": Objects3D.CUBOID, + "C0": C_0, + "x0": x0_rand, + "y0": y0_rand, + "z0": z0_rand, + "a": a_el3, + "b": b_el3, + "c": c_el3, + "phi1": phi1, +} -GROUND_TRUTH = TomoP3D.Object(N_size, [el1,el2,el3]) +GROUND_TRUTH = TomoP3D.Object(N_size, [el1, el2, el3]) GROUND_TRUTH[GROUND_TRUTH > 0.7] = 0.43336788 GROUND_TRUTH[(GROUND_TRUTH > 0.0) & (GROUND_TRUTH < 0.29999998)] = 0.43336788 -sliceSel = (int)(N_size/2) -plt.figure() +sliceSel = (int)(N_size / 2) +plt.figure() plt.subplot(131) -plt.imshow(GROUND_TRUTH[sliceSel,:,:]) -plt.title('Ideal Phantom1') +plt.imshow(GROUND_TRUTH[sliceSel, :, :]) +plt.title("Ideal Phantom1") plt.subplot(132) -plt.imshow(GROUND_TRUTH[:,sliceSel,:]) -plt.title('Ideal Phantom2') +plt.imshow(GROUND_TRUTH[:, sliceSel, :]) +plt.title("Ideal Phantom2") plt.subplot(133) -plt.imshow(GROUND_TRUTH[:,:,sliceSel]) -plt.title('Ideal Phantom3') +plt.imshow(GROUND_TRUTH[:, :, sliceSel]) +plt.title("Ideal Phantom3") plt.show() -#%% +# %% """ print("Generating artificial phantom") el1_add = {'Obj': Objects3D.ELLIPSOID, @@ -390,123 +424,142 @@ Object -= Object_crystal """ # forming dictionaries with artifact types -_noise_ = {'noise_type' : 'Gaussian', - 'noise_amplitude' : 0.02, # noise amplitude - 'noise_seed' : None} +_noise_ = { + "noise_type": "Gaussian", + "noise_amplitude": 0.02, # noise amplitude + "noise_seed": None, +} # adding zingers and stripes -_zingers_ = {'zingers_percentage' : 1.5, - 'zingers_modulus' : 50} +_zingers_ = {"zingers_percentage": 1.5, "zingers_modulus": 50} -_stripes_ = {'stripes_percentage' : 1.2, - 'stripes_maxthickness' : 3.0, - 'stripes_intensity' : 0.25, - 'stripes_type' : 'partial', - 'stripes_variability' : 0.005} +_stripes_ = { + "stripes_percentage": 1.2, + "stripes_maxthickness": 3.0, + "stripes_intensity": 0.25, + "stripes_type": "partial", + "stripes_variability": 0.005, +} Object = artefacts_mix(GROUND_TRUTH, **_noise_, **_zingers_, **_stripes_) - -sliceSel = (int)(N_size/2) -plt.figure() +sliceSel = (int)(N_size / 2) +plt.figure() plt.subplot(131) -plt.imshow(Object[:,sliceSel,:]) -plt.title('Distorted Phantom1') +plt.imshow(Object[:, sliceSel, :]) +plt.title("Distorted Phantom1") plt.subplot(132) -plt.imshow(Object[sliceSel,:,:]) -plt.title('Distorted Phantom2') +plt.imshow(Object[sliceSel, :, :]) +plt.title("Distorted Phantom2") plt.subplot(133) -plt.imshow(Object[:,:,sliceSel]) -plt.title('Distorted Phantom3') +plt.imshow(Object[:, :, sliceSel]) +plt.title("Distorted Phantom3") plt.show() -#%% -print ("Simulate synthetic flat fields, add flat field background to the projections and add noise") +# %% +print( + "Simulate synthetic flat fields, add flat field background to the projections and add noise" +) from tomophantom.flatsgen import synth_flats from tomophantom.artefacts import artefacts_mix -I0 = 75000; # Source intensity -flatsnum = 20 # the number of the flat fields required - -angles_num = int(np.pi*N_size); # angles number -angles = np.linspace(0.0,179.9,angles_num,dtype='float32') -angles_rad = angles*(np.pi/180.0) -P = N_size #detectors - -Rectools = RecToolsDIR(DetectorsDimH = P, # Horizontal detector dimension - DetectorsDimV = N_size, # Vertical detector dimension (3D case) - CenterRotOffset = 0.0, # Center of Rotation scalar - AnglesVec = angles_rad, # A vector of projection angles in radians - ObjSize = N_size, # Reconstructed object dimensions (scalar) - device_projector='gpu') +I0 = 75000 +# Source intensity +flatsnum = 20 # the number of the flat fields required + +angles_num = int(np.pi * N_size) +# angles number +angles = np.linspace(0.0, 179.9, angles_num, dtype="float32") +angles_rad = angles * (np.pi / 180.0) +P = N_size # detectors + +Rectools = RecToolsDIR( + DetectorsDimH=P, # Horizontal detector dimension + DetectorsDimV=N_size, # Vertical detector dimension (3D case) + CenterRotOffset=0.0, # Center of Rotation scalar + AnglesVec=angles_rad, # A vector of projection angles in radians + ObjSize=N_size, # Reconstructed object dimensions (scalar) + device_projector="gpu", +) projection_data3D = Rectools.FORWPROJ(Object) intens_max_clean = np.max(projection_data3D) -_fresnel_propagator_ = {'fresnel_dist_observation' : 40, - 'fresnel_scale_factor' : 10, - 'fresnel_wavelenght' : 0.007} +_fresnel_propagator_ = { + "fresnel_dist_observation": 40, + "fresnel_scale_factor": 10, + "fresnel_wavelenght": 0.007, +} projection_data3D_fresnel = artefacts_mix(projection_data3D, **_fresnel_propagator_) -#%% -[projData3D_noisy, flatsSIM] = synth_flats(projection_data3D_fresnel, - source_intensity = I0, source_variation=0.01,\ - arguments_Bessel = (1,10,10,12),\ - specklesize = 15,\ - kbar = 0.3, - jitter = 0.1, - sigmasmooth = 3, flatsnum=flatsnum) -#del projData3D_analyt - -print ("Normalise projections using ToMoBAR software") +# %% +[projData3D_noisy, flatsSIM] = synth_flats( + projection_data3D_fresnel, + source_intensity=I0, + source_variation=0.01, + arguments_Bessel=(1, 10, 10, 12), + specklesize=15, + kbar=0.3, + jitter=0.1, + sigmasmooth=3, + flatsnum=flatsnum, +) +# del projData3D_analyt + +print("Normalise projections using ToMoBAR software") from tomobar.supp.suppTools import normaliser # normalise the data, the required format is [detectorsX, Projections, detectorsY] -projData3D_norm = normaliser(projData3D_noisy, flatsSIM, darks=None, log='true', method='mean') +projData3D_norm = normaliser( + projData3D_noisy, flatsSIM, darks=None, log="true", method="mean" +) -#del projData3D_noisy -intens_max = 0.3*np.max(projData3D_norm) +# del projData3D_noisy +intens_max = 0.3 * np.max(projData3D_norm) sliceSel = 150 -plt.figure() +plt.figure() plt.subplot(131) -plt.imshow(projData3D_norm[:,sliceSel,:],vmin=0, vmax=intens_max) -plt.title('Normalised 2D Projection (erroneous)') +plt.imshow(projData3D_norm[:, sliceSel, :], vmin=0, vmax=intens_max) +plt.title("Normalised 2D Projection (erroneous)") plt.subplot(132) -plt.imshow(projData3D_norm[sliceSel,:,:],vmin=0, vmax=intens_max) -plt.title('Sinogram view') +plt.imshow(projData3D_norm[sliceSel, :, :], vmin=0, vmax=intens_max) +plt.title("Sinogram view") plt.subplot(133) -plt.imshow(projData3D_norm[:,:,sliceSel],vmin=0, vmax=intens_max) -plt.title('Tangentogram view') +plt.imshow(projData3D_norm[:, :, sliceSel], vmin=0, vmax=intens_max) +plt.title("Tangentogram view") plt.show() -#%% +# %% from tomobar.methodsDIR import RecToolsDIR -RectoolsDIR = RecToolsDIR(DetectorsDimH = P, # DetectorsDimH # detector dimension (horizontal) - DetectorsDimV = N_size, # DetectorsDimV # detector dimension (vertical) for 3D case only - CenterRotOffset = None, # Center of Rotation (CoR) scalar (for 3D case only) - AnglesVec = angles_rad, # array of angles in radians - ObjSize = N_size, # a scalar to define reconstructed object dimensions - device_projector = 'gpu') - -print ("Reconstruction using FBP from tomobar") -recNumerical= RectoolsDIR.FBP(projData3D_norm) # FBP reconstruction + +RectoolsDIR = RecToolsDIR( + DetectorsDimH=P, # DetectorsDimH # detector dimension (horizontal) + DetectorsDimV=N_size, # DetectorsDimV # detector dimension (vertical) for 3D case only + CenterRotOffset=None, # Center of Rotation (CoR) scalar (for 3D case only) + AnglesVec=angles_rad, # array of angles in radians + ObjSize=N_size, # a scalar to define reconstructed object dimensions + device_projector="gpu", +) + +print("Reconstruction using FBP from tomobar") +recNumerical = RectoolsDIR.FBP(projData3D_norm) # FBP reconstruction recNumerical *= intens_max_clean -sliceSel = int(0.5*N_size) +sliceSel = int(0.5 * N_size) max_val = 1 -#plt.gray() -plt.figure() +# plt.gray() +plt.figure() plt.subplot(131) -plt.imshow(recNumerical[sliceSel,:,:],vmin=0, vmax=max_val) -plt.title('3D Reconstruction, axial view') +plt.imshow(recNumerical[sliceSel, :, :], vmin=0, vmax=max_val) +plt.title("3D Reconstruction, axial view") plt.subplot(132) -plt.imshow(recNumerical[:,sliceSel,:],vmin=0, vmax=max_val) -plt.title('3D Reconstruction, coronal view') +plt.imshow(recNumerical[:, sliceSel, :], vmin=0, vmax=max_val) +plt.title("3D Reconstruction, coronal view") plt.subplot(133) -plt.imshow(recNumerical[:,:,sliceSel],vmin=0, vmax=max_val) -plt.title('3D Reconstruction, sagittal view') +plt.imshow(recNumerical[:, :, sliceSel], vmin=0, vmax=max_val) +plt.title("3D Reconstruction, sagittal view") plt.show() -#%% \ No newline at end of file +# %% diff --git a/docs/source/conf.py b/docs/source/conf.py index b31d98e..edd2329 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -94,7 +94,7 @@ html_theme_options = { "logo": { "image_dark": "_static/tomophantom_logo_dark.png", - }, + }, "show_toc_level": 1, "navbar_align": "left", # [left, content, right] For testing that the navbar items align properly } diff --git a/tests/conftest.py b/tests/conftest.py index f6797bd..4d969cb 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -5,9 +5,11 @@ import pytest import tomophantom + @pytest.fixture(scope="session") def sino_model11_64(): from tomophantom.TomoP2D import ModelSino + N_size = 64 # set the desired dimension of the phantom model = 11 angles_num = int(0.5 * np.pi * N_size) @@ -24,9 +26,11 @@ def sino_model11_64(): return (sino_an, angles_num, P) + @pytest.fixture(scope="session") def sino_model11_3d_64(): from tomophantom.TomoP3D import ModelSino + N_size = 64 # set the desired dimension of the phantom # Projection geometry related parameters: Horiz_det = int(np.sqrt(2) * N_size) # detector column count (horizontal) @@ -40,8 +44,6 @@ def sino_model11_3d_64(): path_library3D = os.path.join(path, "phantomlib", "Phantom3DLibrary.dat") # Generatea 3d sino - sino3d = ModelSino( - 11, N_size, Horiz_det, Vert_det, angles, path_library3D - ) + sino3d = ModelSino(11, N_size, Horiz_det, Vert_det, angles, path_library3D) - return (sino3d, Vert_det, angles_num, Horiz_det) \ No newline at end of file + return (sino3d, Vert_det, angles_num, Horiz_det) diff --git a/tests/test_flats_gen.py b/tests/test_flats_gen.py index 83c40aa..ce5cb3d 100644 --- a/tests/test_flats_gen.py +++ b/tests/test_flats_gen.py @@ -2,6 +2,7 @@ from tomophantom.flatsgen import synth_flats + def test_3d_models_sino(sino_model11_3d_64): sino3d, Vert_det, angles_num, Horiz_det = sino_model11_3d_64 @@ -10,20 +11,20 @@ def test_3d_models_sino(sino_model11_3d_64): flatsnum = 20 # the number of the flat fields required [projData3D_raw, flats_combined3D, speckel_map] = synth_flats( - sino3d, - source_intensity=I0, - detectors_miscallibration=0.05, - variations_number=3, - arguments_Bessel=(1, 25), - specklesize=2, - kbar=2, - jitter_projections=0.0, - sigmasmooth=3, - flatsnum=flatsnum, + sino3d, + source_intensity=I0, + detectors_miscallibration=0.05, + variations_number=3, + arguments_Bessel=(1, 25), + specklesize=2, + kbar=2, + jitter_projections=0.0, + sigmasmooth=3, + flatsnum=flatsnum, ) assert 0.0 <= np.max(projData3D_raw) <= 65535 assert projData3D_raw.dtype == np.uint16 assert projData3D_raw.shape == (Vert_det, angles_num, Horiz_det) assert flats_combined3D.dtype == np.uint16 - assert flats_combined3D.shape == (Vert_det, flatsnum, Horiz_det) \ No newline at end of file + assert flats_combined3D.shape == (Vert_det, flatsnum, Horiz_det) diff --git a/tests/test_generators.py b/tests/test_generators.py index aac0858..836189d 100644 --- a/tests/test_generators.py +++ b/tests/test_generators.py @@ -1,14 +1,14 @@ import numpy as np from numpy.testing import assert_allclose -from tomophantom.generator import foam2D,foam3D +from tomophantom.generator import foam2D, foam3D eps = 1e-05 def test_2d_generator(): N_size = 64 # set the desired dimension of the phantom - tot_objects = 20 # the total number of objects to generate + tot_objects = 20 # the total number of objects to generate # define ranges for parameters x0min = -0.9 @@ -21,7 +21,19 @@ def test_2d_generator(): 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", + ) assert 0.0 <= np.max(Objfoam2D) <= 10 assert Objfoam2D.dtype == np.float32 @@ -29,8 +41,8 @@ def test_2d_generator(): def test_3d_generator(): - N_size = 64 # define the grid size - tot_objects = 10 # the total number of objects to generate + N_size = 64 # define the grid size + tot_objects = 10 # the total number of objects to generate # define ranges for parameters x0min = -0.9 @@ -44,7 +56,21 @@ def test_3d_generator(): ab_min = 0.01 ab_max = 0.25 - (Objfoam3D, myObjects) = foam3D(x0min, x0max, y0min, y0max, z0min, z0max, c0min, c0max, ab_min, ab_max, N_size, tot_objects, object_type = 'mix') + (Objfoam3D, myObjects) = foam3D( + x0min, + x0max, + y0min, + y0max, + z0min, + z0max, + c0min, + c0max, + ab_min, + ab_max, + N_size, + tot_objects, + object_type="mix", + ) assert 0.0 <= np.max(Objfoam3D) <= 10 assert Objfoam3D.dtype == np.float32 diff --git a/tests/test_models_sino3d.py b/tests/test_models_sino3d.py index 3ba04b6..5d81c73 100644 --- a/tests/test_models_sino3d.py +++ b/tests/test_models_sino3d.py @@ -4,7 +4,13 @@ import pytest import os import tomophantom -from tomophantom.TomoP3D import _check_params3d, ModelSino, ModelSinoSub, ModelSinoTemporal, ModelSinoTemporalSub +from tomophantom.TomoP3D import ( + _check_params3d, + ModelSino, + ModelSinoSub, + ModelSinoTemporal, + ModelSinoTemporalSub, +) eps = 1e-05 @@ -24,13 +30,12 @@ def test_3d_models_sino(model): path_library3D = os.path.join(path, "phantomlib", "Phantom3DLibrary.dat") # Generatea 3d sino - sino3d = ModelSino( - model, N_size, Horiz_det, Vert_det, angles, path_library3D - ) + sino3d = ModelSino(model, N_size, Horiz_det, Vert_det, angles, path_library3D) assert 0.0 <= np.max(sino3d) <= 1000 assert sino3d.dtype == np.float32 assert sino3d.shape == (Vert_det, angles_num, Horiz_det) + def test_3d_models_sino_sub(): N_size = 64 # set the desired dimension of the phantom # Projection geometry related parameters: @@ -43,16 +48,17 @@ def test_3d_models_sino_sub(): # one can specify an exact path to the parameters file path = os.path.dirname(tomophantom.__file__) path_library3D = os.path.join(path, "phantomlib", "Phantom3DLibrary.dat") - subset_tuple = (15,30) + subset_tuple = (15, 30) # Generate a 3d sino sino3d = ModelSinoSub( - 11, N_size, Horiz_det, Vert_det, subset_tuple, angles, path_library3D + 11, N_size, Horiz_det, Vert_det, subset_tuple, angles, path_library3D ) assert 0.0 <= np.max(sino3d) <= 1000 assert sino3d.dtype == np.float32 assert sino3d.shape == (15, angles_num, Horiz_det) + @pytest.mark.parametrize("model", list(np.arange(100, 103, 1))) def test_3d_models_sino_temporal(model): N_size = 64 # set the desired dimension of the phantom @@ -66,12 +72,12 @@ def test_3d_models_sino_temporal(model): # one can specify an exact path to the parameters file path = os.path.dirname(tomophantom.__file__) path_library3D = os.path.join(path, "phantomlib", "Phantom3DLibrary.dat") - + params = _check_params3d(model, path_library3D) # Generatea 4d sino sino4d = ModelSinoTemporal( - model, N_size, Horiz_det, Vert_det, angles, path_library3D + model, N_size, Horiz_det, Vert_det, angles, path_library3D ) assert 0.0 <= np.max(sino4d) <= 1000 assert sino4d.dtype == np.float32 @@ -87,17 +93,17 @@ def test_3d_models_sino_temporal_sub(): angles_num = int(0.5 * np.pi * N_size) # angles number angles = np.linspace(0.0, 179.9, angles_num, dtype="float32") # in degrees - + # one can specify an exact path to the parameters file path = os.path.dirname(tomophantom.__file__) path_library3D = os.path.join(path, "phantomlib", "Phantom3DLibrary.dat") params = _check_params3d(model, path_library3D) - subset_tuple = (15,30) + subset_tuple = (15, 30) # Generate a 3d sino sino4d = ModelSinoTemporalSub( - model, N_size, Horiz_det, Vert_det, subset_tuple, angles, path_library3D + model, N_size, Horiz_det, Vert_det, subset_tuple, angles, path_library3D ) assert 0.0 <= np.max(sino4d) <= 1000 assert sino4d.dtype == np.float32 - assert sino4d.shape == (params[3], 15, angles_num, Horiz_det) \ No newline at end of file + assert sino4d.shape == (params[3], 15, angles_num, Horiz_det) diff --git a/tests/test_object_sino3d.py b/tests/test_object_sino3d.py index bbd78f2..7421700 100644 --- a/tests/test_object_sino3d.py +++ b/tests/test_object_sino3d.py @@ -3,6 +3,7 @@ from tomophantom.TomoP3D import Objects3D, ObjectSino + def test_3d_object_sino(): N3D_size = 64 # set the desired dimension of the phantom Horiz_det = int(np.sqrt(2) * N3D_size) # detector column count (horizontal) diff --git a/tests/test_sino2d_artefacts.py b/tests/test_sino2d_artefacts.py index 98d5b69..7d676d2 100644 --- a/tests/test_sino2d_artefacts.py +++ b/tests/test_sino2d_artefacts.py @@ -6,6 +6,7 @@ eps = 1e-05 + @pytest.mark.parametrize("noise_type", ["Poisson", "Gaussian"]) def test_2d_sino_noise(noise_type, sino_model11_64): sino_an, angles_num, P = sino_model11_64 diff --git a/tomophantom/TomoP2D.py b/tomophantom/TomoP2D.py index 61ea8b7..9ad53e2 100644 --- a/tomophantom/TomoP2D.py +++ b/tomophantom/TomoP2D.py @@ -31,7 +31,7 @@ "Model", "ModelTemporal", "ModelSino", - "ModelSinoTemporal", + "ModelSinoTemporal", "Object", "ObjectSino", "SinoNum", @@ -138,12 +138,13 @@ def ModelTemporal( ) return phantom + def ModelSino( model_no: int, phantom_size: int, detector_size: int, angles: np.ndarray, - models_library_path: Path, + models_library_path: Path, ) -> np.ndarray: """Generate 2D analytical sinogram for corresponding models in the library file Phantom2DLibrary.dat. @@ -182,6 +183,7 @@ def ModelSino( ) return sino2d.transpose() + def ModelSinoTemporal( model_no: int, phantom_size: int, @@ -228,6 +230,7 @@ def ModelSinoTemporal( ) return sino2d_t + def Object(phantom_size: int, obj_params: Union[list, dict]) -> np.ndarray: """Generate a 2D analytical phantom for the standalone geometrical object that is parametrised in the "obj_params" dictionary. @@ -267,8 +270,12 @@ def Object(phantom_size: int, obj_params: Union[list, dict]) -> np.ndarray: ) return object2d + def ObjectSino( - phantom_size: int, detector_size: int, angles: np.ndarray, obj_params: Union[list, dict] + phantom_size: int, + detector_size: int, + angles: np.ndarray, + obj_params: Union[list, dict], ) -> np.ndarray: """Generate a 2D analytical sinogram for the standalone geometrical object that is parametrised in the "obj_params" dictionary. @@ -311,6 +318,7 @@ def ObjectSino( ) return sino2d.transpose() + def SinoNum(input: np.ndarray, detector_X: int, angles: np.ndarray) -> np.ndarray: """Numerical calculation of 2D sinogram from the 2D input. @@ -342,6 +350,7 @@ def SinoNum(input: np.ndarray, detector_X: int, angles: np.ndarray) -> np.ndarra ) return sinogram.transpose() + def __testParams(obj): """Performs a simple type check of the input parameters and a range check""" if not type(obj) is dict: diff --git a/tomophantom/TomoP3D.py b/tomophantom/TomoP3D.py index 44556fa..fb0ed4e 100644 --- a/tomophantom/TomoP3D.py +++ b/tomophantom/TomoP3D.py @@ -34,13 +34,13 @@ "Model", "ModelSub", "ModelTemporal", - "ModelTemporalSub", + "ModelTemporalSub", "ModelSino", "ModelSinoSub", - 'ModelSinoTemporal', - 'ModelSinoTemporalSub', + "ModelSinoTemporal", + "ModelSinoTemporalSub", "Object", - 'ObjectSino', + "ObjectSino", ] @@ -178,6 +178,7 @@ def ModelSub( ) return phantom + def ModelTemporal( model_no: int, phantom_size: Union[int, Tuple[int, int, int]], @@ -283,6 +284,7 @@ def ModelTemporalSub( ) return phantom + def ModelSino( model_no: int, phantom_size: int, @@ -306,8 +308,10 @@ def ModelSino( np.ndarray: The generated 3D projection data of the following shape [detector_vert, AngTot, detector_horiz]. """ if type(phantom_size) == tuple: - raise ValueError('Please give a scalar for the phantom size (2d slice), if one needs non-cubic data generator please see ModelSinoSub function') - + raise ValueError( + "Please give a scalar for the phantom size (2d slice), if one needs non-cubic data generator please see ModelSinoSub function" + ) + # check the validity of model's parameters params = np.ascontiguousarray(np.zeros([12], dtype=ctypes.c_int)) params = _check_params3d(model_no, models_library_path) @@ -335,7 +339,8 @@ def ModelSino( angles_total, models_library_path, ) - return np.swapaxes(sino3d,0,1) + return np.swapaxes(sino3d, 0, 1) + def ModelSinoSub( model_no: int, @@ -362,7 +367,9 @@ def ModelSinoSub( np.ndarray: The generated 3D projection data of the following shape [detector_vert, AngTot, detector_horiz]. """ if type(phantom_size) == tuple: - raise ValueError('Please give a scalar for phantom size, projection data cannot be obtained for non-cubic phantom') + raise ValueError( + "Please give a scalar for phantom size, projection data cannot be obtained for non-cubic phantom" + ) Z1, Z2 = [int(i) for i in sub_index] @@ -377,7 +384,7 @@ def ModelSinoSub( raise ValueError("Range of the higher index is incorrect") sub_size = Z2 - Z1 # the size of the vertical slab - + # check the validity of model's parameters params = np.ascontiguousarray(np.zeros([12], dtype=ctypes.c_int)) params = _check_params3d(model_no, models_library_path) @@ -405,7 +412,8 @@ def ModelSinoSub( angles_total, models_library_path, ) - return np.swapaxes(sino3d,0,1) + return np.swapaxes(sino3d, 0, 1) + def ModelSinoTemporal( model_no: int, @@ -431,7 +439,9 @@ def ModelSinoTemporal( """ if type(phantom_size) == tuple: - raise ValueError('Please give a scalar for phantom size, projection data cannot be obtained for non-cubic phantom') + raise ValueError( + "Please give a scalar for phantom size, projection data cannot be obtained for non-cubic phantom" + ) # check the validity of model's parameters params = np.ascontiguousarray(np.zeros([12], dtype=ctypes.c_int)) @@ -441,13 +451,15 @@ def ModelSinoTemporal( angles_total = len(angles) if params[3] == 1: raise ValueError( - "The selected model is stationary (3D), use 'ModelSino' function instead" - ) + "The selected model is stationary (3D), use 'ModelSino' function instead" + ) else: sino4d = np.zeros( - [params[3], angles_total, detector_vert, detector_horiz], dtype="float32", order="C" + [params[3], angles_total, detector_vert, detector_horiz], + dtype="float32", + order="C", ) - + external.c_model_sino3d( np.ascontiguousarray(sino4d), model_no, @@ -460,7 +472,8 @@ def ModelSinoTemporal( angles_total, models_library_path, ) - return np.swapaxes(sino4d,1,2) + return np.swapaxes(sino4d, 1, 2) + def ModelSinoTemporalSub( model_no: int, @@ -487,7 +500,9 @@ def ModelSinoTemporalSub( np.ndarray: The generated 3D projection data of the following shape [time, detector_vert, AngTot, detector_horiz]. """ if type(phantom_size) == tuple: - raise ValueError('Please give a scalar for phantom size, projection data cannot be obtained for non-cubic phantom') + raise ValueError( + "Please give a scalar for phantom size, projection data cannot be obtained for non-cubic phantom" + ) Z1, Z2 = [int(i) for i in sub_index] @@ -502,7 +517,7 @@ def ModelSinoTemporalSub( raise ValueError("Range of the higher index is incorrect") sub_size = Z2 - Z1 # the size of the vertical slab - + # check the validity of model's parameters params = np.ascontiguousarray(np.zeros([12], dtype=ctypes.c_int)) params = _check_params3d(model_no, models_library_path) @@ -512,11 +527,13 @@ def ModelSinoTemporalSub( if params[3] == 1: raise ValueError( "The selected model is stationary (3D), use 'ModelSino' function instead" - ) + ) else: sino4d = np.zeros( - [params[3], angles_total, sub_size, detector_horiz], dtype="float32", order="C" - ) + [params[3], angles_total, sub_size, detector_horiz], + dtype="float32", + order="C", + ) external.c_model_sino3d( np.ascontiguousarray(sino4d), @@ -530,7 +547,7 @@ def ModelSinoTemporalSub( angles_total, models_library_path, ) - return np.swapaxes(sino4d,1,2) + return np.swapaxes(sino4d, 1, 2) def Object( @@ -599,7 +616,7 @@ def ObjectSino( detector_horiz: int, detector_vert: int, angles: np.ndarray, - obj_params: Union[list, dict] + obj_params: Union[list, dict], ) -> np.ndarray: """Generates a 3D analytical projection data for the standalone geometrical figure that is parametrised in the "obj_params" dictionary. @@ -607,22 +624,26 @@ def ObjectSino( Args: phantom_size (int): A scalar for phantom dimensions. detector_horiz (int): Size of the horizontal detector in pixels. - detector_vert (int): Size of the vertical detector in pixels. - angles (np.ndarray): Angles vector in degrees. + detector_vert (int): Size of the vertical detector in pixels. + angles (np.ndarray): Angles vector in degrees. obj_params (a list of dicts or dict): A dictionary with parameters of an object, see demos. Returns: np.ndarray: The generated 3D projection data for an object. """ if type(phantom_size) == tuple: - raise ValueError('Please give a scalar for phantom size, projection data cannot be obtained for non-cubic phantom') + raise ValueError( + "Please give a scalar for phantom size, projection data cannot be obtained for non-cubic phantom" + ) if type(obj_params) is dict: obj_params = [obj_params] - + angles_total = len(angles) tt = 0 - sino3d = np.zeros([angles_total, detector_vert, detector_horiz], dtype="float32", order="C") + sino3d = np.zeros( + [angles_total, detector_vert, detector_horiz], dtype="float32", order="C" + ) # unpacking obj_params dictionary for obj in obj_params: if __testParams3d(obj): @@ -660,7 +681,7 @@ def ObjectSino( phi3, tt, ) - return np.swapaxes(sino3d,0,1) + return np.swapaxes(sino3d, 0, 1) def __testParams3d(obj): diff --git a/tomophantom/artefacts.py b/tomophantom/artefacts.py index 7a9e4d2..29b1eea 100644 --- a/tomophantom/artefacts.py +++ b/tomophantom/artefacts.py @@ -8,16 +8,17 @@ import random from typing import Union, Any + def artefacts_mix(data: np.ndarray, **artefacts_dict: Any) -> Union[np.ndarray, list]: - """A module to generate and apply a mix of various typical imaging artefacts to add to the simulated data. - One can build various dictionaries with keywords arguments specified + """A module to generate and apply a mix of various typical imaging artefacts to add to the simulated data. + One can build various dictionaries with keywords arguments specified bellow and pass it to the function as: `artefacts_mix(data, **noise_dict, **zingers, **etc)`. - - DISCLAIMER: Note that most of the features are experimental and do not always reflect the - accurate modelling of the inaccuracies. + + DISCLAIMER: Note that most of the features are experimental and do not always reflect the + accurate modelling of the inaccuracies. Args: - data (np.ndarray): 2D or 3D numpy array (sinogram or 3D projection data). + data (np.ndarray): 2D or 3D numpy array (sinogram or 3D projection data). The input data must be of the following shape: 2D sinogram (anglesDim, DetectorsHoriz), 3D projection data (DetectorsVert, anglesDim, DetectorsHoriz). **artefacts_dict (dict): A dictionary with keyword arguments related to different artefact type. @@ -43,9 +44,9 @@ def artefacts_mix(data: np.ndarray, **artefacts_dict: Any) -> Union[np.ndarray, verbose (bool): Make the output of modules verbose. Returns: - np.ndarray: 2D/3D numpy array with artefacts applied to input data. - [np.ndarray, shifts]: a list of 2D/3D numpy array and simulated shifts. - """ + np.ndarray: 2D/3D numpy array with artefacts applied to input data. + [np.ndarray, shifts]: a list of 2D/3D numpy array and simulated shifts. + """ ####### VERBOSE ######## if "verbose" not in artefacts_dict: verbose = True @@ -205,12 +206,12 @@ def artefacts_mix(data: np.ndarray, **artefacts_dict: Any) -> Union[np.ndarray, # NOISE if _noise_["noise_type"] is not None: sino_artifacts = noise( - data=sino_artifacts, - sigma=_noise_["noise_amplitude"], - noisetype=_noise_["noise_type"], - seed=_noise_["noise_seed"], - prelog=_noise_["noise_prelog"], - ) + data=sino_artifacts, + sigma=_noise_["noise_amplitude"], + noisetype=_noise_["noise_type"], + seed=_noise_["noise_seed"], + prelog=_noise_["noise_prelog"], + ) if verbose is True: print("{} noise has been added to the data.".format(_noise_["noise_type"])) @@ -222,7 +223,14 @@ def artefacts_mix(data: np.ndarray, **artefacts_dict: Any) -> Union[np.ndarray, return np.float32(sino_artifacts) -def stripes(data: np.ndarray, percentage: float, maxthickness: int, intensity_thresh: float, stripe_type: str, variability: float) -> np.ndarray: +def stripes( + data: np.ndarray, + percentage: float, + maxthickness: int, + intensity_thresh: float, + stripe_type: str, + variability: float, +) -> np.ndarray: """Function to add stripes (constant offsets) to sinograms or 3D projection data which results in rings in the reconstructed image. @@ -235,7 +243,7 @@ def stripes(data: np.ndarray, percentage: float, maxthickness: int, intensity_th variability (float): Variability multiplier to incorporate change of intensity in the stripe. Raises: - ValueError: Percentage is out of range. + ValueError: Percentage is out of range. ValueError: Thickness is out of range. Returns: @@ -329,7 +337,7 @@ def stripes(data: np.ndarray, percentage: float, maxthickness: int, intensity_th def zingers(data: np.ndarray, percentage: float, modulus: int) -> np.ndarray: """Adding zingers (zero single pixels or small 4 pixels clusters) - to sinograms or 6 voxels to 3D projection data. + to sinograms or 6 voxels to 3D projection data. Args: data (np.ndarray): 2D sinogram or 3D projection data. The input data must be of the following shape: 2D sinogram (anglesDim, DetectorsHoriz), @@ -338,7 +346,7 @@ def zingers(data: np.ndarray, percentage: float, modulus: int) -> np.ndarray: modulus (int): Modulus to control the amount of 4/6 pixel clusters to be added. Raises: - ValueError: Percentage is out of range. + ValueError: Percentage is out of range. ValueError: Modulus must be positive. Returns: @@ -384,7 +392,13 @@ def zingers(data: np.ndarray, percentage: float, modulus: int) -> np.ndarray: return sino_zingers -def noise(data: np.ndarray, sigma: Union[int, float], noisetype: str, seed: bool=True, prelog: bool=False) -> Union[list, np.ndarray]: +def noise( + data: np.ndarray, + sigma: Union[int, float], + noisetype: str, + seed: bool = True, + prelog: bool = False, +) -> Union[list, np.ndarray]: """Adding random noise to data (adapted from LD-CT simulator) Args: @@ -440,7 +454,7 @@ def datashifts(data: np.ndarray, maxamplitude: int) -> list: maxamplitude (int): The misilighnment ammplitude. Defines the maximal amplitude of each shift in pixels. Returns: - list: 2D or 3d data with misalignment and shifts vectors [data, shifts]. + list: 2D or 3d data with misalignment and shifts vectors [data, shifts]. """ if data.ndim == 2: (anglesDim, DetectorsDimH) = np.shape(data) @@ -491,9 +505,10 @@ def datashifts_subpixel(data: np.ndarray, maxamplitude: float) -> list: maxamplitude (float): The misilighnment ammplitude. Defines the maximal amplitude of each shift. Returns: - list: 2D or 3d data with misalignment and shifts vectors [data, shifts]. + list: 2D or 3d data with misalignment and shifts vectors [data, shifts]. """ from skimage import transform as tf + if data.ndim == 2: shifts = np.zeros([1, 2], dtype="float32") random_shift_x = random.uniform( @@ -531,7 +546,7 @@ def datashifts_subpixel(data: np.ndarray, maxamplitude: float) -> list: def pve(data: np.ndarray, pve_strength: int) -> np.ndarray: - """Applying Partial Volume effect (smoothing) to data. + """Applying Partial Volume effect (smoothing) to data. Args: data (np.ndarray): 2D sinogram or 3D projection data. The input data must be of the following shape: 2D sinogram (anglesDim, DetectorsHoriz), @@ -563,8 +578,10 @@ def pve(data: np.ndarray, pve_strength: int) -> np.ndarray: return data_pve -def fresnel_propagator(data:np.ndarray, dist_observation: int, scale_factor: float, wavelenght:float) -> np.ndarray: - """Fresnel propagator applied to data. +def fresnel_propagator( + data: np.ndarray, dist_observation: int, scale_factor: float, wavelenght: float +) -> np.ndarray: + """Fresnel propagator applied to data. Adapted from the script by Adrián Carbajal-Domínguez, adrian.carbajal@ujat.mx Args: @@ -575,7 +592,7 @@ def fresnel_propagator(data:np.ndarray, dist_observation: int, scale_factor: flo wavelenght (float): Wavelength. Returns: - np.ndarray: 2D or 3D data with propagation. + np.ndarray: 2D or 3D data with propagation. """ data_fresnel = data.copy() diff --git a/tomophantom/ctypes/external.py b/tomophantom/ctypes/external.py index 0e0f1ef..7ed6117 100644 --- a/tomophantom/ctypes/external.py +++ b/tomophantom/ctypes/external.py @@ -199,6 +199,7 @@ def c_model3d( ) return output + def c_object3d( output, N1, @@ -242,6 +243,7 @@ def c_object3d( ) return output + def c_model_sino3d( output, model_no, @@ -259,7 +261,7 @@ def c_model_sino3d( dtype.as_c_float_p(output), dtype.as_c_int(model_no), dtype.as_c_long(detector_horiz), - dtype.as_c_long(detector_vert), + dtype.as_c_long(detector_vert), dtype.as_c_long(Z1), dtype.as_c_long(Z2), dtype.as_c_long(phantom_size), @@ -269,6 +271,7 @@ def c_model_sino3d( ) return output + def c_object_sino3d( output, detector_horiz, @@ -292,7 +295,11 @@ def c_object_sino3d( tt, ): LIB_TOMOPHANTOM.TomoP3DObjectSino_core.restype = dtype.as_c_void_p() - if (("gaussian" in str(objectName)) or ("paraboloid" in str(objectName)) or ("ellipsoid" in str(objectName))): + if ( + ("gaussian" in str(objectName)) + or ("paraboloid" in str(objectName)) + or ("ellipsoid" in str(objectName)) + ): LIB_TOMOPHANTOM.TomoP3DObjectSino_core( dtype.as_c_float_p(output), dtype.as_c_long(detector_horiz), @@ -315,7 +322,7 @@ def c_object_sino3d( dtype.as_c_float(phi1), dtype.as_c_long(tt), ) - elif ("elliptical_cylinder" in str(objectName)): + elif "elliptical_cylinder" in str(objectName): LIB_TOMOPHANTOM.TomoP3DObjectSino_core( dtype.as_c_float_p(output), dtype.as_c_long(detector_horiz), @@ -361,4 +368,4 @@ def c_object_sino3d( dtype.as_c_float(-phi1), dtype.as_c_long(tt), ) - return output \ No newline at end of file + return output diff --git a/tomophantom/flatsgen.py b/tomophantom/flatsgen.py index f04e382..4599b1d 100644 --- a/tomophantom/flatsgen.py +++ b/tomophantom/flatsgen.py @@ -14,13 +14,21 @@ from tomophantom.artefacts import noise from tomophantom.supp.speckle_routines import simulate_speckles_with_shot_noise -def synth_flats(projData3D_clean: np.ndarray, source_intensity: int, - detectors_miscallibration: float = 0.05, - variations_number: int = 3, arguments_Bessel: tuple=(1,25), - specklesize: int = 2, kbar: int = 2, sigmasmooth: int = 3, - jitter_projections: float = 0.0, flatsnum: int =20) -> list: + +def synth_flats( + projData3D_clean: np.ndarray, + source_intensity: int, + detectors_miscallibration: float = 0.05, + variations_number: int = 3, + arguments_Bessel: tuple = (1, 25), + specklesize: int = 2, + kbar: int = 2, + sigmasmooth: int = 3, + jitter_projections: float = 0.0, + flatsnum: int = 20, +) -> list: """Function to generate synthetic flat field images and raw data for projection data normalisation. - This is more realistic modelling of stripes and various normalisation artefacts. + This is more realistic modelling of stripes and various normalisation artefacts. Args: projData3D_clean (np.ndarray): 3D projection data of the following shape: (DetectorsVert, anglesDim, DetectorsHoriz). @@ -36,82 +44,129 @@ def synth_flats(projData3D_clean: np.ndarray, source_intensity: int, Returns: list: List that contains [np.uint16(projData3D_raw), np.uint16(simulated_flats), blurred_speckles_map] - """ + """ [DetectorsDimV, projectionsNo, DetectorsDimH] = np.shape(projData3D_clean) - + # output datasets - flats_combined3D = np.zeros((DetectorsDimV,flatsnum, DetectorsDimH),dtype='uint16') - projData3D_raw = np.zeros(np.shape(projData3D_clean),dtype='float32') - + flats_combined3D = np.zeros( + (DetectorsDimV, flatsnum, DetectorsDimH), dtype="uint16" + ) + projData3D_raw = np.zeros(np.shape(projData3D_clean), dtype="float32") + # normalise the data projData3D_clean /= np.max(projData3D_clean) - + # using spherical Bessel functions to emulate the background (scintillator) variations - func = spherical_yn(1, np.linspace(arguments_Bessel[0], arguments_Bessel[1], DetectorsDimV,dtype='float32')) + func = spherical_yn( + 1, + np.linspace( + arguments_Bessel[0], arguments_Bessel[1], DetectorsDimV, dtype="float32" + ), + ) func += abs(np.min(func)) - - flatfield = np.zeros((DetectorsDimV,DetectorsDimH)) - for i in range(0,DetectorsDimH): - flatfield[:,i] = func - for i in range(0,DetectorsDimH): - flatfield[:,i] += np.flipud(func) - - if (specklesize != 0.0): + + flatfield = np.zeros((DetectorsDimV, DetectorsDimH)) + for i in range(0, DetectorsDimH): + flatfield[:, i] = func + for i in range(0, DetectorsDimH): + flatfield[:, i] += np.flipud(func) + + if specklesize != 0.0: # using speckle generator routines to create a photon count texture in the background - speckle_background = simulate_speckles_with_shot_noise([DetectorsDimV, DetectorsDimH], 1, specklesize, kbar) + speckle_background = simulate_speckles_with_shot_noise( + [DetectorsDimV, DetectorsDimH], 1, specklesize, kbar + ) else: - speckle_background = np.ones((DetectorsDimV,DetectorsDimH)) + speckle_background = np.ones((DetectorsDimV, DetectorsDimH)) # model miscallibrated detectors (a possible path to generate ring artifacts) blurred_speckles_map = np.zeros((DetectorsDimV, DetectorsDimH, variations_number)) - for i in range(0,variations_number): - speckles = simulate_speckles_with_shot_noise([DetectorsDimV, DetectorsDimH], 1, 10, 0.03) - #blur the speckled background + for i in range(0, variations_number): + speckles = simulate_speckles_with_shot_noise( + [DetectorsDimV, DetectorsDimH], 1, 10, 0.03 + ) + # blur the speckled background blurred_speckles = gaussian_filter(speckles.copy(), sigma=sigmasmooth) # threshold the result - blurred_speckles[blurred_speckles < 0.6*np.max(blurred_speckles)] = 0 - blurred_speckles_map[:,:,i] = blurred_speckles + blurred_speckles[blurred_speckles < 0.6 * np.max(blurred_speckles)] = 0 + blurred_speckles_map[:, :, i] = blurred_speckles blurred_speckles_map /= np.max(blurred_speckles_map) - - sinusoidal_response = np.sin(np.linspace(0,1.5*np.pi,projectionsNo)) + np.random.random(projectionsNo) * 0.1 + + sinusoidal_response = ( + np.sin(np.linspace(0, 1.5 * np.pi, projectionsNo)) + + np.random.random(projectionsNo) * 0.1 + ) sinusoidal_response /= np.max(sinusoidal_response) - exponential_response = np.exp(np.linspace(0,np.pi,projectionsNo)) + np.random.random(projectionsNo) * 0.1 + exponential_response = ( + np.exp(np.linspace(0, np.pi, projectionsNo)) + + np.random.random(projectionsNo) * 0.1 + ) exponential_response /= np.max(exponential_response) # prepeare flat fields - for i in range(0,flatsnum): + for i in range(0, flatsnum): # add speckled background to the initial image with the Bessel background - flatfield_combined = flatfield.copy() + 0.5*(speckle_background/np.max(speckle_background)) + flatfield_combined = flatfield.copy() + 0.5 * ( + speckle_background / np.max(speckle_background) + ) flatfield_combined /= np.max(flatfield_combined) - - #adding Poisson noise to flat fields - flatfield_poisson = noise(flatfield_combined*source_intensity, source_intensity, noisetype='Poisson') + + # adding Poisson noise to flat fields + flatfield_poisson = noise( + flatfield_combined * source_intensity, source_intensity, noisetype="Poisson" + ) flatfield_poisson /= np.max(flatfield_poisson) - flats_combined3D[:,i,:] = np.uint16(flatfield_poisson*65535) + flats_combined3D[:, i, :] = np.uint16(flatfield_poisson * 65535) # convert synthetic projections to raw-data like projection ready for normalisation - for i in range(0,projectionsNo): - proj_exp = np.exp(-projData3D_clean[:,i,:])*source_intensity*flatfield_poisson # raw projection - for j in range(0,variations_number): + for i in range(0, projectionsNo): + proj_exp = ( + np.exp(-projData3D_clean[:, i, :]) * source_intensity * flatfield_poisson + ) # raw projection + for j in range(0, variations_number): if j == 0: # adding a consistent offset for certain detectors - proj_exp += blurred_speckles_map[:,:,j]*detectors_miscallibration*source_intensity + proj_exp += ( + blurred_speckles_map[:, :, j] + * detectors_miscallibration + * source_intensity + ) if j == 1: # adding a sinusoidal-like response offset for certain detectors - proj_exp += sinusoidal_response[i]*blurred_speckles_map[:,:,j]*detectors_miscallibration*source_intensity + proj_exp += ( + sinusoidal_response[i] + * blurred_speckles_map[:, :, j] + * detectors_miscallibration + * source_intensity + ) if j == 2: # adding an exponential response offset for certain detectors - proj_exp += exponential_response[i]*blurred_speckles_map[:,:,j]*detectors_miscallibration*source_intensity - - projection_poisson = noise(proj_exp, source_intensity, noisetype='Poisson') + proj_exp += ( + exponential_response[i] + * blurred_speckles_map[:, :, j] + * detectors_miscallibration + * source_intensity + ) + + projection_poisson = noise(proj_exp, source_intensity, noisetype="Poisson") # apply jitter to projections if jitter_projections != 0.0: - horiz_shift = random.uniform(-jitter_projections,jitter_projections) #generate random directional shift - vert_shift = random.uniform(-jitter_projections,jitter_projections) #generate random directional shift - projection_poisson = shift(projection_poisson.copy(),[vert_shift,horiz_shift], mode='reflect') - projData3D_raw[:,i,:] = projection_poisson + horiz_shift = random.uniform( + -jitter_projections, jitter_projections + ) # generate random directional shift + vert_shift = random.uniform( + -jitter_projections, jitter_projections + ) # generate random directional shift + projection_poisson = shift( + projection_poisson.copy(), [vert_shift, horiz_shift], mode="reflect" + ) + projData3D_raw[:, i, :] = projection_poisson projData3D_raw /= np.max(projData3D_raw) - return [np.uint16(projData3D_raw*65535), np.uint16(flats_combined3D), blurred_speckles_map] + return [ + np.uint16(projData3D_raw * 65535), + np.uint16(flats_combined3D), + blurred_speckles_map, + ] diff --git a/tomophantom/generator.py b/tomophantom/generator.py index 703a854..8da77cd 100644 --- a/tomophantom/generator.py +++ b/tomophantom/generator.py @@ -6,167 +6,238 @@ The TomoPhantom package is released under Apache License, Version 2.0 @author: Daniil Kazantsev """ + + def rand_init2D(x0min, x0max, y0min, y0max, c0min, c0max, ab_min, ab_max): import numpy as np - x0 = np.random.uniform(low = x0min, high=x0max) - y0 = np.random.uniform(low = y0min, high=y0max) - c0 = np.random.uniform(low = c0min, high=c0max) - ab = np.random.uniform(low = ab_min, high=ab_max) - return (x0,y0,c0,ab) + + x0 = np.random.uniform(low=x0min, high=x0max) + y0 = np.random.uniform(low=y0min, high=y0max) + c0 = np.random.uniform(low=c0min, high=c0max) + ab = np.random.uniform(low=ab_min, high=ab_max) + return (x0, y0, c0, ab) + def rand_init3D(x0min, x0max, y0min, y0max, z0min, z0max, c0min, c0max, ab_min, ab_max): import numpy as np - x0 = np.random.uniform(low = x0min, high=x0max) - y0 = np.random.uniform(low = y0min, high=y0max) - z0 = np.random.uniform(low = z0min, high=z0max) - c0 = np.random.uniform(low = c0min, high=c0max) - ab = np.random.uniform(low = ab_min, high=ab_max) - return (x0,y0,z0,c0,ab) + + x0 = np.random.uniform(low=x0min, high=x0max) + y0 = np.random.uniform(low=y0min, high=y0max) + z0 = np.random.uniform(low=z0min, high=z0max) + c0 = np.random.uniform(low=c0min, high=c0max) + ab = np.random.uniform(low=ab_min, high=ab_max) + return (x0, y0, z0, c0, ab) + # Function to generate 2D foam-like structures using randomly located circles -def foam2D(x0min, x0max, y0min, y0max, c0min, c0max, ab_min, ab_max, N_size, tot_objects, object_type): +def foam2D( + x0min, + x0max, + y0min, + y0max, + c0min, + c0max, + ab_min, + ab_max, + N_size, + tot_objects, + object_type, +): import numpy as np import math import random - #2D functions - from tomophantom import TomoP2D + + # 2D functions + from tomophantom import TomoP2D from tomophantom.TomoP2D import Objects2D - attemptsNo = 2000 # the number of attempts to fit the object + + attemptsNo = 2000 # the number of attempts to fit the object # objects accepted: 'ellipse', 'parabola', 'gaussian', 'mix' mix_objects = False - if (object_type == 'ellipse'): + if object_type == "ellipse": object_type = Objects2D.ELLIPSE - elif (object_type == 'parabola'): + elif object_type == "parabola": object_type = Objects2D.PARABOLA - elif (object_type == 'gaussian'): + elif object_type == "gaussian": object_type = Objects2D.GAUSSIAN - elif (object_type == 'mix'): + elif object_type == "mix": mix_objects = True - else: - raise TypeError('object_type can be only ellipse, parabola, gaussian or mix') + else: + raise TypeError("object_type can be only ellipse, parabola, gaussian or mix") X0 = np.float32(np.zeros(tot_objects)) Y0 = np.float32(np.zeros(tot_objects)) AB = np.float32(np.zeros(tot_objects)) - C0_var = np.float32(np.zeros(tot_objects)) - for i in range(0,tot_objects): - (x0,y0,c0,ab) = rand_init2D(x0min, x0max, y0min, y0max, c0min, c0max, ab_min, ab_max) - if (i > 0): + C0_var = np.float32(np.zeros(tot_objects)) + for i in range(0, tot_objects): + (x0, y0, c0, ab) = rand_init2D( + x0min, x0max, y0min, y0max, c0min, c0max, ab_min, ab_max + ) + if i > 0: breakj = False - for j in range(0,attemptsNo): - if (breakj == True): - (x0,y0,c0,ab) = rand_init2D(x0min, x0max, y0min, y0max, c0min, c0max, ab_min, ab_max) + for j in range(0, attemptsNo): + if breakj == True: + (x0, y0, c0, ab) = rand_init2D( + x0min, x0max, y0min, y0max, c0min, c0max, ab_min, ab_max + ) breakj = False else: - for l in range(0,i): # checks consistency with previously created objects - dist = math.sqrt((X0[l]-x0)**2 + (Y0[l]-y0)**2) - if (dist < (ab + AB[l])) or ((abs(x0) + ab)**2 + (abs(y0) + ab)**2 > 1.0): + for l in range( + 0, i + ): # checks consistency with previously created objects + dist = math.sqrt((X0[l] - x0) ** 2 + (Y0[l] - y0) ** 2) + if (dist < (ab + AB[l])) or ( + (abs(x0) + ab) ** 2 + (abs(y0) + ab) ** 2 > 1.0 + ): breakj = True - break - if (breakj == False): # re-initialise if doesn't fit the criteria + break + if breakj == False: # re-initialise if doesn't fit the criteria X0[i] = x0 Y0[i] = y0 AB[i] = ab C0_var[i] = c0 break - if (AB[i] == 0.0): + if AB[i] == 0.0: X0[i] = x0 Y0[i] = y0 AB[i] = 0.0001 - C0_var[i] = c0 - - myObjects = [] # dictionary of objects - for obj in range(0,len(X0)): - if (mix_objects == True): - rand_obj = random.randint(0,2) - if (rand_obj == 0): + C0_var[i] = c0 + + myObjects = [] # dictionary of objects + for obj in range(0, len(X0)): + if mix_objects == True: + rand_obj = random.randint(0, 2) + if rand_obj == 0: object_type = Objects2D.ELLIPSE - if (rand_obj == 1): + if rand_obj == 1: object_type = Objects2D.PARABOLA - if (rand_obj == 2): - object_type = Objects2D.GAUSSIAN - curr_obj = {'Obj': object_type, - 'C0' : C0_var[obj], - 'x0' : X0[obj], - 'y0' : Y0[obj], - 'a' : AB[obj], - 'b' : AB[obj], - 'phi': 0.0} + if rand_obj == 2: + object_type = Objects2D.GAUSSIAN + curr_obj = { + "Obj": object_type, + "C0": C0_var[obj], + "x0": X0[obj], + "y0": Y0[obj], + "a": AB[obj], + "b": AB[obj], + "phi": 0.0, + } myObjects.append(curr_obj) Object = TomoP2D.Object(N_size, myObjects) - return (Object,myObjects) + return (Object, myObjects) + # Function to generate 3D foam-like structures using randomly located spheres -def foam3D(x0min, x0max, y0min, y0max, z0min, z0max, c0min, c0max, ab_min, ab_max, N_size, tot_objects, object_type): +def foam3D( + x0min, + x0max, + y0min, + y0max, + z0min, + z0max, + c0min, + c0max, + ab_min, + ab_max, + N_size, + tot_objects, + object_type, +): import numpy as np import math import random - #3D functions - from tomophantom import TomoP3D + + # 3D functions + from tomophantom import TomoP3D from tomophantom.TomoP3D import Objects3D - attemptsNo = 2000 + + attemptsNo = 2000 # objects accepted: 'ellipsoid', 'paraboloid', 'gaussian', 'mix' mix_objects = False - if (object_type == 'ellipsoid'): + if object_type == "ellipsoid": object_type = Objects3D.ELLIPSOID - elif (object_type == 'paraboloid'): + elif object_type == "paraboloid": object_type = Objects3D.PARABOLOID - elif (object_type == 'gaussian'): + elif object_type == "gaussian": object_type = Objects3D.GAUSSIAN - elif (object_type == 'mix'): + elif object_type == "mix": mix_objects = True - else: - raise TypeError('object_type can be only ellipse, parabola, gaussian or mix') + else: + raise TypeError("object_type can be only ellipse, parabola, gaussian or mix") X0 = np.float32(np.zeros(tot_objects)) Y0 = np.float32(np.zeros(tot_objects)) Z0 = np.float32(np.zeros(tot_objects)) AB = np.float32(np.zeros(tot_objects)) - C0_var = np.float32(np.zeros(tot_objects)) - for i in range(0,tot_objects): - (x0,y0,z0,c0,ab) = rand_init3D(x0min, x0max, y0min, y0max, z0min, z0max, c0min, c0max, ab_min, ab_max) - if (i > 0): + C0_var = np.float32(np.zeros(tot_objects)) + for i in range(0, tot_objects): + (x0, y0, z0, c0, ab) = rand_init3D( + x0min, x0max, y0min, y0max, z0min, z0max, c0min, c0max, ab_min, ab_max + ) + if i > 0: breakj = False - for j in range(0,attemptsNo): + for j in range(0, attemptsNo): if breakj: - (x0,y0,z0,c0,ab) = rand_init3D(x0min, x0max, y0min, y0max, z0min, z0max, c0min, c0max, ab_min, ab_max) + (x0, y0, z0, c0, ab) = rand_init3D( + x0min, + x0max, + y0min, + y0max, + z0min, + z0max, + c0min, + c0max, + ab_min, + ab_max, + ) breakj = False else: - for l in range(0,i): # checks consistency with previously created objects - dist = math.sqrt((X0[l]-x0)**2 + (Y0[l]-y0)**2 + (Z0[l]-z0)**2) - if (dist < (ab + AB[l])) or ((abs(x0) + ab)**2 + (abs(y0) + ab)**2 + (abs(z0) + ab)**2 > 1.0): + for l in range( + 0, i + ): # checks consistency with previously created objects + dist = math.sqrt( + (X0[l] - x0) ** 2 + (Y0[l] - y0) ** 2 + (Z0[l] - z0) ** 2 + ) + if (dist < (ab + AB[l])) or ( + (abs(x0) + ab) ** 2 + + (abs(y0) + ab) ** 2 + + (abs(z0) + ab) ** 2 + > 1.0 + ): breakj = True break - if (breakj == False): # re-initialise if doesn't fit the criteria + if breakj == False: # re-initialise if doesn't fit the criteria X0[i] = x0 Y0[i] = y0 Z0[i] = z0 AB[i] = ab C0_var[i] = c0 break - if (AB[i] == 0.0): + if AB[i] == 0.0: X0[i] = x0 Y0[i] = y0 AB[i] = 0.0001 - C0_var[i] = c0 + C0_var[i] = c0 - myObjects = [] # dictionary of objects - for obj in range(0,len(X0)): - if (mix_objects == True): - rand_obj = random.randint(0,2) - if (rand_obj == 0): + myObjects = [] # dictionary of objects + for obj in range(0, len(X0)): + if mix_objects == True: + rand_obj = random.randint(0, 2) + if rand_obj == 0: object_type = Objects3D.ELLIPSOID - if (rand_obj == 1): + if rand_obj == 1: object_type = Objects3D.PARABOLOID - if (rand_obj == 2): - object_type = Objects3D.GAUSSIAN - curr_obj = {'Obj': object_type, - 'C0' : C0_var[obj], - 'x0' : X0[obj], - 'y0' : Y0[obj], - 'z0' : Z0[obj], - 'a' : AB[obj], - 'b' : AB[obj], - 'c' : AB[obj], - 'phi1': 0.0} + if rand_obj == 2: + object_type = Objects3D.GAUSSIAN + curr_obj = { + "Obj": object_type, + "C0": C0_var[obj], + "x0": X0[obj], + "y0": Y0[obj], + "z0": Z0[obj], + "a": AB[obj], + "b": AB[obj], + "c": AB[obj], + "phi1": 0.0, + } myObjects.append(curr_obj) Object3D = TomoP3D.Object(N_size, myObjects) - return (Object3D,myObjects) \ No newline at end of file + return (Object3D, myObjects) diff --git a/tomophantom/supp/speckle_routines.py b/tomophantom/supp/speckle_routines.py index e6a5be3..35d6acc 100755 --- a/tomophantom/supp/speckle_routines.py +++ b/tomophantom/supp/speckle_routines.py @@ -4,119 +4,125 @@ import numpy as np from scipy.ndimage import gaussian_filter -def sample2ddist(asicshape,kbar,dist): + +def sample2ddist(asicshape, kbar, dist): """ samples photons from a 2D distribution """ - #Definitions - nasic = asicshape[0]*asicshape[1] - nphotons = np.round(kbar*nasic).astype(int) - photons = np.zeros(asicshape) - - #sample photon positions from given distribution - sampl = np.random.choice(nasic,nphotons,p=dist) - - #convert to 2d and accumulate photons - for i in range (nphotons): - ix = int(sampl[i]/asicshape[1]) - iy = int(sampl[i]%asicshape[1]) - photons[ix,iy] += 1 - + # Definitions + nasic = asicshape[0] * asicshape[1] + nphotons = np.round(kbar * nasic).astype(int) + photons = np.zeros(asicshape) + + # sample photon positions from given distribution + sampl = np.random.choice(nasic, nphotons, p=dist) + + # convert to 2d and accumulate photons + for i in range(nphotons): + ix = int(sampl[i] / asicshape[1]) + iy = int(sampl[i] % asicshape[1]) + photons[ix, iy] += 1 + return photons -def model_speckles(asicshape,specklesize): + +def model_speckles(asicshape, specklesize): """ models a speckle pattern with a given speckle size (integer) based on random phasors (see JC Goodman - Speckle Phenomena in Optics - Appendix G) """ - #definitions speckles - randphasors = np.zeros([asicshape[0],asicshape[1]],dtype = complex) - nphasorsx = np.round(asicshape[0]/specklesize).astype(int) - nphasorsy = np.round(asicshape[1]/specklesize).astype(int) + # definitions speckles + randphasors = np.zeros([asicshape[0], asicshape[1]], dtype=complex) + nphasorsx = np.round(asicshape[0] / specklesize).astype(int) + nphasorsy = np.round(asicshape[1] / specklesize).astype(int) - #random phasors + # random phasors for i in range(nphasorsx): for j in range(nphasorsy): - randphase = np.random.random_sample()*2.j*np.pi - randphasors[i,j] = np.exp(randphase) + randphase = np.random.random_sample() * 2.0j * np.pi + randphasors[i, j] = np.exp(randphase) - #speckles + # speckles specklefield = np.fft.fft2(randphasors) - speckleint = np.absolute(specklefield)*np.absolute(specklefield) - speckleint/=np.sum(speckleint) #normalise - - return speckleint + speckleint = np.absolute(specklefield) * np.absolute(specklefield) + speckleint /= np.sum(speckleint) # normalise + + return speckleint -def model_speckles_modes(asicshape,specklesize,modes): + +def model_speckles_modes(asicshape, specklesize, modes): """ models a speckle pattern with a given speckle size (integer) and - contrast (1/modes) + contrast (1/modes) """ - - #definitions - nx,ny = asicshape[0],asicshape[1] + + # definitions + nx, ny = asicshape[0], asicshape[1] # add speckle patterns as many times as the number of modes - dist = model_speckles([nx,ny],specklesize) - for m in range(modes-1): - dist += model_speckles([nx,ny],specklesize) - dist/=np.sum(dist) # normalise + dist = model_speckles([nx, ny], specklesize) + for m in range(modes - 1): + dist += model_speckles([nx, ny], specklesize) + dist /= np.sum(dist) # normalise + + return dist - return dist -def simulate_speckles_with_shot_noise(asicshape,modes,specklesize,kbar): +def simulate_speckles_with_shot_noise(asicshape, modes, specklesize, kbar): """ returns a 2D speckle pattern with given speckle size, contrast - (1/modes) and photon density kbar (photons/pixel). + (1/modes) and photon density kbar (photons/pixel). """ - #definitions - nx,ny = asicshape[0],asicshape[1] - dist = model_speckles_modes(asicshape,specklesize,modes).flatten() + # definitions + nx, ny = asicshape[0], asicshape[1] + dist = model_speckles_modes(asicshape, specklesize, modes).flatten() # sample photons from the distribution - sig =sample2ddist([nx,ny],kbar,dist) + sig = sample2ddist([nx, ny], kbar, dist) return sig -def simulate_shot_noise(asicshape,kbar): + +def simulate_shot_noise(asicshape, kbar): """ returns a 2d shot noise pattern with given dimensions and photon density kbar (photons/pixel) """ - #definitions - nx,ny = asicshape[0],asicshape[1] + # definitions + nx, ny = asicshape[0], asicshape[1] # flat 2d distribution - dist = np.ones((nx,ny)).flatten() - dist/=np.sum(dist) + dist = np.ones((nx, ny)).flatten() + dist /= np.sum(dist) # sample photons from the distribution - sig =sample2ddist([nx,ny],kbar,dist) + sig = sample2ddist([nx, ny], kbar, dist) + + return np.array(sig, dtype=float) - return np.array(sig,dtype = float) -def simulate_charge_sharing(asicshape,modes,specklesize,kbar,sigma=1): +def simulate_charge_sharing(asicshape, modes, specklesize, kbar, sigma=1): """ returns a 2D speckle pattern with given speckle size, contrast (1/modes) and photon density kbar (photons/pixel) including charge sharing simulated with a gaussian at each photon position with with sigma """ - #Definitions - nasic = asicshape[0]*asicshape[1] - nphotons = np.round(kbar*nasic).astype(int) + # Definitions + nasic = asicshape[0] * asicshape[1] + nphotons = np.round(kbar * nasic).astype(int) photons = np.zeros(asicshape) photons_blur = np.zeros(asicshape) - dist = model_speckles_modes(asicshape,specklesize,modes).flatten() + dist = model_speckles_modes(asicshape, specklesize, modes).flatten() - #sample photon positions from given distribution - sampl = np.random.choice(nasic,nphotons,p=dist) + # sample photon positions from given distribution + sampl = np.random.choice(nasic, nphotons, p=dist) - #convert to 2d and accumulate photons - for i in range (nphotons): - ix = int(sampl[i]/asicshape[1]) - iy = int(sampl[i]%asicshape[1]) + # convert to 2d and accumulate photons + for i in range(nphotons): + ix = int(sampl[i] / asicshape[1]) + iy = int(sampl[i] % asicshape[1]) photons_blur = np.zeros(asicshape) - photons_blur[ix,iy] = 1 - photons_blur = gaussian_filter(photons_blur,sigma=sigma) + photons_blur[ix, iy] = 1 + photons_blur = gaussian_filter(photons_blur, sigma=sigma) photons += photons_blur - return photons \ No newline at end of file + return photons