Simple demo ip (#79)

* Removed regularizers

* DemoIP running. MC Lipschitz incorrect

* Fixed MC Lipschitz

* Removed old import

5 files changed, 107 insertions, 540 deletions

diff --git a/Wrappers/Python/ccpi/astra/astra_ops.py b/Wrappers/Python/ccpi/astra/astra_ops.py
index d8c2f4b..45b8238 100644
--- a/Wrappers/Python/ccpi/astra/astra_ops.py
+++ b/Wrappers/Python/ccpi/astra/astra_ops.py
@@ -71,7 +71,9 @@ class AstraProjectorSimple(Operator):
     
     def create_image_data(self):
         inputsize = self.size()[1]
-        return DataContainer(numpy.random.randn(inputsize[0],inputsize[1]))
+        return DataContainer(numpy.random.randn(inputsize[0],
+                                                inputsize[1]))
+
 
 class AstraProjectorMC(Operator):
     """ASTRA Multichannel projector"""
@@ -93,7 +95,7 @@ class AstraProjectorMC(Operator):
                                         device=device)
                 
         # Initialise empty for singular value.
-        self.s1 = 50
+        self.s1 = None
     
     def direct(self, IM):
         self.fp.set_input(IM)
@@ -122,80 +124,9 @@ class AstraProjectorMC(Operator):
                  (self.volume_geometry.voxel_num_x, \
                   self.volume_geometry.voxel_num_y) )
     
-
-class AstraProjector(Operator):
-    """A simple 2D/3D parallel/fan beam projection/backprojection class based 
-    on ASTRA toolbox"""
-    def __init__(self, DetWidth, DetectorsDim, SourceOrig, OrigDetec, 
-                 AnglesVec, ObjSize, projtype, device):
-        super(AstraProjector, self).__init__()
-        self.DetectorsDim = DetectorsDim
-        self.AnglesVec = AnglesVec
-        self.ProjNumb = len(AnglesVec)
-        self.ObjSize = ObjSize
-        if projtype == 'parallel':
-            self.proj_geom = astra.create_proj_geom('parallel', DetWidth, 
-                                                    DetectorsDim, AnglesVec)
-        elif projtype == 'fanbeam':
-            self.proj_geom = astra.create_proj_geom('fanflat', DetWidth,
-                                                    DetectorsDim, AnglesVec, 
-                                                    SourceOrig, OrigDetec)
-        else:
-            print ("Please select for projtype between 'parallel' and 'fanbeam'")
-        self.vol_geom = astra.create_vol_geom(ObjSize, ObjSize)
-        if device == 'cpu':
-            self.proj_id = astra.create_projector('line', self.proj_geom, 
-                                                  self.vol_geom) # for CPU
-            self.device = 1
-        elif device == 'gpu':
-            self.proj_id = astra.create_projector('cuda', self.proj_geom, 
-                                                  self.vol_geom) # for GPU
-            self.device = 0
-        else:
-            print ("Select between 'cpu' or 'gpu' for device")
-        self.s1 = None
-    def direct(self, IM):
-        """Applying forward projection to IM [2D or 3D array]"""
-        if numpy.ndim(IM.as_array()) == 3:
-            slices = numpy.size(IM.as_array(),numpy.ndim(IM.as_array())-1)
-            DATA = numpy.zeros((self.ProjNumb,self.DetectorsDim,slices), 
-                               'float32')
-            for i in range(0,slices):
-                sinogram_id, DATA[:,:,i] = \
-                    astra.create_sino(IM[:,:,i].as_array(), self.proj_id)
-                astra.data2d.delete(sinogram_id)
-            astra.data2d.delete(self.proj_id)
-        else:
-            sinogram_id, DATA = astra.create_sino(IM.as_array(), self.proj_id)
-            astra.data2d.delete(sinogram_id)
-            astra.data2d.delete(self.proj_id)
-        return AcquisitionData(DATA)
-    def adjoint(self, DATA):
-        """Applying backprojection to DATA [2D or 3D]"""
-        if numpy.ndim(DATA) == 3:
-           slices = numpy.size(DATA.as_array(),numpy.ndim(DATA.as_array())-1)
-           IM = numpy.zeros((self.ObjSize,self.ObjSize,slices), 'float32')
-           for i in range(0,slices):
-               rec_id, IM[:,:,i] = \
-                   astra.create_backprojection(DATA[:,:,i].as_array(), 
-                                               self.proj_id)
-               astra.data2d.delete(rec_id)
-           astra.data2d.delete(self.proj_id)
-        else:
-            rec_id, IM = astra.create_backprojection(DATA.as_array(), 
-                                                     self.proj_id)        
-            astra.data2d.delete(rec_id)
-            astra.data2d.delete(self.proj_id)
-        return ImageData(IM)
-    
-    def delete(self):
-        astra.data2d.delete(self.proj_id)
-    
-    def get_max_sing_val(self):
-        self.s1, sall, svec = PowerMethodNonsquare(self,10)
-        return self.s1
+    def create_image_data(self):
+        inputsize = self.size()[1]
+        return DataContainer(numpy.random.randn(self.volume_geometry.channels,
+                                                inputsize[0],
+                                                inputsize[1]))
     
