From 62b30291105e8a48633629350fc9820b404da2ff Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Thu, 17 Aug 2017 16:33:09 +0100
Subject: initial revision

---
 src/Python/test/astra_test.py        | 85 ++++++++++++++++++++++++++++++++++++
 src/Python/test/simple_astra_test.py | 25 +++++++++++
 2 files changed, 110 insertions(+)
 create mode 100644 src/Python/test/astra_test.py
 create mode 100644 src/Python/test/simple_astra_test.py

(limited to 'src/Python/test')

diff --git a/src/Python/test/astra_test.py b/src/Python/test/astra_test.py
new file mode 100644
index 0000000..42c375a
--- /dev/null
+++ b/src/Python/test/astra_test.py
@@ -0,0 +1,85 @@
+import astra
+import numpy
+import filefun
+
+
+# read in the same data as the DemoRD2
+angles = filefun.dlmread("DemoRD2/angles.csv")
+darks_ar = filefun.dlmread("DemoRD2/darks_ar.csv", separator=",")
+flats_ar = filefun.dlmread("DemoRD2/flats_ar.csv", separator=",")
+
+if True:
+    Sino3D = numpy.load("DemoRD2/Sino3D.npy")
+else:
+    sino = filefun.dlmread("DemoRD2/sino_01.csv", separator=",")
+    a = map (lambda x:x, numpy.shape(sino))
+    a.append(20)
+
+    Sino3D = numpy.zeros(tuple(a), dtype="float")
+
+    for i in range(1,numpy.shape(Sino3D)[2]+1):
+        print("Read file DemoRD2/sino_%02d.csv" % i)
+        sino = filefun.dlmread("DemoRD2/sino_%02d.csv" % i, separator=",")
+        Sino3D.T[i-1] = sino.T
+    
+Weights3D = numpy.asarray(Sino3D, dtype="float")
+
+##angles_rad = angles*(pi/180); % conversion to radians
+##size_det = size(data_raw3D,1); % detectors dim
+##angSize = size(data_raw3D, 2); % angles dim
+##slices_tot = size(data_raw3D, 3); % no of slices
+##recon_size = 950; % reconstruction size
+
+
+angles_rad = angles * numpy.pi /180.
+size_det, angSize, slices_tot = numpy.shape(Sino3D)
+size_det, angSize, slices_tot = [int(i) for i in numpy.shape(Sino3D)]
+recon_size = 950
+Z_slices = 3;
+det_row_count = Z_slices;
+
+#proj_geom = astra_create_proj_geom('parallel3d', 1, 1,
+# det_row_count, size_det, angles_rad);
+
+detectorSpacingX = 1.0
+detectorSpacingY = detectorSpacingX
+proj_geom = astra.create_proj_geom('parallel3d',
+                                            detectorSpacingX,
+                                            detectorSpacingY,
+                                            det_row_count,
+                                            size_det,
+                                            angles_rad)
+
+#vol_geom = astra_create_vol_geom(recon_size,recon_size,Z_slices);
+vol_geom = astra.create_vol_geom(recon_size,recon_size,Z_slices);
+
+sino = numpy.zeros((size_det, angSize, slices_tot), dtype="float")
+
+#weights = ones(size(sino));
+weights = numpy.ones(numpy.shape(sino))
+
+#####################################################################
+## PowerMethod for Lipschitz constant
+
+N = vol_geom['GridColCount']
+x1 = numpy.random.rand(1,N,N)
+#sqweight = sqrt(weights(:,:,1));
+sqweight = numpy.sqrt(weights.T[0]).T
+##proj_geomT = proj_geom;
+proj_geomT = proj_geom.copy()
+##proj_geomT.DetectorRowCount = 1;
+proj_geomT['DetectorRowCount'] = 1
+##vol_geomT = vol_geom;
+vol_geomT = vol_geom.copy()
+##vol_geomT.GridSliceCount = 1;
+vol_geomT['GridSliceCount'] = 1
+
+##[sino_id, y] = astra_create_sino3d_cuda(x1, proj_geomT, vol_geomT);
+
+#sino_id, y = astra.create_sino3d_gpu(x1, proj_geomT, vol_geomT);
+sino_id, y = astra.create_sino(x1, proj_geomT, vol_geomT);
+
+##y = sqweight.*y;
+##astra_mex_data3d('delete', sino_id);
+        
+
diff --git a/src/Python/test/simple_astra_test.py b/src/Python/test/simple_astra_test.py
new file mode 100644
index 0000000..905eeea
--- /dev/null
+++ b/src/Python/test/simple_astra_test.py
@@ -0,0 +1,25 @@
+import astra
+import numpy
+
+detectorSpacingX = 1.0
+detectorSpacingY = 1.0
+det_row_count = 128
+det_col_count = 128
+
+angles_rad = numpy.asarray([i for i in range(360)], dtype=float) / 180. * numpy.pi
+
+proj_geom = astra.creators.create_proj_geom('parallel3d',
+                                            detectorSpacingX,
+                                            detectorSpacingY,
+                                            det_row_count,
+                                            det_col_count,
+                                            angles_rad)
+
+image_size_x = 64
+image_size_y = 64
+image_size_z = 32
+
+vol_geom = astra.creators.create_vol_geom(image_size_x,image_size_y,image_size_z)
+
+x1 = numpy.random.rand(image_size_z,image_size_y,image_size_x)
+sino_id, y = astra.creators.create_sino3d_gpu(x1, proj_geom, vol_geom)
-- 
cgit v1.2.3


From cc3a464ec587e95ddfd421cd3836a7677dfb9744 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Wed, 23 Aug 2017 16:51:18 +0100
Subject: export/import data from hdf5

Added file to export the data from DemoRD2.m to HDF5 to pass it to Python.
Added file to import the data from DemoRD2.m from HDF5.
---
 src/Python/test/readhd5.py | 28 ++++++++++++++++++++++++++++
 1 file changed, 28 insertions(+)
 create mode 100644 src/Python/test/readhd5.py

(limited to 'src/Python/test')

diff --git a/src/Python/test/readhd5.py b/src/Python/test/readhd5.py
new file mode 100644
index 0000000..1e19e14
--- /dev/null
+++ b/src/Python/test/readhd5.py
@@ -0,0 +1,28 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Aug 23 16:34:49 2017
+
+@author: ofn77899
+"""
+
+import h5py
+import numpy
+
+def getEntry(nx, location):
+    for item in nx[location].keys():
+        print (item)
+        
+filename = r'C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\Demos\DendrData.h5'
+nx = h5py.File(filename, "r")
+#getEntry(nx, '/')
+# I have exported the entries as children of /
+entries = [entry for entry in nx['/'].keys()]
+print (entries)
+
+Sino3D = numpy.asarray(nx.get('/Sino3D'))
+Weights3D = numpy.asarray(nx.get('/Weights3D'))
+angSize = numpy.asarray(nx.get('/angSize'), dtype=int)[0]
+angles_rad = numpy.asarray(nx.get('/angles_rad'))
+recon_size = numpy.asarray(nx.get('/recon_size'), dtype=int)[0]
+size_det = numpy.asarray(nx.get('/size_det'), dtype=int)[0]
+slices_tot = numpy.asarray(nx.get('/slices_tot'), dtype=int)[0]
\ No newline at end of file
-- 
cgit v1.2.3


From 3fffd568589137b17d1fbe44e55a757e3745a3b1 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Wed, 11 Oct 2017 15:42:05 +0100
Subject: added simple_astra_test.py

---
 src/Python/test/simple_astra_test.py | 25 +++++++++++++++++++++++++
 1 file changed, 25 insertions(+)
 create mode 100644 src/Python/test/simple_astra_test.py

(limited to 'src/Python/test')

diff --git a/src/Python/test/simple_astra_test.py b/src/Python/test/simple_astra_test.py
new file mode 100644
index 0000000..905eeea
--- /dev/null
+++ b/src/Python/test/simple_astra_test.py
@@ -0,0 +1,25 @@
+import astra
+import numpy
+
+detectorSpacingX = 1.0
+detectorSpacingY = 1.0
+det_row_count = 128
+det_col_count = 128
+
+angles_rad = numpy.asarray([i for i in range(360)], dtype=float) / 180. * numpy.pi
+
+proj_geom = astra.creators.create_proj_geom('parallel3d',
+                                            detectorSpacingX,
+                                            detectorSpacingY,
+                                            det_row_count,
+                                            det_col_count,
+                                            angles_rad)
+
+image_size_x = 64
+image_size_y = 64
+image_size_z = 32
+
+vol_geom = astra.creators.create_vol_geom(image_size_x,image_size_y,image_size_z)
+
+x1 = numpy.random.rand(image_size_z,image_size_y,image_size_x)
+sino_id, y = astra.creators.create_sino3d_gpu(x1, proj_geom, vol_geom)
-- 
cgit v1.2.3


From 0611d34c31fa1e706c3bcd7e17651f7555469e00 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Thu, 17 Aug 2017 16:33:09 +0100
Subject: initial revision

---
 src/Python/test/simple_astra_test.py | 25 -------------------------
 1 file changed, 25 deletions(-)
 delete mode 100644 src/Python/test/simple_astra_test.py

(limited to 'src/Python/test')

diff --git a/src/Python/test/simple_astra_test.py b/src/Python/test/simple_astra_test.py
deleted file mode 100644
index 905eeea..0000000
--- a/src/Python/test/simple_astra_test.py
+++ /dev/null
@@ -1,25 +0,0 @@
-import astra
-import numpy
-
-detectorSpacingX = 1.0
-detectorSpacingY = 1.0
-det_row_count = 128
-det_col_count = 128
-
-angles_rad = numpy.asarray([i for i in range(360)], dtype=float) / 180. * numpy.pi
-
-proj_geom = astra.creators.create_proj_geom('parallel3d',
-                                            detectorSpacingX,
-                                            detectorSpacingY,
-                                            det_row_count,
-                                            det_col_count,
-                                            angles_rad)
-
-image_size_x = 64
-image_size_y = 64
-image_size_z = 32
-
-vol_geom = astra.creators.create_vol_geom(image_size_x,image_size_y,image_size_z)
-
-x1 = numpy.random.rand(image_size_z,image_size_y,image_size_x)
-sino_id, y = astra.creators.create_sino3d_gpu(x1, proj_geom, vol_geom)
-- 
cgit v1.2.3


From a203949c84484fe2641e39451f033d20d445b1f3 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Wed, 23 Aug 2017 16:51:18 +0100
Subject: export/import data from hdf5

Added file to export the data from DemoRD2.m to HDF5 to pass it to Python.
Added file to import the data from DemoRD2.m from HDF5.
---
 src/Python/test/readhd5.py | 28 ++++++++++++++++++++++++++++
 1 file changed, 28 insertions(+)
 create mode 100644 src/Python/test/readhd5.py

(limited to 'src/Python/test')

diff --git a/src/Python/test/readhd5.py b/src/Python/test/readhd5.py
new file mode 100644
index 0000000..1e19e14
--- /dev/null
+++ b/src/Python/test/readhd5.py
@@ -0,0 +1,28 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Aug 23 16:34:49 2017
+
+@author: ofn77899
+"""
+
+import h5py
+import numpy
+
+def getEntry(nx, location):
+    for item in nx[location].keys():
+        print (item)
+        
+filename = r'C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\Demos\DendrData.h5'
+nx = h5py.File(filename, "r")
+#getEntry(nx, '/')
+# I have exported the entries as children of /
+entries = [entry for entry in nx['/'].keys()]
+print (entries)
+
+Sino3D = numpy.asarray(nx.get('/Sino3D'))
+Weights3D = numpy.asarray(nx.get('/Weights3D'))
+angSize = numpy.asarray(nx.get('/angSize'), dtype=int)[0]
+angles_rad = numpy.asarray(nx.get('/angles_rad'))
+recon_size = numpy.asarray(nx.get('/recon_size'), dtype=int)[0]
+size_det = numpy.asarray(nx.get('/size_det'), dtype=int)[0]
+slices_tot = numpy.asarray(nx.get('/slices_tot'), dtype=int)[0]
\ No newline at end of file
-- 
cgit v1.2.3


