diff options
| -rw-r--r-- | Wrappers/Python/ccpi/reconstruction/ops.py | 102 |
1 files changed, 30 insertions, 72 deletions
diff --git a/Wrappers/Python/ccpi/reconstruction/ops.py b/Wrappers/Python/ccpi/reconstruction/ops.py index 00a4a3c..af8ead6 100644 --- a/Wrappers/Python/ccpi/reconstruction/ops.py +++ b/Wrappers/Python/ccpi/reconstruction/ops.py @@ -1,6 +1,23 @@ +# -*- 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 -import numpy as np -import astra +# 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 from scipy.sparse.linalg import svds from framework import DataSet, VolumeData, SinogramData, DataSetProcessor @@ -8,7 +25,7 @@ from framework import DataSet, VolumeData, SinogramData, DataSetProcessor # to not just use generic DataSet -class Operator: +class Operator(object): def direct(self,x): return x def adjoint(self,x): @@ -33,18 +50,20 @@ class ForwardBackProjector(Operator): def __init__(self): # do nothing i = 1 + super(ForwardBackProjector, self).__init__() class LinearOperatorMatrix(Operator): def __init__(self,A): self.A = A self.s1 = None # Largest singular value, initially unknown + super(LinearOperatorMatrix, self).__init__() def direct(self,x): - return DataSet(np.dot(self.A,x.as_array())) + return DataSet(numpy.dot(self.A,x.as_array())) def adjoint(self,x): - return DataSet(np.dot(self.A.transpose(),x.as_array())) + return DataSet(numpy.dot(self.A.transpose(),x.as_array())) def size(self): return self.A.shape @@ -60,6 +79,7 @@ class LinearOperatorMatrix(Operator): class Identity(Operator): def __init__(self): self.s1 = 1.0 + super(Identity, self).__init__() def direct(self,x): return x @@ -73,81 +93,19 @@ class Identity(Operator): 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) + x0 = DataSet(numpy.random.randn(inputsize[0],inputsize[1])) + s = numpy.zeros(numiters) # Loop - for it in np.arange(numiters): + for it in numpy.arange(numiters): x1 = op.adjoint(op.direct(x0)) - x1norm = np.sqrt((x1**2).sum()) + x1norm = numpy.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 + return numpy.sqrt(s[-1]), numpy.sqrt(s), x0 #def PowerMethod(op,numiters): # # Initialise random |