Merge branch 'demos' into update_fix_test

35 files changed, 2101 insertions, 1284 deletions

diff --git a/Wrappers/Python/ccpi/io/NEXUSDataReader.py b/Wrappers/Python/ccpi/io/NEXUSDataReader.py
new file mode 100644
index 0000000..cf67e27
--- /dev/null
+++ b/Wrappers/Python/ccpi/io/NEXUSDataReader.py
@@ -0,0 +1,143 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Apr  3 10:30:25 2019
+
+@author: evelina
+"""
+
+
+import numpy
+import os
+from ccpi.framework import AcquisitionData, AcquisitionGeometry, ImageData, ImageGeometry
+
+
+h5pyAvailable = True
+try:
+    import h5py 
+except:
+    h5pyAvailable = False
+
+
+class NEXUSDataReader(object):
+    
+    def __init__(self,
+                 **kwargs):
+        
+        '''
+        Constructor
+        
+        Input:
+            
+            nexus_file      full path to NEXUS file
+        '''
+        
+        self.nexus_file = kwargs.get('nexus_file', None)
+        
+        if self.nexus_file is not None:
+            self.set_up(nexus_file = self.nexus_file,
+                        roi = self.roi,
+                        binning = self.binning)
+            
+    def set_up(self, 
+               nexus_file = None):
+        
+        self.nexus_file = nexus_file
+        
+        # check that h5py library is installed
+        if (h5pyAvailable == False):
+            raise Exception('h5py is not available, cannot load NEXUS files.')
+            
+        if self.nexus_file == None:
+            raise Exception('Path to nexus file is required.')
+        
+        # check if nexus file exists
+        if not(os.path.isfile(self.nexus_file)):
+            raise Exception('File\n {}\n does not exist.'.format(self.nexus_file))  
+        
+    def load_data(self):
+        
+        '''
+        Parse NEXUS file and returns either ImageData or Acquisition Data 
+        depending on file content
+        '''
+        
+        try:
+            with h5py.File(self.nexus_file,'r') as file:
+                
+                if (file.attrs['creator'] != 'NEXUSDataWriter.py'):
+                    raise Exception('We can parse only files created by NEXUSDataWriter.py')
+                
+                ds_data = file['entry1/tomo_entry/data/data']
+                data = numpy.array(ds_data, dtype = 'float32')
+                
+                dimension_labels = []
+                
+                for i in range(data.ndim):
+                    dimension_labels.append(ds_data.attrs['dim{}'.format(i)])
+                
+                if ds_data.attrs['data_type'] == 'ImageData':
+                    self._geometry = ImageGeometry(voxel_num_x = ds_data.attrs['voxel_num_x'],
+                                                   voxel_num_y = ds_data.attrs['voxel_num_y'],
+                                                   voxel_num_z = ds_data.attrs['voxel_num_z'],
+                                                   voxel_size_x = ds_data.attrs['voxel_size_x'],
+                                                   voxel_size_y = ds_data.attrs['voxel_size_y'],
+                                                   voxel_size_z = ds_data.attrs['voxel_size_z'],
+                                                   center_x = ds_data.attrs['center_x'],
+                                                   center_y = ds_data.attrs['center_y'],
+                                                   center_z = ds_data.attrs['center_z'],
+                                                   channels = ds_data.attrs['channels'])
+                    
+                    return ImageData(array = data,
+                                     deep_copy = False,
+                                     geometry = self._geometry,
+                                     dimension_labels = dimension_labels)
+                    
+                else:   # AcquisitionData
+                    self._geometry = AcquisitionGeometry(geom_type = ds_data.attrs['geom_type'],
+                                                         dimension = ds_data.attrs['dimension'],
+                                                         dist_source_center = ds_data.attrs['dist_source_center'],
+                                                         dist_center_detector = ds_data.attrs['dist_center_detector'],
+                                                         pixel_num_h = ds_data.attrs['pixel_num_h'],
+                                                         pixel_size_h = ds_data.attrs['pixel_size_h'],
+                                                         pixel_num_v = ds_data.attrs['pixel_num_v'],
+                                                         pixel_size_v = ds_data.attrs['pixel_size_v'],
+                                                         channels = ds_data.attrs['channels'],
+                                                         angles = numpy.array(file['entry1/tomo_entry/data/rotation_angle'], dtype = 'float32'))
+                                                         #angle_unit = file['entry1/tomo_entry/data/rotation_angle'].attrs['units'])
+                                             
+                    return AcquisitionData(array = data,
+                                           deep_copy = False,
+                                           geometry = self._geometry,
+                                           dimension_labels = dimension_labels)
+                    
+        except:
+            print("Error reading nexus file")
+            raise
+                
+    def get_geometry(self):
+        
+        '''
+        Return either ImageGeometry or AcquisitionGeometry 
+        depepnding on file content
+        '''
+        
+        return self._geometry
+
+
+'''
+# usage example
+reader = NEXUSDataReader()
+reader.set_up(nexus_file = '/home/evelina/test_nexus.nxs')
+acquisition_data = reader.load_data()
+print(acquisition_data)
+ag = reader.get_geometry()
+print(ag)
+
+reader = NEXUSDataReader()
+reader.set_up(nexus_file = '/home/evelina/test_nexus_im.nxs')
+image_data = reader.load_data()
+print(image_data)
+ig = reader.get_geometry()
+print(ig)
+'''
diff --git a/Wrappers/Python/ccpi/io/NEXUSDataWriter.py b/Wrappers/Python/ccpi/io/NEXUSDataWriter.py
new file mode 100644
index 0000000..f780f79
--- /dev/null
+++ b/Wrappers/Python/ccpi/io/NEXUSDataWriter.py
@@ -0,0 +1,151 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu May  2 10:11:20 2019
+
+@author: evelina
+"""
+
+
+import numpy
+import os
+from ccpi.framework import AcquisitionData, AcquisitionGeometry, ImageData, ImageGeometry
+import datetime
+
+
+h5pyAvailable = True
+try:
+    import h5py 
+except:
+    h5pyAvailable = False
+
+
+class NEXUSDataWriter(object):
+    
+    def __init__(self,
+                 **kwargs):
+        
+        self.data_container = kwargs.get('data_container', None)
+        self.file_name = kwargs.get('file_name', None)
+        
+        if ((self.data_container is not None) and (self.file_name is not None)):
+            self.set_up(data_container = self.data_container,
+                        file_name = self.file_name)
+        
+    def set_up(self,
+               data_container = None,
+               file_name = None):
+        
+        self.data_container = data_container
+        self.file_name = file_name
+        
+        if not ((isinstance(self.data_container, ImageData)) or 
+                (isinstance(self.data_container, AcquisitionData))):
+            raise Exception('Writer supports only following data types:\n' +
+                            ' - ImageData\n - AcquisitionData')
+        
+        # check that h5py library is installed
+        if (h5pyAvailable == False):
+            raise Exception('h5py is not available, cannot load NEXUS files.')
+    
+    def write_file(self):
+        
+        # if the folder does not exist, create the folder
+        if not os.path.isdir(os.path.dirname(self.file_name)):
+            os.mkdir(os.path.dirname(self.file_name))
+            
+        # create the file
+        with h5py.File(self.file_name, 'w') as f:
+            
+            # give the file some important attributes
+            f.attrs['file_name'] = self.file_name
+            f.attrs['file_time'] = str(datetime.datetime.utcnow())
+            f.attrs['creator'] = 'NEXUSDataWriter.py'
+            f.attrs['NeXus_version'] = '4.3.0'
+            f.attrs['HDF5_Version'] = h5py.version.hdf5_version
+            f.attrs['h5py_version'] = h5py.version.version
+            
+            # create the NXentry group
+            nxentry = f.create_group('entry1/tomo_entry')
+            nxentry.attrs['NX_class'] = 'NXentry'
+            
+            # create dataset to store data
+            ds_data = f.create_dataset('entry1/tomo_entry/data/data', 
+                                       (self.data_container.as_array().shape), 
+                                       dtype = 'float32', 
+                                       data = self.data_container.as_array())
+            
+            # set up dataset attributes
+            if (isinstance(self.data_container, ImageData)):
+                ds_data.attrs['data_type'] = 'ImageData'
+            else:
+                ds_data.attrs['data_type'] = 'AcquisitionData'
+            
+            for i in range(self.data_container.as_array().ndim):
+                ds_data.attrs['dim{}'.format(i)] = self.data_container.dimension_labels[i]
+            
+            if (isinstance(self.data_container, AcquisitionData)):      
+                ds_data.attrs['geom_type'] = self.data_container.geometry.geom_type
+                ds_data.attrs['dimension'] = self.data_container.geometry.dimension
+                ds_data.attrs['dist_source_center'] = self.data_container.geometry.dist_source_center
+                ds_data.attrs['dist_center_detector'] = self.data_container.geometry.dist_center_detector
+                ds_data.attrs['pixel_num_h'] = self.data_container.geometry.pixel_num_h
+                ds_data.attrs['pixel_size_h'] = self.data_container.geometry.pixel_size_h
+                ds_data.attrs['pixel_num_v'] = self.data_container.geometry.pixel_num_v
+                ds_data.attrs['pixel_size_v'] = self.data_container.geometry.pixel_size_v
+                ds_data.attrs['channels'] = self.data_container.geometry.channels
+                ds_data.attrs['n_angles'] = self.data_container.geometry.angles.shape[0]
+                
+                # create the dataset to store rotation angles
+                ds_angles = f.create_dataset('entry1/tomo_entry/data/rotation_angle', 
+                                             (self.data_container.geometry.angles.shape), 
+                                             dtype = 'float32', 
+                                             data = self.data_container.geometry.angles)
+                
+                #ds_angles.attrs['units'] = self.data_container.geometry.angle_unit
+            
+            else:   # ImageData
+                
+                ds_data.attrs['voxel_num_x'] = self.data_container.geometry.voxel_num_x
+                ds_data.attrs['voxel_num_y'] = self.data_container.geometry.voxel_num_y
+                ds_data.attrs['voxel_num_z'] = self.data_container.geometry.voxel_num_z
+                ds_data.attrs['voxel_size_x'] = self.data_container.geometry.voxel_size_x
+                ds_data.attrs['voxel_size_y'] = self.data_container.geometry.voxel_size_y
+                ds_data.attrs['voxel_size_z'] = self.data_container.geometry.voxel_size_z
+                ds_data.attrs['center_x'] = self.data_container.geometry.center_x
+                ds_data.attrs['center_y'] = self.data_container.geometry.center_y
+                ds_data.attrs['center_z'] = self.data_container.geometry.center_z
+                ds_data.attrs['channels'] = self.data_container.geometry.channels
+
+
+'''
+# usage example
+xtek_file = '/home/evelina/nikon_data/SophiaBeads_256_averaged.xtekct'
+reader = NikonDataReader()
+reader.set_up(xtek_file = xtek_file,
+              binning = [3, 1],
+              roi = [200, 500, 1500, 2000],
+              normalize = True)
+
+data = reader.load_projections()
+ag = reader.get_geometry()
+
+writer = NEXUSDataWriter()
+writer.set_up(file_name = '/home/evelina/test_nexus.nxs',
+              data_container = data)
+
+writer.write_file()
+
+ig = ImageGeometry(voxel_num_x = 100,
+                   voxel_num_y = 100)
+
+im = ImageData(array = numpy.zeros((100, 100), dtype = 'float'),
+               geometry = ig)
+
+im_writer = NEXUSDataWriter()
+
+writer.set_up(file_name = '/home/evelina/test_nexus_im.nxs',
+              data_container = im)
+
+writer.write_file()
+'''
diff --git a/Wrappers/Python/ccpi/io/NikonDataReader.py b/Wrappers/Python/ccpi/io/NikonDataReader.py
new file mode 100644
index 0000000..703b65b
--- /dev/null
+++ b/Wrappers/Python/ccpi/io/NikonDataReader.py
@@ -0,0 +1,281 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Apr  3 10:30:25 2019
+
+@author: evelina
+"""
+
+from ccpi.framework import AcquisitionData, AcquisitionGeometry
+import numpy
+import os
+
+
+pilAvailable = True
+try:    
+    from PIL import Image
+except:
+    pilAvailable = False
+    
+        
+class NikonDataReader(object):
+    
+    def __init__(self, 
+                 **kwargs):
+        '''
+        Constructor
+        
+        Input:
+            
+            xtek_file       full path to .xtexct file
+            
+            roi             region-of-interest to load. If roi = -1 (default), 
+                            full projections will be loaded. Otherwise roi is 
+                            given by [(row0, row1), (column0, column1)], where 
+                            row0, column0 are coordinates of top left corner and 
+                            row1, column1 are coordinates of bottom right corner.
+                            
+            binning         number of pixels to bin (combine) along 0 (column) 
+                            and 1 (row) dimension. If binning = [1, 1] (default),
+                            projections in original resolution are loaded. Note, 
+                            if binning[0] != binning[1], then loaded projections
+                            will have anisotropic pixels, which are currently not 
+                            supported by the Framework
+            
+            normalize       normalize loaded projections by detector 
+                            white level (I_0). Default value is False, 
+                            i.e. no normalization.
+                            
+            flip            default = False, flip projections in the left-right direction
+                    
+        '''
+        
+        self.xtek_file = kwargs.get('xtek_file', None)
+        self.roi = kwargs.get('roi', -1)
+        self.binning = kwargs.get('binning', [1, 1])
+        self.normalize = kwargs.get('normalize', False)
+        self.flip = kwargs.get('flip', False)
+        
+        if self.xtek_file is not None:
+            self.set_up(xtek_file = self.xtek_file,
+                        roi = self.roi,
+                        binning = self.binning,
+                        normalize = self.normalize,
+                        flip = self.flip)
+            
+    def set_up(self, 
+               xtek_file = None, 
+               roi = -1, 
+               binning = [1, 1],
+               normalize = False,
+               flip = False):
+        
+        self.xtek_file = xtek_file
+        self.roi = roi
+        self.binning = binning
+        self.normalize = normalize
+        self.flip = flip
+        
+        if self.xtek_file == None:
+            raise Exception('Path to xtek file is required.')
+        
+        # check if xtek file exists
+        if not(os.path.isfile(self.xtek_file)):
+            raise Exception('File\n {}\n does not exist.'.format(self.xtek_file))  
+                
+        # check that PIL library is installed
+        if (pilAvailable == False):
+            raise Exception("PIL (pillow) is not available, cannot load TIFF files.")
+                
+        # parse xtek file
+        with open(self.xtek_file, 'r') as f:
+            content = f.readlines()    
+                
+        content = [x.strip() for x in content]
+        
+        for line in content:
+            # filename of TIFF files
+            if line.startswith("Name"):
+                self._experiment_name = line.split('=')[1]
+            # number of projections
+            elif line.startswith("Projections"):
+                num_projections = int(line.split('=')[1])
+            # white level - used for normalization
+            elif line.startswith("WhiteLevel"):
+                self._white_level = float(line.split('=')[1])
+            # number of pixels along Y axis
+            elif line.startswith("DetectorPixelsY"):
+                pixel_num_v_0 = int(line.split('=')[1])
+            # number of pixels along X axis
+            elif line.startswith("DetectorPixelsX"):
+                pixel_num_h_0 = int(line.split('=')[1])
+            # pixel size along X axis
+            elif line.startswith("DetectorPixelSizeX"):
+                pixel_size_h_0 = float(line.split('=')[1])
+            # pixel size along Y axis
+            elif line.startswith("DetectorPixelSizeY"):
+                pixel_size_v_0 = float(line.split('=')[1])
+            # source to center of rotation distance
+            elif line.startswith("SrcToObject"):
+                source_x = float(line.split('=')[1])
+            # source to detector distance
+            elif line.startswith("SrcToDetector"):
+                detector_x = float(line.split('=')[1])
+            # initial angular position of a rotation stage
+            elif line.startswith("InitialAngle"):
+                initial_angle = float(line.split('=')[1])
+            # angular increment (in degrees)
+            elif line.startswith("AngularStep"):
+                angular_step = float(line.split('=')[1])
+                
+        if self.roi == -1:
+            self._roi_par = [(0, pixel_num_v_0), \
+                              (0, pixel_num_h_0)]
+        else:
+            self._roi_par = self.roi.copy()
+            if self._roi_par[0] == -1:
+                self._roi_par[0] = (0, pixel_num_v_0)
+            if self._roi_par[1] == -1:
+                self._roi_par[1] = (0, pixel_num_h_0)
+                
+        # calculate number of pixels and pixel size
+        if (self.binning == [1, 1]):
+            pixel_num_v = self._roi_par[0][1] - self._roi_par[0][0]
+            pixel_num_h = self._roi_par[1][1] - self._roi_par[1][0]
+            pixel_size_v = pixel_size_v_0
+            pixel_size_h = pixel_size_h_0
+        else:
+            pixel_num_v = (self._roi_par[0][1] - self._roi_par[0][0]) // self.binning[0]
+            pixel_num_h = (self._roi_par[1][1] - self._roi_par[1][0]) // self.binning[1]
+            pixel_size_v = pixel_size_v_0 * self.binning[0]
+            pixel_size_h = pixel_size_h_0 * self.binning[1]
+        
+        '''
+        Parse the angles file .ang or _ctdata.txt file and returns the angles
+        as an numpy array. 
+        '''
+        input_path = os.path.dirname(self.xtek_file)
+        angles_ctdata_file = os.path.join(input_path, '_ctdata.txt')
+        angles_named_file = os.path.join(input_path, self._experiment_name+'.ang')
+        angles = numpy.zeros(num_projections, dtype = 'float')
+        
+        # look for _ctdata.txt
+        if os.path.exists(angles_ctdata_file):
+            # read txt file with angles
+            with open(angles_ctdata_file) as f:
+                content = f.readlines()
+            # skip firt three lines
+            # read the middle value of 3 values in each line as angles in degrees
+            index = 0
+            for line in content[3:]:
+                angles[index] = float(line.split(' ')[1])
+                index += 1
+            angles = angles + initial_angle
+        
+        # look for ang file
+        elif os.path.exists(angles_named_file):
+            # read the angles file which is text with first line as header
+            with open(angles_named_file) as f:
+                content = f.readlines()
+            # skip first line
+            index = 0
+            for line in content[1:]:
+                angles[index] = float(line.split(':')[1])
+                index += 1
+            angles = numpy.flipud(angles + initial_angle) # angles are in the reverse order
+            
+        else:   # calculate angles based on xtek file
+            angles = initial_angle + angular_step * range(num_projections)
+        
+        # fill in metadata
+        self._ag = AcquisitionGeometry(geom_type = 'cone', 
+                                       dimension = '3D', 
+                                       angles = angles, 
+                                       pixel_num_h = pixel_num_h, 
+                                       pixel_size_h = pixel_size_h, 
+                                       pixel_num_v = pixel_num_v, 
+                                       pixel_size_v = pixel_size_v, 
+                                       dist_source_center = source_x, 
+                                       dist_center_detector = detector_x - source_x, 
+                                       channels = 1,
+                                       angle_unit = 'degree')
+
+    def get_geometry(self):
+        
+        '''
+        Return AcquisitionGeometry object
+        '''
+        
+        return self._ag
+        
+    def load_projections(self):
+        
+        '''
+        Load projections and return AcquisitionData container
+        '''
+            
+        # get path to projections
+        path_projection = os.path.dirname(self.xtek_file)
+        
+        # get number of projections
+        num_projections = numpy.shape(self._ag.angles)[0]
+        
+        # allocate array to store projections    
+        data = numpy.zeros((num_projections, self._ag.pixel_num_v, self._ag.pixel_num_h), dtype = float)
+        
+        for i in range(num_projections):
+            
+            filename = (path_projection + '/' + self._experiment_name + '_{:04d}.tif').format(i + 1)
+            
+            try:
+                tmp = numpy.asarray(Image.open(filename), dtype = float)
+            except:
+                print('Error reading\n {}\n file.'.format(filename))
+                raise
+                
+            if (self.binning == [1, 1]):
+                data[i, :, :] = tmp[self._roi_par[0][0]:self._roi_par[0][1], self._roi_par[1][0]:self._roi_par[1][1]]
+            else:
+                shape = (self._ag.pixel_num_v, self.binning[0], 
+                         self._ag.pixel_num_h, self.binning[1])
+                data[i, :, :] = tmp[self._roi_par[0][0]:(self._roi_par[0][0] + (((self._roi_par[0][1] - self._roi_par[0][0]) // self.binning[0]) * self.binning[0])), \
+                                    self._roi_par[1][0]:(self._roi_par[1][0] + (((self._roi_par[1][1] - self._roi_par[1][0]) // self.binning[1]) * self.binning[1]))].reshape(shape).mean(-1).mean(1)
+        
+        if (self.normalize):
+            data /= self._white_level
+            data[data > 1] = 1
+        
+        if self.flip:
+            return AcquisitionData(array = data[:, :, ::-1], 
+                                   deep_copy = False,
+                                   geometry = self._ag,
+                                   dimension_labels = ['angle', \
+                                                       'vertical', \
+                                                       'horizontal'])
+        else:
+            return AcquisitionData(array = data, 
+                                   deep_copy = False,
+                                   geometry = self._ag,
+                                   dimension_labels = ['angle', \
+                                                       'vertical', \
+                                                       'horizontal'])
+
+
+'''
+# usage example
+xtek_file = '/home/evelina/nikon_data/SophiaBeads_256_averaged.xtekct'
+reader = NikonDataReader()
+reader.set_up(xtek_file = xtek_file,
+              binning = [1, 1],
+              roi = -1,
+              normalize = True,
+              flip = True)
+
+data = reader.load_projections()
+print(data)
+ag = reader.get_geometry()
+print(ag)
+
+plt.imshow(data.as_array()[1, :, :])
+plt.show()
+'''
diff --git a/Wrappers/Python/ccpi/io/__init__.py b/Wrappers/Python/ccpi/io/__init__.py
index 9233d7a..455faba 100644
--- a/Wrappers/Python/ccpi/io/__init__.py
+++ b/Wrappers/Python/ccpi/io/__init__.py
@@ -15,4 +15,8 @@
 #   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.
\ No newline at end of file
+#   limitations under the License.

+

+from .NEXUSDataReader import NEXUSDataReader

+from .NEXUSDataWriter import NEXUSDataWriter

+from .NikonDataReader import NikonDataReader

diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py
index 647ae98..c8fd0d4 100755
--- a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py
@@ -70,7 +70,8 @@ class FISTA(Algorithm):
         
         self.t = 0.5*(1 + numpy.sqrt(1 + 4*(self.t_old**2)))
         
-        self.x.subtract(self.x_old, out=self.y)
+#        self.x.subtract(self.x_old, out=self.y)
+        self.y = self.x - self.x_old
         self.y.__imul__ ((self.t_old-1)/self.t)
         self.y.__iadd__( self.x )
         
diff --git a/Wrappers/Python/ccpi/processors/Resizer.py b/Wrappers/Python/ccpi/processors/Resizer.py
new file mode 100755
index 0000000..7509a90
--- /dev/null
+++ b/Wrappers/Python/ccpi/processors/Resizer.py
@@ -0,0 +1,262 @@
+# -*- 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 ccpi.framework import DataProcessor, AcquisitionData, ImageData

+import warnings

+

+

+class Resizer(DataProcessor):

+

+    def __init__(self,

+                 roi = -1,

+                 binning = 1):

+        

+        '''

+        Constructor

+        

+        Input:

+            

+            roi             region-of-interest to crop. If roi = -1 (default), then no crop. 

+                            Otherwise roi is given by a list with ndim elements, 

+                            where each element is either -1 if no crop along this 

+                            dimension or a tuple with beginning and end coodinates to crop to.

+                            Example:

+                                to crop 4D array along 2nd dimension:

+                                roi = [-1, -1, (100, 900), -1]

+                            

+            binning         number of pixels to bin (combine) along each dimension.

+                            If binning = 1, then projections in original resolution are loaded. 

+                            Otherwise, binning is given by a list with ndim integers. 

+                            Example:

+                                to rebin 3D array along 1st direction:

+                                binning = [1, 5, 1]

+        '''

+

+        kwargs = {'roi': roi,

+                  'binning': binning}

+

+        super(Resizer, self).__init__(**kwargs)

+    

+    def check_input(self, data):

+        if not ((isinstance(data, ImageData)) or 

+                (isinstance(data, AcquisitionData))):

+            raise Exception('Processor supports only following data types:\n' +

+                            ' - ImageData\n - AcquisitionData')

+        elif (data.geometry == None):

+            raise Exception('Geometry is not defined.')

+        else:

+            return True

+    

+    def process(self):

+

+        data = self.get_input()

+        ndim = len(data.dimension_labels)

+        

+        geometry_0 = data.geometry

+        geometry = geometry_0.clone()

+        

+        if (self.roi == -1):

+            roi_par = [-1] * ndim

+        else:

+            roi_par = self.roi.copy()

+            if (len(roi_par) != ndim):

+                raise Exception('Number of dimensions and number of elements in roi parameter do not match')

+

+        if (self.binning == 1):

+            binning = [1] * ndim

+        else:

+            binning = self.binning.copy()

+            if (len(binning) != ndim):

+                raise Exception('Number of dimensions and number of elements in binning parameter do not match')

+                

+        if (isinstance(data, ImageData)):

+            if ((all(x == -1 for x in roi_par)) and (all(x == 1 for x in binning))):

+                for key in data.dimension_labels:

+                    if data.dimension_labels[key] == 'channel':

+                        geometry.channels = geometry_0.channels

+                        roi_par[key] = (0, geometry.channels)

+                    elif data.dimension_labels[key] == 'horizontal_y':

+                        geometry.voxel_size_y = geometry_0.voxel_size_y

+                        geometry.voxel_num_y = geometry_0.voxel_num_y

+                        roi_par[key] = (0, geometry.voxel_num_y)

+                    elif data.dimension_labels[key] == 'vertical':

+                        geometry.voxel_size_z = geometry_0.voxel_size_z

+                        geometry.voxel_num_z = geometry_0.voxel_num_z

+                        roi_par[key] = (0, geometry.voxel_num_z)

+                    elif data.dimension_labels[key] == 'horizontal_x':

+                        geometry.voxel_size_x = geometry_0.voxel_size_x

+                        geometry.voxel_num_x = geometry_0.voxel_num_x

+                        roi_par[key] = (0, geometry.voxel_num_x)

+            else:

+                for key in data.dimension_labels:

+                    if data.dimension_labels[key] == 'channel':

+                        if (roi_par[key] != -1):

+                            geometry.channels = (roi_par[key][1] - roi_par[key][0]) // binning[key]

+                            roi_par[key] = (roi_par[key][0], roi_par[key][0] + ((roi_par[key][1] - roi_par[key][0]) // binning[key]) * binning[key])

+                        else:

+                            geometry.channels = geometry_0.channels // binning[key]

+                            roi_par[key] = (0, geometry.channels * binning[key])

+                    elif data.dimension_labels[key] == 'horizontal_y':

+                        if (roi_par[key] != -1):

+                            geometry.voxel_num_y = (roi_par[key][1] - roi_par[key][0]) // binning[key]

+                            geometry.voxel_size_y = geometry_0.voxel_size_y * binning[key]

+                            roi_par[key] = (roi_par[key][0], roi_par[key][0] + ((roi_par[key][1] - roi_par[key][0]) // binning[key]) * binning[key])

+                        else:

+                            geometry.voxel_num_y = geometry_0.voxel_num_y // binning[key]

+                            geometry.voxel_size_y = geometry_0.voxel_size_y * binning[key]

+                            roi_par[key] = (0, geometry.voxel_num_y * binning[key])

+                    elif data.dimension_labels[key] == 'vertical':

+                        if (roi_par[key] != -1):

+                            geometry.voxel_num_z = (roi_par[key][1] - roi_par[key][0]) // binning[key]

+                            geometry.voxel_size_z = geometry_0.voxel_size_z * binning[key]

+                            roi_par[key] = (roi_par[key][0], roi_par[key][0] + ((roi_par[key][1] - roi_par[key][0]) // binning[key]) * binning[key])

+                        else:

+                            geometry.voxel_num_z = geometry_0.voxel_num_z // binning[key]

+                            geometry.voxel_size_z = geometry_0.voxel_size_z * binning[key]

+                            roi_par[key] = (0, geometry.voxel_num_z * binning[key])

+                    elif data.dimension_labels[key] == 'horizontal_x':

+                        if (roi_par[key] != -1):

+                            geometry.voxel_num_x = (roi_par[key][1] - roi_par[key][0]) // binning[key]

+                            geometry.voxel_size_x = geometry_0.voxel_size_x * binning[key]

+                            roi_par[key] = (roi_par[key][0], roi_par[key][0]+ ((roi_par[key][1] - roi_par[key][0]) // binning[key]) * binning[key])

+                        else:

+                            geometry.voxel_num_x = geometry_0.voxel_num_x // binning[key]

+                            geometry.voxel_size_x = geometry_0.voxel_size_x * binning[key]

+                            roi_par[key] = (0, geometry.voxel_num_x * binning[key])

+            

+        else: # AcquisitionData

+            if ((all(x == -1 for x in roi_par)) and (all(x == 1 for x in binning))):

+                for key in data.dimension_labels:

+                    if data.dimension_labels[key] == 'channel':

+                        geometry.channels = geometry_0.channels

+                        roi_par[key] = (0, geometry.channels)

+                    elif data.dimension_labels[key] == 'angle':

+                        geometry.angles = geometry_0.angles

+                        roi_par[key] = (0, len(geometry.angles))

+                    elif data.dimension_labels[key] == 'vertical':

+                        geometry.pixel_size_v = geometry_0.pixel_size_v

+                        geometry.pixel_num_v = geometry_0.pixel_num_v

+                        roi_par[key] = (0, geometry.pixel_num_v)

+                    elif data.dimension_labels[key] == 'horizontal':

+                        geometry.pixel_size_h = geometry_0.pixel_size_h

+                        geometry.pixel_num_h = geometry_0.pixel_num_h

+                        roi_par[key] = (0, geometry.pixel_num_h)

+            else:

+                for key in data.dimension_labels:

+                    if data.dimension_labels[key] == 'channel':

+                        if (roi_par[key] != -1):

+                            geometry.channels = (roi_par[key][1] - roi_par[key][0]) // binning[key]

+                            roi_par[key] = (roi_par[key][0], roi_par[key][0] + ((roi_par[key][1] - roi_par[key][0]) // binning[key]) * binning[key])

+                        else:

+                            geometry.channels = geometry_0.channels // binning[key]

+                            roi_par[key] = (0, geometry.channels * binning[key])

+                    elif data.dimension_labels[key] == 'angle':

+                        if (roi_par[key] != -1):

+                            geometry.angles = geometry_0.angles[roi_par[key][0]:roi_par[key][1]]

+                        else:

+                            geometry.angles = geometry_0.angles

+                            roi_par[key] = (0, len(geometry.angles))

+                        if (binning[key] != 1):

+                            binning[key] = 1

+                            warnings.warn('Rebinning in angular dimensions is not supported: \n binning[{}] is set to 1.'.format(key))

+                    elif data.dimension_labels[key] == 'vertical':

+                        if (roi_par[key] != -1):

+                            geometry.pixel_num_v = (roi_par[key][1] - roi_par[key][0]) // binning[key]

+                            geometry.pixel_size_v = geometry_0.pixel_size_v * binning[key]

+                            roi_par[key] = (roi_par[key][0], roi_par[key][0] + ((roi_par[key][1] - roi_par[key][0]) // binning[key]) * binning[key])

+                        else:

+                            geometry.pixel_num_v = geometry_0.pixel_num_v // binning[key]

+                            geometry.pixel_size_v = geometry_0.pixel_size_v * binning[key]

+                            roi_par[key] = (0, geometry.pixel_num_v * binning[key])

+                    elif data.dimension_labels[key] == 'horizontal':

+                        if (roi_par[key] != -1):

+                            geometry.pixel_num_h = (roi_par[key][1] - roi_par[key][0]) // binning[key]

+                            geometry.pixel_size_h = geometry_0.pixel_size_h * binning[key]

+                            roi_par[key] = (roi_par[key][0], roi_par[key][0] + ((roi_par[key][1] - roi_par[key][0]) // binning[key]) * binning[key])

+                        else:

+                            geometry.pixel_num_h = geometry_0.pixel_num_h // binning[key]

+                            geometry.pixel_size_h = geometry_0.pixel_size_h * binning[key]

+                            roi_par[key] = (0, geometry.pixel_num_h * binning[key])

+                            

+        if ndim == 2:

+            n_pix_0 = (roi_par[0][1] - roi_par[0][0]) // binning[0]

+            n_pix_1 = (roi_par[1][1] - roi_par[1][0]) // binning[1]

+            shape = (n_pix_0, binning[0], 

+                     n_pix_1, binning[1])

+            data_resized = data.as_array()[roi_par[0][0]:(roi_par[0][0] + n_pix_0 * binning[0]), 

+                                           roi_par[1][0]:(roi_par[1][0] + n_pix_1 * binning[1])].reshape(shape).mean(-1).mean(1)

+        if ndim == 3:

+            n_pix_0 = (roi_par[0][1] - roi_par[0][0]) // binning[0]

+            n_pix_1 = (roi_par[1][1] - roi_par[1][0]) // binning[1]

+            n_pix_2 = (roi_par[2][1] - roi_par[2][0]) // binning[2]

+            shape = (n_pix_0, binning[0], 

+                     n_pix_1, binning[1],

+                     n_pix_2, binning[2])

+            data_resized = data.as_array()[roi_par[0][0]:(roi_par[0][0] + n_pix_0 * binning[0]), 

+                                           roi_par[1][0]:(roi_par[1][0] + n_pix_1 * binning[1]), 

+                                           roi_par[2][0]:(roi_par[2][0] + n_pix_2 * binning[2])].reshape(shape).mean(-1).mean(1).mean(2)

+        if ndim == 4:

+            n_pix_0 = (roi_par[0][1] - roi_par[0][0]) // binning[0]

+            n_pix_1 = (roi_par[1][1] - roi_par[1][0]) // binning[1]

+            n_pix_2 = (roi_par[2][1] - roi_par[2][0]) // binning[2]

+            n_pix_3 = (roi_par[3][1] - roi_par[3][0]) // binning[3]

+            shape = (n_pix_0, binning[0], 

+                     n_pix_1, binning[1],

+                     n_pix_2, binning[2],

+                     n_pix_3, binning[3])

+            data_resized = data.as_array()[roi_par[0][0]:(roi_par[0][0] + n_pix_0 * binning[0]), 

+                                           roi_par[1][0]:(roi_par[1][0] + n_pix_1 * binning[1]), 

+                                           roi_par[2][0]:(roi_par[2][0] + n_pix_2 * binning[2]), 

+                                           roi_par[3][0]:(roi_par[3][0] + n_pix_3 * binning[3])].reshape(shape).mean(-1).mean(1).mean(2).mean(3)

+

+        out = type(data)(array = data_resized, 

+                         deep_copy = False,

+                         dimension_labels = data.dimension_labels,

+                         geometry = geometry)

+        

+        return out

+

+

+'''

+#usage exaample

+ig = ImageGeometry(voxel_num_x = 200, 

+                   voxel_num_y = 200, 

+                   voxel_num_z = 200, 

+                   voxel_size_x = 1, 

+                   voxel_size_y = 1, 

+                   voxel_size_z = 1, 

+                   center_x = 0, 

+                   center_y = 0, 

+                   center_z = 0, 

+                   channels = 200)

+

+im = ImageData(array = numpy.zeros((200, 200, 200, 200)),

+               geometry = ig,

+               deep_copy = False, 

+               dimension_labels = ['channel',\

+                                   'vertical',\

+                                   'horizontal_y',\

+                                   'horizontal_x'])

+            

+            

+resizer = Resizer(binning = [1, 1, 7, 1], roi = -1)

+resizer.input = im

+data_resized = resizer.process()

+print(data_resized)

+'''

diff --git a/Wrappers/Python/ccpi/processors/__init__.py b/Wrappers/Python/ccpi/processors/__init__.py
index f8d050e..cba5897 100755
--- a/Wrappers/Python/ccpi/processors/__init__.py
+++ b/Wrappers/Python/ccpi/processors/__init__.py
@@ -7,3 +7,4 @@ Created on Tue Apr 30 13:51:09 2019
 

 from .CenterOfRotationFinder import CenterOfRotationFinder

 from .Normalizer import Normalizer

+from .Resizer import Resizer

diff --git a/Wrappers/Python/data/shapes.png b/Wrappers/Python/data/shapes.png
new file mode 100644
index 0000000..dd4f680
--- /dev/null
+++ b/Wrappers/Python/data/shapes.png
diff --git a/Wrappers/Python/demos/CGLS_examples/CGLS_Tikhonov.py b/Wrappers/Python/demos/CGLS_examples/CGLS_Tikhonov.py
index d1cbe20..653e191 100644
--- a/Wrappers/Python/demos/CGLS_examples/CGLS_Tikhonov.py
+++ b/Wrappers/Python/demos/CGLS_examples/CGLS_Tikhonov.py
@@ -82,18 +82,18 @@ plt.colorbar()
 plt.show()
 
 # Setup and run the CGLS algorithm  
-alpha = 50
-Grad = Gradient(ig)
-
-# Form Tikhonov as a Block CGLS structure
-op_CGLS = BlockOperator( Aop, alpha * Grad, shape=(2,1))
-block_data = BlockDataContainer(noisy_data, Grad.range_geometry().allocate())
-
-x_init = ig.allocate()      
-cgls = CGLS(x_init=x_init, operator=op_CGLS, data=block_data)
-cgls.max_iteration = 1000
-cgls.update_objective_interval = 200
-cgls.run(1000,verbose=False)
+#alpha = 50
+#Grad = Gradient(ig)
+#
+## Form Tikhonov as a Block CGLS structure
+#op_CGLS = BlockOperator( Aop, alpha * Grad, shape=(2,1))
+#block_data = BlockDataContainer(noisy_data, Grad.range_geometry().allocate())
+#
+#x_init = ig.allocate()      
+#cgls = CGLS(x_init=x_init, operator=op_CGLS, data=block_data)
+#cgls.max_iteration = 1000
+#cgls.update_objective_interval = 200
+#cgls.run(1000,verbose=False)
 
 #%%
 # Show results
diff --git a/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py b/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py
index 0875c20..672d4bc 100644
--- a/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py
+++ b/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py
@@ -32,7 +32,7 @@ Problem:     min_x || A x - g ||_{2}^{2}
 """
 
 