From 776070e22bf95491275a023f3a5ac00cea356714 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Wed, 11 Oct 2017 15:12:42 +0100
Subject: read and plot the hdf5

---
 src/Python/test/readhd5.py | 15 ++++++++++++++-
 1 file changed, 14 insertions(+), 1 deletion(-)

(limited to 'src/Python/test')

diff --git a/src/Python/test/readhd5.py b/src/Python/test/readhd5.py
index 1e19e14..b042341 100644
--- a/src/Python/test/readhd5.py
+++ b/src/Python/test/readhd5.py
@@ -25,4 +25,17 @@ angSize = numpy.asarray(nx.get('/angSize'), dtype=int)[0]
 angles_rad = numpy.asarray(nx.get('/angles_rad'))
 recon_size = numpy.asarray(nx.get('/recon_size'), dtype=int)[0]
 size_det = numpy.asarray(nx.get('/size_det'), dtype=int)[0]
-slices_tot = numpy.asarray(nx.get('/slices_tot'), dtype=int)[0]
\ No newline at end of file
+slices_tot = numpy.asarray(nx.get('/slices_tot'), dtype=int)[0]
+
+#from ccpi.viewer.CILViewer2D import CILViewer2D
+#v = CILViewer2D()
+#v.setInputAsNumpy(Weights3D)
+#v.startRenderLoop()
+
+import matplotlib.pyplot as plt
+fig = plt.figure()
+
+a=fig.add_subplot(1,1,1)
+a.set_title('noise')
+imgplot = plt.imshow(Weights3D[0].T)
+plt.show()
-- 
cgit v1.2.3


From 5c978b706192bc5885c7e5001a4bc4626f63d29f Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Wed, 11 Oct 2017 15:49:18 +0100
Subject: initial revision

---
 src/Python/test/simple_astra_test.py | 25 +++++++++++++++++++++++++
 1 file changed, 25 insertions(+)
 create mode 100644 src/Python/test/simple_astra_test.py

(limited to 'src/Python/test')

diff --git a/src/Python/test/simple_astra_test.py b/src/Python/test/simple_astra_test.py
new file mode 100644
index 0000000..905eeea
--- /dev/null
+++ b/src/Python/test/simple_astra_test.py
@@ -0,0 +1,25 @@
+import astra
+import numpy
+
+detectorSpacingX = 1.0
+detectorSpacingY = 1.0
+det_row_count = 128
+det_col_count = 128
+
+angles_rad = numpy.asarray([i for i in range(360)], dtype=float) / 180. * numpy.pi
+
+proj_geom = astra.creators.create_proj_geom('parallel3d',
+                                            detectorSpacingX,
+                                            detectorSpacingY,
+                                            det_row_count,
+                                            det_col_count,
+                                            angles_rad)
+
+image_size_x = 64
+image_size_y = 64
+image_size_z = 32
+
+vol_geom = astra.creators.create_vol_geom(image_size_x,image_size_y,image_size_z)
+
+x1 = numpy.random.rand(image_size_z,image_size_y,image_size_x)
+sino_id, y = astra.creators.create_sino3d_gpu(x1, proj_geom, vol_geom)
-- 
cgit v1.2.3


From 2353624fcb8241222e2044cb9d10ffa7c11c87c6 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Fri, 13 Oct 2017 16:55:15 +0100
Subject: changes for vishighmem

---
 src/Python/test/readhd5.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

(limited to 'src/Python/test')

diff --git a/src/Python/test/readhd5.py b/src/Python/test/readhd5.py
index 1e19e14..406fda9 100644
--- a/src/Python/test/readhd5.py
+++ b/src/Python/test/readhd5.py
@@ -12,7 +12,7 @@ def getEntry(nx, location):
     for item in nx[location].keys():
         print (item)
         
-filename = r'C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\Demos\DendrData.h5'
+filename = r'/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/demos/DendrData.h5'
 nx = h5py.File(filename, "r")
 #getEntry(nx, '/')
 # I have exported the entries as children of /
@@ -25,4 +25,4 @@ angSize = numpy.asarray(nx.get('/angSize'), dtype=int)[0]
 angles_rad = numpy.asarray(nx.get('/angles_rad'))
 recon_size = numpy.asarray(nx.get('/recon_size'), dtype=int)[0]
 size_det = numpy.asarray(nx.get('/size_det'), dtype=int)[0]
-slices_tot = numpy.asarray(nx.get('/slices_tot'), dtype=int)[0]
\ No newline at end of file
+slices_tot = numpy.asarray(nx.get('/slices_tot'), dtype=int)[0]
-- 
cgit v1.2.3


From 4fd4f187a70c0e4f56d5194b09ab4a528d20ee51 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Tue, 24 Oct 2017 10:38:37 +0100
Subject: Builds 2 packages fista and regularizers

---
 src/Python/test/test_reconstructor.py | 301 ++++++++++++++++++++++++++++++++++
 1 file changed, 301 insertions(+)
 create mode 100644 src/Python/test/test_reconstructor.py

(limited to 'src/Python/test')

