diff options
| -rw-r--r-- | Wrappers/Python/ccpi/reconstruction/algs.py | 128 | ||||
| -rw-r--r-- | Wrappers/Python/ccpi/reconstruction/astra_ops.py | 99 | ||||
| -rw-r--r-- | Wrappers/Python/ccpi/reconstruction/funcs.py | 118 | ||||
| -rw-r--r-- | Wrappers/Python/ccpi/reconstruction/ops.py | 162 |
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 |