-from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry
+from ccpi.framework import ImageData, TestData, AcquisitionGeometry
 
 import numpy as np 
 import numpy                          
@@ -42,17 +42,17 @@ from ccpi.optimisation.algorithms import PDHG, CGLS, FISTA
 
 from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, FunctionOperatorComposition
 from ccpi.astra.ops import AstraProjectorSimple
-import astra          
+import astra   
+import os, sys
 
-# Create Ground truth phantom and Sinogram
 
-N = 128
-x = np.zeros((N,N))
-x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
-x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1
-
-data = ImageData(x)
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
+# Load Data 
+loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))  
+                     
+N = 50
+M = 50
+data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,M), scale=(0,1))
+ig = data.geometry
 
 detectors = N
 angles = np.linspace(0, np.pi, N, dtype=np.float32)
@@ -69,13 +69,14 @@ sin = Aop.direct(data)
 
 noisy_data = sin 
 
+###############################################################################
 # Setup and run Astra CGLS algorithm
 vol_geom = astra.create_vol_geom(N, N)
 proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles)
-proj_id = astra.create_projector('strip', proj_geom, vol_geom)
+proj_id = astra.create_projector('linear', proj_geom, vol_geom)
 
 # Create a sinogram from a phantom
-sinogram_id, sinogram = astra.create_sino(x, proj_id)
+sinogram_id, sinogram = astra.create_sino(data.as_array(), proj_id)
 
 # Create a data object for the reconstruction
 rec_id = astra.data2d.create('-vol', vol_geom)