diff --git a/src/Python/test/test_reconstructor.py b/src/Python/test/test_reconstructor.py
new file mode 100644
index 0000000..07668ba
--- /dev/null
+++ b/src/Python/test/test_reconstructor.py
@@ -0,0 +1,301 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Aug 23 16:34:49 2017
+
+@author: ofn77899
+Based on DemoRD2.m
+"""
+
+import h5py
+import numpy
+
+from ccpi.fista.FISTAReconstructor import FISTAReconstructor
+import astra
+import matplotlib.pyplot as plt
+
+def RMSE(signal1, signal2):
+    '''RMSE Root Mean Squared Error'''
+    if numpy.shape(signal1) == numpy.shape(signal2):
+        err = (signal1 - signal2)
+        err = numpy.sum( err * err )/numpy.size(signal1);  # MSE
+        err = sqrt(err);                                   # RMSE
+        return err
+    else:
+        raise Exception('Input signals must have the same shape')
+  
+filename = r'/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/demos/DendrData.h5'
+nx = h5py.File(filename, "r")
+#getEntry(nx, '/')
+# I have exported the entries as children of /
+entries = [entry for entry in nx['/'].keys()]
+print (entries)
+
+Sino3D = numpy.asarray(nx.get('/Sino3D'), dtype="float32")
+Weights3D = numpy.asarray(nx.get('/Weights3D'), dtype="float32")
+angSize = numpy.asarray(nx.get('/angSize'), dtype=int)[0]
+angles_rad = numpy.asarray(nx.get('/angles_rad'), dtype="float32")
+recon_size = numpy.asarray(nx.get('/recon_size'), dtype=int)[0]
+size_det = numpy.asarray(nx.get('/size_det'), dtype=int)[0]
+slices_tot = numpy.asarray(nx.get('/slices_tot'), dtype=int)[0]
+
+Z_slices = 20
+det_row_count = Z_slices
+# next definition is just for consistency of naming
+det_col_count = size_det
+
+detectorSpacingX = 1.0
+detectorSpacingY = detectorSpacingX
+
+
+proj_geom = astra.creators.create_proj_geom('parallel3d',
+                                            detectorSpacingX,
+                                            detectorSpacingY,
+                                            det_row_count,
+                                            det_col_count,
+                                            angles_rad)
+
+#vol_geom = astra_create_vol_geom(recon_size,recon_size,Z_slices);
+image_size_x = recon_size
+image_size_y = recon_size
+image_size_z = Z_slices
+vol_geom = astra.creators.create_vol_geom( image_size_x,
+                                           image_size_y,
+                                           image_size_z)
+
+## First pass the arguments to the FISTAReconstructor and test the
+## Lipschitz constant
+
+fistaRecon = FISTAReconstructor(proj_geom,
+                                vol_geom,
+                                Sino3D ,
+                                weights=Weights3D)
+
+print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant']))
+fistaRecon.setParameter(number_of_iterations = 12)
+fistaRecon.setParameter(Lipschitz_constant = 767893952.0)
+fistaRecon.setParameter(ring_alpha = 21)
+fistaRecon.setParameter(ring_lambda_R_L1 = 0.002)
+
+## Ordered subset
+if False:
+    subsets = 16
+    angles = fistaRecon.getParameter('projector_geometry')['ProjectionAngles']
+    #binEdges = numpy.linspace(angles.min(),
+    #                          angles.max(),
+    #                          subsets + 1)
+    binsDiscr, binEdges = numpy.histogram(angles, bins=subsets)
+    # get rearranged subset indices
+    IndicesReorg = numpy.zeros((numpy.shape(angles)))
+    counterM = 0
+    for ii in range(binsDiscr.max()):
+        counter = 0
+        for jj in range(subsets):
+            curr_index = ii + jj  + counter
+            #print ("{0} {1} {2}".format(binsDiscr[jj] , ii, counterM))
+            if binsDiscr[jj] > ii:
+                if (counterM < numpy.size(IndicesReorg)):
+                    IndicesReorg[counterM] = curr_index
+                counterM = counterM + 1
+                
+            counter = counter + binsDiscr[jj] - 1
+
+
+if False:
+    print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant']))
+    print ("prepare for iteration")
+    fistaRecon.prepareForIteration()
+    
+    
+
+    print("initializing  ...")
+    if False:
+        # if X doesn't exist
+        #N = params.vol_geom.GridColCount
+        N = vol_geom['GridColCount']
+        print ("N " + str(N))
+        X = numpy.zeros((N,N,SlicesZ), dtype=numpy.float)
+    else:
+        #X = fistaRecon.initialize()
+        X = numpy.load("X.npy")
+    
+    print (numpy.shape(X))
+    X_t = X.copy()
+    print ("initialized")
+    proj_geom , vol_geom, sino , \
+                      SlicesZ = fistaRecon.getParameter(['projector_geometry' ,
+                                                            'output_geometry',
+                                                            'input_sinogram',
+                                                         'SlicesZ'])
+
+    #fistaRecon.setParameter(number_of_iterations = 3)
+    iterFISTA = fistaRecon.getParameter('number_of_iterations')
+    # errors vector (if the ground truth is given)
+    Resid_error = numpy.zeros((iterFISTA));
+    # objective function values vector
+    objective = numpy.zeros((iterFISTA)); 
+
+      
+    t = 1
+
+
+    print ("starting iterations")
+##    % Outer FISTA iterations loop
+    for i in range(fistaRecon.getParameter('number_of_iterations')):
+        X_old = X.copy()
+        t_old = t
+        r_old = fistaRecon.r.copy()
+        if fistaRecon.getParameter('projector_geometry')['type'] == 'parallel' or \
+           fistaRecon.getParameter('projector_geometry')['type'] == 'fanflat' or \
+           fistaRecon.getParameter('projector_geometry')['type'] == 'fanflat_vec' :
+            # if the geometry is parallel use slice-by-slice
+            # projection-backprojection routine
+            #sino_updt = zeros(size(sino),'single');
+            proj_geomT = proj_geom.copy()
+            proj_geomT['DetectorRowCount'] = 1
+            vol_geomT = vol_geom.copy()
+            vol_geomT['GridSliceCount'] = 1;
+            sino_updt = numpy.zeros(numpy.shape(sino), dtype=numpy.float)
+            for kkk in range(SlicesZ):
+                sino_id, sino_updt[kkk] = \
+                         astra.creators.create_sino3d_gpu(
+                             X_t[kkk:kkk+1], proj_geom, vol_geom)
+                astra.matlab.data3d('delete', sino_id)
+        else:
+            # for divergent 3D geometry (watch the GPU memory overflow in
+            # ASTRA versions < 1.8)
+            #[sino_id, sino_updt] = astra_create_sino3d_cuda(X_t, proj_geom, vol_geom);
+            sino_id, sino_updt = astra.creators.create_sino3d_gpu(
+                X_t, proj_geom, vol_geom)
+
+        ## RING REMOVAL
+        residual = fistaRecon.residual
+        lambdaR_L1 , alpha_ring , weights , L_const= \
+                   fistaRecon.getParameter(['ring_lambda_R_L1',
+                                           'ring_alpha' , 'weights',
+                                           'Lipschitz_constant'])
+        r_x = fistaRecon.r_x
+        SlicesZ, anglesNumb, Detectors = \
+                    numpy.shape(fistaRecon.getParameter('input_sinogram'))
+        if lambdaR_L1 > 0 :
+             print ("ring removal")
+             for kkk in range(anglesNumb):
+                 
+                 residual[:,kkk,:] = (weights[:,kkk,:]).squeeze() * \
+                                       ((sino_updt[:,kkk,:]).squeeze() - \
+                                        (sino[:,kkk,:]).squeeze() -\
+                                        (alpha_ring * r_x)
+                                        )
+             vec = residual.sum(axis = 1)
+             #if SlicesZ > 1:
+             #    vec = vec[:,1,:].squeeze()
+             fistaRecon.r = (r_x - (1./L_const) * vec).copy()
+             objective[i] = (0.5 * (residual ** 2).sum())
+##            % the ring removal part (Group-Huber fidelity)
+##            for kkk = 1:anglesNumb
+##                residual(:,kkk,:) =  squeeze(weights(:,kkk,:)).*
+##                 (squeeze(sino_updt(:,kkk,:)) -
+##                 (squeeze(sino(:,kkk,:)) - alpha_ring.*r_x));
+##            end
+##            vec = sum(residual,2);
+##            if (SlicesZ > 1)
+##                vec = squeeze(vec(:,1,:));
+##            end
+##            r = r_x - (1./L_const).*vec;
+##            objective(i) = (0.5*sum(residual(:).^2)); % for the objective function output
+
+        
+
+        # Projection/Backprojection Routine
+        if fistaRecon.getParameter('projector_geometry')['type'] == 'parallel' or \
+           fistaRecon.getParameter('projector_geometry')['type'] == 'fanflat' or\
+           fistaRecon.getParameter('projector_geometry')['type'] == 'fanflat_vec':
+            x_temp = numpy.zeros(numpy.shape(X),dtype=numpy.float32)
+            print ("Projection/Backprojection Routine")
+            for kkk in range(SlicesZ):
+                
+                x_id, x_temp[kkk] = \
+                         astra.creators.create_backprojection3d_gpu(
+                             residual[kkk:kkk+1],
+                             proj_geomT, vol_geomT)
+                astra.matlab.data3d('delete', x_id)
+        else:
+            x_id, x_temp = \
+                  astra.creators.create_backprojection3d_gpu(
+                      residual, proj_geom, vol_geom)            
+
+        X = X_t - (1/L_const) * x_temp
+        astra.matlab.data3d('delete', sino_id)
+        astra.matlab.data3d('delete', x_id)
+        
+
+        ## REGULARIZATION
+        ## SKIPPING FOR NOW
+        ## Should be simpli
+        # regularizer = fistaRecon.getParameter('regularizer')
+        # for slices:
+        # out = regularizer(input=X)
+        print ("skipping regularizer")
+
+
+        ## FINAL
+        print ("final")
+        lambdaR_L1 = fistaRecon.getParameter('ring_lambda_R_L1')
+        if lambdaR_L1 > 0:
+            fistaRecon.r = numpy.max(
+                numpy.abs(fistaRecon.r) - lambdaR_L1 , 0) * \
+                numpy.sign(fistaRecon.r)
+        t = (1 + numpy.sqrt(1 + 4 * t**2))/2
+        X_t = X + (((t_old -1)/t) * (X - X_old))
+
+        if lambdaR_L1 > 0:
+            fistaRecon.r_x = fistaRecon.r + \
+                             (((t_old-1)/t) * (fistaRecon.r - r_old))
+
+        if fistaRecon.getParameter('region_of_interest') is None:
+            string = 'Iteration Number {0} | Objective {1} \n'
+            print (string.format( i, objective[i]))
+        else:
+            ROI , X_ideal = fistaRecon.getParameter('region_of_interest',
+                                                    'ideal_image')
+            
+            Resid_error[i] = RMSE(X*ROI, X_ideal*ROI)
+            string = 'Iteration Number {0} | RMS Error {1} | Objective {2} \n'
+            print (string.format(i,Resid_error[i], objective[i]))
+            
+##        if (lambdaR_L1 > 0)
+##            r =  max(abs(r)-lambdaR_L1, 0).*sign(r); % soft-thresholding operator for ring vector
+##        end
+##        
+##        t = (1 + sqrt(1 + 4*t^2))/2; % updating t
+##        X_t = X + ((t_old-1)/t).*(X - X_old); % updating X
+##        
+##        if (lambdaR_L1 > 0)
+##            r_x = r + ((t_old-1)/t).*(r - r_old); % updating r
+##        end
+##        
+##        if (show == 1)
+##            figure(10); imshow(X(:,:,slice), [0 maxvalplot]);
+##            if (lambdaR_L1 > 0)
+##                figure(11); plot(r); title('Rings offset vector')
+##            end
+##            pause(0.01);
+##        end
+##        if (strcmp(X_ideal, 'none' ) == 0)
+##            Resid_error(i) = RMSE(X(ROI), X_ideal(ROI));
+##            fprintf('%s %i %s %s %.4f  %s %s %f \n', 'Iteration Number:', i, '|', 'Error RMSE:', Resid_error(i), '|', 'Objective:', objective(i));
+##        else
+##            fprintf('%s %i  %s %s %f \n', 'Iteration Number:', i, '|', 'Objective:', objective(i));
+##        end
+else:
+    fistaRecon = FISTAReconstructor(proj_geom,
+                                vol_geom,
+                                Sino3D ,
+                                weights=Weights3D)
+
+    print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant']))
+    fistaRecon.setParameter(number_of_iterations = 12)
+    fistaRecon.setParameter(Lipschitz_constant = 767893952.0)
+    fistaRecon.setParameter(ring_alpha = 21)
+    fistaRecon.setParameter(ring_lambda_R_L1 = 0.002)
+    fistaRecon.prepareForIteration()
+    X = fistaRecon.iterate(numpy.load("X.npy"))
-- 
cgit v1.2.3


From 909a7bb4d71bdb14d4e68f42c2297f6154a77ed0 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Tue, 24 Oct 2017 11:27:17 +0100
Subject: use system package

---
 src/Python/test/test_reconstructor.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

(limited to 'src/Python/test')

diff --git a/src/Python/test/test_reconstructor.py b/src/Python/test/test_reconstructor.py
index 07668ba..377f3cf 100644
--- a/src/Python/test/test_reconstructor.py
+++ b/src/Python/test/test_reconstructor.py
@@ -9,7 +9,7 @@ Based on DemoRD2.m
 import h5py
 import numpy
 
-from ccpi.fista.FISTAReconstructor import FISTAReconstructor
+from ccpi.reconstruction.FISTAReconstructor import FISTAReconstructor
 import astra
 import matplotlib.pyplot as plt
 
-- 
cgit v1.2.3


From bb4f7dc7e3a3bf4b4da18e36a2fc69e2195c5a96 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Tue, 24 Oct 2017 15:47:40 +0100
Subject: moved to test directory

---
 src/Python/test/test_reconstructor-os.py | 338 +++++++++++++++++++++++++++++++
 1 file changed, 338 insertions(+)
 create mode 100644 src/Python/test/test_reconstructor-os.py

(limited to 'src/Python/test')

diff --git a/src/Python/test/test_reconstructor-os.py b/src/Python/test/test_reconstructor-os.py
new file mode 100644
index 0000000..a36feda
--- /dev/null
+++ b/src/Python/test/test_reconstructor-os.py
@@ -0,0 +1,338 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Aug 23 16:34:49 2017
+
+@author: ofn77899
+Based on DemoRD2.m
+"""
+
+import h5py
+import numpy
+
+from ccpi.reconstruction.FISTAReconstructor import FISTAReconstructor
+import astra
+import matplotlib.pyplot as plt
+
+def RMSE(signal1, signal2):
+    '''RMSE Root Mean Squared Error'''
+    if numpy.shape(signal1) == numpy.shape(signal2):
+        err = (signal1 - signal2)
+        err = numpy.sum( err * err )/numpy.size(signal1);  # MSE
+        err = sqrt(err);                                   # RMSE
+        return err
+    else:
+        raise Exception('Input signals must have the same shape')
+  
+filename = r'/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/demos/DendrData.h5'
+nx = h5py.File(filename, "r")
+#getEntry(nx, '/')
+# I have exported the entries as children of /
+entries = [entry for entry in nx['/'].keys()]
+print (entries)
+
+Sino3D = numpy.asarray(nx.get('/Sino3D'), dtype="float32")
+Weights3D = numpy.asarray(nx.get('/Weights3D'), dtype="float32")
+angSize = numpy.asarray(nx.get('/angSize'), dtype=int)[0]
+angles_rad = numpy.asarray(nx.get('/angles_rad'), dtype="float32")
+recon_size = numpy.asarray(nx.get('/recon_size'), dtype=int)[0]
+size_det = numpy.asarray(nx.get('/size_det'), dtype=int)[0]
+slices_tot = numpy.asarray(nx.get('/slices_tot'), dtype=int)[0]
+
+Z_slices = 20
+det_row_count = Z_slices
+# next definition is just for consistency of naming
+det_col_count = size_det
+
+detectorSpacingX = 1.0
+detectorSpacingY = detectorSpacingX
+
+
+proj_geom = astra.creators.create_proj_geom('parallel3d',
+                                            detectorSpacingX,
+                                            detectorSpacingY,
+                                            det_row_count,
+                                            det_col_count,
+                                            angles_rad)
+
+#vol_geom = astra_create_vol_geom(recon_size,recon_size,Z_slices);
+image_size_x = recon_size
+image_size_y = recon_size
+image_size_z = Z_slices
+vol_geom = astra.creators.create_vol_geom( image_size_x,
+                                           image_size_y,
+                                           image_size_z)
+
+## First pass the arguments to the FISTAReconstructor and test the
+## Lipschitz constant
+
+fistaRecon = FISTAReconstructor(proj_geom,
+                                vol_geom,
+                                Sino3D ,
+                                weights=Weights3D)
+
+print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant']))
+fistaRecon.setParameter(number_of_iterations = 12)
+fistaRecon.setParameter(Lipschitz_constant = 767893952.0)
+fistaRecon.setParameter(ring_alpha = 21)
+fistaRecon.setParameter(ring_lambda_R_L1 = 0.002)
+
+
+## Ordered subset
+if True:
+    subsets = 16
+    fistaRecon.setParameter(subsets=subsets)
+    fistaRecon.createOrderedSubsets()
+    angles = fistaRecon.getParameter('projector_geometry')['ProjectionAngles']
+    #binEdges = numpy.linspace(angles.min(),
+    #                          angles.max(),
+    #                          subsets + 1)
+    binsDiscr, binEdges = numpy.histogram(angles, bins=subsets)
+    # get rearranged subset indices
+    IndicesReorg = numpy.zeros((numpy.shape(angles)))
+    counterM = 0
+    for ii in range(binsDiscr.max()):
+        counter = 0
+        for jj in range(subsets):
+            curr_index = ii + jj  + counter
+            #print ("{0} {1} {2}".format(binsDiscr[jj] , ii, counterM))
+            if binsDiscr[jj] > ii:
+                if (counterM < numpy.size(IndicesReorg)):
+                    IndicesReorg[counterM] = curr_index
+                counterM = counterM + 1
+                
+            counter = counter + binsDiscr[jj] - 1
+
+
+if True:
+    print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant']))
+    print ("prepare for iteration")
+    fistaRecon.prepareForIteration()
+    
+    
+
+    print("initializing  ...")
+    if False:
+        # if X doesn't exist
+        #N = params.vol_geom.GridColCount
+        N = vol_geom['GridColCount']
+        print ("N " + str(N))
+        X = numpy.zeros((N,N,SlicesZ), dtype=numpy.float)
+    else:
+        #X = fistaRecon.initialize()
+        X = numpy.load("X.npy")
+    
+    print (numpy.shape(X))
+    X_t = X.copy()
+    print ("initialized")
+    proj_geom , vol_geom, sino , \
+        SlicesZ, weights , alpha_ring = fistaRecon.getParameter(
+            ['projector_geometry' , 'output_geometry',
+             'input_sinogram', 'SlicesZ' ,  'weights', 'ring_alpha'])
+    lambdaR_L1 , alpha_ring , weights , L_const= \
+                       fistaRecon.getParameter(['ring_lambda_R_L1',
+                                               'ring_alpha' , 'weights',
+                                               'Lipschitz_constant'])
+
+    #fistaRecon.setParameter(number_of_iterations = 3)
+    iterFISTA = fistaRecon.getParameter('number_of_iterations')
+    # errors vector (if the ground truth is given)
+    Resid_error = numpy.zeros((iterFISTA));
+    # objective function values vector
+    objective = numpy.zeros((iterFISTA)); 
+
+      
+    t = 1
+    
+
+    ## additional for 
+    proj_geomSUB = proj_geom.copy()
+    fistaRecon.residual2 = numpy.zeros(numpy.shape(fistaRecon.pars['input_sinogram']))
+    residual2 = fistaRecon.residual2
+    sino_updt_FULL = fistaRecon.residual.copy()
+    r_x = fistaRecon.r.copy()
+    
+    print ("starting iterations")
+##    % Outer FISTA iterations loop
+    for i in range(fistaRecon.getParameter('number_of_iterations')):
+##        % With OS approach it becomes trickier to correlate independent subsets, hence additional work is required
+##        % one solution is to work with a full sinogram at times
+##        if ((i >= 3) && (lambdaR_L1 > 0))
+##            [sino_id2, sino_updt2] = astra_create_sino3d_cuda(X, proj_geom, vol_geom);
+##            astra_mex_data3d('delete', sino_id2);
+##        end
+        # With OS approach it becomes trickier to correlate independent subsets,
+        # hence additional work is required one solution is to work with a full
+        # sinogram at times
+
+        r_old = fistaRecon.r.copy()
+        t_old = t
+        SlicesZ, anglesNumb, Detectors = \
+                    numpy.shape(fistaRecon.getParameter('input_sinogram'))        ## https://github.com/vais-ral/CCPi-FISTA_Reconstruction/issues/4
+        if (i > 1 and lambdaR_L1 > 0) :
+            for kkk in range(anglesNumb):
+                 
+                 residual2[:,kkk,:] = (weights[:,kkk,:]).squeeze() * \
+                                       ((sino_updt_FULL[:,kkk,:]).squeeze() - \
+                                        (sino[:,kkk,:]).squeeze() -\
+                                        (alpha_ring * r_x)
+                                        )
+            
+            vec = fistaRecon.residual.sum(axis = 1)
+            #if SlicesZ > 1:
+            #    vec = vec[:,1,:] # 1 or 0?
+            r_x = fistaRecon.r_x
+            fistaRecon.r = (r_x - (1./L_const) * vec).copy()
+
+        # subset loop
+        counterInd = 1
+        geometry_type = fistaRecon.getParameter('projector_geometry')['type']
+        if geometry_type == 'parallel' or \
+           geometry_type == 'fanflat' or \
+           geometry_type == 'fanflat_vec' :
+            
+            for kkk in range(SlicesZ):
+                sino_id, sinoT[kkk] = \
+                         astra.creators.create_sino3d_gpu(
+                             X_t[kkk:kkk+1], proj_geomSUB, vol_geom)
+                sino_updt_Sub[kkk] = sinoT.T.copy()
+                
+        else:
+            sino_id, sino_updt_Sub = \
+                astra.creators.create_sino3d_gpu(X_t, proj_geomSUB, vol_geom)
+        
+        astra.matlab.data3d('delete', sino_id)
+  
+        for ss in range(fistaRecon.getParameter('subsets')):
+            print ("Subset {0}".format(ss))
+            X_old = X.copy()
+            t_old = t
+            
+            # the number of projections per subset
+            numProjSub = fistaRecon.getParameter('os_bins')[ss]
+            CurrSubIndices = fistaRecon.getParameter('os_indices')\
+                             [counterInd:counterInd+numProjSub-1]
+            mask = numpy.zeros(numpy.shape(angles), dtype=bool)
+            cc = 0
+            for i in range(len(CurrSubIndices)):
+                mask[int(CurrSubIndices[i])] = True
+            proj_geomSUB['ProjectionAngles'] = angles[mask]
+
+            shape = list(numpy.shape(fistaRecon.getParameter('input_sinogram')))
+            shape[1] = numProjSub
+            sino_updt_Sub = numpy.zeros(shape)
+
+            if geometry_type == 'parallel' or \
+               geometry_type == 'fanflat' or \
+               geometry_type == 'fanflat_vec' :
+
+                for kkk in range(SlicesZ):
+                    sino_id, sinoT = astra.creators.create_sino3d_gpu (
+                        X_t[kkk:kkk+1] , proj_geomSUB, vol_geom)
+                    sino_updt_Sub[kkk] = sinoT.T.copy()
+                    
+            else:
+                # for 3D geometry (watch the GPU memory overflow in ASTRA < 1.8)
+                sino_id, sino_updt_Sub = \
+                     astra.creators.create_sino3d_gpu (X_t, proj_geomSUB, vol_geom)
+                
+            astra.matlab.data3d('delete', sino_id)
+                
+            
+
+
+            ## RING REMOVAL
+            residual = fistaRecon.residual
+            
+            
+            if lambdaR_L1 > 0 :
+                 print ("ring removal")
+                 residualSub = numpy.zeros(shape)
+    ##             for a chosen subset
+    ##                for kkk = 1:numProjSub
+    ##                    indC = CurrSubIndeces(kkk);
+    ##                    residualSub(:,kkk,:) =  squeeze(weights(:,indC,:)).*(squeeze(sino_updt_Sub(:,kkk,:)) - (squeeze(sino(:,indC,:)) - alpha_ring.*r_x));
+    ##                    sino_updt_FULL(:,indC,:) = squeeze(sino_updt_Sub(:,kkk,:)); % filling the full sinogram
+    ##                end
+                 for kkk in range(numProjSub):
+                     print ("ring removal indC ... {0}".format(kkk))
+                     indC = int(CurrSubIndices[kkk])
+                     residualSub[:,kkk,:] = weights[:,indC,:].squeeze() * \
+                            (sino_updt_Sub[:,kkk,:].squeeze() - \
+                             sino[:,indC,:].squeeze() - alpha_ring * r_x)
+                     # filling the full sinogram
+                     sino_updt_FULL[:,indC,:] = sino_updt_Sub[:,kkk,:].squeeze()
+
+            else:
+                #PWLS model
+                residualSub = weights[:,CurrSubIndices,:] * \
+                              ( sino_updt_Sub - sino[:,CurrSubIndices,:].squeeze() )
+                objective[i] = 0.5 * numpy.linalg.norm(residualSub)
+
+            if geometry_type == 'parallel' or \
+               geometry_type == 'fanflat' or \
+               geometry_type == 'fanflat_vec' :
+                # if geometry is 2D use slice-by-slice projection-backprojection
+                # routine
+                x_temp = numpy.zeros(numpy.shape(X), dtype=numpy.float32)
+                for kkk in range(SlicesZ):
+                    
+                    x_id, x_temp[kkk] = \
+                             astra.creators.create_backprojection3d_gpu(
+                                 residualSub[kkk:kkk+1],
+                                 proj_geomSUB, vol_geom)
+                    
+            else:
+                x_id, x_temp = \
+                      astra.creators.create_backprojection3d_gpu(
+                          residualSub, proj_geomSUB, vol_geom)
+
+            astra.matlab.data3d('delete', x_id)
+            X = X_t - (1/L_const) * x_temp
+
+        
+
+        ## REGULARIZATION
+        ## SKIPPING FOR NOW
+        ## Should be simpli
+        # regularizer = fistaRecon.getParameter('regularizer')
+        # for slices:
+        # out = regularizer(input=X)
+        print ("skipping regularizer")
+
+
+        ## FINAL
+        print ("final")
+        lambdaR_L1 = fistaRecon.getParameter('ring_lambda_R_L1')
+        if lambdaR_L1 > 0:
+            fistaRecon.r = numpy.max(
+                numpy.abs(fistaRecon.r) - lambdaR_L1 , 0) * \
+                numpy.sign(fistaRecon.r)
+            # updating r
+            r_x = fistaRecon.r + ((t_old-1)/t) * (fistaRecon.r - r_old)
+        
+
+        if fistaRecon.getParameter('region_of_interest') is None:
+            string = 'Iteration Number {0} | Objective {1} \n'
+            print (string.format( i, objective[i]))
+        else:
+            ROI , X_ideal = fistaRecon.getParameter('region_of_interest',
+                                                    'ideal_image')
+            
+            Resid_error[i] = RMSE(X*ROI, X_ideal*ROI)
+            string = 'Iteration Number {0} | RMS Error {1} | Objective {2} \n'
+            print (string.format(i,Resid_error[i], objective[i]))
+            
+
+else:
+    fistaRecon = FISTAReconstructor(proj_geom,
+                                vol_geom,
+                                Sino3D ,
+                                weights=Weights3D)
+
+    print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant']))
+    fistaRecon.setParameter(number_of_iterations = 12)
+    fistaRecon.setParameter(Lipschitz_constant = 767893952.0)
+    fistaRecon.setParameter(ring_alpha = 21)
+    fistaRecon.setParameter(ring_lambda_R_L1 = 0.002)
+    fistaRecon.prepareForIteration()
+    X = fistaRecon.iterate(numpy.load("X.npy"))
-- 
cgit v1.2.3


