moved modular FISTA here

4 files changed, 507 insertions, 0 deletions

diff --git a/Wrappers/Python/ccpi/reconstruction/algs.py b/Wrappers/Python/ccpi/reconstruction/algs.py
new file mode 100644
index 0000000..bfe9570
--- /dev/null
+++ b/Wrappers/Python/ccpi/reconstruction/algs.py
@@ -0,0 +1,128 @@
+# -*- 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 Jakob Jorgensen, Daniil Kazantsev 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.
+
+import numpy
+import time
+
+from funcs import BaseFunction
+
+def FISTA(x_init, f=None, g=None, opt=None):
+
+    # default inputs
+    if f   is None: f = BaseFunction()
+    if g   is None: g = BaseFunction()
+    if opt is None: opt = {'tol': 1e-4, 'iter': 1000}
+
+    # algorithmic parameters
+    tol      = opt['tol']
+    max_iter = opt['iter']
+    
+    # initialization
+    x_old = x_init
+    y = x_init;
+    
+    timing = numpy.zeros(max_iter)
+    criter = numpy.zeros(max_iter)
+    
+    invL = 1/f.L
+    
+    t_old = 1
+
+    # algorithm loop
+    for it in range(0, max_iter):
+    
+        time0 = time.time()
+        
+        u = y - invL*f.grad(y)
+        
+        x = g.prox(u,invL)
+        
+        t = 0.5*(1 + numpy.sqrt(1 + 4*(t_old**2)))
+        
+        y = x + (t_old-1)/t*(x-x_old)
+        
+        x_old = x
+        t_old = t
+        
+        # time and criterion
+        timing[it] = time.time() - time0
+        criter[it] = f.fun(x) + g.fun(x);
+        
+        # stopping rule
+        #if np.linalg.norm(x - x_old) < tol * np.linalg.norm(x_old) and it > 10:
+        #   break
+    
+        #print(it, 'out of', 10, 'iterations', end='\r');
+
+    criter = criter[0:it+1];
+    timing = numpy.cumsum(timing[0:it+1]);
+    
+    return x, it, timing, criter
+
+def FBPD(x_init, f=None, g=None, h=None, opt=None):
+
+    # default inputs
+    if f   is None: f = BaseFunction()
+    if g   is None: g = BaseFunction()
+    if h   is None: h = BaseFunction()
+    if opt is None: opt = {'tol': 1e-4, 'iter': 1000}
+
+    # algorithmic parameters
+    tol      = opt['tol']
+    max_iter = opt['iter']
+    
+    # step-sizes
+    tau   = 2 / (g.L + 2);
+    sigma = (1/tau - g.L/2) / h.L;
+    
+    inv_sigma = 1/sigma
+
+    # initialization
+    x = x_init
+    y = h.dir_op(x);
+    
+    timing = numpy.zeros(max_iter)
+    criter = numpy.zeros(max_iter)
+
+    # algorithm loop
+    for it in range(0, max_iter):
+    
+        t = time.time()
+    
+        # primal forward-backward step
+        x_old = x;
+        x = x - tau * ( g.grad(x) + h.adj_op(y) );
+        x = f.prox(x, tau);
+    
+        # dual forward-backward step
+        y = y + sigma * h.dir_op(2*x - x_old);
+        y = y - sigma * h.prox(inv_sigma*y, inv_sigma);   
+
+        # time and criterion
+        timing[it] = time.time() - t
+        criter[it] = f.fun(x) + g.fun(x) + h.fun(h.dir_op(x));
+           
+        # stopping rule
+        #if np.linalg.norm(x - x_old) < tol * np.linalg.norm(x_old) and it > 10:
+        #   break
+
+    criter = criter[0:it+1];
+    timing = numpy.cumsum(timing[0:it+1]);
+    
+    return x, it, timing, criter
diff --git a/Wrappers/Python/ccpi/reconstruction/astra_ops.py b/Wrappers/Python/ccpi/reconstruction/astra_ops.py
new file mode 100644
index 0000000..0c95b73
--- /dev/null
+++ b/Wrappers/Python/ccpi/reconstruction/astra_ops.py
@@ -0,0 +1,99 @@
+# -*- coding: utf-8 -*-
+#    This work is independent part of the Core Imaging Library developed by
+#    Visual Analytics and Imaging System Group of the Science Technology
+#    Facilities Council, STFC
+#    This program is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU General Public License as published by
+#    the Free Software Foundation, either version 3 of the License, or
+#    (at your option) any later version.
+#
+#    This program is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#    GNU General Public License for more details.
+#
+#    You should have received a copy of the GNU General Public License
+#    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+from ccpi.reconstruction.ops import Operator
+import astra
+import numpy
+from ccpi.framework import SinogramData, VolumeData
+from ccpi.reconstruction.ops import PowerMethodNonsquare
+
+class AstraProjector(Operator):
+    """A simple 2D/3D parallel/fan beam projection/backprojection class based 
+    on ASTRA toolbox"""
+    def __init__(self, DetWidth, DetectorsDim, SourceOrig, OrigDetec, 
+                 AnglesVec, ObjSize, projtype, device):
+        super(AstraProjector, self).__init__()
+        self.DetectorsDim = DetectorsDim
+        self.AnglesVec = AnglesVec
+        self.ProjNumb = len(AnglesVec)
+        self.ObjSize = ObjSize
+        if projtype == 'parallel':
+            self.proj_geom = astra.create_proj_geom('parallel', DetWidth, 
+                                                    DetectorsDim, AnglesVec)
+        elif projtype == 'fanbeam':
+            self.proj_geom = astra.create_proj_geom('fanflat', DetWidth,
+                                                    DetectorsDim, AnglesVec, 
+                                                    SourceOrig, OrigDetec)
+        else:
+            print ("Please select for projtype between 'parallel' and 'fanbeam'")
+        self.vol_geom = astra.create_vol_geom(ObjSize, ObjSize)
+        if device == 'cpu':
+            self.proj_id = astra.create_projector('line', self.proj_geom, 
+                                                  self.vol_geom) # for CPU
+            self.device = 1
+        elif device == 'gpu':
+            self.proj_id = astra.create_projector('cuda', self.proj_geom, 
+                                                  self.vol_geom) # for GPU
+            self.device = 0
+        else:
+            print ("Select between 'cpu' or 'gpu' for device")
+        self.s1 = None
+    def direct(self, IM):
+        """Applying forward projection to IM [2D or 3D array]"""
+        if numpy.ndim(IM.as_array()) == 3:
+            slices = numpy.size(IM.as_array(),numpy.ndim(IM.as_array())-1)
+            DATA = numpy.zeros((self.ProjNumb,self.DetectorsDim,slices), 
+                               'float32')
+            for i in range(0,slices):
+                sinogram_id, DATA[:,:,i] = \
+                    astra.create_sino(IM[:,:,i].as_array(), self.proj_id)
+                astra.data2d.delete(sinogram_id)
+            astra.data2d.delete(self.proj_id)
+        else:
+            sinogram_id, DATA = astra.create_sino(IM.as_array(), self.proj_id)
+            astra.data2d.delete(sinogram_id)
+            astra.data2d.delete(self.proj_id)
+        return SinogramData(DATA)
+    def adjoint(self, DATA):
+        """Applying backprojection to DATA [2D or 3D]"""
+        if numpy.ndim(DATA) == 3:
+           slices = numpy.size(DATA.as_array(),numpy.ndim(DATA.as_array())-1)
+           IM = numpy.zeros((self.ObjSize,self.ObjSize,slices), 'float32')
+           for i in range(0,slices):
+               rec_id, IM[:,:,i] = \
+                   astra.create_backprojection(DATA[:,:,i].as_array(), 
+                                               self.proj_id)
+               astra.data2d.delete(rec_id)
+           astra.data2d.delete(self.proj_id)
+        else:
+            rec_id, IM = astra.create_backprojection(DATA.as_array(), 
+                                                     self.proj_id)        
+            astra.data2d.delete(rec_id)
+            astra.data2d.delete(self.proj_id)
+        return VolumeData(IM)
+    
+    def delete(self):
+        astra.data2d.delete(self.proj_id)
+    
+    def get_max_sing_val(self):
+        self.s1, sall, svec = PowerMethodNonsquare(self,10)
+        return self.s1
+    
+    def size(self):
+        return ( (self.AnglesVec.size, self.DetectorsDim), \
+                 (self.ObjSize, self.ObjSize) )
+
diff --git a/Wrappers/Python/ccpi/reconstruction/funcs.py b/Wrappers/Python/ccpi/reconstruction/funcs.py
new file mode 100644
index 0000000..c6f29d2
--- /dev/null
+++ b/Wrappers/Python/ccpi/reconstruction/funcs.py
@@ -0,0 +1,118 @@
+# -*- 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 Jakob Jorgensen, Daniil Kazantsev 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.
+
+
+class BaseFunction(object):
+    def __init__(self):   pass
+    def fun(self,x):      return 0
+    def grad(self,x):     return 0
+    def prox(self,x,tau): return x
+    def dir_op(self,x):   return x
+    def adj_op(self,x):   return x
+
+# Define a class for 2-norm
+class Norm2sq(BaseFunction):
+    '''
+    f(x) = c*||A*x-b||_2^2
+    
+    which has 
+    
+    grad[f](x) = 2*c*A^T*(A*x-b)
+    
+    and Lipschitz constant
+    
+    L = 2*c*||A||_2^2 = 2*s1(A)^2
+    
+    where s1(A) is the largest singular value of A.
+    
+    '''
+    
+    def __init__(self,A,b,c=1.0):
+        self.A = A  # Should be an operator, default identity
+        self.b = b  # Default zero DataSet?
+        self.c = c  # Default 1.
+        
+        # Compute the Lipschitz parameter from the operator.
+        # Initialise to None instead and only call when needed.
+        self.L = 2.0*self.c*(self.A.get_max_sing_val()**2)
+        super(Norm2sq, self).__init__()
+    
+    def grad(self,x):
+        #return 2*self.c*self.A.adjoint( self.A.direct(x) - self.b )
+        return 2.0*self.c*self.A.adjoint( self.A.direct(x) - self.b )
+    
+    def fun(self,x):
+        #return self.c* np.sum(np.square((self.A.direct(x) - self.b).ravel()))
+        return self.c*( ( (self.A.direct(x)-self.b)**2).sum() )
+
+## Define a class to represent a least squares data fidelity
+#class LeastSquares(BaseFunction):
+#    
+#    b     = None
+#    A     = None
+#    L     = None
+#    
+#    def __init__(self, A, b):
+#        self.A    = A
+#        #b.shape = (b.shape[0],1)
+#        self.b    = b
+#        
+#        # Compute the Lipschitz parameter from the operator
+#        # Initialise to None instead and only call when needed.
+#        self.L = self.A.get_max_sing_val()**2
+#    
+#    def grad(self, x):
+#        #return np.dot(self.A.transpose(),  (np.dot(self.A,x) - self.b)  )
+#        return self.A.adjoint( self.A.direct(x) - self.b )
+#    
+#    def fun(self, x):
+#        # p = np.dot(self.A, x)
+#        return 0.5*( ( (self.A.direct(x)-self.b)**2).sum() )
+
+# Define a class to represent the zero-function, to test pure least squares 
+# minimization using FISTA
+class ZeroFun(BaseFunction):
+    
+    def __init__(self,gamma=0,L=1):
+        self.gamma = gamma
+        self.L = L
+        super(ZeroFun, self).__init__()
+    
+    def fun(self,x):
+        return 0
+    
+    def prox(self,x,tau):
+        return x
+
+# A more interesting example, least squares plus 1-norm minimization.
+# Define class to represent 1-norm including prox function
+class Norm1(BaseFunction):
+    
+    def __init__(self,gamma):
+        # Do nothing
+        self.gamma = gamma
+        self.L = 1
+        super(Norm1, self).__init__()
+    
+    def fun(self,x):
+        return self.gamma*(x.abs().sum())
+    
+    def prox(self,x,tau):
+        return (x.abs() - tau*self.gamma).maximum(0) * x.sign()
+    
diff --git a/Wrappers/Python/ccpi/reconstruction/ops.py b/Wrappers/Python/ccpi/reconstruction/ops.py
new file mode 100644
index 0000000..00a4a3c
--- /dev/null
+++ b/Wrappers/Python/ccpi/reconstruction/ops.py
@@ -0,0 +1,162 @@
+
+import numpy as np
+import astra
+from scipy.sparse.linalg import svds
+from framework import DataSet, VolumeData, SinogramData, DataSetProcessor
+
+# Maybe operators need to know what types they take as inputs/outputs
+# to not just use generic DataSet
+
+
+class Operator:
+    def direct(self,x):
+        return x
+    def adjoint(self,x):
+        return x
+    def size(self):
+        # To be defined for specific class
+        return None
+
+# Or should we rather have an attribute isLinear instead of separate class?
+
+#class OperatorLinear(Operator):
+#    
+#    def __init__():
+
+class ForwardBackProjector(Operator):
+    
+    # The constructor should set up everything, ie at least hold equivalent of 
+    # projection geometry and volume geometry, so that when calling direct and 
+    # adjoint methods, only the volume/sinogram is needed as input. Quite 
+    # similar to opTomo operator.
+    
+    def __init__(self):
+        # do nothing
+        i  = 1
+    
+
+class LinearOperatorMatrix(Operator):
+    def __init__(self,A):
+        self.A = A
+        self.s1 = None   # Largest singular value, initially unknown
+        
+    def direct(self,x):
+        return DataSet(np.dot(self.A,x.as_array()))
+    
+    def adjoint(self,x):
+        return DataSet(np.dot(self.A.transpose(),x.as_array()))
+    
+    def size(self):
+        return self.A.shape
+    
+    def get_max_sing_val(self):
+        # If unknown, compute and store. If known, simply return it.
+        if self.s1 is None:
+            self.s1 = svds(self.A,1,return_singular_vectors=False)[0]
+            return self.s1
+        else:
+            return self.s1
+
+class Identity(Operator):
+    def __init__(self):
+        self.s1 = 1.0
+        
+    def direct(self,x):
+        return x
+    
+    def adjoint(self,x):
+        return x
+    
+    def size(self):
+        return NotImplemented
+    
+    def get_max_sing_val(self):
+        return self.s1
+
+class AstraProjector:
+    """A simple 2D/3D parallel/fan beam projection/backprojection class based on ASTRA toolbox"""
+    def __init__(self, DetWidth, DetectorsDim, SourceOrig, OrigDetec, AnglesVec, ObjSize, projtype, device):
+        self.DetectorsDim = DetectorsDim
+        self.AnglesVec = AnglesVec
+        self.ProjNumb = len(AnglesVec)
+        self.ObjSize = ObjSize
+        if projtype == 'parallel':
+            self.proj_geom = astra.create_proj_geom('parallel', DetWidth, DetectorsDim, AnglesVec)
+        elif projtype == 'fanbeam':
+            self.proj_geom = astra.create_proj_geom('fanflat', DetWidth, DetectorsDim, AnglesVec, SourceOrig, OrigDetec)
+        else:
+            print ("Please select for projtype between 'parallel' and 'fanbeam'")
+        self.vol_geom = astra.create_vol_geom(ObjSize, ObjSize)
+        if device == 'cpu':
+            self.proj_id = astra.create_projector('line', self.proj_geom, self.vol_geom) # for CPU
+            self.device = 1
+        elif device == 'gpu':
+            self.proj_id = astra.create_projector('cuda', self.proj_geom, self.vol_geom) # for GPU
+            self.device = 0
+        else:
+            print ("Select between 'cpu' or 'gpu' for device")
+        self.s1 = None
+    def direct(self, IM):
+        """Applying forward projection to IM [2D or 3D array]"""
+        if np.ndim(IM.as_array()) == 3:
+            slices = np.size(IM.as_array(),np.ndim(IM.as_array())-1)
+            DATA = np.zeros((self.ProjNumb,self.DetectorsDim,slices), 'float32')
+            for i in range(0,slices):
+                sinogram_id, DATA[:,:,i] = astra.create_sino(IM[:,:,i].as_array(), self.proj_id)
+                astra.data2d.delete(sinogram_id)
+            astra.data2d.delete(self.proj_id)
+        else:
+            sinogram_id, DATA = astra.create_sino(IM.as_array(), self.proj_id)
+            astra.data2d.delete(sinogram_id)
+            astra.data2d.delete(self.proj_id)
+        return SinogramData(DATA)
+    def adjoint(self, DATA):
+        """Applying backprojection to DATA [2D or 3D]"""
+        if np.ndim(DATA) == 3:
+           slices = np.size(DATA.as_array(),np.ndim(DATA.as_array())-1)
+           IM = np.zeros((self.ObjSize,self.ObjSize,slices), 'float32')
+           for i in range(0,slices):
+               rec_id, IM[:,:,i] = astra.create_backprojection(DATA[:,:,i].as_array(), self.proj_id)
+               astra.data2d.delete(rec_id)
+           astra.data2d.delete(self.proj_id)
+        else:
+            rec_id, IM = astra.create_backprojection(DATA.as_array(), self.proj_id)        
+            astra.data2d.delete(rec_id)
+            astra.data2d.delete(self.proj_id)
+        return VolumeData(IM)
+    
+    def delete(self):
+        astra.data2d.delete(self.proj_id)
+    
+    def get_max_sing_val(self):
+        self.s1, sall, svec = PowerMethodNonsquare(self,10)
+        return self.s1
+    
+    def size(self):
+        return ( (self.AnglesVec.size, self.DetectorsDim), \
+                 (self.ObjSize, self.ObjSize) )
+
+def PowerMethodNonsquare(op,numiters):
+    # Initialise random
+    inputsize = op.size()[1]
+    x0 = DataSet(np.random.randn(inputsize[0],inputsize[1]))
+    s = np.zeros(numiters)
+    # Loop
+    for it in np.arange(numiters):
+        x1 = op.adjoint(op.direct(x0))
+        x1norm = np.sqrt((x1**2).sum())
+        s[it] = (x1*x0).sum() / (x0*x0).sum()
+        x0 = (1.0/x1norm)*x1
+    return np.sqrt(s[-1]), np.sqrt(s), x0
+
+#def PowerMethod(op,numiters):
+#    # Initialise random
+#    x0 = np.random.randn(400)
+#    s = np.zeros(numiters)
+#    # Loop
+#    for it in np.arange(numiters):
+#        x1 = np.dot(op.transpose(),np.dot(op,x0))
+#        x1norm = np.sqrt(np.sum(np.dot(x1,x1)))
+#        s[it] = np.dot(x1,x0) / np.dot(x1,x0)
+#        x0 = (1.0/x1norm)*x1
+#    return s, x0
\ No newline at end of file