@@ -92,6 +93,7 @@ astra.algorithm.run(alg_id, 1000)
 
 recon_cgls_astra = astra.data2d.get(rec_id)
 
+###############################################################################
 # Setup and run the CGLS algorithm  
 x_init = ig.allocate()             
 cgls = CGLS(x_init=x_init, operator=Aop, data=noisy_data)
@@ -99,12 +101,15 @@ cgls.max_iteration = 1000
 cgls.update_objective_interval = 200
 cgls.run(1000, verbose=False)
 
+###############################################################################
 # Setup and run the PDHG algorithm 
 operator = Aop
 f = L2NormSquared(b = noisy_data)
 g = ZeroFunction()
+
 ## Compute operator Norm
 normK = operator.norm()
+
 ## Primal & dual stepsizes
 sigma = 1
 tau = 1/(sigma*normK**2)
@@ -112,17 +117,17 @@ tau = 1/(sigma*normK**2)
 pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True)
 pdhg.max_iteration = 1000
 pdhg.update_objective_interval = 200
-pdhg.run(1000, verbose=False)
+pdhg.run(1000, verbose=True)
 
+###############################################################################
 # Setup and run the FISTA algorithm 
 fidelity = FunctionOperatorComposition(L2NormSquared(b=noisy_data), Aop)
 regularizer = ZeroFunction()
 
-opt = {'memopt':True}
-fista = FISTA(x_init=x_init , f=fidelity, g=regularizer, opt=opt)
+fista = FISTA(x_init=x_init , f=fidelity, g=regularizer)
 fista.max_iteration = 1000
 fista.update_objective_interval = 200
-fista.run(1000, verbose=False)
+fista.run(1000, verbose=True)
 
 #%% Show results
 
@@ -131,18 +136,22 @@ plt.suptitle('Reconstructions ', fontsize=16)
 
 plt.subplot(2,2,1)
 plt.imshow(cgls.get_output().as_array())
+plt.colorbar()
 plt.title('CGLS reconstruction')
 
 plt.subplot(2,2,2)
 plt.imshow(fista.get_output().as_array())
+plt.colorbar()
 plt.title('FISTA reconstruction')
 
 plt.subplot(2,2,3)
 plt.imshow(pdhg.get_output().as_array())
+plt.colorbar()
 plt.title('PDHG reconstruction')
 
 plt.subplot(2,2,4)
 plt.imshow(recon_cgls_astra)
+plt.colorbar()
 plt.title('CGLS astra')
 
 diff1 = pdhg.get_output() - cgls.get_output()
@@ -167,28 +176,14 @@ plt.title('Diff CLGS astra vs CGLS')
 plt.colorbar()
 
 
+#%%
 
+print('Primal Objective (FISTA) {} '.format(fista.objective[-1]))
+print('Primal Objective (CGLS) {} '.format(cgls.objective[-1]))
+print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
 
 
+true_obj = (Aop.direct(cglsd.get_output())-noisy_data).squared_norm()
+print('True objective {}'.format(true_obj))
 
 
-
-
-
-
-
-
-
-
-
-
-
-
-#
-#
-#
-#
-#
-#
-#
-#
diff --git a/Wrappers/Python/demos/FISTA_examples/FISTA_CGLS.py b/Wrappers/Python/demos/FISTA_examples/FISTA_CGLS.py
index 68fb649..1a96a2d 100644
--- a/Wrappers/Python/demos/FISTA_examples/FISTA_CGLS.py
+++ b/Wrappers/Python/demos/FISTA_examples/FISTA_CGLS.py
@@ -21,14 +21,15 @@
 
 from ccpi.framework import AcquisitionGeometry
 from ccpi.optimisation.algorithms import FISTA
-from ccpi.optimisation.functions import IndicatorBox, ZeroFunction, L2NormSquared, FunctionOperatorComposition
-from ccpi.astra.ops import AstraProjectorSimple
+from ccpi.optimisation.functions import IndicatorBox, ZeroFunction, \
+                         L2NormSquared, FunctionOperatorComposition
+from ccpi.astra.operators import AstraProjectorSimple
 
 import numpy as np
 import matplotlib.pyplot as plt
 from ccpi.framework import TestData
 import os, sys
-
+from ccpi.optimisation.funcs import Norm2sq
 
 # Load Data 
 loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
@@ -117,16 +118,15 @@ plt.show()
 # demonstrated in the rest of this file. In general all methods need an initial 
 # guess and some algorithm options to be set:
 
-f = FunctionOperatorComposition(0.5 * L2NormSquared(b=sin), Aop)
+f = FunctionOperatorComposition(L2NormSquared(b=sin), Aop)
+#f = Norm2sq(Aop, sin, c=0.5, memopt=True)
 g = ZeroFunction()
 
 x_init = ig.allocate()
-opt = {'memopt':True}
-
-fista = FISTA(x_init=x_init, f=f, g=g, opt=opt)
-fista.max_iteration = 100
-fista.update_objective_interval = 1
-fista.run(100, verbose=True)
+fista = FISTA(x_init=x_init, f=f, g=g)
+fista.max_iteration = 1000
+fista.update_objective_interval = 100
+fista.run(1000, verbose=True)
 
 plt.figure()
 plt.imshow(fista.get_output().as_array())
@@ -134,18 +134,12 @@ plt.title('FISTA unconstrained')
 plt.colorbar()
 plt.show()
 
-plt.figure()
-plt.semilogy(fista.objective)
-plt.title('FISTA objective')
-plt.show()
-
-#%%
 
 # Run FISTA for least squares with lower/upper bound 
-fista0 = FISTA(x_init=x_init, f=f, g=IndicatorBox(lower=0,upper=1), opt=opt)
-fista0.max_iteration = 100
-fista0.update_objective_interval = 1
-fista0.run(100, verbose=True)
+fista0 = FISTA(x_init=x_init, f=f, g=IndicatorBox(lower=0,upper=1))
+fista0.max_iteration = 1000
+fista0.update_objective_interval = 100
+fista0.run(1000, verbose=True)
 
 plt.figure()
 plt.imshow(fista0.get_output().as_array())
@@ -153,10 +147,10 @@ plt.title('FISTA constrained in [0,1]')
 plt.colorbar()
 plt.show()
 
-plt.figure()
-plt.semilogy(fista0.objective)
-plt.title('FISTA constrained in [0,1]')
-plt.show()
+#plt.figure()
+#plt.semilogy(fista0.objective)
+#plt.title('FISTA constrained in [0,1]')
+#plt.show()
 
 #%% Check with CVX solution
 
@@ -178,10 +172,10 @@ if cvx_not_installable:
     vol_geom = astra.create_vol_geom(N, N)
     proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles)
 
-    proj_id = astra.create_projector('strip', proj_geom, vol_geom)
+    proj_id = astra.create_projector('linear', proj_geom, vol_geom)
 
     matrix_id = astra.projector.matrix(proj_id)
-
+    
     ProjMat = astra.matrix.get(matrix_id)
     
     tmp = sin.as_array().ravel()
@@ -216,4 +210,6 @@ if cvx_not_installable:
     plt.legend()
     plt.title('Middle Line Profiles')
     plt.show()