-    def size(self):
-        return ( (self.AnglesVec.size, self.DetectorsDim), \
-                 (self.ObjSize, self.ObjSize) )
-
diff --git a/Wrappers/Python/ccpi/reconstruction/ops.py b/Wrappers/Python/ccpi/reconstruction/ops.py
index fa30b0c..2d8e90d 100644
--- a/Wrappers/Python/ccpi/reconstruction/ops.py
+++ b/Wrappers/Python/ccpi/reconstruction/ops.py
@@ -134,7 +134,7 @@ class FiniteDiff2D(Operator):
         return self.s1
 
 
-def PowerMethodNonsquare(op,numiters):
+def PowerMethodNonsquareOld(op,numiters):
     # Initialise random
     # Jakob's
     #inputsize = op.size()[1]
@@ -173,6 +173,34 @@ def PowerMethodNonsquare(op,numiters):
 #        x0 = (1.0/x1norm)*x1
 #    return s, x0
     
+
+def PowerMethodNonsquare(op,numiters):
+    # Initialise random
+    # Jakob's
+    #inputsize = op.size()[1]
+    #x0 = ImageContainer(numpy.random.randn(*inputsize)
+    # Edo's
+    #vg = ImageGeometry(voxel_num_x=inputsize[0],
+    #                   voxel_num_y=inputsize[1], 
+    #                   voxel_num_z=inputsize[2])
+    #
+    #x0 = ImageData(geometry = vg, dimension_labels=['vertical','horizontal_y','horizontal_x'])
+    #print (x0)
+    #x0.fill(numpy.random.randn(*x0.shape))
+    
+    x0 = op.create_image_data()
+    
+    s = numpy.zeros(numiters)
+    # Loop
+    for it in numpy.arange(numiters):
+        x1 = op.adjoint(op.direct(x0))
+        x1norm = numpy.sqrt((x1**2).sum())
+        #print ("x0 **********" ,x0)
+        #print ("x1 **********" ,x1)
+        s[it] = (x1*x0).sum() / (x0*x0).sum()
+        x0 = (1.0/x1norm)*x1
+    return numpy.sqrt(s[-1]), numpy.sqrt(s), x0
+
 class CCPiProjectorSimple(Operator):
     """ASTRA projector modified to use DataSet and geometry."""
     def __init__(self, geomv, geomp):
@@ -225,4 +253,4 @@ class CCPiProjectorSimple(Operator):
                        #.subset(['horizontal_x','horizontal_y','vertical'])
         print (x0)
         x0.fill(numpy.random.randn(*x0.shape))
-        return x0
+        return x0
\ No newline at end of file
diff --git a/Wrappers/Python/wip/DemoRecIP.py b/Wrappers/Python/wip/DemoRecIP.py
index 483a202..442e40e 100644
--- a/Wrappers/Python/wip/DemoRecIP.py
+++ b/Wrappers/Python/wip/DemoRecIP.py
@@ -7,50 +7,19 @@ Reading multi-channel data and reconstruction using FISTA modular
 import numpy as np
 import matplotlib.pyplot as plt
 