From 7f6e90ed9569e6f935813d8ceb6b3c00feed3bc0 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Tue, 24 Oct 2017 15:48:06 +0100
Subject: saves to file

---
 src/Python/test/test_reconstructor.py | 1 +
 1 file changed, 1 insertion(+)

(limited to 'src/Python/test')

diff --git a/src/Python/test/test_reconstructor.py b/src/Python/test/test_reconstructor.py
index 377f3cf..be28624 100644
--- a/src/Python/test/test_reconstructor.py
+++ b/src/Python/test/test_reconstructor.py
@@ -299,3 +299,4 @@ else:
     fistaRecon.setParameter(ring_lambda_R_L1 = 0.002)
     fistaRecon.prepareForIteration()
     X = fistaRecon.iterate(numpy.load("X.npy"))
+    numpy.save("X_out.npy", X)
-- 
cgit v1.2.3


From 455ca86825c157512f61441d3d27b8148ca795a7 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Tue, 24 Oct 2017 16:37:21 +0100
Subject: Add regularization step

Add regularization step
OS seems to work
---
 src/Python/test/test_reconstructor-os.py | 22 ++++++++++++++++------
 src/Python/test/test_reconstructor.py    |  7 +++++++
 2 files changed, 23 insertions(+), 6 deletions(-)