-            
\ No newline at end of file
+            
+    print('Primal Objective (CVX) {} '.format(obj.value))
+    print('Primal Objective (FISTA) {} '.format(fista0.loss[1]))    
\ No newline at end of file
diff --git a/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py b/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py
index 41219f4..6007990 100644
--- a/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py
+++ b/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py
@@ -21,7 +21,7 @@
 
 """ 
 
-Tikhonov for Poisson denoising using FISTA algorithm:
+"Tikhonov regularization" for Poisson denoising using FISTA algorithm:
 
 Problem:     min_x, x>0  \alpha * ||\nabla x||_{2}^{2} + \int x - g * log(x) 
 
@@ -52,14 +52,14 @@ from skimage.util import random_noise
 loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
 
 # Load Data                      
-N = 150
-M = 150
+N = 100
+M = 100
 data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,M), scale=(0,1))
 
 ig = data.geometry
 ag = ig
 
-# Create Noisy data. Add Gaussian noise
+# Create Noisy data with Poisson noise
 n1 = random_noise(data.as_array(), mode = 'poisson', seed = 10)
 noisy_data = ImageData(n1)
 
@@ -75,7 +75,6 @@ plt.title('Noisy Data')
 plt.colorbar()
 plt.show()
 
-#%%
 # Regularisation Parameter
 alpha = 10
 
@@ -114,8 +113,7 @@ fid.proximal = KL_Prox_PosCone
 reg = FunctionOperatorComposition(alpha * L2NormSquared(), operator)
 
 x_init = ig.allocate()
-opt = {'memopt':True}
-fista = FISTA(x_init=x_init , f=reg, g=fid, opt=opt)
+fista = FISTA(x_init=x_init , f=reg, g=fid)
 fista.max_iteration = 2000
 fista.update_objective_interval = 500
 fista.run(2000, verbose=True)
@@ -196,7 +194,6 @@ if cvx_not_installable:
     plt.title('Middle Line Profiles')
     plt.show()
             
-    #TODO what is the output of fista.objective, fista.loss
     print('Primal Objective (CVX) {} '.format(obj.value))
     print('Primal Objective (FISTA) {} '.format(fista.loss[1]))
 
diff --git a/Wrappers/Python/demos/PDHG_examples/.DS_Store b/Wrappers/Python/demos/PDHG_examples/.DS_Store
deleted file mode 100644
index 141508d..0000000
--- a/Wrappers/Python/demos/PDHG_examples/.DS_Store
+++ /dev/null
diff --git a/Wrappers/Python/demos/PDHG_examples/ColorbayDemo.py b/Wrappers/Python/demos/PDHG_examples/ColorbayDemo.py
index a735323..e69060f 100644
--- a/Wrappers/Python/demos/PDHG_examples/ColorbayDemo.py
+++ b/Wrappers/Python/demos/PDHG_examples/ColorbayDemo.py
@@ -57,9 +57,8 @@ elif phantom == 'powder':
         arrays[k] = numpy.array(v)
     XX = arrays['S']    
     X = numpy.transpose(XX,(0,2,1,3))
-    X = X[0:20]
+    X = X[100:120]
     
-        
     
 #%% Setup Geometry of Colorbay
 
@@ -125,11 +124,16 @@ plt.show()
 
 #%% CGLS
 
+def callback(iteration, objective, x):
+    plt.imshow(x.as_array()[5])
+    plt.colorbar()
+    plt.show()
+
 x_init = ig2d.allocate()      
 cgls1 = CGLS(x_init=x_init, operator=Aall, data=data2d)
 cgls1.max_iteration = 100
 cgls1.update_objective_interval = 1
-cgls1.run(5,verbose=True)
+cgls1.run(5,verbose=True, callback = callback)
 
 plt.imshow(cgls1.get_output().subset(channel=5).array)
 plt.title('CGLS')
@@ -148,7 +152,7 @@ cgls2 = CGLS(x_init=x_init, operator=op_CGLS, data=block_data)
 cgls2.max_iteration = 100
 cgls2.update_objective_interval = 1
 
-cgls2.run(10,verbose=True)
+cgls2.run(10,verbose=True, callback=callback)
 
 plt.imshow(cgls2.get_output().subset(channel=5).array)
 plt.title('Tikhonov')
@@ -174,8 +178,9 @@ f = BlockFunction(f1, f2)
 g = ZeroFunction()
 
 # Compute operator Norm
-normK = 8.70320267279591 # Run one time no need to compute again takes time
-    
+#normK = 8.70320267279591 # For powder Run one time no need to compute again takes time
+normK = 14.60320657253632 # for carbon
+
 # Primal & dual stepsizes
 sigma = 1
 tau = 1/(sigma*normK**2)
@@ -184,11 +189,11 @@ tau = 1/(sigma*normK**2)
 pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
 pdhg.max_iteration = 2000
 pdhg.update_objective_interval = 100
-pdhg.run(1000, verbose =True)
+pdhg.run(1000, verbose =True, callback=callback)
 
 
 #%% Show sinograms
-channel_ind = [10,15,15]
+channel_ind = [10,15,19]
 
 plt.figure(figsize=(15,15))
 
diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py
index 8453d20..9dbcf3e 100755
--- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py
+++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py
@@ -241,7 +241,7 @@ if cvx_not_installable:
         solver = MOSEK
     else:
         solver = SCS      
-    
+        
     # fidelity
     if noise == 's&p':
         fidelity = pnorm( u - noisy_data.as_array(),1)
diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising.py
index 74e7901..d848b9f 100755
--- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising.py
+++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising.py
@@ -239,7 +239,7 @@ if cvx_not_installable:
         solver = MOSEK
     else:
         solver = SCS      
-    
+
     # fidelity
     if noise == 's&p':
         fidelity = pnorm( u - noisy_data.as_array(),1)
@@ -248,8 +248,7 @@ if cvx_not_installable:
         solver = SCS
     elif noise == 'gaussian':
         fidelity = 0.5 * sum_squares(noisy_data.as_array() - u)
-    
-            
+                
     obj =  Minimize( regulariser +  fidelity)
     prob = Problem(obj)
     result = prob.solve(verbose = True, solver = solver)
diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising_2D_time.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising_2D_time.py
new file mode 100644
index 0000000..febe76d
--- /dev/null
+++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising_2D_time.py
@@ -0,0 +1,243 @@
+#========================================================================
+# Copyright 2019 Science Technology Facilities Council
+# Copyright 2019 University of Manchester
+#
+# This work is part of the Core Imaging Library developed by Science Technology
+# Facilities Council and University of Manchester
+#
+# 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.txt
+#
+# 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.
+#
+#=========================================================================
+""" 
+
+Total Variation (Dynamic) Denoising using PDHG algorithm and Tomophantom:
+
+
+Problem:     min_{x} \alpha * ||\nabla x||_{2,1} + \frac{1}{2} * || x - g ||_{2}^{2}
+
+             \alpha: Regularization parameter
+             
+             \nabla: Gradient operator 
+             
+             g: 2D Dynamic noisy data with Gaussian Noise
+                          
+             K = \nabla  
+                                                                
+"""
+
+from ccpi.framework import ImageData, ImageGeometry
+
+import numpy as np 
+import numpy                          
+import matplotlib.pyplot as plt
+
+from ccpi.optimisation.algorithms import PDHG
+
+from ccpi.optimisation.operators import BlockOperator, Gradient, Identity
+from ccpi.optimisation.functions import ZeroFunction, L2NormSquared,  \
+                      MixedL21Norm, BlockFunction
+
+from ccpi.astra.ops import AstraProjectorMC
+
+import os
+import tomophantom
+from tomophantom import TomoP2D
+
+# Create phantom for TV 2D dynamic tomography 
+
+model = 102  # note that the selected model is temporal (2D + time)
+N = 128 # set dimension of the phantom
+
+path = os.path.dirname(tomophantom.__file__)
+path_library2D = os.path.join(path, "Phantom2DLibrary.dat")
+phantom_2Dt = TomoP2D.ModelTemporal(model, N, path_library2D)
+
+plt.close('all')
+plt.figure(1)
+plt.rcParams.update({'font.size': 21})
+plt.title('{}''{}'.format('2D+t phantom using model no.',model))
+for sl in range(0,np.shape(phantom_2Dt)[0]):
+    im = phantom_2Dt[sl,:,:]
+    plt.imshow(im, vmin=0, vmax=1)
+#    plt.pause(.1)
+#    plt.draw
+
+# Setup geometries    
+ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N, channels = np.shape(phantom_2Dt)[0])
+data = ImageData(phantom_2Dt, geometry=ig)
+ag = ig
+                      
+# Create noisy data. Apply Gaussian noise
+np.random.seed(10)
+noisy_data = ImageData( data.as_array() + np.random.normal(0, 0.25, size=ig.shape) )
+
+# time-frames index
+tindex = [8, 16, 24]
+
+fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(10, 10))
+plt.subplot(1,3,1)
+plt.imshow(noisy_data.as_array()[tindex[0],:,:])
+plt.axis('off')
+plt.title('Time {}'.format(tindex[0]))
+plt.subplot(1,3,2)
+plt.imshow(noisy_data.as_array()[tindex[1],:,:])
+plt.axis('off')
+plt.title('Time {}'.format(tindex[1]))
+plt.subplot(1,3,3)
+plt.imshow(noisy_data.as_array()[tindex[2],:,:])
+plt.axis('off')
+plt.title('Time {}'.format(tindex[2]))
+
+fig.subplots_adjust(bottom=0.1, top=0.9, left=0.1, right=0.8,
+                    wspace=0.02, hspace=0.02)
+
+plt.show()   
+
+# Regularisation Parameter
+alpha = 0.3
+
+# Create operators
+op1 = Gradient(ig, correlation='Space')
+op2 = Gradient(ig, correlation='SpaceChannels')
+op3 = Identity(ig, ag)
+
+# Create BlockOperator
+operator1 = BlockOperator(op1, op3, shape=(2,1) ) 
+operator2 = BlockOperator(op2, op3, shape=(2,1) )
+
+# Create functions
+      
+f1 = alpha * MixedL21Norm()
+f2 = 0.5 * L2NormSquared(b = noisy_data)    
+f = BlockFunction(f1, f2)  
+                                                                        
+g = ZeroFunction()
+    
+# Compute operator Norm
+normK1 = operator1.norm()
+normK2 = operator2.norm()
+
+#%%
+# Primal & dual stepsizes
+sigma1 = 1
+tau1 = 1/(sigma1*normK1**2)
+
+sigma2 = 1
+tau2 = 1/(sigma2*normK2**2)
+
+# Setup and run the PDHG algorithm
+pdhg1 = PDHG(f=f,g=g,operator=operator1, tau=tau1, sigma=sigma1)
+pdhg1.max_iteration = 2000
+pdhg1.update_objective_interval = 200
+pdhg1.run(2000)
+
+# Setup and run the PDHG algorithm
+pdhg2 = PDHG(f=f,g=g,operator=operator2, tau=tau2, sigma=sigma2)
+pdhg2.max_iteration = 2000
+pdhg2.update_objective_interval = 200
+pdhg2.run(2000)
+
+
+#%%
+
+tindex = [8, 16, 24]
+fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(10, 8))
+
+plt.subplot(3,3,1)
+plt.imshow(phantom_2Dt[tindex[0],:,:])
+plt.axis('off')
+plt.title('Time {}'.format(tindex[0]))
+
+plt.subplot(3,3,2)
+plt.imshow(phantom_2Dt[tindex[1],:,:])
+plt.axis('off')
+plt.title('Time {}'.format(tindex[1]))
+
+plt.subplot(3,3,3)
+plt.imshow(phantom_2Dt[tindex[2],:,:])
+plt.axis('off')
+plt.title('Time {}'.format(tindex[2]))
+
+plt.subplot(3,3,4)
+plt.imshow(pdhg1.get_output().as_array()[tindex[0],:,:])
+plt.axis('off')
+plt.subplot(3,3,5)
+plt.imshow(pdhg1.get_output().as_array()[tindex[1],:,:])
+plt.axis('off')
+plt.subplot(3,3,6)
+plt.imshow(pdhg1.get_output().as_array()[tindex[2],:,:])
+plt.axis('off')
+
+
+plt.subplot(3,3,7)
+plt.imshow(pdhg2.get_output().as_array()[tindex[0],:,:])
+plt.axis('off')
+plt.subplot(3,3,8)
+plt.imshow(pdhg2.get_output().as_array()[tindex[1],:,:])
+plt.axis('off')
+plt.subplot(3,3,9)
+plt.imshow(pdhg2.get_output().as_array()[tindex[2],:,:])
+plt.axis('off')
+
+
+im = plt.imshow(pdhg1.get_output().as_array()[tindex[0],:,:])
+
+
+fig.subplots_adjust(bottom=0.1, top=0.9, left=0.1, right=0.8,
+                    wspace=0.02, hspace=0.02)
+
+cb_ax = fig.add_axes([0.83, 0.1, 0.02, 0.8])
+cbar = fig.colorbar(im, cax=cb_ax)
+plt.show()
+
+#%%
+import matplotlib.animation as animation
+fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(10, 30))
+ims1 = []
+ims2 = []
+ims3 = []
+for sl in range(0,np.shape(phantom_2Dt)[0]):
+    
+    plt.subplot(1,3,1)
+    im1 = plt.imshow(phantom_2Dt[sl,:,:], animated=True)
+    
+    plt.subplot(1,3,2)
+    im2 = plt.imshow(pdhg1.get_output().as_array()[sl,:,:])
+    
+    plt.subplot(1,3,3)
+    im3 = plt.imshow(pdhg2.get_output().as_array()[sl,:,:])
+    
+    ims1.append([im1])
+    ims2.append([im2])
+    ims3.append([im3])
+    
+    
+ani1 = animation.ArtistAnimation(fig, ims1, interval=500,
+                                repeat_delay=10)
+
+ani2 = animation.ArtistAnimation(fig, ims2, interval=500,
+                                repeat_delay=10)
+
+ani3 = animation.ArtistAnimation(fig, ims3, interval=500,
+                                repeat_delay=10)
+plt.show()        
+#    plt.pause(0.25)
+#    plt.show()
+
+
+
+
+
+
+
+
diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian_3D.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising_Gaussian_3D.py
index c86ddc9..15709cd 100644
--- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian_3D.py
+++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising_Gaussian_3D.py
@@ -1,42 +1,58 @@
-# -*- 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-2019 Evangelos Papoutsellis and 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 ccpi.framework import ImageData, ImageGeometry
-
-import numpy as np 
-import numpy                          
+#========================================================================
+# Copyright 2019 Science Technology Facilities Council
+# Copyright 2019 University of Manchester
+#
+# This work is part of the Core Imaging Library developed by Science Technology
+# Facilities Council and University of Manchester
+#
+# 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.txt
+#
+# 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.
+#
+#=========================================================================
+""" 
+
+Total Variation (3D) Denoising using PDHG algorithm and Tomophantom:
+
+
+Problem:     min_{x} \alpha * ||\nabla x||_{2,1} + \frac{1}{2} * || x - g ||_{2}^{2}
+
+             \alpha: Regularization parameter
+             
+             \nabla: Gradient operator 
+             
+             g: 3D Noisy Data with Gaussian Noise
+                          
+             Method = 0 ( PDHG - split ) :  K = [ \nabla,
+                                                 Identity]
+                          
+                                                                    
+             Method = 1 (PDHG - explicit ):  K = \nabla  
+                                                                
+"""
+
+from ccpi.framework import ImageData, ImageGeometry                       
 import matplotlib.pyplot as plt
-
 from ccpi.optimisation.algorithms import PDHG
-
-from ccpi.optimisation.operators import BlockOperator, Identity, Gradient
-from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \
-                      MixedL21Norm, BlockFunction
+from ccpi.optimisation.operators import Gradient
+from ccpi.optimisation.functions import L2NormSquared, MixedL21Norm
 
 from skimage.util import random_noise
 
-# Create phantom for TV Gaussian denoising
 import timeit
 import os
 from tomophantom import TomoP3D
 import tomophantom
 
+# Create a phantom from Tomophantom
 print ("Building 3D phantom using TomoPhantom software")
 tic=timeit.default_timer()
 model = 13 # select a model number from the library
@@ -47,12 +63,13 @@ path_library3D = os.path.join(path, "Phantom3DLibrary.dat")
 #This will generate a N x N x N phantom (3D)
 phantom_tm = TomoP3D.Model(model, N, path_library3D)
 
-# Create noisy data. Add Gaussian noise
+# Create noisy data. Apply Gaussian noise
 ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N, voxel_num_z=N)
 ag = ig
 n1 = random_noise(phantom_tm, mode = 'gaussian', mean=0, var = 0.001, seed=10)
 noisy_data = ImageData(n1)
 
+# Show results
 sliceSel = int(0.5*N)
 plt.figure(figsize=(15,15))
 plt.subplot(3,1,1)
@@ -69,52 +86,27 @@ plt.title('Sagittal View')
 plt.colorbar()
 plt.show()   
 
-
 # Regularisation Parameter
 alpha = 0.05
-
-method = '0'
-
-if method == '0':
-
-    # Create operators
-    op1 = Gradient(ig)
-    op2 = Identity(ig, ag)
-
-    # Create BlockOperator
-    operator = BlockOperator(op1, op2, shape=(2,1) ) 
-
-    # Create functions
-      
-    f1 = alpha * MixedL21Norm()
-    f2 = 0.5 * L2NormSquared(b = noisy_data)    
-    f = BlockFunction(f1, f2)  
-                                      
-    g = ZeroFunction()
-    
-else:
     
-    # Without the "Block Framework"
-    operator = Gradient(ig)
-    f =  alpha * MixedL21Norm()
-    g =  0.5 * L2NormSquared(b = noisy_data)
-        
-    
-# Compute operator Norm
+# Setup and run the PDHG algorithm
+operator = Gradient(ig)
+f =  alpha * MixedL21Norm()
+g =  0.5 * L2NormSquared(b = noisy_data)
+            
 normK = operator.norm()
 
-# Primal & dual stepsizes
 sigma = 1
 tau = 1/(sigma*normK**2)
 
-
-# Setup and run the PDHG algorithm
 pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True)
-pdhg.max_iteration = 2000
+pdhg.max_iteration = 1000
 pdhg.update_objective_interval = 200
-pdhg.run(2000)
+pdhg.run(1000, verbose = True)
 
+# Show results
 fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(10, 8))
+fig.suptitle('TV Reconstruction',fontsize=20)
 
 plt.subplot(2,3,1)
 plt.imshow(noisy_data.as_array()[sliceSel,:,:],vmin=0, vmax=1)
diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Tomo2D.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Tomo2D.py
new file mode 100644
index 0000000..f179e95
--- /dev/null
+++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Tomo2D.py
@@ -0,0 +1,173 @@
+#========================================================================
+# Copyright 2019 Science Technology Facilities Council
+# Copyright 2019 University of Manchester
+#
+# This work is part of the Core Imaging Library developed by Science Technology
+# Facilities Council and University of Manchester
+#
+# 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.txt
+#
+# 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.
+#
+#=========================================================================
+
+""" 
+
+Total Variation 2D Tomography Reconstruction using PDHG algorithm:
+
+
+Problem:     min_u  \alpha * ||\nabla u||_{2,1} + \frac{1}{2}||Au - g||^{2}
+             min_u, u>0  \alpha * ||\nabla u||_{2,1} + \int A u  - g log (Au + \eta)
+
+             \nabla: Gradient operator              
+             A: System Matrix
+             g: Noisy sinogram 
+             \eta: Background noise
+             
+             \alpha: Regularization parameter
+ 
+"""
+
+from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData
+
+import numpy as np 
+import numpy                          
+import matplotlib.pyplot as plt
+
+from ccpi.optimisation.algorithms import PDHG
+
+from ccpi.optimisation.operators import BlockOperator, Gradient
+from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \
+                      MixedL21Norm, BlockFunction, KullbackLeibler, IndicatorBox
+                      
+from ccpi.astra.ops import AstraProjectorSimple
+from ccpi.framework import TestData
+from PIL import Image
+import os, sys
+if int(numpy.version.version.split('.')[1]) > 12:
+    from skimage.util import random_noise
+else:
+    from demoutil import random_noise
+
+import scipy.io
+
+# user supplied input
+if len(sys.argv) > 1:
+    which_noise = int(sys.argv[1])
+else:
+    which_noise = 1
+
+# Load 256 shepp-logan
+data256 = scipy.io.loadmat('phantom.mat')['phantom256']
+data = ImageData(numpy.array(Image.fromarray(data256).resize((256,256))))
+N, M = data.shape
+ig = ImageGeometry(voxel_num_x=N, voxel_num_y=M)
+
+# Add it to testdata or use tomophantom
+#loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
+#data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(50, 50))
+#ig = data.geometry
+
+# Create acquisition data and geometry
+detectors = N
+angles = np.linspace(0, np.pi, 180)
+ag = AcquisitionGeometry('parallel','2D',angles, detectors)
+
+# Select device
+device = '0'
+#device = input('Available device: GPU==1 / CPU==0 ')
+if device=='1':
+    dev = 'gpu'
+else:
+    dev = 'cpu'
+    
+Aop = AstraProjectorSimple(ig, ag, dev)
+sin = Aop.direct(data)
+
+# Create noisy data. Apply Gaussian noise
+noises = ['gaussian', 'poisson']
+noise = noises[which_noise]
+
+if noise == 'poisson':
+    scale = 5
+    eta = 0
+    noisy_data = AcquisitionData(np.random.poisson( scale * (eta + sin.as_array()))/scale, ag)
+elif noise == 'gaussian':
+    n1 = np.random.normal(0, 1, size = ag.shape)
+    noisy_data = AcquisitionData(n1 + sin.as_array(), ag)
+else:
+    raise ValueError('Unsupported Noise ', noise)
+
+# Show Ground Truth and Noisy Data
+plt.figure(figsize=(10,10))
+plt.subplot(1,2,2)
+plt.imshow(data.as_array())
+plt.title('Ground Truth')
+plt.colorbar()
+plt.subplot(1,2,1)
+plt.imshow(noisy_data.as_array())
+plt.title('Noisy Data')
+plt.colorbar()
+plt.show()
+
+# Create operators
+op1 = Gradient(ig)
+op2 = Aop
+
+# Create BlockOperator
+operator = BlockOperator(op1, op2, shape=(2,1) ) 
+
+# Compute operator Norm
+normK = operator.norm()
+
+# Create functions
+if noise == 'poisson':
+    alpha = 3
+    f2 = KullbackLeibler(noisy_data)  
+    g =  IndicatorBox(lower=0) 
+    sigma = 1
+    tau = 1/(sigma*normK**2)    
+    
+elif noise == 'gaussian':   
+    alpha = 20
+    f2 = 0.5 * L2NormSquared(b=noisy_data)                                         
+    g = ZeroFunction()
+    sigma = 10
+    tau = 1/(sigma*normK**2)     
+    
+f1 = alpha * MixedL21Norm() 
+f = BlockFunction(f1, f2)    
+
+# Setup and run the PDHG algorithm
+pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
+pdhg.max_iteration = 2000
+pdhg.update_objective_interval = 200
+pdhg.run(2000)
+
+plt.figure(figsize=(15,15))
+plt.subplot(3,1,1)
+plt.imshow(data.as_array())
+plt.title('Ground Truth')
+plt.colorbar()
+plt.subplot(3,1,2)
+plt.imshow(noisy_data.as_array())
+plt.title('Noisy Data')
+plt.colorbar()
+plt.subplot(3,1,3)
+plt.imshow(pdhg.get_output().as_array())
+plt.title('TV Reconstruction')
+plt.colorbar()
+plt.show()
+plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), data.as_array()[int(N/2),:], label = 'GTruth')
+plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction')
+plt.legend()
+plt.title('Middle Line Profiles')
+plt.show()
diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py
index e16c5a6..78d4980 100644
--- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py
+++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py
@@ -195,7 +195,6 @@ try:
 except ImportError:
     cvx_not_installable = False
 
-
 if cvx_not_installable:
 
     ##Construct problem    
diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/phantom.mat b/Wrappers/Python/demos/PDHG_examples/GatherAll/phantom.mat
new file mode 100755
index 0000000..c465bbe
--- /dev/null
+++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/phantom.mat
diff --git a/Wrappers/Python/demos/PDHG_examples/IMATDemo.py b/Wrappers/Python/demos/PDHG_examples/IMATDemo.py
new file mode 100644
index 0000000..2051860
--- /dev/null
+++ b/Wrappers/Python/demos/PDHG_examples/IMATDemo.py
@@ -0,0 +1,339 @@
+
+
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Mar 25 12:50:27 2019
+
+@author: vaggelis
+"""
+
+from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer, AcquisitionGeometry, AcquisitionData
+from astropy.io import fits
+import numpy
+import matplotlib.pyplot as plt
+
+from ccpi.optimisation.algorithms import CGLS, PDHG
+from ccpi.optimisation.functions import MixedL21Norm, L2NormSquared, BlockFunction, ZeroFunction, KullbackLeibler, IndicatorBox
+from ccpi.optimisation.operators import Gradient, BlockOperator
+
+from ccpi.astra.operators import  AstraProjectorMC, AstraProjectorSimple
+
+import pickle
+
+
+# load file
+
+#filename_sino = '/media/newhd/shared/DataProcessed/IMAT_beamtime_Feb_2019/preprocessed_test_flat/sino/rebin_slice_350/sino_log_rebin_282.fits' 
+#filename_sino = '/media/newhd/shared/DataProcessed/IMAT_beamtime_Feb_2019/preprocessed_test_flat/sino/rebin_slice_350/sino_log_rebin_564.fits'
+#filename_sino = '/media/newhd/shared/DataProcessed/IMAT_beamtime_Feb_2019/preprocessed_test_flat/sino/rebin_slice_350/sino_log_rebin_141.fits'
+filename_sino = '/media/newhd/shared/DataProcessed/IMAT_beamtime_Feb_2019/preprocessed_test_flat/sino/rebin_slice_350/sino_log_rebin_80_channels.fits'
+
+sino_handler = fits.open(filename_sino)
+sino = numpy.array(sino_handler[0].data, dtype=float)
+
+# change axis order: channels, angles, detectors
+sino_new = numpy.rollaxis(sino, 2)
+sino_handler.close()
+
+
+sino_shape = sino_new.shape
+
+num_channels = sino_shape[0] # channelss
+num_pixels_h = sino_shape[2] # detectors
+num_pixels_v = sino_shape[2] # detectors
+num_angles = sino_shape[1] # angles
+
+
+ig = ImageGeometry(voxel_num_x = num_pixels_h, voxel_num_y = num_pixels_v, channels = num_channels)
+
+with open("/media/newhd/vaggelis/CCPi/IMAT_reconstruction/CCPi-Framework/Wrappers/Python/ccpi/optimisation/IMAT_data/golden_angles_new.txt") as f:
+    angles_string = [line.rstrip() for line in f]
+    angles = numpy.array(angles_string).astype(float)
+                                                                                 
+
+ag = AcquisitionGeometry('parallel', '2D',  angles * numpy.pi / 180, pixel_num_h = num_pixels_h, channels = num_channels)
+op_MC = AstraProjectorMC(ig, ag, 'gpu')
+
+sino_aqdata = AcquisitionData(sino_new, ag)
+result_bp = op_MC.adjoint(sino_aqdata)
+
+#%%
+
+channel = [40, 60]
+for j in range(2):        
+    z4 = sino_aqdata.as_array()[channel[j]] 
+    plt.figure(figsize=(10,6))
+    plt.imshow(z4, cmap='viridis')
+    plt.axis('off')
+    plt.savefig('Sino_141/Sinogram_ch_{}_.png'.format(channel[j]), bbox_inches='tight', transparent=True)
+    plt.show() 
+
+#%%
+
+def callback(iteration, objective, x):
+    plt.imshow(x.as_array()[40])
+    plt.colorbar()
+    plt.show()
+
+#%%
+# CGLS 
+
+x_init = ig.allocate()      
+cgls1 = CGLS(x_init=x_init, operator=op_MC, data=sino_aqdata)
+cgls1.max_iteration = 100
+cgls1.update_objective_interval = 2
+cgls1.run(20,verbose=True, callback=callback)
+
+plt.imshow(cgls1.get_output().subset(channel=20).array)
+plt.title('CGLS')
+plt.colorbar()
+plt.show()
+
+#%%
+with open('Sino_141/CGLS/CGLS_{}_iter.pkl'.format(20), 'wb') as f:
+    z = cgls1.get_output()
+    pickle.dump(z, f)  
+    
+#%% 
+#% Tikhonov Space
+    
+x_init = ig.allocate() 
+alpha = [1,3,5,10,20,50]
+
+for a in alpha:
+    
+    Grad = Gradient(ig, correlation = Gradient.CORRELATION_SPACE)
+    operator = BlockOperator(op_MC, a * Grad, shape=(2,1))
+    blockData = BlockDataContainer(sino_aqdata, \
+                                   Grad.range_geometry().allocate())
+    cgls2 = CGLS()
+    cgls2.max_iteration = 500
+    cgls2.set_up(x_init, operator, blockData)
+    cgls2.update_objective_interval = 50
+    cgls2.run(100,verbose=True)
+        
+    with open('Sino_141/CGLS_Space/CGLS_a_{}.pkl'.format(a), 'wb') as f:
+        z = cgls2.get_output()
+        pickle.dump(z, f)
+   
+#% Tikhonov SpaceChannels
+
+for a1 in alpha:
+    
+    Grad1 = Gradient(ig, correlation = Gradient.CORRELATION_SPACECHANNEL)
+    operator1 = BlockOperator(op_MC, a1 * Grad1, shape=(2,1))
+    blockData1 = BlockDataContainer(sino_aqdata, \
+                                   Grad1.range_geometry().allocate())
+    cgls3 = CGLS()
+    cgls3.max_iteration = 500
+    cgls3.set_up(x_init, operator1, blockData1)
+    cgls3.update_objective_interval = 10
+    cgls3.run(100, verbose=True)
+    
+    with open('Sino_141/CGLS_SpaceChannels/CGLS_a_{}.pkl'.format(a1), 'wb') as f1:
+        z1 = cgls3.get_output()
+        pickle.dump(z1, f1)
+    
+  
+    
+#%%
+#
+ig_tmp = ImageGeometry(voxel_num_x = num_pixels_h, voxel_num_y = num_pixels_v)
+ag_tmp = AcquisitionGeometry('parallel', '2D',  angles * numpy.pi / 180, pixel_num_h = num_pixels_h)
+op_tmp = AstraProjectorSimple(ig_tmp, ag_tmp, 'gpu')
+normK1 = op_tmp.norm()
+
+alpha_TV = [2, 5, 10] # for powder
+
+# Create operators
+op1 = Gradient(ig, correlation=Gradient.CORRELATION_SPACECHANNEL)
+op2 = op_MC
+
+# Create BlockOperator
+operator = BlockOperator(op1, op2, shape=(2,1) ) 
+
+
+for alpha in alpha_TV:
+# Create functions      
+    f1 = alpha * MixedL21Norm()
+    
+    f2 = KullbackLeibler(sino_aqdata)    
+    f = BlockFunction(f1, f2)  
+    g = IndicatorBox(lower=0)
+    
+    # Compute operator Norm
+    normK = numpy.sqrt(8 + normK1**2) 
+    
+    # Primal & dual stepsizes
+    sigma = 1
+    tau = 1/(sigma*normK**2)
+    
+    # Setup and run the PDHG algorithm
+    pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
+    pdhg.max_iteration = 5000
+    pdhg.update_objective_interval = 500
+#    pdhg.run(2000, verbose=True, callback=callback)   
+    pdhg.run(5000, verbose=True, callback=callback)  
+#    
+    with open('Sino_141/TV_SpaceChannels/TV_a = {}.pkl'.format(alpha), 'wb') as f3:
+        z3 = pdhg.get_output()
+        pickle.dump(z3, f3)    
+#    
+#
+#
+#
+##%%
+#        
+#ig_tmp = ImageGeometry(voxel_num_x = num_pixels_h, voxel_num_y = num_pixels_v)
+#ag_tmp = AcquisitionGeometry('parallel', '2D',  angles * numpy.pi / 180, pixel_num_h = num_pixels_h)
+#op_tmp = AstraProjectorSimple(ig_tmp, ag_tmp, 'gpu')
+#normK1 = op_tmp.norm()
+#
+#alpha_TV = 10 # for powder
+#
+## Create operators
+#op1 = Gradient(ig, correlation=Gradient.CORRELATION_SPACECHANNEL)
+#op2 = op_MC
+#
+## Create BlockOperator
+#operator = BlockOperator(op1, op2, shape=(2,1) ) 
+#
+#
+## Create functions      
+#f1 = alpha_TV * MixedL21Norm()
+#f2 = 0.5 * L2NormSquared(b=sino_aqdata)    
+#f = BlockFunction(f1, f2)  
+#g = ZeroFunction()
+#
+## Compute operator Norm
+##normK = 8.70320267279591 # For powder Run one time no need to compute again takes time
+#normK = numpy.sqrt(8 + normK1**2) # for carbon
+#
+## Primal & dual stepsizes
+#sigma = 0.1
+#tau = 1/(sigma*normK**2)
+#
+#def callback(iteration, objective, x):
+#    plt.imshow(x.as_array()[100])
+#    plt.colorbar()
+#    plt.show()
+#
+## Setup and run the PDHG algorithm
+#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
+#pdhg.max_iteration = 2000
+#pdhg.update_objective_interval = 100
+#pdhg.run(2000, verbose=True)       
+#        
+#        
+#        
+#        
+#        
+#        
+        
+        
+        
+        
+        
+        
+        
+        
+        
+        
+#%%
+    
+#with open('/media/newhd/vaggelis/CCPi/IMAT_reconstruction/CCPi-Framework/Wrappers/Python/ccpi/optimisation/CGLS_Tikhonov/CGLS_Space/CGLS_Space_a = 50.pkl', 'wb') as f:
+#    z = cgls2.get_output()
+#    pickle.dump(z, f)
+#      
+        
+        #%%
+with open('Sino_141/CGLS_Space/CGLS_Space_a_20.pkl', 'rb') as f1:
+    x = pickle.load(f1)
+    
+with open('Sino_141/CGLS_SpaceChannels/CGLS_SpaceChannels_a_20.pkl', 'rb') as f1:
+    x1 = pickle.load(f1)    
+    
+   
+    
+#    
+plt.imshow(x.as_array()[40]*mask)
+plt.colorbar()
+plt.show()
+
+plt.imshow(x1.as_array()[40]*mask)
+plt.colorbar()
+plt.show()
+
+plt.plot(x.as_array()[40,100,:])
+plt.plot(x1.as_array()[40,100,:])
+plt.show()
+
+#%%
+
+# Show results
+
+def circ_mask(h, w, center=None, radius=None):
+
+    if center is None: # use the middle of the image
+        center = [int(w/2), int(h/2)]
+    if radius is None: # use the smallest distance between the center and image walls
+        radius = min(center[0], center[1], w-center[0], h-center[1])
+
+    Y, X = numpy.ogrid[:h, :w]
+    dist_from_center = numpy.sqrt((X - center[0])**2 + (Y-center[1])**2)
+
+    mask = dist_from_center <= radius
+    return mask
+
+mask = circ_mask(141, 141, center=None, radius = 55)
+plt.imshow(numpy.multiply(x.as_array()[40],mask))
+plt.show()
+#%%
+#channel = [100, 200, 300]
+#
+#for i in range(3):     
+#    tmp = cgls1.get_output().as_array()[channel[i]] 
+#    
+#    z =  tmp * mask
+#    plt.figure(figsize=(10,6))
+#    plt.imshow(z, vmin=0, cmap='viridis')
+#    plt.axis('off')
+##    plt.clim(0, 0.02)
+##    plt.colorbar()
+##    del z
+#    plt.savefig('CGLS_282/CGLS_Chan_{}.png'.format(channel[i]), bbox_inches='tight', transparent=True)
+#    plt.show()
+#    
+#    
+##%% Line Profiles
+#    
+#n1, n2, n3 =  cgs.get_output().as_array().shape
+#mask = circ_mask(564, 564, center=None, radius = 220)
+#material = ['Cu', 'Fe', 'Ni']
+#ycoords = [200, 300, 380]
+#    
+#for i in range(3):
+#    z = cgs.get_output().as_array()[channel[i]]  * mask
+#    
+#    for k1 in range(len(ycoords)):
+#        plt.plot(numpy.arange(0,n2), z[ycoords[k1],:])
+#        plt.title('Channel {}: {}'.format(channel[i], material[k1]))
+#        plt.savefig('CGLS/line_profile_chan_{}_material_{}.png'.\
+#                        format(channel[i], material[k1]), bbox_inches='tight')
+#        plt.show()
+#    
+#    
+#
+#
+#
+##%%
+#
+#%%
+        
+
+
+#%%
+
+#plt.imshow(pdhg.get_output().subset(channel=100).as_array())
+#plt.show()
diff --git a/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Color_Denoising.py b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Color_Denoising.py
new file mode 100644
index 0000000..ddf5ace
--- /dev/null
+++ b/Wrappers/Python/demos/PDHG_examples/PDHG_TV_Color_Denoising.py
@@ -0,0 +1,115 @@
+#========================================================================
+# Copyright 2019 Science Technology Facilities Council
+# Copyright 2019 University of Manchester
+#
+# This work is part of the Core Imaging Library developed by Science Technology
+# Facilities Council and University of Manchester
+#
+# 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.txt
+#
+# 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.
+#
+#=========================================================================
+
+""" 
+Total Variation Denoising using PDHG algorithm:
+Problem:     min_x, x>0  \alpha * ||\nabla x||_{2,1} + ||x-g||_{1}
+             \alpha: Regularization parameter
+             
+             \nabla: Gradient operator 
+             
+             g: Noisy Data with Salt & Pepper Noise
+             
+             
+             Method = 0 ( PDHG - split ) :  K = [ \nabla,
+                                                 Identity]
+                          
+                                                                    
+             Method = 1 (PDHG - explicit ):  K = \nabla    
+             
+             
+"""
+
+import numpy                          
+import matplotlib.pyplot as plt
+
+from ccpi.optimisation.algorithms import PDHG
+
+from ccpi.optimisation.operators import Gradient, BlockOperator, FiniteDiff
+from ccpi.optimisation.functions import MixedL21Norm, MixedL11Norm, L2NormSquared, BlockFunction, L1Norm                      
+from ccpi.framework import TestData, ImageGeometry
+import os, sys
+if int(numpy.version.version.split('.')[1]) > 12:
+    from skimage.util import random_noise
+else:
+    from demoutil import random_noise
+
+loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
+data = loader.load(TestData.PEPPERS, size=(256,256))
+ig = data.geometry
+ag = ig
+
+# Create noisy data. 
+n1 = random_noise(data.as_array(), mode = 'gaussian', var = 0.15, seed = 50)
+noisy_data = ig.allocate()
+noisy_data.fill(n1)
+
+# Show Ground Truth and Noisy Data
+plt.figure(figsize=(10,5))
+plt.subplot(1,2,1)
+plt.imshow(data.as_array())
+plt.title('Ground Truth')
+plt.colorbar()
+plt.subplot(1,2,2)
+plt.imshow(noisy_data.as_array())
+plt.title('Noisy Data')
+plt.colorbar()
+plt.show()
+
+# Regularisation Parameter
+operator = Gradient(ig, correlation=Gradient.CORRELATION_SPACE)
+f1 =  5 * MixedL21Norm()
+g = 0.5 * L2NormSquared(b=noisy_data)
+            
+# Compute operator Norm
+normK = operator.norm()
+
+# Primal & dual stepsizes
+sigma = 1
+tau = 1/(sigma*normK**2)
+
+# Setup and run the PDHG algorithm
+pdhg1 = PDHG(f=f1,g=g,operator=operator, tau=tau, sigma=sigma)
+pdhg1.max_iteration = 2000
+pdhg1.update_objective_interval = 200
+pdhg1.run(1000)
+
+
+# Show results
+plt.figure(figsize=(10,10))
+plt.subplot(2,2,1)
+plt.imshow(data.as_array())
+plt.title('Ground Truth')
+plt.colorbar()
+plt.subplot(2,2,2)
+plt.imshow(noisy_data.as_array())
+plt.title('Noisy Data')
+plt.colorbar()
+plt.subplot(2,2,3)
+plt.imshow(pdhg1.get_output().as_array())
+plt.title('TV Reconstruction')
+plt.colorbar()
+plt.subplot(2,2,4)
+plt.imshow(pdhg2.get_output().as_array())
+plt.title('TV Reconstruction')
+plt.colorbar()
+plt.show()
+
diff --git a/Wrappers/Python/demos/PDHG_examples/TV_Denoising/.DS_Store b/Wrappers/Python/demos/PDHG_examples/TV_Denoising/.DS_Store
deleted file mode 100644
index 5008ddf..0000000
--- a/Wrappers/Python/demos/PDHG_examples/TV_Denoising/.DS_Store
+++ /dev/null
diff --git a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py
index bc0081f..e74e1c6 100644
--- a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py
+++ b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py
@@ -23,18 +23,20 @@
 Total Generalised Variation (TGV) Tomography 2D using PDHG algorithm:
 
 
