Adding tests with pytest

.github/workflows/build.yml

@@ -39,7 +39,7 @@ jobs:
     - name: test
       run: |
         conda activate "${{ steps.reqs.outputs.envname }}"
-        PYTHONPATH=./src/Python python -m unittest discover ./test
+        PYTHONPATH=./src/Python python -m pytest -s ./test
     - if: always()
       name: Post Run conda-incubator/setup-miniconda@v3
       shell: bash
@@ -71,7 +71,7 @@ jobs:
         cmake --build ./build_proj --target install
         pip install ./src/Python
     - name: test
-      run: PYTHONPATH=./src/Python python -m unittest discover ./test
+      run: PYTHONPATH=./src/Python python -m pytest -s ./test
   conda:
     defaults: {run: {shell: 'bash -el {0}'}}
     runs-on: ubuntu-latest

Installation.md

@@ -56,7 +56,7 @@ conda install ccpi-regulariser -c ccpi -c conda-forge
 #### Python (conda-build)
 
 ```sh
-conda build recipe/ --numpy 1.23 --python 3.10
+conda build recipe/ --numpy 1.26 --python 3.10 -c conda-forge
 conda install ccpi-regulariser --use-local --force-reinstall # doesn't work?
 conda install -c file://${CONDA_PREFIX}/conda-bld/ ccpi-regulariser --force-reinstall # try this one
 cd demos/

demos/Matlab_demos/demoMatlab_3Ddenoise.m

@@ -14,7 +14,7 @@
 slices = 15;
 vol3D = zeros(N,N,slices, 'single');
 Ideal3D = zeros(N,N,slices, 'single');
-Im = double(imread('lena_gray_512.tif'))/255;  % loading image
+Im = double(imread('peppers.tif'))/255;  % loading image
 for i = 1:slices
 vol3D(:,:,i) = Im + .05*randn(size(Im)); 
 Ideal3D(:,:,i) = Im;

demos/Matlab_demos/demoMatlab_denoise.m

@@ -9,7 +9,7 @@
 addpath(Path2);
 addpath(Path3);
 
-Im = double(imread('lena_gray_512.tif'))/255;  % loading image
+Im = double(imread('peppers.tif'))/255;  % loading image
 u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
 figure; imshow(u0, [0 1]); title('Noisy image');
 %%

demos/SoftwareX_supp/Demo_RealData_Recon_SX.py

@@ -28,62 +28,66 @@
 import time
 
 # load dendritic projection data
-h5f = h5py.File('data/DendrData_3D.h5','r')
-dataRaw = h5f['dataRaw'][:]
-flats = h5f['flats'][:]
-darks = h5f['darks'][:]
-angles_rad = h5f['angles_rad'][:]
+h5f = h5py.File("data/DendrData_3D.h5", "r")
+dataRaw = h5f["dataRaw"][:]
+flats = h5f["flats"][:]
+darks = h5f["darks"][:]
+angles_rad = h5f["angles_rad"][:]
 h5f.close()
-#%%
+# %%
 # normalise the data [detectorsVert, Projections, detectorsHoriz]
-data_norm = normaliser(dataRaw, flats, darks, log='log')
+data_norm = normaliser(dataRaw, flats, darks, log="log")
 del dataRaw, darks, flats
 
 intens_max = 2.3
 plt.figure()
 plt.subplot(131)
-plt.imshow(data_norm[:,150,:],vmin=0, vmax=intens_max)
-plt.title('2D Projection (analytical)')
+plt.imshow(data_norm[:, 150, :], vmin=0, vmax=intens_max)
+plt.title("2D Projection (analytical)")
 plt.subplot(132)
-plt.imshow(data_norm[300,:,:],vmin=0, vmax=intens_max)
-plt.title('Sinogram view')
+plt.imshow(data_norm[300, :, :], vmin=0, vmax=intens_max)
+plt.title("Sinogram view")
 plt.subplot(133)
-plt.imshow(data_norm[:,:,600],vmin=0, vmax=intens_max)
-plt.title('Tangentogram view')
+plt.imshow(data_norm[:, :, 600], vmin=0, vmax=intens_max)
+plt.title("Tangentogram view")
 plt.show()
 
-detectorHoriz = np.size(data_norm,2)
-det_y_crop = [i for i in range(0,detectorHoriz-22)]
-N_size = 950 # reconstruction domain
+detectorHoriz = np.size(data_norm, 2)
+det_y_crop = [i for i in range(0, detectorHoriz - 22)]
+N_size = 950  # reconstruction domain
 time_label = int(time.time())
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("%%%%%%%%%%%%Reconstructing with FBP method %%%%%%%%%%%%%%%%%")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+# %%
+print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+print("%%%%%%%%%%%%Reconstructing with FBP method %%%%%%%%%%%%%%%%%")
+print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
 from tomobar.methodsDIR import RecToolsDIR
 
-RectoolsDIR = RecToolsDIR(DetectorsDimH = np.size(det_y_crop),  # DetectorsDimH # detector dimension (horizontal)
-                    DetectorsDimV = 100,  # DetectorsDimV # detector dimension (vertical) for 3D case only
-                    AnglesVec = angles_rad, # array of angles in radians
-                    ObjSize = N_size, # a scalar to define reconstructed object dimensions
-                    device='gpu')
+RectoolsDIR = RecToolsDIR(
+    DetectorsDimH=np.size(
+        det_y_crop
+    ),  # DetectorsDimH # detector dimension (horizontal)
+    DetectorsDimV=100,  # DetectorsDimV # detector dimension (vertical) for 3D case only
+    AnglesVec=angles_rad,  # array of angles in radians
+    ObjSize=N_size,  # a scalar to define reconstructed object dimensions
+    device="gpu",
+)
 
-FBPrec = RectoolsDIR.FBP(data_norm[0:100,:,det_y_crop])
+FBPrec = RectoolsDIR.FBP(data_norm[0:100, :, det_y_crop])
 
 sliceSel = 50
 max_val = 0.003
 plt.figure()
 plt.subplot(131)
-plt.imshow(FBPrec[sliceSel,:,:],vmin=0, vmax=max_val, cmap="gray")
-plt.title('FBP Reconstruction, axial view')
+plt.imshow(FBPrec[sliceSel, :, :], vmin=0, vmax=max_val, cmap="gray")
+plt.title("FBP Reconstruction, axial view")
 
 plt.subplot(132)
-plt.imshow(FBPrec[:,sliceSel,:],vmin=0, vmax=max_val, cmap="gray")
-plt.title('FBP Reconstruction, coronal view')
+plt.imshow(FBPrec[:, sliceSel, :], vmin=0, vmax=max_val, cmap="gray")
+plt.title("FBP Reconstruction, coronal view")
 
 plt.subplot(133)
-plt.imshow(FBPrec[:,:,sliceSel],vmin=0, vmax=max_val, cmap="gray")
-plt.title('FBP Reconstruction, sagittal view')
+plt.imshow(FBPrec[:, :, sliceSel], vmin=0, vmax=max_val, cmap="gray")
+plt.title("FBP Reconstruction, sagittal view")
 plt.show()
 
 # saving to tiffs (16bit)
@@ -98,44 +102,51 @@
     tiff.write_image(np.uint16(FBPrec[i,:,:]*multiplier))
     tiff.close()
 """
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("Reconstructing with ADMM method using tomobar software")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+# %%
+print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+print("Reconstructing with ADMM method using tomobar software")
+print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
 # initialise tomobar ITERATIVE reconstruction class ONCE
 from tomobar.methodsIR import RecToolsIR
-RectoolsIR = RecToolsIR(DetectorsDimH =  np.size(det_y_crop),  # DetectorsDimH # detector dimension (horizontal)
-                    DetectorsDimV = 100,  # DetectorsDimV # detector dimension (vertical) 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, GH (wip), Students t (wip)
-                    nonnegativity='ENABLE', # enable nonnegativity constraint (set to 'ENABLE')
-                    OS_number = None, # the number of subsets, NONE/(or > 1) ~ classical / ordered subsets
-                    tolerance = 0.0, # tolerance to stop inner (regularisation) iterations earlier
-                    device='gpu')
-#%%
-print ("Reconstructing with ADMM method using SB-TV penalty")
-RecADMM_reg_sbtv = RectoolsIR.ADMM(data_norm[0:100,:,det_y_crop],
-                              rho_const = 2000.0, \
-                              iterationsADMM = 15, \
-                              regularisation = 'SB_TV', \
-                              regularisation_parameter = 0.00085,\
-                              regularisation_iterations = 50)
+
+RectoolsIR = RecToolsIR(
+    DetectorsDimH=np.size(
+        det_y_crop
+    ),  # DetectorsDimH # detector dimension (horizontal)
+    DetectorsDimV=100,  # DetectorsDimV # detector dimension (vertical) 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, GH (wip), Students t (wip)
+    nonnegativity="ENABLE",  # enable nonnegativity constraint (set to 'ENABLE')
+    OS_number=None,  # the number of subsets, NONE/(or > 1) ~ classical / ordered subsets
+    tolerance=0.0,  # tolerance to stop inner (regularisation) iterations earlier
+    device="gpu",
+)
+# %%
+print("Reconstructing with ADMM method using SB-TV penalty")
+RecADMM_reg_sbtv = RectoolsIR.ADMM(
+    data_norm[0:100, :, det_y_crop],
+    rho_const=2000.0,
+    iterationsADMM=15,
+    regularisation="SB_TV",
+    regularisation_parameter=0.00085,
+    regularisation_iterations=50,
+)
 
 sliceSel = 50
 max_val = 0.003
 plt.figure()
 plt.subplot(131)
-plt.imshow(RecADMM_reg_sbtv[sliceSel,:,:],vmin=0, vmax=max_val, cmap="gray")
-plt.title('3D ADMM-SB-TV Reconstruction, axial view')
+plt.imshow(RecADMM_reg_sbtv[sliceSel, :, :], vmin=0, vmax=max_val, cmap="gray")
+plt.title("3D ADMM-SB-TV Reconstruction, axial view")
 
 plt.subplot(132)
-plt.imshow(RecADMM_reg_sbtv[:,sliceSel,:],vmin=0, vmax=max_val, cmap="gray")
-plt.title('3D ADMM-SB-TV Reconstruction, coronal view')
+plt.imshow(RecADMM_reg_sbtv[:, sliceSel, :], vmin=0, vmax=max_val, cmap="gray")
+plt.title("3D ADMM-SB-TV Reconstruction, coronal view")
 
 plt.subplot(133)
-plt.imshow(RecADMM_reg_sbtv[:,:,sliceSel],vmin=0, vmax=max_val, cmap="gray")
-plt.title('3D ADMM-SB-TV Reconstruction, sagittal view')
+plt.imshow(RecADMM_reg_sbtv[:, :, sliceSel], vmin=0, vmax=max_val, cmap="gray")
+plt.title("3D ADMM-SB-TV Reconstruction, sagittal view")
 plt.show()
 
 
@@ -149,33 +160,35 @@
     tiff.close()
 """
 # Saving recpnstructed data with a unique time label
-np.save('Dendr_ADMM_SBTV'+str(time_label)+'.npy', RecADMM_reg_sbtv)
+np.save("Dendr_ADMM_SBTV" + str(time_label) + ".npy", RecADMM_reg_sbtv)
 del RecADMM_reg_sbtv
-#%%
-print ("Reconstructing with ADMM method using ROF-LLT penalty")
-RecADMM_reg_rofllt = RectoolsIR.ADMM(data_norm[0:100,:,det_y_crop],
-                              rho_const = 2000.0, \
-                              iterationsADMM = 15, \
-                              regularisation = 'LLT_ROF', \
-                              regularisation_parameter = 0.0009,\
-                              regularisation_parameter2 = 0.0007,\
-                              time_marching_parameter = 0.001,\
-                              regularisation_iterations = 550)
+# %%
+print("Reconstructing with ADMM method using ROF-LLT penalty")
+RecADMM_reg_rofllt = RectoolsIR.ADMM(
+    data_norm[0:100, :, det_y_crop],
+    rho_const=2000.0,
+    iterationsADMM=15,
+    regularisation="LLT_ROF",
+    regularisation_parameter=0.0009,
+    regularisation_parameter2=0.0007,
+    time_marching_parameter=0.001,
+    regularisation_iterations=550,
+)
 
 sliceSel = 50
 max_val = 0.003
 plt.figure()
 plt.subplot(131)
-plt.imshow(RecADMM_reg_rofllt[sliceSel,:,:],vmin=0, vmax=max_val)
-plt.title('3D ADMM-ROFLLT Reconstruction, axial view')
+plt.imshow(RecADMM_reg_rofllt[sliceSel, :, :], vmin=0, vmax=max_val)
+plt.title("3D ADMM-ROFLLT Reconstruction, axial view")
 
 plt.subplot(132)
-plt.imshow(RecADMM_reg_rofllt[:,sliceSel,:],vmin=0, vmax=max_val)
-plt.title('3D ADMM-ROFLLT Reconstruction, coronal view')
+plt.imshow(RecADMM_reg_rofllt[:, sliceSel, :], vmin=0, vmax=max_val)
+plt.title("3D ADMM-ROFLLT Reconstruction, coronal view")
 
 plt.subplot(133)
-plt.imshow(RecADMM_reg_rofllt[:,:,sliceSel],vmin=0, vmax=max_val)
-plt.title('3D ADMM-ROFLLT Reconstruction, sagittal view')
+plt.imshow(RecADMM_reg_rofllt[:, :, sliceSel], vmin=0, vmax=max_val)
+plt.title("3D ADMM-ROFLLT Reconstruction, sagittal view")
 plt.show()
 
 # saving to tiffs (16bit)
@@ -189,31 +202,33 @@
 """
 
 # Saving recpnstructed data with a unique time label
-np.save('Dendr_ADMM_ROFLLT'+str(time_label)+'.npy', RecADMM_reg_rofllt)
+np.save("Dendr_ADMM_ROFLLT" + str(time_label) + ".npy", RecADMM_reg_rofllt)
 del RecADMM_reg_rofllt
-#%%
-print ("Reconstructing with ADMM method using TGV penalty")
-RecADMM_reg_tgv = RectoolsIR.ADMM(data_norm[0:100,:,det_y_crop],
-                              rho_const = 2000.0, \
-                              iterationsADMM = 15, \
-                              regularisation = 'TGV', \
-                              regularisation_parameter = 0.01,\
-                              regularisation_iterations = 500)
+# %%
+print("Reconstructing with ADMM method using TGV penalty")
+RecADMM_reg_tgv = RectoolsIR.ADMM(
+    data_norm[0:100, :, det_y_crop],
+    rho_const=2000.0,
+    iterationsADMM=15,
+    regularisation="TGV",
+    regularisation_parameter=0.01,
+    regularisation_iterations=500,
+)
 
 sliceSel = 50
 max_val = 0.003
 plt.figure()
 plt.subplot(131)
-plt.imshow(RecADMM_reg_tgv[sliceSel,:,:],vmin=0, vmax=max_val)
-plt.title('3D ADMM-TGV Reconstruction, axial view')
+plt.imshow(RecADMM_reg_tgv[sliceSel, :, :], vmin=0, vmax=max_val)
+plt.title("3D ADMM-TGV Reconstruction, axial view")
 
 plt.subplot(132)
-plt.imshow(RecADMM_reg_tgv[:,sliceSel,:],vmin=0, vmax=max_val)
-plt.title('3D ADMM-TGV Reconstruction, coronal view')
+plt.imshow(RecADMM_reg_tgv[:, sliceSel, :], vmin=0, vmax=max_val)
+plt.title("3D ADMM-TGV Reconstruction, coronal view")
 
 plt.subplot(133)
-plt.imshow(RecADMM_reg_tgv[:,:,sliceSel],vmin=0, vmax=max_val)
-plt.title('3D ADMM-TGV Reconstruction, sagittal view')
+plt.imshow(RecADMM_reg_tgv[:, :, sliceSel], vmin=0, vmax=max_val)
+plt.title("3D ADMM-TGV Reconstruction, sagittal view")
 plt.show()
 
 # saving to tiffs (16bit)
@@ -226,6 +241,6 @@
     tiff.close()
 """
 # Saving recpnstructed data with a unique time label
-np.save('Dendr_ADMM_TGV'+str(time_label)+'.npy', RecADMM_reg_tgv)
+np.save("Dendr_ADMM_TGV" + str(time_label) + ".npy", RecADMM_reg_tgv)
 del RecADMM_reg_tgv
-#%%
+# %%