-import sys
+#import sys
 #sys.path.append('../../../data/')
 from read_IPdata import read_IPdata
 
-from ccpi.reconstruction.astra_ops import AstraProjector
-from ccpi.reconstruction.funcs import Norm2sq , BaseFunction
+from ccpi.astra.astra_ops import AstraProjectorSimple, AstraProjectorMC
+from ccpi.reconstruction.funcs import Norm2sq, Norm1, BaseFunction
 from ccpi.reconstruction.algs import FISTA
 #from ccpi.reconstruction.funcs import BaseFunction
 
-from ccpi.framework import VolumeData, SinogramData
-
-
-from ccpi.filters.cpu_regularizers_boost import SplitBregman_TV , FGP_TV ,\
-                                                 LLT_model, PatchBased_Regul ,\
-                                                 TGV_PD
-
-
-# TV regularization class for FGP_TV method
-class TV_FGP(BaseFunction):
-    def __init__(self,lambdaReg,iterationsTV):
-        # set parameters
-        self.lambdaReg = lambdaReg
-        self.iterationsTV = iterationsTV
-        super(TV_FGP, self).__init__()
-    def fun(self,x):
-        # function to calculate energy from utils can be used here
-        return 0
-    def prox(self,x,Lipshitz):
-        pars = {'algorithm' : FGP_TV , \
-                'input' : x.as_array(),
-                'regularization_parameter':self.lambdaReg*Lipshitz, \
-                'number_of_iterations' :self.iterationsTV ,\
-                'tolerance_constant':1e-4,\
-                'TV_penalty': 0}
-        
-        out = FGP_TV (pars['input'], 
-                      pars['regularization_parameter'],
-                      pars['number_of_iterations'],
-                      pars['tolerance_constant'], 
-                      pars['TV_penalty'])
-        return out[0]
+from ccpi.framework import ImageData, AcquisitionData, AcquisitionGeometry, ImageGeometry
 
 # read IP paper data into a dictionary
-dataDICT = read_IPdata()
+dataDICT = read_IPdata('..\..\..\data\IP_data70channels.mat')
 
 # Set ASTRA Projection-backprojection class (fan-beam geometry)
 DetWidth = dataDICT.get('im_size')[0] * dataDICT.get('det_width')[0] / \
@@ -61,43 +30,81 @@ OrigDetec = dataDICT.get('im_size')[0] * \
             (dataDICT.get('src_to_det')[0] - dataDICT.get('src_to_rotc')[0]) /\
             dataDICT.get('dom_width')[0]
 
-# Set up ASTRA projector
-Aop = AstraProjector(DetWidth, dataDICT.get('detectors_numb')[0], SourceOrig, 
-                     OrigDetec, (np.pi/180)*dataDICT.get('theta')[0], 
-                     dataDICT.get('im_size')[0],'fanbeam','gpu') 
-# initiate a class object
+N = dataDICT.get('im_size')[0]
 
-sino = dataDICT.get('data_norm')[0][:,:,34] # select mid-channel 
-b = SinogramData(sino)
-# Try forward and backprojection
-#backprj = Aop.adjoint(b)
+vg = ImageGeometry(voxel_num_x=dataDICT.get('im_size')[0],
+                   voxel_num_y=dataDICT.get('im_size')[0],
+                   channels=1)
+
+pg = AcquisitionGeometry('cone',
+                         '2D',
+                         angles=(np.pi/180)*dataDICT.get('theta')[0],
+                         pixel_num_h=dataDICT.get('detectors_numb')[0],
+                         pixel_size_h=DetWidth,
+                         dist_source_center=SourceOrig, 
+                         dist_center_detector=OrigDetec,
+                         channels=1)
 
-# Create least squares object instance with projector and data.
-f = Norm2sq(Aop,b,c=0.5)
+
+sino = dataDICT.get('data_norm')[0][:,:,34] # select mid-channel 
+b = AcquisitionData(sino,geometry=pg)
 
 # Initial guess