-Problem:     min_{x>0} \alpha * ||\nabla x - w||_{2,1} +
-                     \beta *  || E w ||_{2,1} +
-                     int A x - g log(Ax + \eta)
+Problem:     min_{x>0} \alpha * ||\nabla x - w||_{2,1} + \beta *  || E w ||_{2,1} +
+                       \frac{1}{2}||Au - g||^{2}
+                     
+            min_{u>0} \alpha * ||\nabla u - w||_{2,1} + \beta *  || E w ||_{2,1} +
+                      int A u - g log(Au + \eta)                     
 
              \alpha: Regularization parameter
              \beta: Regularization parameter
              
              \nabla: Gradient operator 
               E: Symmetrized Gradient operator
-              A: Projection Matrix
+              A: System Matrix
              
-             g: Noisy Data with Poisson Noise
+             g: Noisy Sinogram 
                           
               K = [ \nabla, - Identity
                    ZeroOperator, E 
@@ -53,54 +55,41 @@ from ccpi.optimisation.algorithms import PDHG
 from ccpi.optimisation.operators import BlockOperator, Gradient, Identity, \
                                     SymmetrizedGradient, ZeroOperator
 from ccpi.optimisation.functions import IndicatorBox, KullbackLeibler, ZeroFunction,\
-                      MixedL21Norm, BlockFunction
+                      MixedL21Norm, BlockFunction, L2NormSquared
 
 from ccpi.astra.ops import AstraProjectorSimple
-from ccpi.framework import TestData
 import os, sys
-from skimage.util import random_noise
 
-loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
 
-# Load Data                      
-#N = 50
-#M = 50
-#data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,M), scale=(0,1))
-#
-#ig = data.geometry
-#ag = ig
-N = 100
-
-data = np.zeros((N,N))
-
-x1 = np.linspace(0, int(N/2), N)
-x2 = np.linspace(int(N/2), 0., N)
-xv, yv = np.meshgrid(x1, x2)
+import tomophantom
+from tomophantom import TomoP2D
 
-xv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1] = yv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1].T
-
-data = xv
-data = ImageData(data/data.max())
+# user supplied input
+if len(sys.argv) > 1:
+    which_noise = int(sys.argv[1])
+else:
+    which_noise = 1 
+    
+# Load Piecewise smooth Shepp-Logan phantom 
+model = 2 # select a model number from the library
+N = 128 # set dimension of the phantom
+# one can specify an exact path to the parameters file
+# path_library2D = '../../../PhantomLibrary/models/Phantom2DLibrary.dat'
+path = os.path.dirname(tomophantom.__file__)
+path_library2D = os.path.join(path, "Phantom2DLibrary.dat")
+#This will generate a N_size x N_size phantom (2D)
+phantom_2D = TomoP2D.Model(model, N, path_library2D)
 
 ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
