Merge pull request #115 from dkazanc/flats

flat field generator has been re-written

Demos/Python/2D/ReconASTRA.py

@@ -62,6 +62,7 @@
             'noise_amplitude' : 1e5, # noise amplitude
             'noise_seed' : 0}
 
+
 noisy_sino = _Artifacts_(sino_an, **_noise_)
 
 plt.figure()

Demos/Python/2D/Recon_artifacts_ASTRA.py

@@ -25,9 +25,7 @@
 from tomophantom.supp.qualitymetrics import QualityTools
 
 model = 15 # select a model
-N_size = 512 # set dimension of the phantom
-# one can specify an exact path to the parameters file
-# path_library2D = '../../../PhantomLibrary/models/Phantom2DLibrary.dat'
+N_size = 256 # set dimension of the phantom
 path = os.path.dirname(tomophantom.__file__)
 path_library2D = os.path.join(path, "Phantom2DLibrary.dat")
 phantom_2D = TomoP2D.Model(model, N_size, path_library2D)
@@ -59,11 +57,10 @@
 plt.close('all')
 # forming dictionaries with artifact types
 _noise_ =  {'noise_type' : 'Poisson',
-            'noise_amplitude' : 10000, # noise amplitude
-            'noise_seed' : 0}
+            'noise_amplitude' : 10000}
 # misalignment dictionary
-_sinoshifts_ = {'sinoshifts_maxamplitude' : 10}
-[noisy_sino_misalign,shifts] = _Artifacts_(sino_an, **_noise_, **_sinoshifts_)
+_sinoshifts_ = {'datashifts_maxamplitude_pixel' : 10}
+[noisy_sino_misalign, shifts] = _Artifacts_(sino_an, **_noise_, **_sinoshifts_)
 
 # adding zingers and stripes
 _zingers_ = {'zingers_percentage' : 2,

Demos/Python/3D/ReconASTRA3D_misalignment.py

@@ -140,7 +140,7 @@
     # reconstruct with the estimated shifts
     RectoolsDIR = RecToolsDIR(DetectorsDimH = Horiz_det,     # Horizontal detector dimension
                         DetectorsDimV = Vert_det,            # Vertical detector dimension (3D case)
-                        CenterRotOffset = -shifts_estimated, # Center of Rotation scalar + shifts passed
+                        CenterRotOffset = shifts_estimated, # Center of Rotation scalar + shifts passed
                         AnglesVec = angles_rad,              # A vector of projection angles in radians
                         ObjSize = N_size,                    # Reconstructed object dimensions (scalar)
                         device_projector='gpu')

Demos/Python/3D/ReconASTRA3D_realistic.py

@@ -28,7 +28,7 @@
 
 print ("Building 3D phantom using TomoPhantom software")
 tic=timeit.default_timer()
-model = 17 # select a model number from the library
+model = 13 # select a model number from the library
 N_size = 256 # Define phantom dimensions using a scalar value (cubic phantom)
 path = os.path.dirname(tomophantom.__file__)
 path_library3D = os.path.join(path, "Phantom3DLibrary.dat")
@@ -65,7 +65,7 @@
 projData3D_analyt= TomoP3D.ModelSino(model, N_size, Horiz_det, Vert_det, angles, path_library3D)
 
 intens_max_clean = np.max(projData3D_analyt)
-sliceSel = 150
+sliceSel = (int)(N_size*0.5)
 plt.figure() 
 plt.subplot(131)
 plt.imshow(projData3D_analyt[:,sliceSel,:],vmin=0, vmax=intens_max_clean)
@@ -79,44 +79,47 @@
 plt.show()
 #%%
 print ("Simulate synthetic flat fields, add flat field background to the projections and add noise")
-I0  = 15000; # Source intensity
+I0  = 50000; # full-beam photon flux intensity
 flatsnum = 20 # the number of the flat fields required
 
-[projData3D_noisy, flatsSIM] = synth_flats(projData3D_analyt, 
-                                           source_intensity = I0, source_variation=0.02,\
-                                           arguments_Bessel = (1,10,10,12),\
-                                           specklesize = 15,\
-                                           kbar = 0.3,
-                                           jitter = 1.0,
-                                           sigmasmooth = 3, flatsnum=flatsnum)
+[projData3D_raw, flats_combined3D, speckel_map] = synth_flats(projData3D_analyt, 
+                                           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)
 #del projData3D_analyt
-plt.figure() 
+plt.figure()
 plt.subplot(121)
-plt.imshow(projData3D_noisy[:,0,:])
+plt.imshow(projData3D_raw[:,0,:])
 plt.title('2D Projection (before normalisation)')
 plt.subplot(122)
-plt.imshow(flatsSIM[:,0,:])
+plt.imshow(flats_combined3D[:,0,:])
 plt.title('A selected simulated flat-field')
 plt.show()
 #%%
 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_raw, flats_combined3D, darks=None, log='true', method='mean')
