diff --git a/.travis.yml b/.travis.yml index 4f01b8d20..e9cd31321 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,13 +1,6 @@ language: python matrix: - include: - - python: 2.7 - env: NUMPY=1.16 - - - python: 2.7 - env: NUMPY=1.15 - - python: 3.5 env: NUMPY=1.15 @@ -16,7 +9,13 @@ matrix: - python: 3.6 env: NUMPY=1.16 - + + - python: 3.7 + env: NUMPY=1.15 + + - python: 3.7 + env: NUMPY=1.16 + os: - linux @@ -41,10 +40,10 @@ install: - export VERSION=`date +%Y.%m` - conda build conda-recipe --numpy=$NUMPY --python=$TRAVIS_PYTHON_VERSION - conda install --channel /home/travis/miniconda/envs/test-environment/conda-bld/ tomobar --offline --override-channels - + after_success: - chmod +x src/Python/conda-recipe/conda_upload.sh - test $TRAVIS_BRANCH = "master" && bash conda-recipe/conda_upload.sh script: - - python tests/test.py + - python test/test_tomobarCPU_DIR.py diff --git a/Demos/Python/DemoFISTA_artifacts2D.py b/Demos/Python/DemoFISTA_artifacts2D.py index c89b9d4ec..6c6ad5f13 100755 --- a/Demos/Python/DemoFISTA_artifacts2D.py +++ b/Demos/Python/DemoFISTA_artifacts2D.py @@ -117,7 +117,7 @@ plt.title('Misaligned noisy FBP Reconstruction') plt.show() -plt.figure() +plt.figure() plt.imshow(abs(FBPrec_ideal-FBPrec_error), vmin=0, vmax=2, cmap="gray") plt.colorbar(ticks=[0, 0.5, 2], orientation='vertical') plt.title('FBP reconstruction differences') @@ -133,7 +133,7 @@ FBPrec_misalign = Rectools.FBP(sino_misalign) # reconstruction with misalignment -plt.figure() +plt.figure() plt.imshow(FBPrec_misalign, vmin=0, vmax=1, cmap="gray") plt.title('FBP reconstruction of misaligned data using known exact shifts') #%% @@ -153,7 +153,7 @@ # data dictionary _data_ = {'projection_norm_data' : sino_artifacts, 'projection_raw_data' : sino_artifacts_raw/np.max(sino_artifacts_raw), - 'OS_number' : 10} + 'OS_number' : 10} lc = RectoolsIR.powermethod(_data_) # calculate Lipschitz constant (run once to initialise) _algorithm_ = {'iterations' : 30, @@ -169,7 +169,7 @@ print("Run FISTA reconstrucion algorithm with regularisation...") RecFISTA_LS_reg = RectoolsIR.FISTA(_data_, _algorithm_, _regularisation_) -# adding Huber data fidelity threshold +# adding Huber data fidelity threshold _data_.update({'huber_threshold' : 7.0}) print(" Run FISTA reconstrucion algorithm with regularisation and Huber data...") RecFISTA_Huber_reg = RectoolsIR.FISTA(_data_, _algorithm_, _regularisation_) @@ -196,7 +196,7 @@ plt.title('FISTA-HuberRing-TV reconstruction') plt.show() -# calculate errors +# calculate errors Qtools = QualityTools(phantom_2D[indicesROI], RecFISTA_LS_reg[indicesROI]) RMSE_FISTA_LS_TV = Qtools.rmse() Qtools = QualityTools(phantom_2D[indicesROI], RecFISTA_Huber_reg[indicesROI]) @@ -256,7 +256,7 @@ plt.title('FISTA-SWLS-Huber-TV reconstruction') plt.show() -# calculate errors +# calculate errors Qtools = QualityTools(phantom_2D[indicesROI], RecFISTA_SWLS[indicesROI]) RMSE_FISTA_SWLS_TV = Qtools.rmse() Qtools = QualityTools(phantom_2D[indicesROI], RecFISTA_SWLS_Huber[indicesROI]) @@ -279,7 +279,7 @@ 'mask_diameter' : 0.9, 'lipschitz_const' : lc} -# Run FISTA reconstrucion algorithm with regularisation +# Run FISTA reconstrucion algorithm with regularisation RecFISTA_LS_GH_reg = RectoolsIR.FISTA(_data_, _algorithm_, _regularisation_) plt.figure() @@ -288,7 +288,7 @@ plt.title('FISTA-OS-GH-TV reconstruction') plt.show() -# calculate errors +# calculate errors Qtools = QualityTools(phantom_2D[indicesROI], RecFISTA_LS_GH_reg[indicesROI]) RMSE_FISTA_LS_GH_TV = Qtools.rmse() print("RMSE for FISTA-PWLS-GH-TV reconstruction is {}".format(RMSE_FISTA_LS_GH_TV)) @@ -312,7 +312,7 @@ 'iterations' : 120, 'device_regulariser': 'gpu'} -# Run FISTA reconstrucion algorithm with regularisation +# Run FISTA reconstrucion algorithm with regularisation RecFISTA_LS_stud_reg = RectoolsIR.FISTA(_data_, _algorithm_, _regularisation_) plt.figure() @@ -321,8 +321,8 @@ plt.title('FISTA-OS-Stidentst-TV reconstruction') plt.show() -# calculate errors +# calculate errors Qtools = QualityTools(phantom_2D[indicesROI], RecFISTA_LS_stud_reg[indicesROI]) RMSE_FISTA_LS_studentst_TV = Qtools.rmse() print("RMSE for FISTA-PWLS-Studentst-TV reconstruction is {}".format(RMSE_FISTA_LS_studentst_TV)) -#%% \ No newline at end of file +#%% diff --git a/Demos/Python/DemoFISTA_artifacts3D.py b/Demos/Python/DemoFISTA_artifacts3D.py index 6f6c6831c..a8bcdafb4 100644 --- a/Demos/Python/DemoFISTA_artifacts3D.py +++ b/Demos/Python/DemoFISTA_artifacts3D.py @@ -50,7 +50,7 @@ # Projection geometry related parameters: Horiz_det = int(N_size) # detector column count (horizontal) -Vert_det = N_size+20 # detector row count (vertical) (no reason for it to be > N) +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) diff --git a/conda-recipe/build.sh b/conda-recipe/build.sh index 84be09e2e..0f2041b20 100755 --- a/conda-recipe/build.sh +++ b/conda-recipe/build.sh @@ -1,4 +1,5 @@ -set -xe +set -xe +cp -rv "$RECIPE_DIR/../test" "$SRC_DIR/" cd $SRC_DIR diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml index 498b87dbf..cee96baae 100755 --- a/conda-recipe/meta.yaml +++ b/conda-recipe/meta.yaml @@ -1,13 +1,20 @@ package: name: tomobar version: {{ environ['VERSION'] }} - + build: preserve_egg_dir: False number: 0 script_env: - VERSION +test: + source_files: + - ./test/ + commands: + - python -c "import os; print (os.getcwd())" + - python -m unittest discover -s test + requirements: build: - python @@ -15,8 +22,8 @@ requirements: - setuptools - cython - cmake - run: + - scipy - python - numpy - libgcc-ng # [unix] diff --git a/src/Python/tomobar/methodsIR.py b/src/Python/tomobar/methodsIR.py index 0bde03b54..ee92d1769 100644 --- a/src/Python/tomobar/methodsIR.py +++ b/src/Python/tomobar/methodsIR.py @@ -22,17 +22,28 @@ import numpy as np from numpy import linalg as LA -import scipy.sparse.linalg try: - from tomobar.supp.addmodules import RING_WEIGHTS -except: - print('____! RING_WEIGHTS C-module failed on import !____') + from ccpi.filters.regularisers import ROF_TV,FGP_TV,PD_TV,SB_TV,LLT_ROF,TGV,NDF,Diff4th,NLTV +except: + print('____! CCPi-regularisation package is missing, please install !____') try: - from ccpi.filters.regularisers import ROF_TV,FGP_TV,PD_TV,SB_TV,LLT_ROF,TGV,NDF,Diff4th,NLTV + import astra +except: + print('____! Astra-toolbox package is missing, please install !____') + +try: + import scipy.sparse.linalg +except: + print('____! Scipy toolbox package is missing, please install !____') + +try: + from tomobar.supp.addmodules import RING_WEIGHTS except: - print('____! CCPi regularisation package is missing, please install !____') + print('____! RING_WEIGHTS C-module failed on import !____') + + def smooth(y, box_pts): # a function to smooth 1D signal @@ -264,7 +275,7 @@ class RecToolsIR: --huber_threshold # threshold for Huber function to apply to data model (supress outliers) --studentst_threshold # threshold for Students't function to apply to data model (supress outliers) --ring_weights_threshold # threshold to produce additional weights to supress ring artifacts - --ring_huber_power # defines the strength of Huber penalty to supress artifacts 1 = Huber, > 1 more penalising + --ring_huber_power # defines the strength of Huber penalty to supress artifacts 1 = Huber, > 1 more penalising --ring_tuple_halfsizes # a tuple for half window sizes as [detector, angles, num of projections] --ringGH_lambda # a parameter for Group-Huber data model to supress full rings of the same intensity --ringGH_accelerate # Group-Huber data model acceleration factor (use carefully to avoid divergence, 50 default) @@ -347,7 +358,7 @@ def __init__(self, from tomobar.supp.astraOP import AstraTools3D self.Atools = AstraTools3D(self.DetectorsDimH, self.DetectorsDimV, self.AnglesVec, self.CenterRotOffset, self.ObjSize) # initiate 3D ASTRA class object return None - + def SIRT(self, _data_, _algorithm_): ###################################################################### @@ -372,7 +383,7 @@ def CGLS(self, _data_, _algorithm_): if (self.geom == '3D'): CGLS_rec = self.Atools.cgls3D(_data_['projection_norm_data'], _algorithm_['iterations']) return CGLS_rec - + def powermethod(self, _data_): # power iteration algorithm to calculate the eigenvalue of the operator (projection matrix) # projection_raw_data is required for PWLS fidelity (self.datafidelity = PWLS), otherwise will be ignored @@ -389,8 +400,8 @@ def powermethod(self, _data_): self.AtoolsOS = AstraToolsOS3D(self.DetectorsDimH, self.DetectorsDimV, self.AnglesVec, self.CenterRotOffset, self.ObjSize, _data_['OS_number']) # initiate 3D ASTRA class OS object niter = 15 # number of power method iterations s = 1.0 - - # classical approach + + # classical approach if (self.geom == '2D'): x1 = np.float32(np.random.randn(self.Atools.ObjSize,self.Atools.ObjSize)) else: @@ -489,9 +500,9 @@ def FISTA(self, _data_, _algorithm_, _regularisation_): r[:,0] = r_x[:,0] - np.multiply(L_const_inv,vec) else: r = r_x - np.multiply(L_const_inv,vec) - + if ((_data_['OS_number'] != 1) and (_data_['ring_weights_threshold'] is not None) and (iter > 0)): - # Ordered subset approach for a better ring model + # Ordered subset approach for a better ring model res_full = self.Atools.forwproj(X_t) - _data_['projection_norm_data'] rings_weights = RING_WEIGHTS(res_full, _data_['ring_tuple_halfsizes'][0], _data_['ring_tuple_halfsizes'][1], _data_['ring_tuple_halfsizes'][2]) ring_function_weight = np.ones(np.shape(res_full)) @@ -703,7 +714,7 @@ def ADMM_Atb(b): # update u variable u = u + (x_hat - z) if (_algorithm_['verbose'] == 'on'): - if (np.mod(iter,(round)(_algorithm_['iterations']/5)+1) == 0): + if (np.mod(iter,(round)(_algorithm_['iterations']/5)+1) == 0): print('ADMM iteration (',iter+1,') using', _regularisation_['method'], 'regularisation for (',(int)(info_vec[0]),') iterations') if (iter == _algorithm_['iterations']-1): print('ADMM stopped at iteration (', iter+1, ')') diff --git a/src/Python/tomobar/supp/astraOP.py b/src/Python/tomobar/supp/astraOP.py index c81545d37..b1b3b5e48 100644 --- a/src/Python/tomobar/supp/astraOP.py +++ b/src/Python/tomobar/supp/astraOP.py @@ -2,17 +2,19 @@ # -*- coding: utf-8 -*- """ A class based on using ASTRA toolbox to perform projection/bakprojection of 2D/3D -data using parallel beam geometry -- SIRT algorithm from ASTRA -- CGLS algorithm from ASTRA +data using parallel beam geometry +- SIRT algorithm from ASTRA +- CGLS algorithm from ASTRA GPLv3 license (ASTRA toolbox) @author: Daniil Kazantsev: https://github.com/dkazanc """ - -import astra import numpy as np +try: + import astra +except: + print('____! Astra-toolbox package is missing, please install !____') #define 2D rotation matrix def rotation_matrix2D(theta): @@ -49,7 +51,7 @@ def vec_geom_init3D(angles_rad, DetectorSpacingX, DetectorSpacingY, CenterRotOff s0 = [0.0, -1.0, 0.0] # source u0 = [DetectorSpacingX, 0.0, 0.0] # detector coordinates v0 = [0.0, 0.0, DetectorSpacingY] # detector coordinates - + vectors = np.zeros([angles_rad.size,12]) for i in range(0,angles_rad.size): if np.ndim(CenterRotOffset) == 0: @@ -91,7 +93,7 @@ def __init__(self, DetectorsDim, AnglesVec, CenterRotOffset, ObjSize, device): print ("Select between 'cpu' or 'gpu' for device") # add optomo operator self.A_optomo = astra.OpTomo(self.proj_id) - + def forwproj(self, image): """Applying forward projection""" sinogram_id, sinogram = astra.create_sino(image, self.proj_id) @@ -105,11 +107,11 @@ def backproj(self, sinogram): astra.data2d.delete(rec_id) return image def fbp2D(self, sinogram): - """perform FBP reconstruction""" + """perform FBP reconstruction""" rec_id = astra.data2d.create( '-vol', self.vol_geom) # Create a data object to hold the sinogram data sinogram_id = astra.data2d.create('-sino', self.proj_geom, sinogram) - + if self.device == 1: cfg = astra.astra_dict('FBP') cfg['ProjectorId'] = self.proj_id @@ -131,11 +133,11 @@ def fbp2D(self, sinogram): astra.data2d.delete(self.proj_id) return recFBP def sirt2D(self, sinogram, iterations): - """perform SIRT reconstruction""" + """perform SIRT reconstruction""" rec_id = astra.data2d.create( '-vol', self.vol_geom) # Create a data object to hold the sinogram data sinogram_id = astra.data2d.create('-sino', self.proj_geom, sinogram) - + if self.device == 1: cfg = astra.astra_dict('SIRT') cfg['ProjectorId'] = self.proj_id @@ -156,11 +158,11 @@ def sirt2D(self, sinogram, iterations): astra.data2d.delete(self.proj_id) return recSIRT def cgls2D(self, sinogram, iterations): - """perform CGLS reconstruction""" + """perform CGLS reconstruction""" rec_id = astra.data2d.create( '-vol', self.vol_geom) # Create a data object to hold the sinogram data sinogram_id = astra.data2d.create('-sino', self.proj_geom, sinogram) - + if self.device == 1: cfg = astra.astra_dict('CGLS') cfg['ProjectorId'] = self.proj_id @@ -180,17 +182,17 @@ def cgls2D(self, sinogram, iterations): astra.data2d.delete(sinogram_id) astra.data2d.delete(self.proj_id) return recCGLS - + class AstraToolsOS: """ - 2D ordered subset parallel beam projection/backprojection class based + 2D ordered subset parallel beam projection/backprojection class based on ASTRA toolbox """ def __init__(self, DetectorsDim, AnglesVec, CenterRotOffset, ObjSize, OS, device): self.DetectorsDim = DetectorsDim self.AnglesVec = AnglesVec self.ObjSize = ObjSize - + ################ arrange ordered-subsets ################ import numpy as np AnglesTot = np.size(AnglesVec) # total number of angles @@ -203,7 +205,7 @@ def __init__(self, DetectorsDim, AnglesVec, CenterRotOffset, ObjSize, OS, device if (indexS < AnglesTot): self.newInd_Vec[sub_ind,proj_ind] = indexS ind_sel += OS - + # create full ASTRA geometry (to calculate Lipshitz constant) vectors = vec_geom_init2D(AnglesVec, 1.0, CenterRotOffset) self.proj_geom = astra.create_proj_geom('parallel_vec', DetectorsDim, vectors) @@ -224,7 +226,7 @@ def __init__(self, DetectorsDim, AnglesVec, CenterRotOffset, ObjSize, OS, device if (self.indVec[self.NumbProjBins-1] == 0): self.indVec = self.indVec[:-1] #shrink vector size anglesOS = self.AnglesVec[self.indVec] # OS-specific angles - + if np.ndim(CenterRotOffset) == 0: # CenterRotOffset is a _scalar_ vectorsOS = vec_geom_init2D(anglesOS, 1.0, CenterRotOffset) #self.proj_geom_OS[sub_ind] = astra.create_proj_geom('parallel', 1.0, DetectorsDim, anglesOS) @@ -279,7 +281,7 @@ def __init__(self, DetColumnCount, DetRowCount, AnglesVec, CenterRotOffset, ObjS self.vol_geom = astra.create_vol_geom(Y,X,Z) self.proj_id = astra.create_projector('cuda3d', self.proj_geom, self.vol_geom) # for GPU self.A_optomo = astra.OpTomo(self.proj_id) - + def forwproj(self, object3D): """Applying forward projection""" proj_id, proj_data = astra.create_sino3d_gpu(object3D, self.proj_geom, self.vol_geom) @@ -291,7 +293,7 @@ def backproj(self, proj_data): astra.data3d.delete(rec_id) return object3D def sirt3D(self, sinogram, iterations): - """perform SIRT reconstruction""" + """perform SIRT reconstruction""" sinogram_id = astra.data3d.create("-sino", self.proj_geom, sinogram) # Create a data object for the reconstruction rec_id = astra.data3d.create('-vol', self.vol_geom) @@ -304,7 +306,7 @@ def sirt3D(self, sinogram, iterations): alg_id = astra.algorithm.create(cfg) # This will have a runtime in the order of 10 seconds. astra.algorithm.run(alg_id, iterations) - + # Get the result recSIRT = astra.data3d.get(rec_id) @@ -313,7 +315,7 @@ def sirt3D(self, sinogram, iterations): astra.data3d.delete(sinogram_id) return recSIRT def cgls3D(self, sinogram, iterations): - """perform CGLS reconstruction""" + """perform CGLS reconstruction""" sinogram_id = astra.data3d.create("-sino", self.proj_geom, sinogram) # Create a data object for the reconstruction rec_id = astra.data3d.create('-vol', self.vol_geom) @@ -326,7 +328,7 @@ def cgls3D(self, sinogram, iterations): alg_id = astra.algorithm.create(cfg) # This will have a runtime in the order of 10 seconds. astra.algorithm.run(alg_id, iterations) - + # Get the result recCGLS = astra.data3d.get(rec_id) @@ -337,7 +339,7 @@ def cgls3D(self, sinogram, iterations): class AstraToolsOS3D: """ - 3D ordered subset parallel beam projection/backprojection class based + 3D ordered subset parallel beam projection/backprojection class based on ASTRA toolbox """ def __init__(self, DetColumnCount, DetRowCount, AnglesVec, CenterRotOffset, ObjSize, OS): @@ -349,7 +351,7 @@ def __init__(self, DetColumnCount, DetRowCount, AnglesVec, CenterRotOffset, ObjS Y=X=ObjSize Z=DetRowCount self.vol_geom = astra.create_vol_geom(Y,X,Z) - + ################ arrange ordered-subsets ################ import numpy as np AnglesTot = np.size(AnglesVec) # total number of angles @@ -362,7 +364,7 @@ def __init__(self, DetColumnCount, DetRowCount, AnglesVec, CenterRotOffset, ObjS if (indexS < AnglesTot): self.newInd_Vec[sub_ind,proj_ind] = indexS ind_sel += OS - + # create full ASTRA geometry (to calculate Lipshitz constant) vectors = vec_geom_init3D(AnglesVec, 1.0, 1.0, CenterRotOffset) self.proj_geom = astra.create_proj_geom('parallel3d_vec', DetRowCount, DetColumnCount, vectors) @@ -373,7 +375,7 @@ def __init__(self, DetColumnCount, DetRowCount, AnglesVec, CenterRotOffset, ObjS if (self.indVec[self.NumbProjBins-1] == 0): self.indVec = self.indVec[:-1] #shrink vector size anglesOS = AnglesVec[self.indVec] # OS-specific angles - + if np.ndim(CenterRotOffset) == 0: # CenterRotOffset is a _scalar_ vectors = vec_geom_init3D(anglesOS, 1.0, 1.0, CenterRotOffset) else: # CenterRotOffset is a _vector_ @@ -398,4 +400,4 @@ def backprojOS(self, proj_data, no_os): """Applying back-projection to a specific subset""" rec_id, object3D = astra.create_backprojection3d_gpu(proj_data, self.proj_geom_OS[no_os], self.vol_geom) astra.data3d.delete(rec_id) - return object3D \ No newline at end of file + return object3D diff --git a/test/__init__.py b/test/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/test/test_tomobarCPU_DIR.py b/test/test_tomobarCPU_DIR.py new file mode 100644 index 000000000..4b4977fc9 --- /dev/null +++ b/test/test_tomobarCPU_DIR.py @@ -0,0 +1,39 @@ +import unittest +import numpy as np +from tomobar.methodsDIR import RecToolsDIR +from tomobar.methodsIR import RecToolsIR + +############################################################################### +class TestTomobar(unittest.TestCase): + + def test_tomobarDIR(self): + model = 4 # select a model + N_size = 64 # set dimension of the phantom + + # create sinogram analytically + angles_num = int(0.5*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 = int(np.sqrt(2)*N_size) # detectors + sino_num = np.ones((P, angles_num)) + + RectoolsDirect = RecToolsDIR(DetectorsDimH = P, # DetectorsDimH # detector dimension (horizontal) + DetectorsDimV = None, # DetectorsDimV # detector dimension (vertical) for 3D case only + CenterRotOffset = 0.0, # Center of Rotation (CoR) scalar + AnglesVec = angles_rad, # array of angles in radians + ObjSize = N_size, # a scalar to define reconstructed object dimensions + device_projector = 'cpu') + """ + RectoolsIterative = RecToolsIR(DetectorsDimH = P, # DetectorsDimH # detector dimension (horizontal) + DetectorsDimV = None, # DetectorsDimV # detector dimension (vertical) for 3D case only + CenterRotOffset = 0.0, # Center of Rotation (CoR) scalar + AnglesVec = angles_rad, # array of angles in radians + ObjSize = N_size, # a scalar to define reconstructed object dimensions + datafidelity = 'LS', + device_projector = 'cpu') + """ + RecFourier = RectoolsDirect.FOURIER(sino_num,'linear') + +############################################################################### +if __name__ == '__main__': + unittest.main() diff --git a/tests/test.py b/tests/test.py deleted file mode 100644 index cc56ac472..000000000 --- a/tests/test.py +++ /dev/null @@ -1,9 +0,0 @@ -import unittest -import numpy as np -from tomobar.methodsDIR import RecToolsDIR -from tomobar.methodsIR import RecToolsIR - -############################################################################### - -if __name__ == '__main__': - unittest.main()