-ag = ig
-
-# Create noisy data. Add Gaussian noise
-scale = 0.5
-eta = 0 
-n1 = scale * np.random.poisson(eta + sin.as_array()/scale)
-
-noisy_data = AcquisitionData(n1, ag)
-
-#Create Acquisition Data and apply poisson noise
+data = ImageData(phantom_2D)
 
+#Create Acquisition Data 
 detectors = N
 angles = np.linspace(0, np.pi, N)
-
 ag = AcquisitionGeometry('parallel','2D',angles, detectors)
 
 #device = input('Available device: GPU==1 / CPU==0 ')
-device = '0'
+device = '1'
 if device=='1':
     dev = 'gpu'
 else:
@@ -109,28 +98,31 @@ else:
 Aop = AstraProjectorSimple(ig, ag, 'cpu')
 sin = Aop.direct(data)
 
-# Create noisy data. Apply Poisson noise
-scale = 0.5
-eta = 0 
-n1 = scale * np.random.poisson(eta + sin.as_array()/scale)
-
-noisy_data = AcquisitionData(n1, ag)
+# Create noisy sinogram.
+noises = ['gaussian', 'poisson']
+noise = noises[which_noise]
+
+if noise == 'poisson':
+    scale = 5
+    eta = 0
+    noisy_data = AcquisitionData(np.random.poisson( scale * (eta + sin.as_array()))/scale, ag)
+elif noise == 'gaussian':
+    n1 = np.random.normal(0, 1, size = ag.shape)
+    noisy_data = AcquisitionData(n1 + sin.as_array(), ag)
+else:
+    raise ValueError('Unsupported Noise ', noise)
 
 # Show Ground Truth and Noisy Data
 plt.figure(figsize=(10,10))
-plt.subplot(2,1,1)
+plt.subplot(1,2,2)
 plt.imshow(data.as_array())
 plt.title('Ground Truth')
 plt.colorbar()
-plt.subplot(2,1,2)
+plt.subplot(1,2,1)
 plt.imshow(noisy_data.as_array())
 plt.title('Noisy Data')
 plt.colorbar()
 plt.show()
-#%%
-# Regularisation Parameters
-alpha = 1
-beta =  5
 
 # Create Operators
 op11 = Gradient(ig)
@@ -143,28 +135,42 @@ op31 = Aop
 op32 = ZeroOperator(op22.domain_geometry(), ag)
 
 operator = BlockOperator(op11, -1*op12, op21, op22, op31, op32, shape=(3,2) ) 
+normK = operator.norm()
+
+# Create functions
+if noise == 'poisson':
+    alpha = 3
+    beta = 6
+    f3 = KullbackLeibler(noisy_data)    
+    g =  BlockFunction(IndicatorBox(lower=0), ZeroFunction()) 
+    
+    # Primal & dual stepsizes
+    sigma = 1
+    tau = 1/(sigma*normK**2)    
+    
+elif noise == 'gaussian':   
+    alpha = 20
+    beta = 50
+    f3 = 0.5 * L2NormSquared(b=noisy_data)                                         
+    g = BlockFunction(ZeroFunction(), ZeroFunction())
+    
+    # Primal & dual stepsizes
+    sigma = 10
+    tau = 1/(sigma*normK**2)     
     
 f1 = alpha * MixedL21Norm()
-f2 = beta * MixedL21Norm() 
-f3 = KullbackLeibler(noisy_data)    
+f2 = beta * MixedL21Norm()  
 f = BlockFunction(f1, f2, f3)         
-
-g =  BlockFunction(IndicatorBox(lower=0), ZeroFunction())
     
 # Compute operator Norm
 normK = operator.norm()
 
-# Primal & dual stepsizes
-sigma = 1
-tau = 1/(sigma*normK**2)
-
-
 # Setup and run the PDHG algorithm
 pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
 pdhg.max_iteration = 3000
 pdhg.update_objective_interval = 500
 pdhg.run(3000)
-
+#%%
 plt.figure(figsize=(15,15))
 plt.subplot(3,1,1)
 plt.imshow(data.as_array())
@@ -179,9 +185,8 @@ plt.imshow(pdhg.get_output()[0].as_array())
 plt.title('TGV Reconstruction')
 plt.colorbar()
 plt.show()
-## 
-plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth')
-plt.plot(np.linspace(0,N,N), pdhg.get_output()[0].as_array()[int(N/2),:], label = 'TGV reconstruction')
+plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), data.as_array()[:, int(N/2)], label = 'GTruth')
+plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), pdhg.get_output()[0].as_array()[:, int(N/2)], label = 'TGV reconstruction')
 plt.legend()
 plt.title('Middle Line Profiles')
 plt.show()