+projData3D_norm *= intens_max_clean
 
 #del projData3D_noisy
-intens_max = 0.3*np.max(projData3D_norm)
-sliceSel = 150
+sliceSel = (int)(N_size*0.5)
 plt.figure() 
 plt.subplot(131)
-plt.imshow(projData3D_norm[:,sliceSel,:],vmin=0, vmax=intens_max)
+plt.imshow(projData3D_norm[:,sliceSel,:],vmin=0, vmax=intens_max_clean)
 plt.title('Normalised 2D Projection (erroneous)')
 plt.subplot(132)
-plt.imshow(projData3D_norm[sliceSel,:,:],vmin=0, vmax=intens_max)
+plt.imshow(projData3D_norm[sliceSel,:,:],vmin=0, vmax=intens_max_clean)
 plt.title('Sinogram view')
 plt.subplot(133)
-plt.imshow(projData3D_norm[:,:,sliceSel],vmin=0, vmax=intens_max)
+plt.imshow(projData3D_norm[:,:,sliceSel],vmin=0, vmax=intens_max_clean)
 plt.title('Tangentogram view')
 plt.show()
 #%%
@@ -131,7 +134,6 @@
 
 print ("Reconstruction using FBP from tomobar")
 recNumerical= RectoolsDIR.FBP(projData3D_norm) # FBP reconstruction
-recNumerical *= intens_max_clean
 
 sliceSel = int(0.5*N_size)
 max_val = 1
@@ -184,7 +186,6 @@
                     'device_regulariser': 'gpu'}
 
 RecFISTA_os_reg = Rectools.FISTA(_data_, _algorithm_, _regularisation_)
-RecFISTA_os_reg *= intens_max_clean
 
 sliceSel = int(0.5*N_size)
 max_val = 1

Demos/Python/3D/ReconASTRA_object3D.py

@@ -20,8 +20,10 @@
 from tomophantom.supp.qualitymetrics import QualityTools
 from tomophantom import TomoP3D
 from tomophantom.TomoP3D import Objects3D
+import timeit
 