(limited to 'src/Python/test')

diff --git a/src/Python/test/test_reconstructor-os.py b/src/Python/test/test_reconstructor-os.py
index a36feda..6c82ae0 100644
--- a/src/Python/test/test_reconstructor-os.py
+++ b/src/Python/test/test_reconstructor-os.py
@@ -12,6 +12,7 @@ import numpy
 from ccpi.reconstruction.FISTAReconstructor import FISTAReconstructor
 import astra
 import matplotlib.pyplot as plt
+from ccpi.imaging.Regularizer import Regularizer
 
 def RMSE(signal1, signal2):
     '''RMSE Root Mean Squared Error'''
@@ -77,6 +78,12 @@ fistaRecon.setParameter(ring_alpha = 21)
 fistaRecon.setParameter(ring_lambda_R_L1 = 0.002)
 
 
+reg = Regularizer(Regularizer.Algorithm.LLT_model)
+reg.setParameter(regularization_parameter=25,
+                          time_step=0.0003,
+                          tolerance_constant=0.0001,
+                          number_of_iterations=300)
+
 ## Ordered subset
 if True:
     subsets = 16
@@ -210,11 +217,12 @@ if True:
             # the number of projections per subset
             numProjSub = fistaRecon.getParameter('os_bins')[ss]
             CurrSubIndices = fistaRecon.getParameter('os_indices')\
-                             [counterInd:counterInd+numProjSub-1]
+                             [counterInd:counterInd+numProjSub]
+            #print ("Len CurrSubIndices {0}".format(numProjSub))
             mask = numpy.zeros(numpy.shape(angles), dtype=bool)
             cc = 0
-            for i in range(len(CurrSubIndices)):
-                mask[int(CurrSubIndices[i])] = True
+            for j in range(len(CurrSubIndices)):
+                mask[int(CurrSubIndices[j])] = True
             proj_geomSUB['ProjectionAngles'] = angles[mask]
 
             shape = list(numpy.shape(fistaRecon.getParameter('input_sinogram')))
@@ -254,7 +262,7 @@ if True:
     ##                    sino_updt_FULL(:,indC,:) = squeeze(sino_updt_Sub(:,kkk,:)); % filling the full sinogram
     ##                end
                  for kkk in range(numProjSub):