diff --git a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D.py b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D.py
new file mode 100644
index 0000000..f179e95
--- /dev/null
+++ b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D.py
@@ -0,0 +1,173 @@
+#========================================================================
+# Copyright 2019 Science Technology Facilities Council
+# Copyright 2019 University of Manchester
+#
+# This work is part of the Core Imaging Library developed by Science Technology
+# Facilities Council and University of Manchester
+#
+# 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.txt
+#
+# 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.
+#
+#=========================================================================
+
+""" 
+
+Total Variation 2D Tomography Reconstruction using PDHG algorithm:
+
+
+Problem:     min_u  \alpha * ||\nabla u||_{2,1} + \frac{1}{2}||Au - g||^{2}
+             min_u, u>0  \alpha * ||\nabla u||_{2,1} + \int A u  - g log (Au + \eta)
+
+             \nabla: Gradient operator              
+             A: System Matrix
+             g: Noisy sinogram 
+             \eta: Background noise
+             
+             \alpha: Regularization parameter
+ 
+"""
+
+from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData
+
+import numpy as np 
+import numpy                          
+import matplotlib.pyplot as plt
+
+from ccpi.optimisation.algorithms import PDHG
+
+from ccpi.optimisation.operators import BlockOperator, Gradient
+from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \
+                      MixedL21Norm, BlockFunction, KullbackLeibler, IndicatorBox
+                      
+from ccpi.astra.ops import AstraProjectorSimple
+from ccpi.framework import TestData
+from PIL import Image
+import os, sys
+if int(numpy.version.version.split('.')[1]) > 12:
+    from skimage.util import random_noise
+else:
+    from demoutil import random_noise
+
+import scipy.io
+
+# user supplied input
+if len(sys.argv) > 1:
+    which_noise = int(sys.argv[1])
+else:
+    which_noise = 1
+
+# Load 256 shepp-logan
+data256 = scipy.io.loadmat('phantom.mat')['phantom256']
+data = ImageData(numpy.array(Image.fromarray(data256).resize((256,256))))
+N, M = data.shape
+ig = ImageGeometry(voxel_num_x=N, voxel_num_y=M)
+
+# Add it to testdata or use tomophantom
+#loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
+#data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(50, 50))
+#ig = data.geometry
+
+# Create acquisition data and geometry
+detectors = N
+angles = np.linspace(0, np.pi, 180)
+ag = AcquisitionGeometry('parallel','2D',angles, detectors)
+
+# Select device
+device = '0'
+#device = input('Available device: GPU==1 / CPU==0 ')
+if device=='1':
+    dev = 'gpu'
+else:
+    dev = 'cpu'
+    
+Aop = AstraProjectorSimple(ig, ag, dev)
+sin = Aop.direct(data)
+
+# Create noisy data. Apply Gaussian noise
+noises = ['gaussian', 'poisson']
+noise = noises[which_noise]
+
+if noise == 'poisson':
+    scale = 5
+    eta = 0
+    noisy_data = AcquisitionData(np.random.poisson( scale * (eta + sin.as_array()))/scale, ag)
+elif noise == 'gaussian':
+    n1 = np.random.normal(0, 1, size = ag.shape)
+    noisy_data = AcquisitionData(n1 + sin.as_array(), ag)
+else:
+    raise ValueError('Unsupported Noise ', noise)
+
+# Show Ground Truth and Noisy Data
+plt.figure(figsize=(10,10))
+plt.subplot(1,2,2)
+plt.imshow(data.as_array())
+plt.title('Ground Truth')
+plt.colorbar()
+plt.subplot(1,2,1)
+plt.imshow(noisy_data.as_array())
+plt.title('Noisy Data')
+plt.colorbar()
+plt.show()
+
+# Create operators
+op1 = Gradient(ig)
+op2 = Aop
+
+# Create BlockOperator
+operator = BlockOperator(op1, op2, shape=(2,1) ) 
+
+# Compute operator Norm
+normK = operator.norm()
+
+# Create functions
+if noise == 'poisson':
+    alpha = 3
+    f2 = KullbackLeibler(noisy_data)  
+    g =  IndicatorBox(lower=0) 
+    sigma = 1
+    tau = 1/(sigma*normK**2)    
+    
+elif noise == 'gaussian':   
+    alpha = 20
+    f2 = 0.5 * L2NormSquared(b=noisy_data)                                         
+    g = ZeroFunction()
+    sigma = 10
+    tau = 1/(sigma*normK**2)     
+    
+f1 = alpha * MixedL21Norm() 
+f = BlockFunction(f1, f2)    
+
+# Setup and run the PDHG algorithm
+pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
+pdhg.max_iteration = 2000
+pdhg.update_objective_interval = 200
+pdhg.run(2000)
+
+plt.figure(figsize=(15,15))
+plt.subplot(3,1,1)
+plt.imshow(data.as_array())
+plt.title('Ground Truth')
+plt.colorbar()
+plt.subplot(3,1,2)
+plt.imshow(noisy_data.as_array())
+plt.title('Noisy Data')
+plt.colorbar()
+plt.subplot(3,1,3)
+plt.imshow(pdhg.get_output().as_array())
+plt.title('TV Reconstruction')
+plt.colorbar()
+plt.show()
+plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), data.as_array()[int(N/2),:], label = 'GTruth')
+plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction')
+plt.legend()
+plt.title('Middle Line Profiles')
+plt.show()
diff --git a/Wrappers/Python/demos/PDHG_examples/Tomo/phantom.mat b/Wrappers/Python/demos/PDHG_examples/Tomo/phantom.mat
new file mode 100755
index 0000000..c465bbe
--- /dev/null
+++ b/Wrappers/Python/demos/PDHG_examples/Tomo/phantom.mat
diff --git a/Wrappers/Python/wip/Demos/.DS_Store b/Wrappers/Python/wip/Demos/.DS_Store
deleted file mode 100644
index 5008ddf..0000000
--- a/Wrappers/Python/wip/Demos/.DS_Store
+++ /dev/null
diff --git a/Wrappers/Python/wip/Demos/PDHG_TGV_Tomo2D.py b/Wrappers/Python/wip/Demos/PDHG_TGV_Tomo2D.py
deleted file mode 100644
index 49d4db6..0000000
--- a/Wrappers/Python/wip/Demos/PDHG_TGV_Tomo2D.py
+++ /dev/null
@@ -1,124 +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-2019 Evangelos Papoutsellis and 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 ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG
-
-from ccpi.optimisation.operators import BlockOperator, Gradient, Identity, \
-                                    SymmetrizedGradient, ZeroOperator
-from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \
-                      MixedL21Norm, BlockFunction
-
-from ccpi.astra.ops import AstraProjectorSimple
-
-# Create phantom for TV 2D tomography 
-N = 75
-
-data = np.zeros((N,N))
-
-x1 = np.linspace(0, int(N/2), N)
-x2 = np.linspace(int(N/2), 0., N)
-xv, yv = np.meshgrid(x1, x2)
-
-xv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1] = yv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1].T
-data = xv
-data = ImageData(data/data.max())
-
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
-
-detectors = N
-angles = np.linspace(0, np.pi, N, dtype=np.float32)
-
-ag = AcquisitionGeometry('parallel','2D',angles, detectors)
-Aop = AstraProjectorSimple(ig, ag, 'gpu')
-sin = Aop.direct(data)
-
-# Create noisy data. Apply Poisson noise
-scale = 0.1
-np.random.seed(5)
-n1 = scale * np.random.poisson(sin.as_array()/scale)
-noisy_data = AcquisitionData(n1, ag)
-
-
-plt.imshow(noisy_data.as_array())
-plt.show()
-#%%
-# Regularisation Parameters
-alpha = 0.7
-beta =  2
-
-# Create Operators
-op11 = Gradient(ig)
-op12 = Identity(op11.range_geometry())
-
-op22 = SymmetrizedGradient(op11.domain_geometry())    
-op21 = ZeroOperator(ig, op22.range_geometry())
-    
-op31 = Aop
-op32 = ZeroOperator(op22.domain_geometry(), ag)
-
-operator = BlockOperator(op11, -1*op12, op21, op22, op31, op32, shape=(3,2) ) 
-    
-f1 = alpha * MixedL21Norm()
-f2 = beta * MixedL21Norm() 
-f3 = KullbackLeibler(noisy_data)    
-f = BlockFunction(f1, f2, f3)         
-g = ZeroFunction()
-    
-# Compute operator Norm
-normK = operator.norm()
-
-# Primal & dual stepsizes
-sigma = 1
-tau = 1/(sigma*normK**2)
-
-
-# Setup and run the PDHG algorithm
-pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True)
-pdhg.max_iteration = 2000
-pdhg.update_objective_interval = 50
-pdhg.run(2000)
-
-plt.figure(figsize=(15,15))
-plt.subplot(3,1,1)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.subplot(3,1,2)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.subplot(3,1,3)
-plt.imshow(pdhg.get_output()[0].as_array())
-plt.title('TGV Reconstruction')
-plt.colorbar()
-plt.show()
-## 
-plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth')
-plt.plot(np.linspace(0,N,N), pdhg.get_output()[0].as_array()[int(N/2),:], label = 'TGV reconstruction')
-plt.legend()
-plt.title('Middle Line Profiles')
-plt.show()
-
-
diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py
deleted file mode 100644
index 5df02b1..0000000
--- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py
+++ /dev/null
@@ -1,225 +0,0 @@
-#========================================================================
-# Copyright 2019 Science Technology Facilities Council
-# Copyright 2019 University of Manchester
-#
-# This work is part of the Core Imaging Library developed by Science Technology
-# Facilities Council and University of Manchester
-#
-# 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.txt
-#
-# 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.
-#
-#=========================================================================
-
-
-""" 
-
-Total Variation Denoising using PDHG algorithm:
-
-             min_{x} max_{y} < K x, y > + g(x) - f^{*}(y) 
-
-
-Problem:     min_{x} \alpha * ||\nabla x||_{2,1} + \frac{1}{2} * || x - g ||_{2}^{2}
-
-             \alpha: Regularization parameter
-             
-             \nabla: Gradient operator 
-             
-             g: Noisy Data with Gaussian Noise
-                          
-             Method = 0:  K = [ \nabla,
-                                 Identity]
-                                                                    
-             Method = 1:  K = \nabla    
-             
-             
-"""
-
-from ccpi.framework import ImageData, ImageGeometry
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG
-
-from ccpi.optimisation.operators import BlockOperator, Identity, Gradient
-from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \
-                      MixedL21Norm, BlockFunction
-                      
-  
-
-
-from Data import *
-
-#%%
-
-
-data = ImageData(plt.imread('camera.png'))
-
-#
-## Create phantom for TV Gaussian denoising
-#N = 200
-#
-#data = np.zeros((N,N))
-#data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
-#data[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1
-#data = ImageData(data)
-#ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
-#ag = ig
-#
-#
-#
-## Replace with http://sipi.usc.edu/database/database.php?volume=misc&image=36#top
-
-
-
-# Create noisy data. Add Gaussian noise
-np.random.seed(10)
-noisy_data = ImageData( data.as_array() + np.random.normal(0, 0.05, size=ig.shape) )
-
-# Show Ground Truth and Noisy Data
-plt.figure(figsize=(15,15))
-plt.subplot(2,1,1)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.subplot(2,1,2)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.show()
-
-# Regularisation Parameter
-alpha = 2
-
-method = '0'
-
-if method == '0':
-
-    # Create operators
-    op1 = Gradient(ig)
-    op2 = Identity(ig, ag)
-
-    # Create BlockOperator
-    operator = BlockOperator(op1, op2, shape=(2,1) ) 
-
-    # Create functions      
-    f1 = alpha * MixedL21Norm()
-    f2 = 0.5 * L2NormSquared(b = noisy_data)    
-    f = BlockFunction(f1, f2)  
-                                      
-    g = ZeroFunction()
-    
-else:
-    
-    # Without the "Block Framework"
-    operator =  Gradient(ig)
-    f =  alpha * MixedL21Norm()
-    g =  0.5 * L2NormSquared(b = noisy_data)
-        
-# Compute Operator Norm
-normK = operator.norm()
-
-# Primal & Dual stepsizes
-sigma = 1
-tau = 1/(sigma*normK**2)
-
-# Setup and Run the PDHG algorithm
-pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
-pdhg.max_iteration = 3000
-pdhg.update_objective_interval = 200
-pdhg.run(3000, verbose=False)
-
-# Show Results
-plt.figure(figsize=(15,15))
-plt.subplot(3,1,1)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.subplot(3,1,2)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.subplot(3,1,3)
-plt.imshow(pdhg.get_output().as_array())
-plt.title('TV Reconstruction')
-plt.colorbar()
-plt.show()
-
-plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth')
-plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction')
-plt.legend()
-plt.title('Middle Line Profiles')
-plt.show()
-
-
-#%% Check with CVX solution
-
-from ccpi.optimisation.operators import SparseFiniteDiff
-
-try:
-    from cvxpy import *
-    cvx_not_installable = True
-except ImportError:
-    cvx_not_installable = False
-
-
-if cvx_not_installable:
-
-    ##Construct problem    
-    u = Variable(ig.shape)
-    
-    DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
-    DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
-    
-    # Define Total Variation as a regulariser
-    regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0))
-    fidelity = 0.5 * sum_squares(u - noisy_data.as_array())
-    
-    # choose solver
-    if 'MOSEK' in installed_solvers():
-        solver = MOSEK
-    else:
-        solver = SCS  
-        
-    obj =  Minimize( regulariser +  fidelity)
-    prob = Problem(obj)
-    result = prob.solve(verbose = True, solver = MOSEK)
-    
-    diff_cvx = numpy.abs( pdhg.get_output().as_array() - u.value )
-        
-    plt.figure(figsize=(15,15))
-    plt.subplot(3,1,1)
-    plt.imshow(pdhg.get_output().as_array())
-    plt.title('PDHG solution')
-    plt.colorbar()
-    plt.subplot(3,1,2)
-    plt.imshow(u.value)
-    plt.title('CVX solution')
-    plt.colorbar()
-    plt.subplot(3,1,3)
-    plt.imshow(diff_cvx)
-    plt.title('Difference')
-    plt.colorbar()
-    plt.show()    
-    
-    plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG')
-    plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX')
-    plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'Truth')
-    
-    plt.legend()
-    plt.title('Middle Line Profiles')
-    plt.show()
-            
-    print('Primal Objective (CVX) {} '.format(obj.value))
-    print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
-
diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py
deleted file mode 100644
index 70f6b9b..0000000
--- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py
+++ /dev/null
@@ -1,207 +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-2019 STFC,  University of Manchester
-
-#   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.
-
-""" 
-
-Total Variation Denoising using PDHG algorithm:
-
-             min_{x} max_{y} < K x, y > + g(x) - f^{*}(y) 
-
-
-Problem:     min_x, x>0  \alpha * ||\nabla x||_{1} + \int x - g * log(x)
-
-             \nabla: Gradient operator              
-             g: Noisy Data with Poisson Noise
-             \alpha: Regularization parameter
-             
-             Method = 0:  K = [ \nabla,
-                                 Identity]
-                                                                    
-             Method = 1:  K = \nabla    
-             
-             
-"""
-
-from ccpi.framework import ImageData, ImageGeometry
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG
-
-from ccpi.optimisation.operators import BlockOperator, Identity, Gradient
-from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \
-                      MixedL21Norm, BlockFunction
-
-from skimage.util import random_noise
-
-# Create phantom for TV Poisson denoising
-N = 100
-
-data = np.zeros((N,N))
-data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
-data[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1
-data = ImageData(data)
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
-ag = ig
-
-# Create noisy data. Apply Poisson noise
-n1 = random_noise(data.as_array(), mode = 'poisson', seed = 10)
-noisy_data = ImageData(n1)
-
-# Regularisation Parameter
-alpha = 2
-
-method = '1'
-
-if method == '0':
-
-    # Create operators
-    op1 = Gradient(ig)
-    op2 = Identity(ig, ag)
-
-    # Create BlockOperator
-    operator = BlockOperator(op1, op2, shape=(2,1) ) 
-
-    # Create functions
-      
-    f1 = alpha * MixedL21Norm()
-    f2 = KullbackLeibler(noisy_data)    
-    f = BlockFunction(f1, f2)  
-                                      
-    g = ZeroFunction()
-    
-else:
-    
-    # Without the "Block Framework"
-    operator = Gradient(ig)
-    f =  alpha * MixedL21Norm()
-    g =  KullbackLeibler(noisy_data)   
-        
-    
-# Compute operator Norm
-normK = operator.norm()
-
-# Primal & dual stepsizes
-sigma = 1
-tau = 1/(sigma*normK**2)
-opt = {'niter':2000, 'memopt': True}
-
-# Setup and run the PDHG algorithm
-pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True)
-pdhg.max_iteration = 2000
-pdhg.update_objective_interval = 50
-
-def pdgap_objectives(niter, objective, solution):
-    
-
-     print( "{:04}/{:04} {:<5} {:.4f} {:<5} {:.4f} {:<5} {:.4f}".\
-                      format(niter, pdhg.max_iteration,'', \
-                             objective[0],'',\
-                             objective[1],'',\
-                             objective[2]))
-
-pdhg.run(2000, callback = pdgap_objectives)
-
-
-plt.figure(figsize=(15,15))
-plt.subplot(3,1,1)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.subplot(3,1,2)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.subplot(3,1,3)
-plt.imshow(pdhg.get_output().as_array())
-plt.title('TV Reconstruction')
-plt.colorbar()
-plt.show()
-## 
-plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth')
-plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction')
-plt.legend()
-plt.title('Middle Line Profiles')
-plt.show()
-
-
-#%% Check with CVX solution
-
-from ccpi.optimisation.operators import SparseFiniteDiff
-
-try:
-    from cvxpy import *
-    cvx_not_installable = True
-except ImportError:
-    cvx_not_installable = False
-
-
-if cvx_not_installable:
-
-    ##Construct problem    
-    u1 = Variable(ig.shape)
-    q = Variable()
-    
-    DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
-    DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
-    
-    # Define Total Variation as a regulariser
-    regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u1), DY.matrix() * vec(u1)]), 2, axis = 0))
-    
-    fidelity = sum( u1 - multiply(noisy_data.as_array(), log(u1)) )
-    constraints = [q>= fidelity, u1>=0]
-        
-    solver = ECOS
-    obj =  Minimize( regulariser +  q)
-    prob = Problem(obj, constraints)
-    result = prob.solve(verbose = True, solver = solver)
-
-    
-    diff_cvx = numpy.abs( pdhg.get_output().as_array() - u1.value )
-        
-    plt.figure(figsize=(15,15))
-    plt.subplot(3,1,1)
-    plt.imshow(pdhg.get_output().as_array())
-    plt.title('PDHG solution')
-    plt.colorbar()
-    plt.subplot(3,1,2)
-    plt.imshow(u1.value)
-    plt.title('CVX solution')
-    plt.colorbar()
-    plt.subplot(3,1,3)
-    plt.imshow(diff_cvx)
-    plt.title('Difference')
-    plt.colorbar()
-    plt.show()    
-    
-    plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG')
-    plt.plot(np.linspace(0,N,N), u1.value[int(N/2),:], label = 'CVX')
-    plt.legend()
-    plt.title('Middle Line Profiles')
-    plt.show()
-            
-    print('Primal Objective (CVX) {} '.format(obj.value))
-    print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
-
-
-
-
-
diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py
deleted file mode 100644
index f5d4ce4..0000000
--- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py
+++ /dev/null
@@ -1,198 +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-2019 Evangelos Papoutsellis and 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.
-
-""" 
-
-Total Variation Denoising using PDHG algorithm:
-
-             min_{x} max_{y} < K x, y > + g(x) - f^{*}(y) 
-
-
-Problem:     min_x, x>0  \alpha * ||\nabla x||_{1} + ||x-g||_{1}
-
-             \nabla: Gradient operator 
-             g: Noisy Data with Salt & Pepper Noise
-             \alpha: Regularization parameter
-             
-             Method = 0:  K = [ \nabla,
-                                 Identity]
-                                                                    
-             Method = 1:  K = \nabla    
-             
-             
-"""
-
-from ccpi.framework import ImageData, ImageGeometry
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG
-
-from ccpi.optimisation.operators import BlockOperator, Identity, Gradient
-from ccpi.optimisation.functions import ZeroFunction, L1Norm, \
-                      MixedL21Norm, BlockFunction
-
-from skimage.util import random_noise
-
-# Create phantom for TV Salt & Pepper denoising
-N = 100
-
-data = np.zeros((N,N))
-data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
-data[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1
-data = ImageData(data)
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
-ag = ig
-
-# Create noisy data. Apply Salt & Pepper noise
-n1 = random_noise(data.as_array(), mode = 's&p', salt_vs_pepper = 0.9, amount=0.2)
-noisy_data = ImageData(n1)
-
-# Regularisation Parameter
-alpha = 2
-
-method = '0'
-
-if method == '0':
-
-    # Create operators
-    op1 = Gradient(ig)
-    op2 = Identity(ig, ag)
-
-    # Create BlockOperator
-    operator = BlockOperator(op1, op2, shape=(2,1) ) 
-
-    # Create functions
-      
-    f1 = alpha * MixedL21Norm()
-    f2 = L1Norm(b = noisy_data)    
-    f = BlockFunction(f1, f2)  
-                                      
-    g = ZeroFunction()
-    
-else:
-    
-    # Without the "Block Framework"
-    operator = Gradient(ig)
-    f =  alpha * MixedL21Norm()
-    g =  L1Norm(b = noisy_data)
-        
-    
-# Compute operator Norm
-normK = operator.norm()
-
-# Primal & dual stepsizes
-sigma = 1
-tau = 1/(sigma*normK**2)
-opt = {'niter':2000, 'memopt': True}
-
-# Setup and run the PDHG algorithm
-pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True)
-pdhg.max_iteration = 2000
-pdhg.update_objective_interval = 50
-pdhg.run(2000)
-
-
-plt.figure(figsize=(15,15))
-plt.subplot(3,1,1)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.subplot(3,1,2)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.subplot(3,1,3)
-plt.imshow(pdhg.get_output().as_array())
-plt.title('TV Reconstruction')
-plt.colorbar()
-plt.show()
-## 
-plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth')
-plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction')
-plt.legend()
-plt.title('Middle Line Profiles')
-plt.show()
-
-
-##%% Check with CVX solution
-
-from ccpi.optimisation.operators import SparseFiniteDiff
-
-try:
-    from cvxpy import *
-    cvx_not_installable = True
-except ImportError:
-    cvx_not_installable = False
-
-
-if cvx_not_installable:
-
-    ##Construct problem    
-    u = Variable(ig.shape)
-    
-    DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
-    DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
-    
-    # Define Total Variation as a regulariser
-    regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0))
-    fidelity = pnorm( u - noisy_data.as_array(),1)
-    
-    # choose solver
-    if 'MOSEK' in installed_solvers():
-        solver = MOSEK
-    else:
-        solver = SCS  
-        
-    obj =  Minimize( regulariser +  fidelity)
-    prob = Problem(obj)
-    result = prob.solve(verbose = True, solver = solver)
-    
-    diff_cvx = numpy.abs( pdhg.get_output().as_array() - u.value )
-        
-    plt.figure(figsize=(15,15))
-    plt.subplot(3,1,1)
-    plt.imshow(pdhg.get_output().as_array())
-    plt.title('PDHG solution')
-    plt.colorbar()
-    plt.subplot(3,1,2)
-    plt.imshow(u.value)
-    plt.title('CVX solution')
-    plt.colorbar()
-    plt.subplot(3,1,3)
-    plt.imshow(diff_cvx)
-    plt.title('Difference')
-    plt.colorbar()
-    plt.show()    
-    
-    plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG')
-    plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX')
-    plt.legend()
-    plt.title('Middle Line Profiles')
-    plt.show()
-            
-    print('Primal Objective (CVX) {} '.format(obj.value))
-    print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
-
-
-
-
-
diff --git a/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Denoising.py b/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Denoising.py
deleted file mode 100644
index 041d4ee..0000000
--- a/Wrappers/Python/wip/Demos/PDHG_Tikhonov_Denoising.py
+++ /dev/null
@@ -1,176 +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-2019 Evangelos Papoutsellis and 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 ccpi.framework import ImageData, ImageGeometry
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG
-
-from ccpi.optimisation.operators import BlockOperator, Identity, Gradient
-from ccpi.optimisation.functions import ZeroFunction, L2NormSquared,  BlockFunction
-
-from skimage.util import random_noise
-
-# Create phantom for TV Salt & Pepper denoising
-N = 100
-
-data = np.zeros((N,N))
-data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
-data[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1
-data = ImageData(data)
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
-ag = ig
-
-# Create noisy data. Apply Salt & Pepper noise
-n1 = random_noise(data.as_array(), mode = 'gaussian', mean=0, var = 0.05, seed=10)
-noisy_data = ImageData(n1)
-
-# Regularisation Parameter
-alpha = 4
-
-method = '0'
-
-if method == '0':
-
-    # Create operators
-    op1 = Gradient(ig)
-    op2 = Identity(ig, ag)
-
-    # Create BlockOperator
-    operator = BlockOperator(op1, op2, shape=(2,1) ) 
-
-    # Create functions
-      
-    f1 = alpha * L2NormSquared()
-    f2 = 0.5 * L2NormSquared(b = noisy_data)    
-    f = BlockFunction(f1, f2)                                       
-    g = ZeroFunction()
-    
-else:
-    
-    # Without the "Block Framework"
-    operator = Gradient(ig)
-    f =  alpha * L2NormSquared()
-    g =  0.5 * L2NormSquared(b = noisy_data)
-        
-    
-# Compute operator Norm
-normK = operator.norm()
-
-# Primal & dual stepsizes
-sigma = 1
-tau = 1/(sigma*normK**2)
-opt = {'niter':2000, 'memopt': True}
-
-# Setup and run the PDHG algorithm
-pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True)
-pdhg.max_iteration = 2000
-pdhg.update_objective_interval = 50
-pdhg.run(2000)
-
-
-plt.figure(figsize=(15,15))
-plt.subplot(3,1,1)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.subplot(3,1,2)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.subplot(3,1,3)
-plt.imshow(pdhg.get_output().as_array())
-plt.title('Tikhonov Reconstruction')
-plt.colorbar()
-plt.show()
-## 
-plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth')
-plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'Tikhonov reconstruction')
-plt.legend()
-plt.title('Middle Line Profiles')
-plt.show()
-
-
-##%% Check with CVX solution
-
-from ccpi.optimisation.operators import SparseFiniteDiff
-
-try:
-    from cvxpy import *
-    cvx_not_installable = True
-except ImportError:
-    cvx_not_installable = False
-
-
-if cvx_not_installable:
-
-    ##Construct problem    
-    u = Variable(ig.shape)
-    
-    DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
-    DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
-    
-    # Define Total Variation as a regulariser
-    
-    regulariser = alpha * sum_squares(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0))
-    fidelity = 0.5 * sum_squares(u - noisy_data.as_array())    
-    
-    # choose solver
-    if 'MOSEK' in installed_solvers():
-        solver = MOSEK
-    else:
-        solver = SCS  
-        
-    obj =  Minimize( regulariser +  fidelity)
-    prob = Problem(obj)
-    result = prob.solve(verbose = True, solver = solver)
-    
-    diff_cvx = numpy.abs( pdhg.get_output().as_array() - u.value )
-        
-    plt.figure(figsize=(15,15))
-    plt.subplot(3,1,1)
-    plt.imshow(pdhg.get_output().as_array())
-    plt.title('PDHG solution')
-    plt.colorbar()
-    plt.subplot(3,1,2)
-    plt.imshow(u.value)
-    plt.title('CVX solution')
-    plt.colorbar()
-    plt.subplot(3,1,3)
-    plt.imshow(diff_cvx)
-    plt.title('Difference')
-    plt.colorbar()
-    plt.show()    
-    
-    plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG')
-    plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX')
-    plt.legend()
-    plt.title('Middle Line Profiles')
-    plt.show()
-            
-    print('Primal Objective (CVX) {} '.format(obj.value))
-    print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
-
-
-
-
-
diff --git a/Wrappers/Python/wip/Demos/fista_test.py b/Wrappers/Python/wip/Demos/fista_test.py
deleted file mode 100644
index dd1f6fa..0000000
--- a/Wrappers/Python/wip/Demos/fista_test.py
+++ /dev/null
@@ -1,127 +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-2019 Evangelos Papoutsellis and 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 ccpi.framework import ImageData, ImageGeometry
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import FISTA, PDHG
-
-from ccpi.optimisation.operators import BlockOperator, Gradient, Identity
-from ccpi.optimisation.functions import L2NormSquared, L1Norm, \
-                                        MixedL21Norm, FunctionOperatorComposition, BlockFunction, ZeroFunction
-                                                
-from skimage.util import random_noise
-
-# Create phantom for TV Gaussian denoising
-N = 100
-
-data = np.zeros((N,N))
-data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
-data[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1
-data = ImageData(data)
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
-ag = ig
-
-# Create noisy data. Add Gaussian noise
-n1 = random_noise(data.as_array(), mode = 's&p', salt_vs_pepper = 0.9, amount=0.2)
-noisy_data = ImageData(n1)
-
-# Regularisation Parameter
-alpha = 5
-
-operator = Gradient(ig)
-
-#fidelity = L1Norm(b=noisy_data)
-#regulariser = FunctionOperatorComposition(alpha * L2NormSquared(), operator)
-
-fidelity = FunctionOperatorComposition(alpha * MixedL21Norm(), operator)
-regulariser = 0.5 * L2NormSquared(b = noisy_data)
-
-x_init = ig.allocate()
-
-## Setup and run the PDHG algorithm
-opt = {'tol': 1e-4, 'memopt':True}
-fista = FISTA(x_init=x_init , f=regulariser, g=fidelity, opt=opt)
-fista.max_iteration = 2000
-fista.update_objective_interval = 50
-fista.run(2000, verbose=True)
-
-plt.figure(figsize=(15,15))
-plt.subplot(3,1,1)
-plt.imshow(data.as_array())
-plt.title('Ground Truth')
-plt.colorbar()
-plt.subplot(3,1,2)
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy Data')
-plt.colorbar()
-plt.subplot(3,1,3)
-plt.imshow(fista.get_output().as_array())
-plt.title('TV Reconstruction')
-plt.colorbar()
-plt.show()
-
-# Compare with PDHG
-method = '0'
-#
-if method == '0':
-#
-#    # Create operators
-    op1 = Gradient(ig)
-    op2 = Identity(ig, ag)
-#
-#    # Create BlockOperator
-    operator = BlockOperator(op1, op2, shape=(2,1) )   
-    f = BlockFunction(alpha * L2NormSquared(), fidelity)                                        
-    g = ZeroFunction()
-  
-## Compute operator Norm
-normK = operator.norm()
-#
-## Primal & dual stepsizes
-sigma = 1
-tau = 1/(sigma*normK**2)
-#
-#
-## Setup and run the PDHG algorithm
-pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True)
-pdhg.max_iteration = 2000
-pdhg.update_objective_interval = 50
-pdhg.run(2000)
-#
-#%%
-plt.figure(figsize=(15,15))
-plt.subplot(3,1,1)
-plt.imshow(fista.get_output().as_array())
-plt.title('FISTA')
-plt.colorbar()
-plt.subplot(3,1,2)
-plt.imshow(pdhg.get_output().as_array())
-plt.title('PDHG')
-plt.colorbar()
-plt.subplot(3,1,3)
-plt.imshow(np.abs(pdhg.get_output().as_array()-fista.get_output().as_array()))
-plt.title('Diff FISTA-PDHG')
-plt.colorbar()
-plt.show()
-
-
diff --git a/Wrappers/Python/wip/demo_colourbay.py b/Wrappers/Python/wip/demo_colourbay.py
index 5dbf2e1..0536b07 100644
--- a/Wrappers/Python/wip/demo_colourbay.py
+++ b/Wrappers/Python/wip/demo_colourbay.py
@@ -18,7 +18,7 @@ from ccpi.optimisation.funcs import Norm2sq, Norm1
 # Permute (numpy.transpose) puts into our default ordering which is 
 # (channel, angle, vertical, horizontal).
 
-pathname = '/media/jakob/050d8d45-fab3-4285-935f-260e6c5f162c1/Data/ColourBay/spectral_data_sets/CarbonPd/'
+pathname = '/media/newhd/shared/Data/ColourBay/spectral_data_sets/CarbonPd/'
 filename = 'carbonPd_full_sinogram_stripes_removed.mat'
 
 X = loadmat(pathname + filename)