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