-                     print ("ring removal indC ... {0}".format(kkk))
+                     #print ("ring removal indC ... {0}".format(kkk))
                      indC = int(CurrSubIndices[kkk])
                      residualSub[:,kkk,:] = weights[:,indC,:].squeeze() * \
                             (sino_updt_Sub[:,kkk,:].squeeze() - \
@@ -297,7 +305,8 @@ if True:
         # regularizer = fistaRecon.getParameter('regularizer')
         # for slices:
         # out = regularizer(input=X)
-        print ("skipping regularizer")
+        print ("regularizer")
+        #X = reg(input=X)
 
 
         ## FINAL
@@ -321,7 +330,8 @@ if True:
             Resid_error[i] = RMSE(X*ROI, X_ideal*ROI)
             string = 'Iteration Number {0} | RMS Error {1} | Objective {2} \n'
             print (string.format(i,Resid_error[i], objective[i]))
-            
+
+    numpy.save("X_out_os.npy", X)
 
 else:
     fistaRecon = FISTAReconstructor(proj_geom,
diff --git a/src/Python/test/test_reconstructor.py b/src/Python/test/test_reconstructor.py
index be28624..3342301 100644
--- a/src/Python/test/test_reconstructor.py
+++ b/src/Python/test/test_reconstructor.py
@@ -76,6 +76,13 @@ fistaRecon.setParameter(Lipschitz_constant = 767893952.0)
 fistaRecon.setParameter(ring_alpha = 21)
 fistaRecon.setParameter(ring_lambda_R_L1 = 0.002)
 
+reg = Regularizer(Regularizer.Algorithm.LLT_model)
+reg.setParameter(regularization_parameter=25,
+                          time_step=0.0003,
+                          tolerance_constant=0.0001,
+                          number_of_iterations=300)
+fistaRecon.setParameter(regularizer = reg)
+
 ## Ordered subset
 if False:
     subsets = 16
-- 
cgit v1.2.3


From ea490be883d6933c481383a42e0d50edb19ece73 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Wed, 25 Oct 2017 16:25:50 +0100
Subject: added to repository

---
 src/Python/test/test_regularizers_3d.py | 380 ++++++++++++++++++++++++++++++++
 1 file changed, 380 insertions(+)
 create mode 100644 src/Python/test/test_regularizers_3d.py

(limited to 'src/Python/test')

diff --git a/src/Python/test/test_regularizers_3d.py b/src/Python/test/test_regularizers_3d.py
new file mode 100644
index 0000000..a2e3027
--- /dev/null
+++ b/src/Python/test/test_regularizers_3d.py
@@ -0,0 +1,380 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Aug  4 11:10:05 2017
+
+@author: ofn77899
+"""
+
+from ccpi.viewer.CILViewer2D import Converter
+import vtk
+
+import regularizers
+import matplotlib.pyplot as plt
+import numpy as np
+import os    
+from enum import Enum
+import timeit
+
+from Regularizer import Regularizer
+
+###############################################################################
+#https://stackoverflow.com/questions/13875989/comparing-image-in-url-to-image-in-filesystem-in-python/13884956#13884956
+#NRMSE a normalization of the root of the mean squared error
+#NRMSE is simply 1 - [RMSE / (maxval - minval)]. Where maxval is the maximum
+# intensity from the two images being compared, and respectively the same for
+# minval. RMSE is given by the square root of MSE: 
+# sqrt[(sum(A - B) ** 2) / |A|],
+# where |A| means the number of elements in A. By doing this, the maximum value
+# given by RMSE is maxval.
+
+def nrmse(im1, im2):
+    a, b = im1.shape
+    rmse = np.sqrt(np.sum((im2 - im1) ** 2) / float(a * b))
+    max_val = max(np.max(im1), np.max(im2))
+    min_val = min(np.min(im1), np.min(im2))
+    return 1 - (rmse / (max_val - min_val))
+###############################################################################
+
+###############################################################################
+#
+#  2D Regularizers
+#
+###############################################################################
+#Example:
+# figure;
+# Im = double(imread('lena_gray_256.tif'))/255;  % loading image
+# u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
+# u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+
+#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\lena_gray_512.tif"
+#reader = vtk.vtkTIFFReader()
+#reader.SetFileName(os.path.normpath(filename))
+#reader.Update()
+##vtk returns 3D images, let's take just the one slice there is as 2D
+#Im = Converter.vtk2numpy(reader.GetOutput()).T[0]/255
+#
+##imgplot = plt.imshow(Im)
+#perc = 0.05
+#u0 = Im + (perc* np.random.normal(size=np.shape(Im)))
+## map the u0 u0->u0>0
+#f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
+#u0 = f(u0).astype('float32')
+#
+### plot 
+#fig = plt.figure()
+##a=fig.add_subplot(3,3,1)
+##a.set_title('Original')
+##imgplot = plt.imshow(Im)
+#
+#a=fig.add_subplot(2,3,1)
+#a.set_title('noise')
+#imgplot = plt.imshow(u0)
+#
+#reg_output = []
+###############################################################################
+## Call regularizer
+#
+######################## SplitBregman_TV #####################################
+## u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+##reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV)
+##
+##out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30,
+##          #tolerance_constant=1e-4, 
+##          TV_Penalty=Regularizer.TotalVariationPenalty.l1)
+#
+##out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30,
+##          tolerance_constant=1e-4, 
+##          TV_Penalty=Regularizer.TotalVariationPenalty.l1)
+#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10. )
+#pars = out2[2]
+#reg_output.append(out2)
+#
+#a=fig.add_subplot(2,3,2)
+#
+#textstr = out2[-1]
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+#        verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output[-1][0])
+#
+####################### FGP_TV #########################################
+## u = FGP_TV(single(u0), 0.05, 100, 1e-04);
+#out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.005,
+#                          number_of_iterations=200)
+#pars = out2[-2]
+#
+#reg_output.append(out2)
+#
+#a=fig.add_subplot(2,3,3)
+#
+#textstr = out2[-1]
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+#        verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output[-1][0])
+#
+####################### LLT_model #########################################
+## * u0 = Im + .03*randn(size(Im)); % adding noise
+## [Den] = LLT_model(single(u0), 10, 0.1, 1);
+##Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0); 
+##input, regularization_parameter , time_step, number_of_iterations,
+##                  tolerance_constant, restrictive_Z_smoothing=0
+#out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25,
+#                          time_step=0.0003,
+#                          tolerance_constant=0.0001,
+#                          number_of_iterations=300)
+#pars = out2[-2]
+#
+#reg_output.append(out2)
+#
+#a=fig.add_subplot(2,3,4)
+#textstr = out2[-1]
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+#        verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output[-1][0])
+#
+####################### PatchBased_Regul #########################################
+## Quick 2D denoising example in Matlab:   
+##   Im = double(imread('lena_gray_256.tif'))/255;  % loading image
+##   u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+##   ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05); 
+#
+#out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05,
+#                          searching_window_ratio=3,
+#                          similarity_window_ratio=1,
+#                          PB_filtering_parameter=0.08)
+#pars = out2[-2]
+#reg_output.append(out2)
+#
+#a=fig.add_subplot(2,3,5)
+#
+#
+#textstr = out2[-1]
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+#        verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output[-1][0])
+#
+#
+####################### TGV_PD #########################################
+## Quick 2D denoising example in Matlab:   
+##   Im = double(imread('lena_gray_256.tif'))/255;  % loading image
+##   u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+##   u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550);
+#
+#
+#out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05,
+#                          first_order_term=1.3,
+#                          second_order_term=1,
+#                          number_of_iterations=550)
+#pars = out2[-2]
+#reg_output.append(out2)
+#
+#a=fig.add_subplot(2,3,6)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+#        verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output[-1][0])
+#
+
+###############################################################################
+#
+#  3D Regularizers
+#
+###############################################################################
+#Example:
+# figure;
+# Im = double(imread('lena_gray_256.tif'))/255;  % loading image
+# u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
+# u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+
+#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Reconstruction\python\test\reconstruction_example.mha"
+filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Simpleflex\data\head.mha"
+
+reader = vtk.vtkMetaImageReader()
+reader.SetFileName(os.path.normpath(filename))
+reader.Update()
+#vtk returns 3D images, let's take just the one slice there is as 2D
+Im = Converter.vtk2numpy(reader.GetOutput())
+Im = Im.astype('float32')
+#imgplot = plt.imshow(Im)
+perc = 0.05
+u0 = Im + (perc* np.random.normal(size=np.shape(Im)))
+# map the u0 u0->u0>0
+f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
+u0 = f(u0).astype('float32')
+converter = Converter.numpy2vtkImporter(u0, reader.GetOutput().GetSpacing(),
+                                        reader.GetOutput().GetOrigin())
+converter.Update()
+writer = vtk.vtkMetaImageWriter()
+writer.SetInputData(converter.GetOutput())
+writer.SetFileName(r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\noisy_head.mha")
+#writer.Write()
+
+
+## plot 
+fig3D = plt.figure(figsize=(20,16))
+
+#a=fig.add_subplot(3,3,1)
+#a.set_title('Original')
+#imgplot = plt.imshow(Im)
+sliceNo = 32
+
+a=fig3D.add_subplot(2,3,1)
+a.set_title('noise')
+imgplot = plt.imshow(u0.T[sliceNo])
+
+reg_output3d = []
+
+##############################################################################
+# Call regularizer
+
+####################### SplitBregman_TV #####################################
+# u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+
+#reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV)
+
+#out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30,
+#          #tolerance_constant=1e-4, 
+#          TV_Penalty=Regularizer.TotalVariationPenalty.l1)
+
+out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30,
+          tolerance_constant=1e-4, 
+          TV_Penalty=Regularizer.TotalVariationPenalty.l1)
+
+
+pars = out2[-2]
+reg_output3d.append(out2)
+
+a=fig3D.add_subplot(2,3,2)
+
+
+textstr = out2[-1]
+
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+        verticalalignment='top', bbox=props)
+imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+
+###################### FGP_TV #########################################
+# u = FGP_TV(single(u0), 0.05, 100, 1e-04);
+#out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.005,
+#                          number_of_iterations=200)
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,3,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+#        verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+
+###################### LLT_model #########################################
+# * u0 = Im + .03*randn(size(Im)); % adding noise
+# [Den] = LLT_model(single(u0), 10, 0.1, 1);
+#Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0); 
+#input, regularization_parameter , time_step, number_of_iterations,
+#                  tolerance_constant, restrictive_Z_smoothing=0
+out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25,
+                          time_step=0.0003,
+                          tolerance_constant=0.0001,
+                          number_of_iterations=300)
+pars = out2[-2]
+reg_output3d.append(out2)
+
+a=fig3D.add_subplot(2,3,3)
+
+
+textstr = out2[-1]
+
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+        verticalalignment='top', bbox=props)
+imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+
+###################### PatchBased_Regul #########################################
+# Quick 2D denoising example in Matlab:   
+#   Im = double(imread('lena_gray_256.tif'))/255;  % loading image
+#   u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+#   ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05); 
+
+out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05,
+                          searching_window_ratio=3,
+                          similarity_window_ratio=1,
+                          PB_filtering_parameter=0.08)
+pars = out2[-2]
+reg_output3d.append(out2)
+
+a=fig3D.add_subplot(2,3,4)
+
+
+textstr = out2[-1]
+
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+        verticalalignment='top', bbox=props)
+imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+
+
+####################### TGV_PD #########################################
+## Quick 2D denoising example in Matlab:   
+##   Im = double(imread('lena_gray_256.tif'))/255;  % loading image
+##   u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+##   u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550);
+#
+#
+#out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05,
+#                          first_order_term=1.3,
+#                          second_order_term=1,
+#                          number_of_iterations=550)
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,3,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+#        verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+fig3D.savefig('test\\3d.png')
+plt.close(fig3D)
\ No newline at end of file
-- 
cgit v1.2.3


From 6b24ef4e1e0780dc1eade61df025f886712339bc Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Wed, 25 Oct 2017 16:35:12 +0100
Subject: development

---
 src/Python/test/test_reconstructor-os.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

(limited to 'src/Python/test')

diff --git a/src/Python/test/test_reconstructor-os.py b/src/Python/test/test_reconstructor-os.py
index 6c82ae0..3f419cf 100644
--- a/src/Python/test/test_reconstructor-os.py
+++ b/src/Python/test/test_reconstructor-os.py
@@ -78,12 +78,16 @@ fistaRecon.setParameter(ring_alpha = 21)
 fistaRecon.setParameter(ring_lambda_R_L1 = 0.002)
 
 
+#reg = Regularizer(Regularizer.Algorithm.FGP_TV)
+#reg.setParameter(regularization_parameter=0.005,
+#                          number_of_iterations=50)
 reg = Regularizer(Regularizer.Algorithm.LLT_model)
 reg.setParameter(regularization_parameter=25,
                           time_step=0.0003,
                           tolerance_constant=0.0001,
                           number_of_iterations=300)
 
+
 ## Ordered subset
 if True:
     subsets = 16
@@ -306,7 +310,7 @@ if True:
         # for slices:
         # out = regularizer(input=X)
         print ("regularizer")
-        #X = reg(input=X)
+        X = reg(input=X)[0]
 
 
         ## FINAL
-- 
cgit v1.2.3


From ff9cc12694172e1e8720f7ea7f5b22e647722e21 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Wed, 25 Oct 2017 16:56:17 +0100
Subject: doing work

---
 src/Python/test/test_regularizers_3d.py | 471 +++++++++++++++++---------------
 1 file changed, 258 insertions(+), 213 deletions(-)

(limited to 'src/Python/test')

diff --git a/src/Python/test/test_regularizers_3d.py b/src/Python/test/test_regularizers_3d.py
index a2e3027..2d11a7e 100644
--- a/src/Python/test/test_regularizers_3d.py
+++ b/src/Python/test/test_regularizers_3d.py
@@ -5,17 +5,17 @@ Created on Fri Aug  4 11:10:05 2017
 @author: ofn77899
 """
 
-from ccpi.viewer.CILViewer2D import Converter
-import vtk
+#from ccpi.viewer.CILViewer2D import Converter
+#import vtk
 
-import regularizers
 import matplotlib.pyplot as plt
 import numpy as np
 import os    
 from enum import Enum
 import timeit
-
-from Regularizer import Regularizer
+#from PIL import Image
+#from Regularizer import Regularizer
+from ccpi.imaging.Regularizer import Regularizer
 
 ###############################################################################
 #https://stackoverflow.com/questions/13875989/comparing-image-in-url-to-image-in-filesystem-in-python/13884956#13884956
@@ -46,77 +46,303 @@ def nrmse(im1, im2):
 # u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
 # u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
 
+
 #filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\lena_gray_512.tif"
+filename = r"/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/data/lena_gray_512.tif"
+#filename = r'/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif'
+
 #reader = vtk.vtkTIFFReader()
 #reader.SetFileName(os.path.normpath(filename))
 #reader.Update()
-##vtk returns 3D images, let's take just the one slice there is as 2D
-#Im = Converter.vtk2numpy(reader.GetOutput()).T[0]/255
+Im = plt.imread(filename)                     
+#Im = Image.open('/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif')/255
+#img.show()
+Im = np.asarray(Im, dtype='float32')
+
+# create a 3D image by stacking N of this images
+
+
+#imgplot = plt.imshow(Im)
+perc = 0.05
+u_n = Im + (perc* np.random.normal(size=np.shape(Im)))
+y,z = np.shape(u_n)
+u_n = np.reshape(u_n , (1,y,z))
+
+u0 = u_n.copy()
+for i in range (19):
+    u_n = Im + (perc* np.random.normal(size=np.shape(Im)))
+    u_n = np.reshape(u_n , (1,y,z))
+
+    u0 = np.vstack ( (u0, u_n) )
+
+# map the u0 u0->u0>0
+f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
+u0 = f(u0).astype('float32')
+
+print ("Passed image shape {0}".format(np.shape(u0)))
+
+## plot 
+fig = plt.figure()
+#a=fig.add_subplot(3,3,1)
+#a.set_title('Original')
+#imgplot = plt.imshow(Im)
+sliceno = 10
+
+a=fig.add_subplot(2,3,1)
+a.set_title('noise')
+imgplot = plt.imshow(u0[sliceno],cmap="gray")
+
+reg_output = []
+##############################################################################
+# Call regularizer
+
+####################### SplitBregman_TV #####################################
+# u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+
+use_object = True
+if use_object:
+    reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV)
+    print (reg.pars)
+    reg.setParameter(input=u0)    
+    reg.setParameter(regularization_parameter=10.)
+    # or 
+    # reg.setParameter(input=u0, regularization_parameter=10., #number_of_iterations=30,
+              #tolerance_constant=1e-4, 
+              #TV_Penalty=Regularizer.TotalVariationPenalty.l1)
+    plotme = reg() [0]
+    pars = reg.pars
+    textstr = reg.printParametersToString() 
+    
+    #out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30,
+              #tolerance_constant=1e-4, 
+    #          TV_Penalty=Regularizer.TotalVariationPenalty.l1)
+    
+#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30,
+#          tolerance_constant=1e-4, 
+#          TV_Penalty=Regularizer.TotalVariationPenalty.l1)
+
+else:
+    out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10. )
+    pars = out2[2]
+    reg_output.append(out2)
+    plotme = reg_output[-1][0]
+    textstr = out2[-1]
+
+a=fig.add_subplot(2,3,2)
+
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+        verticalalignment='top', bbox=props)
+imgplot = plt.imshow(plotme[sliceno],cmap="gray")
+
+###################### FGP_TV #########################################
+# u = FGP_TV(single(u0), 0.05, 100, 1e-04);
+out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.0005,
+                          number_of_iterations=50)
+pars = out2[-2]
+
+reg_output.append(out2)
+
+a=fig.add_subplot(2,3,3)
+
+textstr = out2[-1]
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+         verticalalignment='top', bbox=props)
+imgplot = plt.imshow(reg_output[-1][0][sliceno])
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+        verticalalignment='top', bbox=props)
+imgplot = plt.imshow(reg_output[-1][0][sliceno],cmap="gray")
+
+###################### LLT_model #########################################
+# * u0 = Im + .03*randn(size(Im)); % adding noise
+# [Den] = LLT_model(single(u0), 10, 0.1, 1);
+#Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0); 
+#input, regularization_parameter , time_step, number_of_iterations,
+#                  tolerance_constant, restrictive_Z_smoothing=0
+out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25,
+                          time_step=0.0003,
+                          tolerance_constant=0.0001,
+                          number_of_iterations=300)
+pars = out2[-2]
+
+reg_output.append(out2)
+
+a=fig.add_subplot(2,3,4)
+
+textstr = out2[-1]
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+         verticalalignment='top', bbox=props)
+imgplot = plt.imshow(reg_output[-1][0][sliceno],cmap="gray")
+
+
+# ###################### PatchBased_Regul #########################################
+# # Quick 2D denoising example in Matlab:   
+# #   Im = double(imread('lena_gray_256.tif'))/255;  % loading image
+# #   u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+# #   ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05); 
+
+out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05,
+                       searching_window_ratio=3,
+                       similarity_window_ratio=1,
+                       PB_filtering_parameter=0.08)
+pars = out2[-2]
+reg_output.append(out2)
+
+a=fig.add_subplot(2,3,5)
+
+
+textstr = out2[-1]
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+        verticalalignment='top', bbox=props)
+imgplot = plt.imshow(reg_output[-1][0][sliceno],cmap="gray")
+
+
+# ###################### TGV_PD #########################################
+# # Quick 2D denoising example in Matlab:   
+# #   Im = double(imread('lena_gray_256.tif'))/255;  % loading image
+# #   u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+# #   u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550);
+
+
+out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05,
+                           first_order_term=1.3,
+                           second_order_term=1,
+                           number_of_iterations=550)
+pars = out2[-2]
+reg_output.append(out2)
+
+a=fig.add_subplot(2,3,6)
+
+
+textstr = out2[-1]
+
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+         verticalalignment='top', bbox=props)
+imgplot = plt.imshow(reg_output[-1][0][sliceno],cmap="gray")
+
+
+plt.show()
+
+################################################################################
+##
+##  3D Regularizers
+##
+################################################################################
+##Example:
+## figure;
+## Im = double(imread('lena_gray_256.tif'))/255;  % loading image
+## u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
+## u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+#
+##filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Reconstruction\python\test\reconstruction_example.mha"
+#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Simpleflex\data\head.mha"
 #