-x_init = VolumeData(np.zeros((dataDICT.get('im_size')[0],
-                              dataDICT.get('im_size')[0])))
+x_init = ImageData(np.zeros((N, N)),geometry=vg)
+
+
+
+
+
+Aop = AstraProjectorSimple(vg,pg,'gpu')
+f = Norm2sq(Aop,b,c=0.5)
 
 # Run FISTA for least squares without regularization
-opt = {'tol': 1e-4, 'iter': 50}
+opt = {'tol': 1e-4, 'iter': 10}
 x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt)
 
 plt.imshow(x_fista0.array)
 plt.show()
 
-"""
 # Now least squares plus 1-norm regularization
-#lam = 1
-#g0 = Norm1(lam)
-# using FGP_TV regularizer
-
-g0 = TV_FGP(lambdaReg = 0.01,iterationsTV=50)
+g1 = Norm1(10)
 
 # Run FISTA for least squares plus 1-norm function.
-opt = {'tol': 1e-4, 'iter': 50}
-x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0, opt)
+x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g1, opt)
 
 plt.imshow(x_fista1.array)
 plt.show()
-"""
\ No newline at end of file
+
+# Multiple channels
+sino_mc = dataDICT.get('data_norm')[0][:,:,32:37] # select mid-channel 
+
+vg_mc = ImageGeometry(voxel_num_x=dataDICT.get('im_size')[0],
+                   voxel_num_y=dataDICT.get('im_size')[0],
+                   channels=5)
+
+pg_mc = AcquisitionGeometry('cone',
+                         '2D',
+                         angles=(np.pi/180)*dataDICT.get('theta')[0],
+                         pixel_num_h=dataDICT.get('detectors_numb')[0],
+                         pixel_size_h=DetWidth,
+                         dist_source_center=SourceOrig, 
+                         dist_center_detector=OrigDetec,
+                         channels=5)
+
+b_mc = AcquisitionData(np.transpose(sino_mc,(2,0,1)),
+                       geometry=pg_mc, 
+                       dimension_labels=("channel","angle","horizontal"))
+
+# ASTRA operator using volume and sinogram geometries
+Aop_mc = AstraProjectorMC(vg_mc, pg_mc, 'gpu')
+
+f_mc = Norm2sq(Aop_mc,b_mc,c=0.5)
+
+# Initial guess
+x_init_mc = ImageData(np.zeros((5, N, N)),geometry=vg_mc)
+
+
+x_fista0_mc, it0_mc, timing0_mc, criter0_mc = FISTA(x_init_mc, f_mc, None, opt)
+
+plt.imshow(x_fista0_mc.as_array()[4,:,:])
+plt.show()
\ No newline at end of file
diff --git a/Wrappers/Python/wip/regularizers.py b/Wrappers/Python/wip/regularizers.py
deleted file mode 100644
index 04ac3aa..0000000
--- a/Wrappers/Python/wip/regularizers.py
+++ /dev/null
@@ -1,399 +0,0 @@
-# -*- coding: utf-8 -*-
-#   This work is part of the Core Imaging Library developed by
-#   Visual Analytics and Imaging System Group of the Science Technology
-#   Facilities Council, STFC
-
-#   Copyright 2018 Edoardo Pasca
-
-#   Licensed under the Apache License, Version 2.0 (the "License");
-#   you may not use this file except in compliance with the License.
-#   You may obtain a copy of the License at
-
-#       http://www.apache.org/licenses/LICENSE-2.0
-
-#   Unless required by applicable law or agreed to in writing, software
-#   distributed under the License is distributed on an "AS IS" BASIS,
-#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-#   See the License for the specific language governing permissions and
-#   limitations under the License.
-
-from __future__ import print_function
-import matplotlib.pyplot as plt
-import numpy as np
-import os    
-from ccpi.framework import DataSetProcessor, DataSetProcessor23D , DataSet
-from ccpi.filters.cpu_regularizers_boost import SplitBregman_TV , FGP_TV ,\
-                                                LLT_model, PatchBased_Regul ,\
-                                                TGV_PD
-#from ccpi.filters.cpu_regularizers_cython import some
-
-try:
-    from ccpi.filters import gpu_regularizers as gpu
-    class PatchBasedRegGPU(DataSetProcessor23D):
-        '''Regularizers DataSetProcessor for PatchBasedReg
-        
-        
-        '''
-        
-        def __init__(self):
-            attributes = {'regularization_parameter':None, 
-                          'searching_window_ratio': None, 
-                          'similarity_window_ratio': None, 
-                          'PB_filtering_parameter': None
-                      }
-            super(PatchBasedRegGPU, self).__init__(**attributes)
-                
-            
-        def process(self):
-            '''Executes the processor
-                
-            '''
-            dsi = self.getInput()
-            out = gpu.NML (dsi.as_array(), 
-                             self.searching_window_ratio,
-                             self.similarity_window_ratio,
-                             self.regularization_parameter,
-                             self.PB_filtering_parameter)  
-            y = DataSet( out , False )
-            return y
-    class Diff4thHajiaboli(DataSetProcessor23D):
-        '''Regularizers DataSetProcessor for PatchBasedReg
-        
-        
-        '''
-        
-        def __init__(self):
-            attributes = {'regularization_parameter':None, 
-                          'searching_window_ratio': None, 
-                          'similarity_window_ratio': None, 
-                          'PB_filtering_parameter': None
-                      }
-            super(Diff4thHajiaboli, self).__init__(self, **attributes)
-                
-            
-        def process(self):
-            '''Executes the processor
-                
-            '''
-            dsi = self.getInput()
-            out = gpu.Diff4thHajiaboli (dsi.as_array(), 
-                                        self.regularization_parameter, 
-                                        self.number_of_iterations,
-                                        self.edge_preserving_parameter)  
-            y = DataSet( out , False )
-            return y
-
-except ImportError as ie:
-    print (ie)
-    
-
-
-class SBTV(DataSetProcessor23D):
-    '''Regularizers DataSetProcessor
-    '''
-    
-    def __init__(self):
-        attributes = {'regularization_parameter':None, 
-                  'number_of_iterations': 35, 
-                  'tolerance_constant': 0.0001, 
-                  'TV_penalty':0
-                  }
-        super(SBTV , self).__init__(**attributes)
-            
-            
-    def process(self):
-        '''Executes the processor
-        
-        '''
-    
-        dsi = self.getInput()
-        out = SplitBregman_TV (dsi.as_array(), 
-                               self.regularization_parameter,
-                               self.number_of_iterations,
-                               self.tolerance_constant,
-                               self.TV_penalty)  
-        y = DataSet( out[0] , False )
-        return y
-    
-class FGPTV(DataSetProcessor23D):
-    '''Regularizers DataSetProcessor
-    '''
-    
-    def __init__(self):
-        attributes = {'regularization_parameter':None, 
-                  'number_of_iterations': 35, 
-                  'tolerance_constant': 0.0001, 
-                  'TV_penalty':0
-                  }
-        super(FGPTV, self).__init__(**attributes)
-            
-        
-    def process(self):
-        '''Executes the processor
-            
-        '''
-        dsi = self.getInput()
-        out = FGP_TV (dsi.as_array(), 
-                      self.regularization_parameter,
-                      self.number_of_iterations,
-                      self.tolerance_constant,
-                      self.TV_penalty)  
-        y = DataSet( out[0] , False )
-        return y
-    
-class LLT(DataSetProcessor23D):
-    '''Regularizers DataSetProcessor for LLT_model
-    
-    
-    '''
-    
-    def __init__(self):
-        attributes = {'regularization_parameter':None, 
-                      'time_step': 0.0001, 
-                      'number_of_iterations': 35, 
-                      'tolerance_constant': 0, 
-                      'restrictive_Z_smoothing': None 
-                  }
-        super(LLT, self).__init__(**attributes)
-            
-        
-    def process(self):
-        '''Executes the processor
-            
-        '''
-        dsi = self.getInput()
-        out = LLT_model (dsi.as_array(), 
-                         self.regularization_parameter,
-                         self.time_step,
-                         self.number_of_iterations,
-                         self.tolerance_constant,
-                         self.restrictive_Z_smoothing)  
-        y = DataSet( out[0] , False )
-        return y
-    
-class PatchBasedReg(DataSetProcessor23D):
-    '''Regularizers DataSetProcessor for PatchBasedReg
-    
-    
-    '''
-    
-    def __init__(self):
-        attributes = {'regularization_parameter':None, 
-                      'searching_window_ratio': None, 
-                      'similarity_window_ratio': None, 
-                      'PB_filtering_parameter': None
-                  }
-        super(PatchBasedReg, self).__init__(**attributes)
-            
-        
-    def process(self):
-        '''Executes the processor
-            
-        '''
-        dsi = self.getInput()
-        out = PatchBased_Regul (dsi.as_array(), 
-                         self.regularization_parameter,
-                         self.searching_window_ratio,
-                         self.similarity_window_ratio,
-                         self.PB_filtering_parameter)  
-        y = DataSet( out[0] , False )
-        return y
-    
-class TGVPD(DataSetProcessor23D):
-    '''Regularizers DataSetProcessor for PatchBasedReg
-    
-    
-    '''
-    
-    def __init__(self,**kwargs):
-        attributes = {'regularization_parameter':None, 
-                      'first_order_term': None, 
-                      'second_order_term': None, 
-                      'number_of_iterations': None
-                  }
-        for key, value in kwargs.items():
-            if key in attributes.keys():
-                attributes[key] = value
-                
-        super(TGVPD, self).__init__(**attributes)
-            
-        
-    def process(self):
-        '''Executes the processor
-            
-        '''
-        dsi = self.getInput()
-        if dsi.number_of_dimensions == 2:
-            out = TGV_PD(dsi.as_array(), 
-                         self.regularization_parameter,
-                         self.first_order_term, 
-                         self.second_order_term , 
-                         self.number_of_iterations)
-            y = DataSet( out[0] , False )
-        elif len(np.shape(input)) == 3:
-            #assuming it's 3D
-            # run independent calls on each slice
-            out3d = dsi.as_array().copy()
-            for i in range(np.shape(dsi.as_array())[0]):
-                out = TGV_PD(dsi.as_array()[i], 
-                                     self.regularization_parameter,
-                                     self.first_order_term, 
-                                     self.second_order_term , 
-                                     self.number_of_iterations)
-                # copy the result in the 3D image
-                out3d[i] = out[0].copy()
-            y = DataSet (out3d , False)
-        return y
-    
-## self contained test
-if __name__ == '__main__':
-    filename = os.path.join(".." , ".." , ".." , ".." , 
-                            "CCPi-FISTA_Reconstruction", "data" ,
-                            "lena_gray_512.tif")
-    Im = plt.imread(filename)                     
-    Im = np.asarray(Im, dtype='float32')
-
-    Im = Im/255
-
-    perc = 0.075
-    u0 = Im + np.random.normal(loc = Im ,
-                                  scale = perc * Im , 
-                                  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')
-
-    lena = DataSet(u0, False, ['X','Y'])
-    
-    ## plot 
-    fig = plt.figure()
-    
-    a=fig.add_subplot(2,3,1)
-    a.set_title('noise')
-    imgplot = plt.imshow(u0#,cmap="gray"
-                         )
-    
-    reg_output = []
-    ##############################################################################
-    # Call regularizer
-    
-  
-    reg3 = SBTV()
-    reg3.number_of_iterations = 40
-    reg3.tolerance_constant = 0.0001
-    reg3.regularization_parameter = 15
-    reg3.TV_penalty = 0
-    reg3.setInput(lena)
-    dataprocessoroutput = reg3.getOutput()
-    
-    
-    # plot
-    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, 'SBTV', transform=a.transAxes, fontsize=14,
-            verticalalignment='top', bbox=props)
-    imgplot = plt.imshow(dataprocessoroutput.as_array(),\
-                         #cmap="gray"
-                         )
-    ##########################################################################
-    
-    reg4 = FGPTV()
-    reg4.number_of_iterations = 200
-    reg4.tolerance_constant = 1e-4
-    reg4.regularization_parameter = 0.05
-    reg4.TV_penalty = 0
-    reg4.setInput(lena)
-    dataprocessoroutput2 = reg4.getOutput()
-   
-    # plot
-    a=fig.add_subplot(2,3,3)    
-    # 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, 'FGPTV', transform=a.transAxes, fontsize=14,
-            verticalalignment='top', bbox=props)
-    imgplot = plt.imshow(dataprocessoroutput2.as_array(),\
-                         #cmap="gray"
-                         )
-    
-    ###########################################################################
-    reg6 = LLT()
-    reg6.regularization_parameter = 5
-    reg6.time_step = 0.00035
-    reg6.number_of_iterations = 350
-    reg6.tolerance_constant = 0.0001
-    reg6.restrictive_Z_smoothing = 0
-    reg6.setInput(lena)
-    llt = reg6.getOutput()
-    # plot
-    a=fig.add_subplot(2,3,4)    
-    # 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, 'LLT', transform=a.transAxes, fontsize=14,
-            verticalalignment='top', bbox=props)
-    imgplot = plt.imshow(llt.as_array(),\
-                         #cmap="gray"
-                         )
-    ###########################################################################
-    
-    reg7 = PatchBasedReg()
-    reg7.regularization_parameter = 0.05
-    reg7.searching_window_ratio = 3
-    reg7.similarity_window_ratio = 1
-    reg7.PB_filtering_parameter = 0.06
-    reg7.setInput(lena)
-    pbr = reg7.getOutput()
-    # plot
-    a=fig.add_subplot(2,3,5)    
-    # 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, 'PatchBasedReg', transform=a.transAxes, fontsize=14,
-            verticalalignment='top', bbox=props)
-    imgplot = plt.imshow(pbr.as_array(),\
-                         #cmap="gray"
-                         )
-    ###########################################################################
-    
-    reg5 = TGVPD()
-    reg5.regularization_parameter = 0.07
-    reg5.first_order_term = 1.3
-    reg5.second_order_term = 1
-    reg5.number_of_iterations = 550
-    reg5.setInput(lena)
-    tgvpd = reg5.getOutput()
-    # plot
-    a=fig.add_subplot(2,3,6)    
-    # 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, 'TGVPD', transform=a.transAxes, fontsize=14,
-            verticalalignment='top', bbox=props)
-    imgplot = plt.imshow(tgvpd.as_array(),\
-                         #cmap="gray"
-                         )
-    if False:
-        #reg4.input = None
-        reg5 = FGPTV()
-        reg5.number_of_iterations = 350
-        reg5.tolerance_constant = 0.01
-        reg5.regularization_parameter = 40
-        reg5.TV_penalty = 0
-        reg5.setInputProcessor(reg3)
-        chain = reg5.process()
-       
-        a=fig.add_subplot(2,3,6)
-        
-        
-        # 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, 'SBTV + FGPTV', transform=a.transAxes, fontsize=14,
-                verticalalignment='top', bbox=props)
-        imgplot = plt.imshow(chain.as_array(),\
-                             #cmap="gray"
-                             )
-    plt.show()
\ No newline at end of file
diff --git a/data/read_IPdata.py b/data/read_IPdata.py
index 2e6a525..a7565d7 100644
--- a/data/read_IPdata.py
+++ b/data/read_IPdata.py
@@ -7,9 +7,9 @@ from scipy import io
 import numpy as np
 from collections import defaultdict
 
-def read_IPdata():
+def read_IPdata(datafile):
     # read data from mat file (specify the location)
-    alldata = io.loadmat('IP_data70channels.mat')
+    alldata = io.loadmat(datafile)
     data_raw = alldata.get('Data_raw') # here is raw projection data
     Phantom_ideal = alldata.get('Phantom_ideal') # here is 70 channels ideal phantom
     Photon_flux = alldata.get('Photon_flux') # photon flux for normalization