-N3D_size = 256
+
+N3D_size = 512
 
 # specify object parameters, here we replicate model 
 obj3D_1 = {'Obj': Objects3D.GAUSSIAN, 
@@ -92,13 +94,17 @@
 from tomobar.methodsDIR import RecToolsDIR
 RectoolsDIR = RecToolsDIR(DetectorsDimH = Horiz_det,  # DetectorsDimH # detector dimension (horizontal)
                     DetectorsDimV = Vert_det,  # DetectorsDimV # detector dimension (vertical) for 3D case only
-                    CenterRotOffset = None, # Center of Rotation (CoR) scalar (for 3D case only)
+                    CenterRotOffset = 0.0, # Center of Rotation (CoR) scalar (for 3D case only)
                     AnglesVec = angles_rad, # array of angles in radians
                     ObjSize = N3D_size, # a scalar to define reconstructed object dimensions
                     device_projector = 'gpu')
 
 print ("Reconstruction using FBP from tomobar")
+tic=timeit.default_timer()
 recNumerical= RectoolsDIR.FBP(ProjData3D) # FBP reconstruction
+toc=timeit.default_timer()
+Run_time = toc - tic
+print("Reconstruction in {} seconds".format(Run_time))
 
 sliceSel = int(0.5*N3D_size)
 max_val = 1

Wrappers/Python/tomophantom/supp/artifacts.py

@@ -2,7 +2,7 @@
 # -*- coding: utf-8 -*-
 #    Note that the TomoPhantom package is released under Apache License, Version 2.0
 #    @author: Daniil Kazantsev
-#    latest update: 29.11.2021
+#    latest update: 07.10.2022
 
 import numpy as np
 import random
@@ -279,37 +279,38 @@ def zingers(data, percentage, modulus):
     sino_zingers[:] = sino_zingers_fl.reshape(sino_zingers.shape)
     return sino_zingers
 
-def noise(data, sigma, noisetype, seed, prelog):
+def noise(data, sigma, noisetype, seed=True, prelog=False):
     # Adding random noise to data (adapted from LD-CT simulator)
-    # noisetype = None, # 'Gaussian', 'Poisson' or None
-    # sigma = 10000, # photon flux (Poisson) or variance for Gaussian
-    # seed = 0, # seeds for noise
-    # prelog: None or True (get the raw pre-log data)
-    sino_noisy = data.copy()
+    # noisetype = 'Gaussian' or 'Poisson'
+    # sigma = 10000, # photon flux (Poisson) or the variance for Gaussian
+    # seed = 0, # initiate random seed with True
+    # prelog: True or False
+    data_noisy = data.copy()
+    maxData = np.max(data)
+    data_noisy=data/maxData #normalising
+    if seed is True:
+        np.random.seed(int(seed))
     if noisetype == 'Gaussian':
         # add normal Gaussian noise
-        if seed is not None:
-            np.random.seed(int(seed))
-        sino_noisy += np.random.normal(loc = 0.0, scale = sigma, size = np.shape(sino_noisy))
-        sino_noisy[sino_noisy<0] = 0
+        data_noisy += np.random.normal(loc = 0.0, scale = sigma, size = np.shape(data_noisy))
+        data_noisy[data_noisy<0] = 0
+        data_noisy *= maxData
     elif noisetype == 'Poisson':
         # add Poisson noise
-        maxSino = np.max(data)
-        if maxSino > 0:
+        if maxData > 0:
             ri = 1.0
-            sig = np.sqrt(11) #standard variance of electronic noise, a characteristic of CT scanner
-            sino_noisy=data/maxSino
-            yb = sigma*np.exp(-sino_noisy) + ri # exponential transform to incident flux
-            sino_raw = np.random.poisson(yb) + np.sqrt(sig)*np.reshape(np.random.randn(np.size(yb)),np.shape(yb))
-            li_hat = -np.log((sino_raw-ri)/sigma)*maxSino # log corrected data
-            li_hat[(sino_raw-ri)<=0] = 0.0
-            sino_noisy = li_hat.copy()
-    else:
-        print ("Select 'Gaussian' or 'Poisson' for noise type")
+            sig = np.sqrt(11) #standard variance of electronic noise, a characteristic of a CT scanner
+            yb = sigma*np.exp(-data_noisy) + ri # exponential transform to incident flux
+            data_raw = np.random.poisson(yb) + np.sqrt(sig)*np.reshape(np.random.randn(np.size(yb)),np.shape(yb))
+            li_hat = -np.log((data_raw-ri)/sigma)*maxData # log corrected data
+            li_hat[(data_raw-ri)<=0] = 0.0
+            data_noisy = li_hat.copy()
+    else:
+        print ("Select 'Gaussian' or 'Poisson' for the noise type")
     if prelog is True:
-        return [sino_noisy,sino_raw]
+        return [data_noisy,data_raw]
     else:
-        return sino_noisy
+        return data_noisy
 
 def datashifts(data, maxamplitude):
     # A function to add random shifts to sinogram rows (an offset for each angular position)