+#reader = vtk.vtkMetaImageReader()
+#reader.SetFileName(os.path.normpath(filename))
+#reader.Update()
+##vtk returns 3D images, let's take just the one slice there is as 2D
+#Im = Converter.vtk2numpy(reader.GetOutput())
+#Im = Im.astype('float32')
 ##imgplot = plt.imshow(Im)
 #perc = 0.05
 #u0 = Im + (perc* np.random.normal(size=np.shape(Im)))
 ## map the u0 u0->u0>0
 #f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
 #u0 = f(u0).astype('float32')
+#converter = Converter.numpy2vtkImporter(u0, reader.GetOutput().GetSpacing(),
+#                                        reader.GetOutput().GetOrigin())
+#converter.Update()
+#writer = vtk.vtkMetaImageWriter()
+#writer.SetInputData(converter.GetOutput())
+#writer.SetFileName(r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\noisy_head.mha")
+##writer.Write()
+#
 #
 ### plot 
-#fig = plt.figure()
+#fig3D = plt.figure()
 ##a=fig.add_subplot(3,3,1)
 ##a.set_title('Original')
 ##imgplot = plt.imshow(Im)
+#sliceNo = 32
 #
-#a=fig.add_subplot(2,3,1)
+#a=fig3D.add_subplot(2,3,1)
 #a.set_title('noise')
-#imgplot = plt.imshow(u0)
+#imgplot = plt.imshow(u0.T[sliceNo])
+#
+#reg_output3d = []
 #
-#reg_output = []
 ###############################################################################
 ## Call regularizer
 #
 ######################## SplitBregman_TV #####################################
 ## u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+#
 ##reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV)
-##
+#
 ##out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30,
 ##          #tolerance_constant=1e-4, 
 ##          TV_Penalty=Regularizer.TotalVariationPenalty.l1)
 #
-##out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30,
-##          tolerance_constant=1e-4, 
-##          TV_Penalty=Regularizer.TotalVariationPenalty.l1)
-#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10. )
-#pars = out2[2]
-#reg_output.append(out2)
+#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30,
+#          tolerance_constant=1e-4, 
+#          TV_Penalty=Regularizer.TotalVariationPenalty.l1)
+#
+#
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,3,2)
 #
-#a=fig.add_subplot(2,3,2)
 #
 #textstr = out2[-1]
+#
+#
 ## these are matplotlib.patch.Patch properties
 #props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
 ## place a text box in upper left in axes coords
 #a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
 #        verticalalignment='top', bbox=props)
-#imgplot = plt.imshow(reg_output[-1][0])
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
 #
 ####################### FGP_TV #########################################
 ## u = FGP_TV(single(u0), 0.05, 100, 1e-04);
 #out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.005,
 #                          number_of_iterations=200)
 #pars = out2[-2]
+#reg_output3d.append(out2)
 #
-#reg_output.append(out2)
+#a=fig3D.add_subplot(2,3,2)
 #
-#a=fig.add_subplot(2,3,3)
 #
 #textstr = out2[-1]
 #
+#
 ## these are matplotlib.patch.Patch properties
 #props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
 ## place a text box in upper left in axes coords
 #a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
 #        verticalalignment='top', bbox=props)
-#imgplot = plt.imshow(reg_output[-1][0])
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
 #
 ####################### LLT_model #########################################
 ## * u0 = Im + .03*randn(size(Im)); % adding noise
@@ -129,17 +355,20 @@ def nrmse(im1, im2):
 #                          tolerance_constant=0.0001,
 #                          number_of_iterations=300)
 #pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,3,2)
 #
-#reg_output.append(out2)
 #
-#a=fig.add_subplot(2,3,4)
 #textstr = out2[-1]
+#
+#
 ## these are matplotlib.patch.Patch properties
 #props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
 ## place a text box in upper left in axes coords
 #a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
 #        verticalalignment='top', bbox=props)
-#imgplot = plt.imshow(reg_output[-1][0])
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
 #
 ####################### PatchBased_Regul #########################################
 ## Quick 2D denoising example in Matlab:   
@@ -152,136 +381,6 @@ def nrmse(im1, im2):
 #                          similarity_window_ratio=1,
 #                          PB_filtering_parameter=0.08)
 #pars = out2[-2]
-#reg_output.append(out2)
-#
-#a=fig.add_subplot(2,3,5)
-#
-#
-#textstr = out2[-1]
-#
-## these are matplotlib.patch.Patch properties
-#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-## place a text box in upper left in axes coords
-#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
-#        verticalalignment='top', bbox=props)
-#imgplot = plt.imshow(reg_output[-1][0])
-#
-#
-####################### TGV_PD #########################################
-## Quick 2D denoising example in Matlab:   
-##   Im = double(imread('lena_gray_256.tif'))/255;  % loading image
-##   u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
-##   u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550);
-#
-#
-#out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05,
-#                          first_order_term=1.3,
-#                          second_order_term=1,
-#                          number_of_iterations=550)
-#pars = out2[-2]
-#reg_output.append(out2)
-#
-#a=fig.add_subplot(2,3,6)
-#
-#
-#textstr = out2[-1]
-#
-#
-## these are matplotlib.patch.Patch properties
-#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-## place a text box in upper left in axes coords
-#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
-#        verticalalignment='top', bbox=props)
-#imgplot = plt.imshow(reg_output[-1][0])
-#
-
-###############################################################################
-#
-#  3D Regularizers
-#
-###############################################################################
-#Example:
-# figure;
-# Im = double(imread('lena_gray_256.tif'))/255;  % loading image
-# u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
-# u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
-
-#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Reconstruction\python\test\reconstruction_example.mha"
-filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Simpleflex\data\head.mha"
-
-reader = vtk.vtkMetaImageReader()
-reader.SetFileName(os.path.normpath(filename))
-reader.Update()
-#vtk returns 3D images, let's take just the one slice there is as 2D
-Im = Converter.vtk2numpy(reader.GetOutput())
-Im = Im.astype('float32')
-#imgplot = plt.imshow(Im)
-perc = 0.05
-u0 = Im + (perc* np.random.normal(size=np.shape(Im)))
-# map the u0 u0->u0>0
-f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
-u0 = f(u0).astype('float32')
-converter = Converter.numpy2vtkImporter(u0, reader.GetOutput().GetSpacing(),
-                                        reader.GetOutput().GetOrigin())
-converter.Update()
-writer = vtk.vtkMetaImageWriter()
-writer.SetInputData(converter.GetOutput())
-writer.SetFileName(r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\noisy_head.mha")
-#writer.Write()
-
-
-## plot 
-fig3D = plt.figure(figsize=(20,16))
-
-#a=fig.add_subplot(3,3,1)
-#a.set_title('Original')
-#imgplot = plt.imshow(Im)
-sliceNo = 32
-
-a=fig3D.add_subplot(2,3,1)
-a.set_title('noise')
-imgplot = plt.imshow(u0.T[sliceNo])
-
-reg_output3d = []
-
-##############################################################################
-# Call regularizer
-
-####################### SplitBregman_TV #####################################
-# u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
-
-#reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV)
-
-#out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30,
-#          #tolerance_constant=1e-4, 
-#          TV_Penalty=Regularizer.TotalVariationPenalty.l1)
-
-out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30,
-          tolerance_constant=1e-4, 
-          TV_Penalty=Regularizer.TotalVariationPenalty.l1)
-
-
-pars = out2[-2]
-reg_output3d.append(out2)
-
-a=fig3D.add_subplot(2,3,2)
-
-
-textstr = out2[-1]
-
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
-        verticalalignment='top', bbox=props)
-imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
-
-###################### FGP_TV #########################################
-# u = FGP_TV(single(u0), 0.05, 100, 1e-04);
-#out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.005,
-#                          number_of_iterations=200)
-#pars = out2[-2]
 #reg_output3d.append(out2)
 #
 #a=fig3D.add_subplot(2,3,2)
@@ -296,67 +395,15 @@ imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
 #a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
 #        verticalalignment='top', bbox=props)
 #imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+#
 
-###################### LLT_model #########################################
-# * u0 = Im + .03*randn(size(Im)); % adding noise
-# [Den] = LLT_model(single(u0), 10, 0.1, 1);
-#Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0); 
-#input, regularization_parameter , time_step, number_of_iterations,
-#                  tolerance_constant, restrictive_Z_smoothing=0
-out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25,
-                          time_step=0.0003,
-                          tolerance_constant=0.0001,
-                          number_of_iterations=300)
-pars = out2[-2]
-reg_output3d.append(out2)
-
-a=fig3D.add_subplot(2,3,3)
-
-
-textstr = out2[-1]
-
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
-        verticalalignment='top', bbox=props)
-imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
-
-###################### PatchBased_Regul #########################################
+###################### TGV_PD #########################################
 # Quick 2D denoising example in Matlab:   
 #   Im = double(imread('lena_gray_256.tif'))/255;  % loading image
 #   u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
-#   ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05); 
-
-out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05,
-                          searching_window_ratio=3,
-                          similarity_window_ratio=1,
-                          PB_filtering_parameter=0.08)
-pars = out2[-2]
-reg_output3d.append(out2)
-
-a=fig3D.add_subplot(2,3,4)
-
-
-textstr = out2[-1]
+#   u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550);
 
 
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
-# place a text box in upper left in axes coords
-a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
-        verticalalignment='top', bbox=props)
-imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
-
-
-####################### TGV_PD #########################################
-## Quick 2D denoising example in Matlab:   
-##   Im = double(imread('lena_gray_256.tif'))/255;  % loading image
-##   u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
-##   u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550);
-#
-#
 #out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05,
 #                          first_order_term=1.3,
 #                          second_order_term=1,
@@ -376,5 +423,3 @@ imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
 #a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
 #        verticalalignment='top', bbox=props)
 #imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
-fig3D.savefig('test\\3d.png')
-plt.close(fig3D)
\ No newline at end of file
-- 
cgit v1.2.3


From 01861a7022cb7855bc1a8cd7f8cfd6282690a4f1 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Wed, 25 Oct 2017 16:57:19 +0100
Subject: added to repository

---
 src/Python/test/test_regularizers.py | 412 +++++++++++++++++++++++++++++++++++
 1 file changed, 412 insertions(+)
 create mode 100644 src/Python/test/test_regularizers.py

(limited to 'src/Python/test')

diff --git a/src/Python/test/test_regularizers.py b/src/Python/test/test_regularizers.py
new file mode 100644
index 0000000..27e4ed3
--- /dev/null
+++ b/src/Python/test/test_regularizers.py
@@ -0,0 +1,412 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Aug  4 11:10:05 2017
+
+@author: ofn77899
+"""
+
+#from ccpi.viewer.CILViewer2D import Converter
+#import vtk
+
+import matplotlib.pyplot as plt
+import numpy as np
+import os    
+from enum import Enum
+import timeit
+#from PIL import Image
+#from Regularizer import Regularizer
+from ccpi.imaging.Regularizer import Regularizer
+
+###############################################################################
+#https://stackoverflow.com/questions/13875989/comparing-image-in-url-to-image-in-filesystem-in-python/13884956#13884956
+#NRMSE a normalization of the root of the mean squared error
+#NRMSE is simply 1 - [RMSE / (maxval - minval)]. Where maxval is the maximum
+# intensity from the two images being compared, and respectively the same for
+# minval. RMSE is given by the square root of MSE: 
+# sqrt[(sum(A - B) ** 2) / |A|],
+# where |A| means the number of elements in A. By doing this, the maximum value
+# given by RMSE is maxval.
+
+def nrmse(im1, im2):
+    a, b = im1.shape
+    rmse = np.sqrt(np.sum((im2 - im1) ** 2) / float(a * b))
+    max_val = max(np.max(im1), np.max(im2))
+    min_val = min(np.min(im1), np.min(im2))
+    return 1 - (rmse / (max_val - min_val))
+###############################################################################
+
+###############################################################################
+#
+#  2D Regularizers
+#
+###############################################################################
+#Example:
+# figure;
+# Im = double(imread('lena_gray_256.tif'))/255;  % loading image
+# u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
+# u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+
+
+#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\lena_gray_512.tif"
+filename = r"/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/data/lena_gray_512.tif"
+#filename = r'/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif'
+
+#reader = vtk.vtkTIFFReader()
+#reader.SetFileName(os.path.normpath(filename))
+#reader.Update()
+Im = plt.imread(filename)                     
+#Im = Image.open('/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif')/255
+#img.show()
+Im = np.asarray(Im, dtype='float32')
+
+
+
+
+#imgplot = plt.imshow(Im)
+perc = 0.05
+u0 = Im + (perc* np.random.normal(size=np.shape(Im)))
+# map the u0 u0->u0>0
+f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
+u0 = f(u0).astype('float32')
+
+## plot 
+fig = plt.figure()
+#a=fig.add_subplot(3,3,1)
+#a.set_title('Original')
+#imgplot = plt.imshow(Im)
+
+a=fig.add_subplot(2,3,1)
+a.set_title('noise')
+imgplot = plt.imshow(u0,cmap="gray")
+
+reg_output = []
+##############################################################################
+# Call regularizer
+
+####################### SplitBregman_TV #####################################
+# u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+
+use_object = True
+if use_object:
+    reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV)
+    print (reg.pars)
+    reg.setParameter(input=u0)    
+    reg.setParameter(regularization_parameter=10.)
+    # or 
+    # reg.setParameter(input=u0, regularization_parameter=10., #number_of_iterations=30,
+              #tolerance_constant=1e-4, 
+              #TV_Penalty=Regularizer.TotalVariationPenalty.l1)
+    plotme = reg() [0]
+    pars = reg.pars
+    textstr = reg.printParametersToString() 
+    
+    #out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30,
+              #tolerance_constant=1e-4, 
+    #          TV_Penalty=Regularizer.TotalVariationPenalty.l1)
+    
+#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30,
+#          tolerance_constant=1e-4, 
+#          TV_Penalty=Regularizer.TotalVariationPenalty.l1)
+
+else:
+    out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10. )
+    pars = out2[2]
+    reg_output.append(out2)
+    plotme = reg_output[-1][0]
+    textstr = out2[-1]
+
+a=fig.add_subplot(2,3,2)
+
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+        verticalalignment='top', bbox=props)
+imgplot = plt.imshow(plotme,cmap="gray")
+
+###################### FGP_TV #########################################
+# u = FGP_TV(single(u0), 0.05, 100, 1e-04);
+out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.0005,
+                          number_of_iterations=50)
+pars = out2[-2]
+
+reg_output.append(out2)
+
+a=fig.add_subplot(2,3,3)
+
+textstr = out2[-1]
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+         verticalalignment='top', bbox=props)
+imgplot = plt.imshow(reg_output[-1][0])
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+        verticalalignment='top', bbox=props)
+imgplot = plt.imshow(reg_output[-1][0],cmap="gray")
+
+###################### LLT_model #########################################
+# * u0 = Im + .03*randn(size(Im)); % adding noise
+# [Den] = LLT_model(single(u0), 10, 0.1, 1);
+#Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0); 
+#input, regularization_parameter , time_step, number_of_iterations,
+#                  tolerance_constant, restrictive_Z_smoothing=0
+out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25,
+                          time_step=0.0003,
+                          tolerance_constant=0.0001,
+                          number_of_iterations=300)
+pars = out2[-2]
+
+reg_output.append(out2)
+
+a=fig.add_subplot(2,3,4)
+
+textstr = out2[-1]
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+         verticalalignment='top', bbox=props)
+imgplot = plt.imshow(reg_output[-1][0],cmap="gray")
+
+
+# ###################### PatchBased_Regul #########################################
+# # Quick 2D denoising example in Matlab:   
+# #   Im = double(imread('lena_gray_256.tif'))/255;  % loading image
+# #   u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+# #   ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05); 
+
+out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05,
+                       searching_window_ratio=3,
+                       similarity_window_ratio=1,
+                       PB_filtering_parameter=0.08)
+pars = out2[-2]
+reg_output.append(out2)
+
+a=fig.add_subplot(2,3,5)
+
+
+textstr = out2[-1]
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+        verticalalignment='top', bbox=props)
+imgplot = plt.imshow(reg_output[-1][0],cmap="gray")
+
+
+# ###################### TGV_PD #########################################
+# # Quick 2D denoising example in Matlab:   
+# #   Im = double(imread('lena_gray_256.tif'))/255;  % loading image
+# #   u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+# #   u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550);
+
+
+out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05,
+                           first_order_term=1.3,
+                           second_order_term=1,
+                           number_of_iterations=550)
+pars = out2[-2]
+reg_output.append(out2)
+
+a=fig.add_subplot(2,3,6)
+
+
+textstr = out2[-1]
+
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+         verticalalignment='top', bbox=props)
+imgplot = plt.imshow(reg_output[-1][0],cmap="gray")
+
+
+plt.show()
+
+################################################################################
+##
+##  3D Regularizers
+##
+################################################################################
+##Example:
+## figure;
+## Im = double(imread('lena_gray_256.tif'))/255;  % loading image
+## u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
+## u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+#
+##filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Reconstruction\python\test\reconstruction_example.mha"
+#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Simpleflex\data\head.mha"
+#
+#reader = vtk.vtkMetaImageReader()
+#reader.SetFileName(os.path.normpath(filename))
+#reader.Update()
+##vtk returns 3D images, let's take just the one slice there is as 2D
+#Im = Converter.vtk2numpy(reader.GetOutput())
+#Im = Im.astype('float32')
+##imgplot = plt.imshow(Im)
+#perc = 0.05
+#u0 = Im + (perc* np.random.normal(size=np.shape(Im)))
+## map the u0 u0->u0>0
+#f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
+#u0 = f(u0).astype('float32')
+#converter = Converter.numpy2vtkImporter(u0, reader.GetOutput().GetSpacing(),
+#                                        reader.GetOutput().GetOrigin())
+#converter.Update()
+#writer = vtk.vtkMetaImageWriter()
+#writer.SetInputData(converter.GetOutput())
+#writer.SetFileName(r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\noisy_head.mha")
+##writer.Write()
+#
+#
+### plot 
+#fig3D = plt.figure()
+##a=fig.add_subplot(3,3,1)
+##a.set_title('Original')
+##imgplot = plt.imshow(Im)
+#sliceNo = 32
+#
+#a=fig3D.add_subplot(2,3,1)
+#a.set_title('noise')
+#imgplot = plt.imshow(u0.T[sliceNo])
+#
+#reg_output3d = []
+#
+###############################################################################
+## Call regularizer
+#
+######################## SplitBregman_TV #####################################
+## u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+#
+##reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV)
+#
+##out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30,
+##          #tolerance_constant=1e-4, 
+##          TV_Penalty=Regularizer.TotalVariationPenalty.l1)
+#
+#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30,
+#          tolerance_constant=1e-4, 
+#          TV_Penalty=Regularizer.TotalVariationPenalty.l1)
+#
+#
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,3,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+#        verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+#
+####################### FGP_TV #########################################
+## u = FGP_TV(single(u0), 0.05, 100, 1e-04);
+#out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.005,
+#                          number_of_iterations=200)
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,3,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+#        verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+#
+####################### LLT_model #########################################
+## * u0 = Im + .03*randn(size(Im)); % adding noise
+## [Den] = LLT_model(single(u0), 10, 0.1, 1);
+##Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0); 
+##input, regularization_parameter , time_step, number_of_iterations,
+##                  tolerance_constant, restrictive_Z_smoothing=0
+#out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25,
+#                          time_step=0.0003,
+#                          tolerance_constant=0.0001,
+#                          number_of_iterations=300)
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,3,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+#        verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+#
+####################### PatchBased_Regul #########################################
+## Quick 2D denoising example in Matlab:   
+##   Im = double(imread('lena_gray_256.tif'))/255;  % loading image
+##   u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+##   ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05); 
+#
+#out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05,
+#                          searching_window_ratio=3,
+#                          similarity_window_ratio=1,
+#                          PB_filtering_parameter=0.08)
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,3,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+#        verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+#
+
+###################### TGV_PD #########################################
+# Quick 2D denoising example in Matlab:   
+#   Im = double(imread('lena_gray_256.tif'))/255;  % loading image
+#   u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+#   u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550);
+
+
+#out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05,
+#                          first_order_term=1.3,
+#                          second_order_term=1,
+#                          number_of_iterations=550)
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,3,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+#        verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
-- 
cgit v1.2.3