From 00781c2c2b1f51bed23db8878f628d72b296ac0c Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 6 Mar 2018 15:25:02 +0000 Subject: moved work-in-progress files to wip --- Wrappers/Python/test/DemoRecIP.py | 103 --------- Wrappers/Python/test/regularizers.py | 399 ----------------------------------- Wrappers/Python/test/simple_demo.py | 103 --------- Wrappers/Python/wip/DemoRecIP.py | 103 +++++++++ Wrappers/Python/wip/regularizers.py | 399 +++++++++++++++++++++++++++++++++++ Wrappers/Python/wip/simple_demo.py | 103 +++++++++ 6 files changed, 605 insertions(+), 605 deletions(-) delete mode 100644 Wrappers/Python/test/DemoRecIP.py delete mode 100644 Wrappers/Python/test/regularizers.py delete mode 100644 Wrappers/Python/test/simple_demo.py create mode 100644 Wrappers/Python/wip/DemoRecIP.py create mode 100644 Wrappers/Python/wip/regularizers.py create mode 100644 Wrappers/Python/wip/simple_demo.py diff --git a/Wrappers/Python/test/DemoRecIP.py b/Wrappers/Python/test/DemoRecIP.py deleted file mode 100644 index 483a202..0000000 --- a/Wrappers/Python/test/DemoRecIP.py +++ /dev/null @@ -1,103 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Reading multi-channel data and reconstruction using FISTA modular -""" - -import numpy as np -import matplotlib.pyplot as plt - -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.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] - -# read IP paper data into a dictionary -dataDICT = read_IPdata() - -# Set ASTRA Projection-backprojection class (fan-beam geometry) -DetWidth = dataDICT.get('im_size')[0] * dataDICT.get('det_width')[0] / \ - dataDICT.get('detectors_numb')[0] -SourceOrig = dataDICT.get('im_size')[0] * dataDICT.get('src_to_rotc')[0] / \ - dataDICT.get('dom_width')[0] -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 - -sino = dataDICT.get('data_norm')[0][:,:,34] # select mid-channel -b = SinogramData(sino) -# Try forward and backprojection -#backprj = Aop.adjoint(b) - -# Create least squares object instance with projector and data. -f = Norm2sq(Aop,b,c=0.5) - -# Initial guess -x_init = VolumeData(np.zeros((dataDICT.get('im_size')[0], - dataDICT.get('im_size')[0]))) - -# Run FISTA for least squares without regularization -opt = {'tol': 1e-4, 'iter': 50} -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) - -# 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) - -plt.imshow(x_fista1.array) -plt.show() -""" \ No newline at end of file diff --git a/Wrappers/Python/test/regularizers.py b/Wrappers/Python/test/regularizers.py deleted file mode 100644 index 04ac3aa..0000000 --- a/Wrappers/Python/test/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/Wrappers/Python/test/simple_demo.py b/Wrappers/Python/test/simple_demo.py deleted file mode 100644 index 1046e7b..0000000 --- a/Wrappers/Python/test/simple_demo.py +++ /dev/null @@ -1,103 +0,0 @@ - -import sys - -sys.path.append("..") - -from ccpi.framework import * -from ccpi.reconstruction.algs import * -from ccpi.reconstruction.funcs import * -from ccpi.reconstruction.ops import * -from ccpi.reconstruction.astra_ops import * -from ccpi.reconstruction.geoms import * - -import numpy as np -import matplotlib.pyplot as plt - -test_case = 2 # 1=parallel2D, 2=cone2D - -# Set up phantom -N = 128 - -x = np.zeros((N,N)) -x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 1.0 -x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 2.0 - -plt.imshow(x) -plt.show() - -vg = VolumeGeometry(N,N,None, 1,1,None) - -Phantom = VolumeData(x,geometry=vg) - -# Set up measurement geometry -angles_num = 20; # angles number - -if test_case==1: - angles = np.linspace(0,np.pi,angles_num,endpoint=False) -elif test_case==2: - angles = np.linspace(0,2*np.pi,angles_num,endpoint=False) -else: - NotImplemented - -det_w = 1.0 -det_num = N -SourceOrig = 200 -OrigDetec = 0 - -# Parallelbeam geometry test -if test_case==1: - pg = SinogramGeometry('parallel', - '2D', - angles, - det_num,det_w) -elif test_case==2: - pg = SinogramGeometry('cone', - '2D', - angles, - det_num, - det_w, - dist_source_center=SourceOrig, - dist_center_detector=OrigDetec) - -# ASTRA operator using volume and sinogram geometries -Aop = AstraProjectorSimple(vg, pg, 'gpu') - -# Unused old astra projector without geometry -# Aop_old = AstraProjector(det_w, det_num, SourceOrig, -# OrigDetec, angles, -# N,'fanbeam','gpu') - -# Try forward and backprojection -b = Aop.direct(Phantom) -out2 = Aop.adjoint(b) - -plt.imshow(b.array) -plt.show() - -plt.imshow(out2.array) -plt.show() - -# Create least squares object instance with projector and data. -f = Norm2sq(Aop,b,c=0.5) - -# Initial guess -x_init = VolumeData(np.zeros(x.shape),geometry=vg) - -# Run FISTA for least squares without regularization -x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None) - -plt.imshow(x_fista0.array) -plt.show() - -# Now least squares plus 1-norm regularization -lam = 0.1 -g0 = Norm1(lam) - -# Run FISTA for least squares plus 1-norm function. -x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0) - -plt.imshow(x_fista1.array) -plt.show() - -# Delete projector -#Aop.delete() diff --git a/Wrappers/Python/wip/DemoRecIP.py b/Wrappers/Python/wip/DemoRecIP.py new file mode 100644 index 0000000..483a202 --- /dev/null +++ b/Wrappers/Python/wip/DemoRecIP.py @@ -0,0 +1,103 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Reading multi-channel data and reconstruction using FISTA modular +""" + +import numpy as np +import matplotlib.pyplot as plt + +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.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] + +# read IP paper data into a dictionary +dataDICT = read_IPdata() + +# Set ASTRA Projection-backprojection class (fan-beam geometry) +DetWidth = dataDICT.get('im_size')[0] * dataDICT.get('det_width')[0] / \ + dataDICT.get('detectors_numb')[0] +SourceOrig = dataDICT.get('im_size')[0] * dataDICT.get('src_to_rotc')[0] / \ + dataDICT.get('dom_width')[0] +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 + +sino = dataDICT.get('data_norm')[0][:,:,34] # select mid-channel +b = SinogramData(sino) +# Try forward and backprojection +#backprj = Aop.adjoint(b) + +# Create least squares object instance with projector and data. +f = Norm2sq(Aop,b,c=0.5) + +# Initial guess +x_init = VolumeData(np.zeros((dataDICT.get('im_size')[0], + dataDICT.get('im_size')[0]))) + +# Run FISTA for least squares without regularization +opt = {'tol': 1e-4, 'iter': 50} +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) + +# 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) + +plt.imshow(x_fista1.array) +plt.show() +""" \ No newline at end of file diff --git a/Wrappers/Python/wip/regularizers.py b/Wrappers/Python/wip/regularizers.py new file mode 100644 index 0000000..04ac3aa --- /dev/null +++ b/Wrappers/Python/wip/regularizers.py @@ -0,0 +1,399 @@ +# -*- 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/Wrappers/Python/wip/simple_demo.py b/Wrappers/Python/wip/simple_demo.py new file mode 100644 index 0000000..1046e7b --- /dev/null +++ b/Wrappers/Python/wip/simple_demo.py @@ -0,0 +1,103 @@ + +import sys + +sys.path.append("..") + +from ccpi.framework import * +from ccpi.reconstruction.algs import * +from ccpi.reconstruction.funcs import * +from ccpi.reconstruction.ops import * +from ccpi.reconstruction.astra_ops import * +from ccpi.reconstruction.geoms import * + +import numpy as np +import matplotlib.pyplot as plt + +test_case = 2 # 1=parallel2D, 2=cone2D + +# Set up phantom +N = 128 + +x = np.zeros((N,N)) +x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 1.0 +x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 2.0 + +plt.imshow(x) +plt.show() + +vg = VolumeGeometry(N,N,None, 1,1,None) + +Phantom = VolumeData(x,geometry=vg) + +# Set up measurement geometry +angles_num = 20; # angles number + +if test_case==1: + angles = np.linspace(0,np.pi,angles_num,endpoint=False) +elif test_case==2: + angles = np.linspace(0,2*np.pi,angles_num,endpoint=False) +else: + NotImplemented + +det_w = 1.0 +det_num = N +SourceOrig = 200 +OrigDetec = 0 + +# Parallelbeam geometry test +if test_case==1: + pg = SinogramGeometry('parallel', + '2D', + angles, + det_num,det_w) +elif test_case==2: + pg = SinogramGeometry('cone', + '2D', + angles, + det_num, + det_w, + dist_source_center=SourceOrig, + dist_center_detector=OrigDetec) + +# ASTRA operator using volume and sinogram geometries +Aop = AstraProjectorSimple(vg, pg, 'gpu') + +# Unused old astra projector without geometry +# Aop_old = AstraProjector(det_w, det_num, SourceOrig, +# OrigDetec, angles, +# N,'fanbeam','gpu') + +# Try forward and backprojection +b = Aop.direct(Phantom) +out2 = Aop.adjoint(b) + +plt.imshow(b.array) +plt.show() + +plt.imshow(out2.array) +plt.show() + +# Create least squares object instance with projector and data. +f = Norm2sq(Aop,b,c=0.5) + +# Initial guess +x_init = VolumeData(np.zeros(x.shape),geometry=vg) + +# Run FISTA for least squares without regularization +x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None) + +plt.imshow(x_fista0.array) +plt.show() + +# Now least squares plus 1-norm regularization +lam = 0.1 +g0 = Norm1(lam) + +# Run FISTA for least squares plus 1-norm function. +x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0) + +plt.imshow(x_fista1.array) +plt.show() + +# Delete projector +#Aop.delete() -- cgit v1.2.3 From 1d0e10be49f157d4f68bba54c0395d84fafd98a3 Mon Sep 17 00:00:00 2001 From: Jakob Jorgensen Date: Wed, 7 Mar 2018 10:02:03 +0000 Subject: Extended geoms with channels attribute, added experimental mc demo script --- Wrappers/Python/ccpi/reconstruction/geoms.py | 45 +++++++++++++++------------- Wrappers/Python/test/simple_mc_demo.py | 43 ++++++++++++++++++++++++++ 2 files changed, 68 insertions(+), 20 deletions(-) create mode 100644 Wrappers/Python/test/simple_mc_demo.py diff --git a/Wrappers/Python/ccpi/reconstruction/geoms.py b/Wrappers/Python/ccpi/reconstruction/geoms.py index edce3b3..9f81fa0 100644 --- a/Wrappers/Python/ccpi/reconstruction/geoms.py +++ b/Wrappers/Python/ccpi/reconstruction/geoms.py @@ -1,16 +1,17 @@ class VolumeGeometry: - def __init__(self, \ - voxel_num_x=None, \ - voxel_num_y=None, \ - voxel_num_z=None, \ - voxel_size_x=None, \ - voxel_size_y=None, \ - voxel_size_z=None, \ - center_x=0, \ - center_y=0, \ - center_z=0): + def __init__(self, + voxel_num_x=None, + voxel_num_y=None, + voxel_num_z=None, + voxel_size_x=None, + voxel_size_y=None, + voxel_size_z=None, + center_x=0, + center_y=0, + center_z=0, + channels=1): self.voxel_num_x = voxel_num_x self.voxel_num_y = voxel_num_y @@ -21,6 +22,7 @@ class VolumeGeometry: self.center_x = center_x self.center_y = center_y self.center_z = center_z + self.channels = channels def getMinX(self): return self.center_x - 0.5*self.voxel_num_x*self.voxel_size_x @@ -43,16 +45,17 @@ class VolumeGeometry: class SinogramGeometry: - def __init__(self, \ - geom_type, \ - dimension, \ - angles, \ - pixel_num_h=None, \ - pixel_size_h=1, \ - pixel_num_v=None, \ - pixel_size_v=1, \ - dist_source_center=None, \ - dist_center_detector=None, \ + def __init__(self, + geom_type, + dimension, + angles, + pixel_num_h=None, + pixel_size_h=1, + pixel_num_v=None, + pixel_size_v=1, + dist_source_center=None, + dist_center_detector=None, + channels=1 ): """ General inputs for standard type projection geometries @@ -91,6 +94,8 @@ class SinogramGeometry: self.pixel_size_h = pixel_size_h self.pixel_num_v = pixel_num_v self.pixel_size_v = pixel_size_v + + self.channels = channels diff --git a/Wrappers/Python/test/simple_mc_demo.py b/Wrappers/Python/test/simple_mc_demo.py new file mode 100644 index 0000000..ca6a89e --- /dev/null +++ b/Wrappers/Python/test/simple_mc_demo.py @@ -0,0 +1,43 @@ + + +import sys + +sys.path.append("..") + +from ccpi.framework import * +from ccpi.reconstruction.algs import * +from ccpi.reconstruction.funcs import * +from ccpi.reconstruction.ops import * +from ccpi.reconstruction.astra_ops import * +from ccpi.reconstruction.geoms import * + +import numpy as np +import matplotlib.pyplot as plt + +test_case = 2 # 1=parallel2D, 2=cone2D + +# Set up phantom +N = 128 + +numchannels = 3 + +x = np.zeros((N,N,numchannels)) + +x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),0] = 1.0 +x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),0] = 2.0 + +x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),1] = 0.7 +x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),1] = 1.2 + +x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),2] = 1.5 +x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),2] = 2.2 + +f, axarr = plt.subplots(1,numchannels) +for k in numpy.arange(3): + axarr[k].imshow(x[:,:,k],vmin=0,vmax=2.5) +plt.show() + +vg = VolumeGeometry(N,N,None, 1,1,None,channels=numchannels) + + +Phantom = VolumeData(x,geometry=vg) \ No newline at end of file -- cgit v1.2.3 From 2f074e52ddbcb69176ad76310c5c51a15658dfa4 Mon Sep 17 00:00:00 2001 From: Jakob Jorgensen Date: Wed, 7 Mar 2018 11:02:10 +0000 Subject: MC FP and BP working --- Wrappers/Python/ccpi/astra/astra_processors.py | 143 +++++++++++++++++++++++ Wrappers/Python/ccpi/framework.py | 30 ++--- Wrappers/Python/ccpi/reconstruction/astra_ops.py | 48 +++++++- Wrappers/Python/test/simple_mc_demo.py | 71 +++++++++-- 4 files changed, 267 insertions(+), 25 deletions(-) diff --git a/Wrappers/Python/ccpi/astra/astra_processors.py b/Wrappers/Python/ccpi/astra/astra_processors.py index 91612b1..650e11b 100644 --- a/Wrappers/Python/ccpi/astra/astra_processors.py +++ b/Wrappers/Python/ccpi/astra/astra_processors.py @@ -1,6 +1,7 @@ from ccpi.framework import DataSetProcessor, DataSet, VolumeData, SinogramData from ccpi.astra.astra_utils import convert_geometry_to_astra import astra +import numpy class AstraForwardProjector(DataSetProcessor): @@ -127,4 +128,146 @@ class AstraBackProjector(DataSetProcessor): DATA = self.getInput() rec_id, IM = astra.create_backprojection(DATA.as_array(), self.proj_id) astra.data2d.delete(rec_id) + return VolumeData(IM,geometry=self.volume_geometry) + +class AstraForwardProjectorMC(DataSetProcessor): + '''AstraForwardProjector + + Forward project VolumeDataSet to SinogramDataSet using ASTRA proj_id. + + Input: VolumeDataSet + Parameter: proj_id + Output: SinogramDataSet + ''' + + def __init__(self, + volume_geometry=None, + sinogram_geometry=None, + proj_id=None, + device='cpu'): + kwargs = { + 'volume_geometry' : volume_geometry, + 'sinogram_geometry' : sinogram_geometry, + 'proj_id' : proj_id, + 'device' : device + } + + #DataSetProcessor.__init__(self, **kwargs) + super(AstraForwardProjectorMC, self).__init__(**kwargs) + + self.setVolumeGeometry(volume_geometry) + self.setSinogramGeometry(sinogram_geometry) + + # Set up ASTRA Volume and projection geometry, not to be stored in self + vol_geom, proj_geom = convert_geometry_to_astra(self.volume_geometry, + self.sinogram_geometry) + + # ASTRA projector, to be stored + if device == 'cpu': + # Note that 'line' is only for parallel (2D) and only one option + self.setProjector(astra.create_projector('line', proj_geom, vol_geom) ) + elif device == 'gpu': + self.setProjector(astra.create_projector('cuda', proj_geom, vol_geom) ) + else: + NotImplemented + + def checkInput(self, dataset): + if dataset.number_of_dimensions == 3 or dataset.number_of_dimensions == 2: + return True + else: + return True + #raise ValueError("Expected input dimensions is 2 or 3, got {0}"\ + # .format(dataset.number_of_dimensions)) + + def setProjector(self, proj_id): + self.proj_id = proj_id + + def setVolumeGeometry(self, volume_geometry): + self.volume_geometry = volume_geometry + + def setSinogramGeometry(self, sinogram_geometry): + self.sinogram_geometry = sinogram_geometry + + def process(self): + IM = self.getInput() + DATA = numpy.zeros((self.sinogram_geometry.angles.size, + self.sinogram_geometry.pixel_num_h, + 1, + self.sinogram_geometry.channels), + 'float32') + for k in range(self.sinogram_geometry.channels): + sinogram_id, DATA[:,:,0,k] = astra.create_sino(IM.as_array()[:,:,0,k], + self.proj_id) + astra.data2d.delete(sinogram_id) + return SinogramData(DATA,geometry=self.sinogram_geometry) + +class AstraBackProjectorMC(DataSetProcessor): + '''AstraBackProjector + + Back project SinogramDataSet to VolumeDataSet using ASTRA proj_id. + + Input: SinogramDataSet + Parameter: proj_id + Output: VolumeDataSet + ''' + + def __init__(self, + volume_geometry=None, + sinogram_geometry=None, + proj_id=None, + device='cpu'): + kwargs = { + 'volume_geometry' : volume_geometry, + 'sinogram_geometry' : sinogram_geometry, + 'proj_id' : proj_id, + 'device' : device + } + + #DataSetProcessor.__init__(self, **kwargs) + super(AstraBackProjectorMC, self).__init__(**kwargs) + + self.setVolumeGeometry(volume_geometry) + self.setSinogramGeometry(sinogram_geometry) + + # Set up ASTRA Volume and projection geometry, not to be stored in self + vol_geom, proj_geom = convert_geometry_to_astra(self.volume_geometry, + self.sinogram_geometry) + + # ASTRA projector, to be stored + if device == 'cpu': + # Note that 'line' is only for parallel (2D) and only one option + self.setProjector(astra.create_projector('line', proj_geom, vol_geom) ) + elif device == 'gpu': + self.setProjector(astra.create_projector('cuda', proj_geom, vol_geom) ) + else: + NotImplemented + + def checkInput(self, dataset): + if dataset.number_of_dimensions == 3 or dataset.number_of_dimensions == 2: + return True + else: + return True + #raise ValueError("Expected input dimensions is 2 or 3, got {0}"\ + # .format(dataset.number_of_dimensions)) + + def setProjector(self, proj_id): + self.proj_id = proj_id + + def setVolumeGeometry(self, volume_geometry): + self.volume_geometry = volume_geometry + + def setSinogramGeometry(self, sinogram_geometry): + self.sinogram_geometry = sinogram_geometry + + def process(self): + DATA = self.getInput() + IM = numpy.zeros((self.volume_geometry.voxel_num_x, + self.volume_geometry.voxel_num_y, + 1, + self.volume_geometry.channels), + 'float32') + for k in range(self.volume_geometry.channels): + rec_id, IM[:,:,0,k] = astra.create_backprojection(DATA.as_array()[:,:,0,k], + self.proj_id) + astra.data2d.delete(rec_id) return VolumeData(IM,geometry=self.volume_geometry) \ No newline at end of file diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py index 5ac17ee..5ee8279 100644 --- a/Wrappers/Python/ccpi/framework.py +++ b/Wrappers/Python/ccpi/framework.py @@ -428,20 +428,20 @@ class VolumeData(DataSet): if type(array) == DataSet: # if the array is a DataSet get the info from there - if not ( array.number_of_dimensions == 2 or \ - array.number_of_dimensions == 3 ): - raise ValueError('Number of dimensions are not 2 or 3: {0}'\ - .format(array.number_of_dimensions)) + #if not ( array.number_of_dimensions == 2 or \ + # array.number_of_dimensions == 3 ): + # raise ValueError('Number of dimensions are not 2 or 3: {0}'\ + # .format(array.number_of_dimensions)) #DataSet.__init__(self, array.as_array(), deep_copy, # array.dimension_labels, **kwargs) super(VolumeData, self).__init__(array.as_array(), deep_copy, array.dimension_labels, **kwargs) elif type(array) == numpy.ndarray: - if not ( array.ndim == 3 or array.ndim == 2 ): - raise ValueError( - 'Number of dimensions are not 3 or 2 : {0}'\ - .format(array.ndim)) + #if not ( array.ndim == 3 or array.ndim == 2 ): + # raise ValueError( + # 'Number of dimensions are not 3 or 2 : {0}'\ + # .format(array.ndim)) if dimension_labels is None: if array.ndim == 3: @@ -476,17 +476,17 @@ class SinogramData(DataSet): if type(array) == DataSet: # if the array is a DataSet get the info from there - if not ( array.number_of_dimensions == 2 or \ - array.number_of_dimensions == 3 ): - raise ValueError('Number of dimensions are not 2 or 3: {0}'\ - .format(array.number_of_dimensions)) + #if not ( array.number_of_dimensions == 2 or \ + # array.number_of_dimensions == 3 ): + # raise ValueError('Number of dimensions are not 2 or 3: {0}'\ + # .format(array.number_of_dimensions)) DataSet.__init__(self, array.as_array(), deep_copy, array.dimension_labels, **kwargs) elif type(array) == numpy.ndarray: - if not ( array.ndim == 3 or array.ndim == 2 ): - raise ValueError('Number of dimensions are != 3: {0}'\ - .format(array.ndim)) + #if not ( array.ndim == 3 or array.ndim == 2 ): + # raise ValueError('Number of dimensions are != 3: {0}'\ + # .format(array.ndim)) if dimension_labels is None: if array.ndim == 3: dimension_labels = ['angle' , diff --git a/Wrappers/Python/ccpi/reconstruction/astra_ops.py b/Wrappers/Python/ccpi/reconstruction/astra_ops.py index 6af7af7..cf138b4 100644 --- a/Wrappers/Python/ccpi/reconstruction/astra_ops.py +++ b/Wrappers/Python/ccpi/reconstruction/astra_ops.py @@ -19,7 +19,7 @@ from ccpi.reconstruction.ops import Operator import numpy from ccpi.framework import SinogramData, VolumeData from ccpi.reconstruction.ops import PowerMethodNonsquare -from ccpi.astra.astra_processors import AstraForwardProjector, AstraBackProjector +from ccpi.astra.astra_processors import * class AstraProjectorSimple(Operator): """ASTRA projector modified to use DataSet and geometry.""" @@ -66,6 +66,52 @@ class AstraProjectorSimple(Operator): self.sinogram_geometry.pixel_num_h), \ (self.volume_geometry.voxel_num_x, \ self.volume_geometry.voxel_num_y) ) + +class AstraProjectorMC(Operator): + """ASTRA Multichannel projector""" + def __init__(self, geomv, geomp, device): + super(AstraProjectorMC, self).__init__() + + # Store volume and sinogram geometries. + self.sinogram_geometry = geomp + self.volume_geometry = geomv + + self.fp = AstraForwardProjectorMC(volume_geometry=geomv, + sinogram_geometry=geomp, + proj_id=None, + device=device) + + self.bp = AstraBackProjectorMC(volume_geometry=geomv, + sinogram_geometry=geomp, + proj_id=None, + device=device) + + # Initialise empty for singular value. + self.s1 = None + + def direct(self, IM): + self.fp.setInput(IM) + out = self.fp.getOutput() + return out + + def adjoint(self, DATA): + self.bp.setInput(DATA) + out = self.bp.getOutput() + return out + + #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): + # Only implemented for 2D + return ( (self.sinogram_geometry.angles.size, \ + self.sinogram_geometry.pixel_num_h), \ + (self.volume_geometry.voxel_num_x, \ + self.volume_geometry.voxel_num_y) ) class AstraProjector(Operator): diff --git a/Wrappers/Python/test/simple_mc_demo.py b/Wrappers/Python/test/simple_mc_demo.py index ca6a89e..1888740 100644 --- a/Wrappers/Python/test/simple_mc_demo.py +++ b/Wrappers/Python/test/simple_mc_demo.py @@ -21,23 +21,76 @@ N = 128 numchannels = 3 -x = np.zeros((N,N,numchannels)) +x = np.zeros((N,N,1,numchannels)) -x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),0] = 1.0 -x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),0] = 2.0 +x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,0] = 1.0 +x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,0] = 2.0 -x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),1] = 0.7 -x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),1] = 1.2 +x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,1] = 0.7 +x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,1] = 1.2 -x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),2] = 1.5 -x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),2] = 2.2 +x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,2] = 1.5 +x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,2] = 2.2 f, axarr = plt.subplots(1,numchannels) for k in numpy.arange(3): - axarr[k].imshow(x[:,:,k],vmin=0,vmax=2.5) + axarr[k].imshow(x[:,:,0,k],vmin=0,vmax=2.5) plt.show() vg = VolumeGeometry(N,N,None, 1,1,None,channels=numchannels) -Phantom = VolumeData(x,geometry=vg) \ No newline at end of file +Phantom = VolumeData(x,geometry=vg) + +# Set up measurement geometry +angles_num = 20; # angles number + +if test_case==1: + angles = np.linspace(0,np.pi,angles_num,endpoint=False) +elif test_case==2: + angles = np.linspace(0,2*np.pi,angles_num,endpoint=False) +else: + NotImplemented + +det_w = 1.0 +det_num = N +SourceOrig = 200 +OrigDetec = 0 + +# Parallelbeam geometry test +if test_case==1: + pg = SinogramGeometry('parallel', + '2D', + angles, + det_num, + det_w, + channels=numchannels) +elif test_case==2: + pg = SinogramGeometry('cone', + '2D', + angles, + det_num, + det_w, + dist_source_center=SourceOrig, + dist_center_detector=OrigDetec, + channels=numchannels) + +# ASTRA operator using volume and sinogram geometries +Aop = AstraProjectorMC(vg, pg, 'gpu') + + +# Try forward and backprojection +b = Aop.direct(Phantom) + +fb, axarrb = plt.subplots(1,numchannels) +for k in numpy.arange(3): + axarrb[k].imshow(b.as_array()[:,:,0,k],vmin=0,vmax=250) +plt.show() + +out2 = Aop.adjoint(b) + +fo, axarro = plt.subplots(1,numchannels) +for k in range(3): + axarro[k].imshow(out2.as_array()[:,:,0,k],vmin=0,vmax=3500) +plt.show() + -- cgit v1.2.3 From c13b4508ec82e6eca27d90c7d72cd4eafe6fbe3e Mon Sep 17 00:00:00 2001 From: Jakob Jorgensen Date: Wed, 7 Mar 2018 11:39:29 +0000 Subject: FISTA seems working with MC --- Wrappers/Python/ccpi/reconstruction/astra_ops.py | 9 ++++-- Wrappers/Python/test/simple_mc_demo.py | 38 ++++++++++++++++++++++++ 2 files changed, 44 insertions(+), 3 deletions(-) diff --git a/Wrappers/Python/ccpi/reconstruction/astra_ops.py b/Wrappers/Python/ccpi/reconstruction/astra_ops.py index cf138b4..1346f22 100644 --- a/Wrappers/Python/ccpi/reconstruction/astra_ops.py +++ b/Wrappers/Python/ccpi/reconstruction/astra_ops.py @@ -87,7 +87,7 @@ class AstraProjectorMC(Operator): device=device) # Initialise empty for singular value. - self.s1 = None + self.s1 = 50 def direct(self, IM): self.fp.setInput(IM) @@ -103,8 +103,11 @@ class AstraProjectorMC(Operator): # astra.data2d.delete(self.proj_id) def get_max_sing_val(self): - self.s1, sall, svec = PowerMethodNonsquare(self,10) - return self.s1 + if self.s1 is None: + self.s1, sall, svec = PowerMethodNonsquare(self,10) + return self.s1 + else: + return self.s1 def size(self): # Only implemented for 2D diff --git a/Wrappers/Python/test/simple_mc_demo.py b/Wrappers/Python/test/simple_mc_demo.py index 1888740..a6959f9 100644 --- a/Wrappers/Python/test/simple_mc_demo.py +++ b/Wrappers/Python/test/simple_mc_demo.py @@ -94,3 +94,41 @@ for k in range(3): axarro[k].imshow(out2.as_array()[:,:,0,k],vmin=0,vmax=3500) plt.show() +# Create least squares object instance with projector and data. +f = Norm2sq(Aop,b,c=0.5) + +# Initial guess +x_init = VolumeData(np.zeros(x.shape),geometry=vg) + +# FISTA options +opt = {'tol': 1e-4, 'iter': 200} + +# Run FISTA for least squares without regularization +x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt) + + +ff0, axarrf0 = plt.subplots(1,numchannels) +for k in numpy.arange(3): + axarrf0[k].imshow(x_fista0.as_array()[:,:,0,k],vmin=0,vmax=2.5) +plt.show() + +plt.semilogy(criter0) +plt.title('Criterion vs iterations, least squares') +plt.show() + +# Now least squares plus 1-norm regularization +lam = 0.1 +g0 = Norm1(lam) + + +# Run FISTA for least squares plus 1-norm function. +x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0, opt) + +ff1, axarrf1 = plt.subplots(1,numchannels) +for k in numpy.arange(3): + axarrf1[k].imshow(x_fista1.as_array()[:,:,0,k],vmin=0,vmax=2.5) +plt.show() + +plt.semilogy(criter1) +plt.title('Criterion vs iterations, least squares plus 1-norm regu') +plt.show() \ No newline at end of file -- cgit v1.2.3 From c8aea54ea281f5990d5ce2966fff30db48e87716 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 7 Mar 2018 14:37:19 +0000 Subject: bugfix in subset, initialize processor with variables --- Wrappers/Python/ccpi/framework.py | 42 ++++++++++++++++++++++++++++++-------- Wrappers/Python/ccpi/processors.py | 14 ++++++++----- 2 files changed, 43 insertions(+), 13 deletions(-) diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py index 5ac17ee..754e0ee 100644 --- a/Wrappers/Python/ccpi/framework.py +++ b/Wrappers/Python/ccpi/framework.py @@ -141,8 +141,9 @@ class DataSet(object): if dimensions is not None: return self.subset(dimensions).as_array() return self.array - - def subset(self, dimensions=None): + + + def subset(self, dimensions=None, **kw): '''Creates a DataSet containing a subset of self according to the labels in dimensions''' if dimensions is None: @@ -164,7 +165,7 @@ class DataSet(object): else: axis_order.append(find_key(self.dimension_labels, dl)) if not proceed: - raise KeyError('Unknown key specified {0}'.format(dl)) + raise KeyError('Subset error: Unknown key specified {0}'.format(dl)) # slice away the unwanted data from the array unwanted_dimensions = self.dimension_labels.copy() @@ -176,13 +177,21 @@ class DataSet(object): #print ("left_dimensions {0}".format(left_dimensions)) #new_shape = [self.shape[ax] for ax in axis_order] #print ("new_shape {0}".format(new_shape)) - command = "self.array" + command = "self.array[" for i in range(self.number_of_dimensions): if self.dimension_labels[i] in unwanted_dimensions.values(): - command = command + "[0]" + value = 0 + for k,v in kw.items(): + if k == self.dimension_labels[i]: + value = v + + command = command + str(value) else: - command = command + "[:]" - #print ("command {0}".format(command)) + command = command + ":" + if i < self.number_of_dimensions -1: + command = command + ',' + command = command + ']' + cleaned = eval(command) # cleaned has collapsed dimensions in the same order of # self.array, but we want it in the order stated in the @@ -413,7 +422,7 @@ class DataSet(object): repres += "Number of dimensions: {0}\n".format(self.number_of_dimensions) repres += "Shape: {0}\n".format(self.shape) repres += "Axis labels: {0}\n".format(self.dimension_labels) - repres += "Representation: {0}\n".format(self.array) + repres += "Representation: \n{0}\n".format(self.array) return repres @@ -772,5 +781,22 @@ if __name__ == '__main__': s = numpy.reshape(numpy.asarray(s), (3,4,4)) sino = SinogramData( s ) + shape = (4,3,2) + a = [i for i in range(2*3*4)] + a = numpy.asarray(a) + a = numpy.reshape(a, shape) + print (numpy.shape(a)) + ds = DataSet(a, True, ['X', 'Y','Z']) + # this means that I expect the X to be of length 2 , + # y of length 3 and z of length 4 + subset = ['Y' ,'Z'] + b0 = ds.subset( subset ) + print ("shape b 3,2? {0}".format(numpy.shape(b0.as_array()))) + # expectation on b is that it is + # 3x2 cut at z = 0 + + subset = ['X' ,'Y'] + b1 = ds.subset( subset , Z=1) + print ("shape b 2,3? {0}".format(numpy.shape(b1.as_array()))) \ No newline at end of file diff --git a/Wrappers/Python/ccpi/processors.py b/Wrappers/Python/ccpi/processors.py index 09e7229..c509a4e 100755 --- a/Wrappers/Python/ccpi/processors.py +++ b/Wrappers/Python/ccpi/processors.py @@ -32,19 +32,23 @@ class Normalizer(DataSetProcessor): Input: SinogramDataSet Parameter: 2D projection with flat field (or stack) 2D projection with dark field (or stack) - Output: SinogramDataSet + Output: SinogramDataSetn ''' - def __init__(self): + def __init__(self, flat_field = None, dark_field = None, tolerance = 1e-5): kwargs = { - 'flat_field' :None, - 'dark_field' :None, + 'flat_field' : None, + 'dark_field' : None, # very small number. Used when there is a division by zero - 'tolerance' : 1e-5 + 'tolerance' : tolerance } #DataSetProcessor.__init__(self, **kwargs) super(Normalizer, self).__init__(**kwargs) + if not flat_field is None: + self.setFlatField(flat_field) + if not dark_field is None: + self.setDarkField(dark_field) def checkInput(self, dataset): if dataset.number_of_dimensions == 3: -- cgit v1.2.3 From c2d2b9f4d17f0f67b205f3f53a0fb8524451dbc1 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 7 Mar 2018 14:39:06 +0000 Subject: added subdir astra (TO BE REMOVED) --- Wrappers/Python/setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Wrappers/Python/setup.py b/Wrappers/Python/setup.py index bbd210f..3f83081 100644 --- a/Wrappers/Python/setup.py +++ b/Wrappers/Python/setup.py @@ -35,7 +35,7 @@ if cil_version == '': setup( name="ccpi-common", version=cil_version, - packages=['ccpi' , 'ccpi.reconstruction'], + packages=['ccpi' , 'ccpi.reconstruction' , 'ccpi.astra'], # Project uses reStructuredText, so ensure that the docutils get # installed or upgraded on the target machine -- cgit v1.2.3 From 406b04e44435dee2b5980512ae484dd816e9332f Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 8 Mar 2018 11:09:10 +0000 Subject: changed default values --- Wrappers/Python/ccpi/reconstruction/geoms.py | 26 +++-- Wrappers/Python/test/simple_mc_demo.py | 134 ------------------------ Wrappers/Python/wip/simple_mc_demo.py | 147 +++++++++++++++++++++++++++ 3 files changed, 163 insertions(+), 144 deletions(-) delete mode 100644 Wrappers/Python/test/simple_mc_demo.py create mode 100755 Wrappers/Python/wip/simple_mc_demo.py diff --git a/Wrappers/Python/ccpi/reconstruction/geoms.py b/Wrappers/Python/ccpi/reconstruction/geoms.py index 9f81fa0..a887c0c 100644 --- a/Wrappers/Python/ccpi/reconstruction/geoms.py +++ b/Wrappers/Python/ccpi/reconstruction/geoms.py @@ -2,12 +2,12 @@ class VolumeGeometry: def __init__(self, - voxel_num_x=None, - voxel_num_y=None, - voxel_num_z=None, - voxel_size_x=None, - voxel_size_y=None, - voxel_size_z=None, + voxel_num_x=0, + voxel_num_y=0, + voxel_num_z=0, + voxel_size_x=1, + voxel_size_y=1, + voxel_size_z=1, center_x=0, center_y=0, center_z=0, @@ -37,10 +37,16 @@ class VolumeGeometry: return self.center_y + 0.5*self.voxel_num_y*self.voxel_size_y def getMinZ(self): - return self.center_z - 0.5*self.voxel_num_z*self.voxel_size_z + if not voxel_num_z == 0: + return self.center_z - 0.5*self.voxel_num_z*self.voxel_size_z + else: + return 0 def getMaxZ(self): - return self.center_z + 0.5*self.voxel_num_z*self.voxel_size_z + if not voxel_num_z == 0: + return self.center_z + 0.5*self.voxel_num_z*self.voxel_size_z + else: + return 0 class SinogramGeometry: @@ -49,9 +55,9 @@ class SinogramGeometry: geom_type, dimension, angles, - pixel_num_h=None, + pixel_num_h=0, pixel_size_h=1, - pixel_num_v=None, + pixel_num_v=0, pixel_size_v=1, dist_source_center=None, dist_center_detector=None, diff --git a/Wrappers/Python/test/simple_mc_demo.py b/Wrappers/Python/test/simple_mc_demo.py deleted file mode 100644 index a6959f9..0000000 --- a/Wrappers/Python/test/simple_mc_demo.py +++ /dev/null @@ -1,134 +0,0 @@ - - -import sys - -sys.path.append("..") - -from ccpi.framework import * -from ccpi.reconstruction.algs import * -from ccpi.reconstruction.funcs import * -from ccpi.reconstruction.ops import * -from ccpi.reconstruction.astra_ops import * -from ccpi.reconstruction.geoms import * - -import numpy as np -import matplotlib.pyplot as plt - -test_case = 2 # 1=parallel2D, 2=cone2D - -# Set up phantom -N = 128 - -numchannels = 3 - -x = np.zeros((N,N,1,numchannels)) - -x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,0] = 1.0 -x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,0] = 2.0 - -x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,1] = 0.7 -x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,1] = 1.2 - -x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,2] = 1.5 -x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,2] = 2.2 - -f, axarr = plt.subplots(1,numchannels) -for k in numpy.arange(3): - axarr[k].imshow(x[:,:,0,k],vmin=0,vmax=2.5) -plt.show() - -vg = VolumeGeometry(N,N,None, 1,1,None,channels=numchannels) - - -Phantom = VolumeData(x,geometry=vg) - -# Set up measurement geometry -angles_num = 20; # angles number - -if test_case==1: - angles = np.linspace(0,np.pi,angles_num,endpoint=False) -elif test_case==2: - angles = np.linspace(0,2*np.pi,angles_num,endpoint=False) -else: - NotImplemented - -det_w = 1.0 -det_num = N -SourceOrig = 200 -OrigDetec = 0 - -# Parallelbeam geometry test -if test_case==1: - pg = SinogramGeometry('parallel', - '2D', - angles, - det_num, - det_w, - channels=numchannels) -elif test_case==2: - pg = SinogramGeometry('cone', - '2D', - angles, - det_num, - det_w, - dist_source_center=SourceOrig, - dist_center_detector=OrigDetec, - channels=numchannels) - -# ASTRA operator using volume and sinogram geometries -Aop = AstraProjectorMC(vg, pg, 'gpu') - - -# Try forward and backprojection -b = Aop.direct(Phantom) - -fb, axarrb = plt.subplots(1,numchannels) -for k in numpy.arange(3): - axarrb[k].imshow(b.as_array()[:,:,0,k],vmin=0,vmax=250) -plt.show() - -out2 = Aop.adjoint(b) - -fo, axarro = plt.subplots(1,numchannels) -for k in range(3): - axarro[k].imshow(out2.as_array()[:,:,0,k],vmin=0,vmax=3500) -plt.show() - -# Create least squares object instance with projector and data. -f = Norm2sq(Aop,b,c=0.5) - -# Initial guess -x_init = VolumeData(np.zeros(x.shape),geometry=vg) - -# FISTA options -opt = {'tol': 1e-4, 'iter': 200} - -# Run FISTA for least squares without regularization -x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt) - - -ff0, axarrf0 = plt.subplots(1,numchannels) -for k in numpy.arange(3): - axarrf0[k].imshow(x_fista0.as_array()[:,:,0,k],vmin=0,vmax=2.5) -plt.show() - -plt.semilogy(criter0) -plt.title('Criterion vs iterations, least squares') -plt.show() - -# Now least squares plus 1-norm regularization -lam = 0.1 -g0 = Norm1(lam) - - -# Run FISTA for least squares plus 1-norm function. -x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0, opt) - -ff1, axarrf1 = plt.subplots(1,numchannels) -for k in numpy.arange(3): - axarrf1[k].imshow(x_fista1.as_array()[:,:,0,k],vmin=0,vmax=2.5) -plt.show() - -plt.semilogy(criter1) -plt.title('Criterion vs iterations, least squares plus 1-norm regu') -plt.show() \ No newline at end of file diff --git a/Wrappers/Python/wip/simple_mc_demo.py b/Wrappers/Python/wip/simple_mc_demo.py new file mode 100755 index 0000000..54c9134 --- /dev/null +++ b/Wrappers/Python/wip/simple_mc_demo.py @@ -0,0 +1,147 @@ + + +import sys + +sys.path.append("..") + +from ccpi.framework import * +from ccpi.reconstruction.algs import * +from ccpi.reconstruction.funcs import * +from ccpi.reconstruction.ops import * +from ccpi.reconstruction.astra_ops import * +from ccpi.reconstruction.geoms import * + +import numpy +import matplotlib.pyplot as plt + +test_case = 1 # 1=parallel2D, 2=cone2D + +# Set up phantom +N = 128 +M = 100 +numchannels = 3 + +vg = VolumeGeometry(voxel_num_x=N,voxel_num_y=N,channels=numchannels) +Phantom = VolumeData(geometry=vg) + +x = Phantom.as_array() +x[0 , round(N/4):round(3*N/4) , round(N/4):round(3*N/4) ] = 1.0 +x[0 , round(N/8):round(7*N/8) , round(3*N/8):round(5*N/8)] = 2.0 + +x[1 , round(N/4):round(3*N/4) , round(N/4):round(3*N/4) ] = 0.7 +x[1 , round(N/8):round(7*N/8) , round(3*N/8):round(5*N/8)] = 1.2 + +x[2 , round(N/4):round(3*N/4) , round(N/4):round(3*N/4) ] = 1.5 +x[2 , round(N/8):round(7*N/8) , round(3*N/8):round(5*N/8)] = 2.2 + +#x = numpy.zeros((N,N,1,numchannels)) + +#x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,0] = 1.0 +#x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,0] = 2.0 + +#x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,1] = 0.7 +#x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,1] = 1.2 + +#x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,2] = 1.5 +#x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,2] = 2.2 + +f, axarr = plt.subplots(1,numchannels) +for k in numpy.arange(3): + axarr[k].imshow(x[k],vmin=0,vmax=2.5) +plt.show() + +#vg = VolumeGeometry(N,N,None, 1,1,None,channels=numchannels) + + +#Phantom = VolumeData(x,geometry=vg) + +# Set up measurement geometry +angles_num = 20; # angles number + +if test_case==1: + angles = numpy.linspace(0,numpy.pi,angles_num,endpoint=False) +elif test_case==2: + angles = numpy.linspace(0,2*numpy.pi,angles_num,endpoint=False) +else: + NotImplemented + +det_w = 1.0 +det_num = N +SourceOrig = 200 +OrigDetec = 0 + +# Parallelbeam geometry test +if test_case==1: + pg = SinogramGeometry('parallel', + '2D', + angles, + det_num, + det_w, + channels=numchannels) +elif test_case==2: + pg = SinogramGeometry('cone', + '2D', + angles, + det_num, + det_w, + dist_source_center=SourceOrig, + dist_center_detector=OrigDetec, + channels=numchannels) + +# ASTRA operator using volume and sinogram geometries +Aop = AstraProjectorMC(vg, pg, 'cpu') + + +# Try forward and backprojection +b = Aop.direct(Phantom) + +fb, axarrb = plt.subplots(1,numchannels) +for k in numpy.arange(3): + axarrb[k].imshow(b.as_array()[k],vmin=0,vmax=250) +plt.show() + +out2 = Aop.adjoint(b) + +fo, axarro = plt.subplots(1,numchannels) +for k in range(3): + axarro[k].imshow(out2.as_array()[k],vmin=0,vmax=3500) +plt.show() + +# Create least squares object instance with projector and data. +f = Norm2sq(Aop,b,c=0.5) + +# Initial guess +x_init = VolumeData(numpy.zeros(x.shape),geometry=vg) + +# FISTA options +opt = {'tol': 1e-4, 'iter': 200} + +# Run FISTA for least squares without regularization +x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt) + + +ff0, axarrf0 = plt.subplots(1,numchannels) +for k in numpy.arange(3): + axarrf0[k].imshow(x_fista0.as_array()[k],vmin=0,vmax=2.5) +plt.show() + +plt.semilogy(criter0) +plt.title('Criterion vs iterations, least squares') +plt.show() + +# Now least squares plus 1-norm regularization +lam = 0.1 +g0 = Norm1(lam) + + +# Run FISTA for least squares plus 1-norm function. +x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0, opt) + +ff1, axarrf1 = plt.subplots(1,numchannels) +for k in numpy.arange(3): + axarrf1[k].imshow(x_fista1.as_array()[k],vmin=0,vmax=2.5) +plt.show() + +plt.semilogy(criter1) +plt.title('Criterion vs iterations, least squares plus 1-norm regu') +plt.show() \ No newline at end of file -- cgit v1.2.3 From 8b3c127facf22d99280257ab81428ce3d72cbd01 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 8 Mar 2018 11:11:38 +0000 Subject: moved files and use SinogramData/VolumeData and Geometries --- Wrappers/Python/wip/simple_demo.py | 4 ++-- Wrappers/Python/wip/simple_mc_demo.py | 5 +++-- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/Wrappers/Python/wip/simple_demo.py b/Wrappers/Python/wip/simple_demo.py index 1046e7b..a970b3e 100644 --- a/Wrappers/Python/wip/simple_demo.py +++ b/Wrappers/Python/wip/simple_demo.py @@ -13,7 +13,7 @@ from ccpi.reconstruction.geoms import * import numpy as np import matplotlib.pyplot as plt -test_case = 2 # 1=parallel2D, 2=cone2D +test_case = 1 # 1=parallel2D, 2=cone2D # Set up phantom N = 128 @@ -60,7 +60,7 @@ elif test_case==2: dist_center_detector=OrigDetec) # ASTRA operator using volume and sinogram geometries -Aop = AstraProjectorSimple(vg, pg, 'gpu') +Aop = AstraProjectorSimple(vg, pg, 'cpu') # Unused old astra projector without geometry # Aop_old = AstraProjector(det_w, det_num, SourceOrig, diff --git a/Wrappers/Python/wip/simple_mc_demo.py b/Wrappers/Python/wip/simple_mc_demo.py index 54c9134..beec717 100755 --- a/Wrappers/Python/wip/simple_mc_demo.py +++ b/Wrappers/Python/wip/simple_mc_demo.py @@ -4,12 +4,12 @@ import sys sys.path.append("..") -from ccpi.framework import * +from ccpi.framework import VolumeData, SinogramData from ccpi.reconstruction.algs import * from ccpi.reconstruction.funcs import * from ccpi.reconstruction.ops import * from ccpi.reconstruction.astra_ops import * -from ccpi.reconstruction.geoms import * +from ccpi.reconstruction.geoms import VolumeGeometry, SinogramGeometry import numpy import matplotlib.pyplot as plt @@ -89,6 +89,7 @@ elif test_case==2: channels=numchannels) # ASTRA operator using volume and sinogram geometries +print("create Projector") Aop = AstraProjectorMC(vg, pg, 'cpu') -- cgit v1.2.3 From 87b2b8a08034eccd676cd49430220cf511f1f8b1 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 8 Mar 2018 11:12:04 +0000 Subject: Multichannel processor works with SinogramData and VolumeData Internally the processor allocates the memory for sinograms and volume data, passing SinogramGeometry and VolumeGeometry to SinogramData and VolumeData. --- Wrappers/Python/ccpi/astra/astra_processors.py | 36 ++++++++++++++------------ 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/Wrappers/Python/ccpi/astra/astra_processors.py b/Wrappers/Python/ccpi/astra/astra_processors.py index 650e11b..8cc0df4 100644 --- a/Wrappers/Python/ccpi/astra/astra_processors.py +++ b/Wrappers/Python/ccpi/astra/astra_processors.py @@ -190,16 +190,18 @@ class AstraForwardProjectorMC(DataSetProcessor): def process(self): IM = self.getInput() - DATA = numpy.zeros((self.sinogram_geometry.angles.size, - self.sinogram_geometry.pixel_num_h, - 1, - self.sinogram_geometry.channels), - 'float32') - for k in range(self.sinogram_geometry.channels): - sinogram_id, DATA[:,:,0,k] = astra.create_sino(IM.as_array()[:,:,0,k], + #create the output Sinogram + DATA = SinogramData(geometry=self.sinogram_geometry) + #DATA = numpy.zeros((self.sinogram_geometry.angles.size, + # self.sinogram_geometry.pixel_num_h, + # 1, + # self.sinogram_geometry.channels), + # 'float32') + for k in range(DATA.geometry.channels): + sinogram_id, DATA.as_array()[k] = astra.create_sino(IM.as_array()[k], self.proj_id) astra.data2d.delete(sinogram_id) - return SinogramData(DATA,geometry=self.sinogram_geometry) + return DATA class AstraBackProjectorMC(DataSetProcessor): '''AstraBackProjector @@ -261,13 +263,15 @@ class AstraBackProjectorMC(DataSetProcessor): def process(self): DATA = self.getInput() - IM = numpy.zeros((self.volume_geometry.voxel_num_x, - self.volume_geometry.voxel_num_y, - 1, - self.volume_geometry.channels), - 'float32') - for k in range(self.volume_geometry.channels): - rec_id, IM[:,:,0,k] = astra.create_backprojection(DATA.as_array()[:,:,0,k], + #IM = numpy.zeros((self.volume_geometry.voxel_num_x, + # self.volume_geometry.voxel_num_y, + # 1, + # self.volume_geometry.channels), + # 'float32') + IM = VolumeData(geometry=self.volume_geometry) + + for k in range(IM.geometry.channels): + rec_id, IM.as_array()[k] = astra.create_backprojection(DATA.as_array()[k], self.proj_id) astra.data2d.delete(rec_id) - return VolumeData(IM,geometry=self.volume_geometry) \ No newline at end of file + return IM \ No newline at end of file -- cgit v1.2.3 From 997e2b22710cba68be023a8d85c0f844710c7e9b Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 8 Mar 2018 11:13:54 +0000 Subject: SinogramData and VolumeData use geometry for initialization One can now pass the sinogram geometry or volume geometry to initialize the SinogramData or VolumeData. It handles multichannel data. --- Wrappers/Python/ccpi/framework.py | 200 ++++++++++++++++++++++++++++---------- 1 file changed, 148 insertions(+), 52 deletions(-) diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py index 82ad65c..7324a5c 100644 --- a/Wrappers/Python/ccpi/framework.py +++ b/Wrappers/Python/ccpi/framework.py @@ -23,6 +23,7 @@ import numpy import sys from datetime import timedelta, datetime import warnings +from ccpi.reconstruction import geoms if sys.version_info[0] >= 3 and sys.version_info[1] >= 4: ABC = abc.ABC @@ -119,7 +120,7 @@ class DataSet(object): if deep_copy: self.array = array[:] else: - self.array = array + self.array = array else: raise TypeError('Array must be NumpyArray, passed {0}'\ .format(type(array))) @@ -430,40 +431,81 @@ class DataSet(object): class VolumeData(DataSet): '''DataSet for holding 2D or 3D dataset''' def __init__(self, - array, + array = None, deep_copy=True, dimension_labels=None, **kwargs): - if type(array) == DataSet: - # if the array is a DataSet get the info from there - #if not ( array.number_of_dimensions == 2 or \ - # array.number_of_dimensions == 3 ): - # raise ValueError('Number of dimensions are not 2 or 3: {0}'\ - # .format(array.number_of_dimensions)) - - #DataSet.__init__(self, array.as_array(), deep_copy, - # array.dimension_labels, **kwargs) - super(VolumeData, self).__init__(array.as_array(), deep_copy, - array.dimension_labels, **kwargs) - elif type(array) == numpy.ndarray: - #if not ( array.ndim == 3 or array.ndim == 2 ): - # raise ValueError( - # 'Number of dimensions are not 3 or 2 : {0}'\ - # .format(array.ndim)) + self.geometry = None + if array is None: + if 'geometry' in kwargs.keys(): + geometry = kwargs['geometry'] + self.geometry = geometry + channels = geometry.channels + horiz_x = geometry.voxel_num_x + horiz_y = geometry.voxel_num_y + vert = 1 if geometry.voxel_num_z is None\ + else geometry.voxel_num_z # this should be 1 for 2D - if dimension_labels is None: - if array.ndim == 3: - dimension_labels = ['horizontal_x' , - 'horizontal_y' , - 'vertical'] + if channels > 1: + if vert > 1: + shape = (channels, vert, horiz_y, horiz_x) + dim_labels = ['channel' ,'vertical' , 'horizontal_y' , + 'horizontal_x'] + else: + shape = (channels , horiz_y, horiz_x) + dim_labels = ['channel' , 'horizontal_y' , + 'horizontal_x'] else: - dimension_labels = ['horizontal' , - 'vertical'] - - #DataSet.__init__(self, array, deep_copy, dimension_labels, **kwargs) - super(VolumeData, self).__init__(array, deep_copy, - dimension_labels, **kwargs) + if vert > 1: + shape = (vert, horiz_y, horiz_x) + dim_labels = ['vertical' , 'horizontal_y' , + 'horizontal_x'] + else: + shape = (horiz_y, horiz_x) + dim_labels = ['horizontal_y' , + 'horizontal_x'] + + array = numpy.zeros( shape , dtype=numpy.float32) + super(VolumeData, self).__init__(array, deep_copy, + dim_labels, **kwargs) + + else: + raise ValueError('Please pass either a DataSet, ' +\ + 'a numpy array or a geometry') + else: + if type(array) == DataSet: + # if the array is a DataSet get the info from there + if not ( array.number_of_dimensions == 2 or \ + array.number_of_dimensions == 3 or \ + array.number_of_dimensions == 4): + raise ValueError('Number of dimensions are not 2 or 3 or 4: {0}'\ + .format(array.number_of_dimensions)) + + #DataSet.__init__(self, array.as_array(), deep_copy, + # array.dimension_labels, **kwargs) + super(VolumeData, self).__init__(array.as_array(), deep_copy, + array.dimension_labels, **kwargs) + elif type(array) == numpy.ndarray: + if not ( array.ndim == 2 or array.ndim == 3 or array.ndim == 4 ): + raise ValueError( + 'Number of dimensions are not 2 or 3 or 4 : {0}'\ + .format(array.ndim)) + + if dimension_labels is None: + if array.ndim == 4: + dimension_labels = ['channel' ,'vertical' , 'horizontal_y' , + 'horizontal_x'] + elif array.ndim == 3: + dimension_labels = ['vertical' , 'horizontal_y' , + 'horizontal_x'] + else: + dimension_labels = ['horizontal_y' , + 'horizontal_x'] + + #DataSet.__init__(self, array, deep_copy, dimension_labels, **kwargs) + super(VolumeData, self).__init__(array, deep_copy, + dimension_labels, **kwargs) # load metadata from kwargs if present for key, value in kwargs.items(): @@ -478,34 +520,79 @@ class VolumeData(DataSet): class SinogramData(DataSet): '''DataSet for holding 2D or 3D sinogram''' def __init__(self, - array, + array = None, deep_copy=True, dimension_labels=None, **kwargs): - - if type(array) == DataSet: - # if the array is a DataSet get the info from there - #if not ( array.number_of_dimensions == 2 or \ - # array.number_of_dimensions == 3 ): - # raise ValueError('Number of dimensions are not 2 or 3: {0}'\ - # .format(array.number_of_dimensions)) - - DataSet.__init__(self, array.as_array(), deep_copy, - array.dimension_labels, **kwargs) - elif type(array) == numpy.ndarray: - #if not ( array.ndim == 3 or array.ndim == 2 ): - # raise ValueError('Number of dimensions are != 3: {0}'\ - # .format(array.ndim)) - if dimension_labels is None: - if array.ndim == 3: - dimension_labels = ['angle' , - 'horizontal' , - 'vertical'] + self.geometry = None + if array is None: + if 'geometry' in kwargs.keys(): + geometry = kwargs['geometry'] + self.geometry = geometry + channels = geometry.channels + horiz = geometry.pixel_num_h + vert = geometry.pixel_num_v + angles = geometry.angles + num_of_angles = numpy.shape(angles)[0] + + + if channels > 1: + if vert > 1: + shape = (channels, num_of_angles , vert, horiz) + dim_labels = ['channel' , ' angle' , + 'vertical' , 'horizontal'] + else: + shape = (channels , num_of_angles, horiz) + dim_labels = ['angle' , 'angle' , + 'horizontal'] else: - dimension_labels = ['angle' , - 'horizontal'] - DataSet.__init__(self, array, deep_copy, dimension_labels, **kwargs) + if vert > 1: + shape = (num_of_angles, vert, horiz) + dim_labels = ['angles' , 'vertical' , + 'horizontal'] + else: + shape = (num_of_angles, horiz) + dim_labels = ['angles' , + 'horizontal'] + + array = numpy.zeros( shape , dtype=numpy.float32) + super(SinogramData, self).__init__(array, deep_copy, + dim_labels, **kwargs) + else: + if type(array) == DataSet: + # if the array is a DataSet get the info from there + if not ( array.number_of_dimensions == 2 or \ + array.number_of_dimensions == 3 or \ + array.number_of_dimensions == 4): + raise ValueError('Number of dimensions are not 2 or 3 or 4: {0}'\ + .format(array.number_of_dimensions)) + + #DataSet.__init__(self, array.as_array(), deep_copy, + # array.dimension_labels, **kwargs) + super(SinogramData, self).__init__(array.as_array(), deep_copy, + array.dimension_labels, **kwargs) + elif type(array) == numpy.ndarray: + if not ( array.ndim == 2 or array.ndim == 3 or array.ndim == 4 ): + raise ValueError( + 'Number of dimensions are not 2 or 3 or 4 : {0}'\ + .format(array.ndim)) + + if dimension_labels is None: + if array.ndim == 4: + dimension_labels = ['channel' ,'angle' , 'vertical' , + 'horizontal'] + elif array.ndim == 3: + dimension_labels = ['angle' , 'vertical' , + 'horizontal'] + else: + dimension_labels = ['angle' , + 'horizontal'] + + #DataSet.__init__(self, array, deep_copy, dimension_labels, **kwargs) + super(SinogramData, self).__init__(array, deep_copy, + dimension_labels, **kwargs) + class DataSetProcessor(object): '''Defines a generic DataSet processor @@ -799,4 +886,13 @@ if __name__ == '__main__': b1 = ds.subset( subset , Z=1) print ("shape b 2,3? {0}".format(numpy.shape(b1.as_array()))) + + # create VolumeData from geometry + vgeometry = geoms.VolumeGeometry(voxel_num_x=2, voxel_num_y=3, channels=2) + vol = VolumeData(geometry=vgeometry) + + sgeometry = geoms.SinogramGeometry(dimension=2, angles=10, + geom_type='parallel', pixel_num_v=3, + pixel_num_h=5 , channels=2) + sino = SinogramData(geometry=sgeometry) \ No newline at end of file -- cgit v1.2.3 From b51914301ba3f7675a600415630a91d81d24074b Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 8 Mar 2018 11:43:05 +0000 Subject: Multichannel projectors inherit from single channel ones --- Wrappers/Python/ccpi/astra/astra_processors.py | 140 +++++-------------------- 1 file changed, 29 insertions(+), 111 deletions(-) diff --git a/Wrappers/Python/ccpi/astra/astra_processors.py b/Wrappers/Python/ccpi/astra/astra_processors.py index 8cc0df4..94a8745 100644 --- a/Wrappers/Python/ccpi/astra/astra_processors.py +++ b/Wrappers/Python/ccpi/astra/astra_processors.py @@ -46,7 +46,8 @@ class AstraForwardProjector(DataSetProcessor): NotImplemented def checkInput(self, dataset): - if dataset.number_of_dimensions == 3 or dataset.number_of_dimensions == 2: + if dataset.number_of_dimensions == 3 or\ + dataset.number_of_dimensions == 2: return True else: raise ValueError("Expected input dimensions is 2 or 3, got {0}"\ @@ -63,9 +64,14 @@ class AstraForwardProjector(DataSetProcessor): def process(self): IM = self.getInput() - sinogram_id, DATA = astra.create_sino(IM.as_array(), self.proj_id) + DATA = SinogramData(geometry=self.sinogram_geometry) + #sinogram_id, DATA = astra.create_sino( IM.as_array(), + # self.proj_id) + sinogram_id, DATA.array = astra.create_sino(IM.as_array(), + self.proj_id) astra.data2d.delete(sinogram_id) - return SinogramData(DATA,geometry=self.sinogram_geometry) + #return SinogramData(array=DATA, geometry=self.sinogram_geometry) + return DATA class AstraBackProjector(DataSetProcessor): '''AstraBackProjector @@ -126,12 +132,14 @@ class AstraBackProjector(DataSetProcessor): def process(self): DATA = self.getInput() - rec_id, IM = astra.create_backprojection(DATA.as_array(), self.proj_id) + IM = VolumeData(geometry=self.volume_geometry) + rec_id, IM.array = astra.create_backprojection(DATA.as_array(), + self.proj_id) astra.data2d.delete(rec_id) - return VolumeData(IM,geometry=self.volume_geometry) + return IM -class AstraForwardProjectorMC(DataSetProcessor): - '''AstraForwardProjector +class AstraForwardProjectorMC(AstraForwardProjector): + '''AstraForwardProjector Multi channel Forward project VolumeDataSet to SinogramDataSet using ASTRA proj_id. @@ -139,72 +147,27 @@ class AstraForwardProjectorMC(DataSetProcessor): Parameter: proj_id Output: SinogramDataSet ''' - - def __init__(self, - volume_geometry=None, - sinogram_geometry=None, - proj_id=None, - device='cpu'): - kwargs = { - 'volume_geometry' : volume_geometry, - 'sinogram_geometry' : sinogram_geometry, - 'proj_id' : proj_id, - 'device' : device - } - - #DataSetProcessor.__init__(self, **kwargs) - super(AstraForwardProjectorMC, self).__init__(**kwargs) - - self.setVolumeGeometry(volume_geometry) - self.setSinogramGeometry(sinogram_geometry) - - # Set up ASTRA Volume and projection geometry, not to be stored in self - vol_geom, proj_geom = convert_geometry_to_astra(self.volume_geometry, - self.sinogram_geometry) - - # ASTRA projector, to be stored - if device == 'cpu': - # Note that 'line' is only for parallel (2D) and only one option - self.setProjector(astra.create_projector('line', proj_geom, vol_geom) ) - elif device == 'gpu': - self.setProjector(astra.create_projector('cuda', proj_geom, vol_geom) ) - else: - NotImplemented - def checkInput(self, dataset): - if dataset.number_of_dimensions == 3 or dataset.number_of_dimensions == 2: + if dataset.number_of_dimensions == 2 or \ + dataset.number_of_dimensions == 3 or \ + dataset.number_of_dimensions == 4: return True else: - return True - #raise ValueError("Expected input dimensions is 2 or 3, got {0}"\ - # .format(dataset.number_of_dimensions)) - - def setProjector(self, proj_id): - self.proj_id = proj_id - - def setVolumeGeometry(self, volume_geometry): - self.volume_geometry = volume_geometry - - def setSinogramGeometry(self, sinogram_geometry): - self.sinogram_geometry = sinogram_geometry - + raise ValueError("Expected input dimensions is 2 or 3, got {0}"\ + .format(dataset.number_of_dimensions)) def process(self): IM = self.getInput() #create the output Sinogram DATA = SinogramData(geometry=self.sinogram_geometry) - #DATA = numpy.zeros((self.sinogram_geometry.angles.size, - # self.sinogram_geometry.pixel_num_h, - # 1, - # self.sinogram_geometry.channels), - # 'float32') + for k in range(DATA.geometry.channels): sinogram_id, DATA.as_array()[k] = astra.create_sino(IM.as_array()[k], self.proj_id) astra.data2d.delete(sinogram_id) return DATA -class AstraBackProjectorMC(DataSetProcessor): - '''AstraBackProjector +class AstraBackProjectorMC(AstraBackProjector): + '''AstraBackProjector Multi channel Back project SinogramDataSet to VolumeDataSet using ASTRA proj_id. @@ -212,62 +175,17 @@ class AstraBackProjectorMC(DataSetProcessor): Parameter: proj_id Output: VolumeDataSet ''' - - def __init__(self, - volume_geometry=None, - sinogram_geometry=None, - proj_id=None, - device='cpu'): - kwargs = { - 'volume_geometry' : volume_geometry, - 'sinogram_geometry' : sinogram_geometry, - 'proj_id' : proj_id, - 'device' : device - } - - #DataSetProcessor.__init__(self, **kwargs) - super(AstraBackProjectorMC, self).__init__(**kwargs) - - self.setVolumeGeometry(volume_geometry) - self.setSinogramGeometry(sinogram_geometry) - - # Set up ASTRA Volume and projection geometry, not to be stored in self - vol_geom, proj_geom = convert_geometry_to_astra(self.volume_geometry, - self.sinogram_geometry) - - # ASTRA projector, to be stored - if device == 'cpu': - # Note that 'line' is only for parallel (2D) and only one option - self.setProjector(astra.create_projector('line', proj_geom, vol_geom) ) - elif device == 'gpu': - self.setProjector(astra.create_projector('cuda', proj_geom, vol_geom) ) - else: - NotImplemented - def checkInput(self, dataset): - if dataset.number_of_dimensions == 3 or dataset.number_of_dimensions == 2: + if dataset.number_of_dimensions == 2 or \ + dataset.number_of_dimensions == 3 or \ + dataset.number_of_dimensions == 4: return True else: - return True - #raise ValueError("Expected input dimensions is 2 or 3, got {0}"\ - # .format(dataset.number_of_dimensions)) - - def setProjector(self, proj_id): - self.proj_id = proj_id - - def setVolumeGeometry(self, volume_geometry): - self.volume_geometry = volume_geometry - - def setSinogramGeometry(self, sinogram_geometry): - self.sinogram_geometry = sinogram_geometry - + raise ValueError("Expected input dimensions is 2 or 3, got {0}"\ + .format(dataset.number_of_dimensions)) def process(self): DATA = self.getInput() - #IM = numpy.zeros((self.volume_geometry.voxel_num_x, - # self.volume_geometry.voxel_num_y, - # 1, - # self.volume_geometry.channels), - # 'float32') + IM = VolumeData(geometry=self.volume_geometry) for k in range(IM.geometry.channels): -- cgit v1.2.3 From 79eeb48ca1856f093cf5a72f26d437b11f4a0c15 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 8 Mar 2018 11:50:19 +0000 Subject: Simple Demo uses the updated VolumeData and SinogramData --- Wrappers/Python/wip/simple_demo.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/Wrappers/Python/wip/simple_demo.py b/Wrappers/Python/wip/simple_demo.py index a970b3e..de2cc99 100644 --- a/Wrappers/Python/wip/simple_demo.py +++ b/Wrappers/Python/wip/simple_demo.py @@ -18,17 +18,16 @@ test_case = 1 # 1=parallel2D, 2=cone2D # Set up phantom N = 128 -x = np.zeros((N,N)) +vg = VolumeGeometry(voxel_num_x=N,voxel_num_y=N) +Phantom = VolumeData(geometry=vg) + +x = Phantom.as_array() x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 1.0 x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 2.0 plt.imshow(x) plt.show() -vg = VolumeGeometry(N,N,None, 1,1,None) - -Phantom = VolumeData(x,geometry=vg) - # Set up measurement geometry angles_num = 20; # angles number -- cgit v1.2.3 From fc7475329efea7234371e96d98fe880f6896b66e Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 8 Mar 2018 11:57:01 +0000 Subject: cleaned imports --- Wrappers/Python/wip/simple_demo.py | 19 ++++++++----------- Wrappers/Python/wip/simple_mc_demo.py | 15 +++++---------- 2 files changed, 13 insertions(+), 21 deletions(-) diff --git a/Wrappers/Python/wip/simple_demo.py b/Wrappers/Python/wip/simple_demo.py index de2cc99..7c5f43e 100644 --- a/Wrappers/Python/wip/simple_demo.py +++ b/Wrappers/Python/wip/simple_demo.py @@ -1,14 +1,11 @@ - -import sys - -sys.path.append("..") - -from ccpi.framework import * -from ccpi.reconstruction.algs import * -from ccpi.reconstruction.funcs import * -from ccpi.reconstruction.ops import * -from ccpi.reconstruction.astra_ops import * -from ccpi.reconstruction.geoms import * +#import sys +#sys.path.append("..") + +from ccpi.framework import VolumeData +from ccpi.reconstruction.algs import FISTA +from ccpi.reconstruction.funcs import Norm2sq, Norm1 +from ccpi.reconstruction.astra_ops import AstraProjectorSimple +from ccpi.reconstruction.geoms import VolumeGeometry, SinogramGeometry import numpy as np import matplotlib.pyplot as plt diff --git a/Wrappers/Python/wip/simple_mc_demo.py b/Wrappers/Python/wip/simple_mc_demo.py index beec717..0d976d7 100755 --- a/Wrappers/Python/wip/simple_mc_demo.py +++ b/Wrappers/Python/wip/simple_mc_demo.py @@ -1,14 +1,10 @@ - - -import sys - -sys.path.append("..") +#import sys +#sys.path.append("..") from ccpi.framework import VolumeData, SinogramData -from ccpi.reconstruction.algs import * -from ccpi.reconstruction.funcs import * -from ccpi.reconstruction.ops import * -from ccpi.reconstruction.astra_ops import * +from ccpi.reconstruction.algs import FISTA +from ccpi.reconstruction.funcs import Norm2sq, Norm1 +from ccpi.reconstruction.astra_ops import AstraProjectorMC from ccpi.reconstruction.geoms import VolumeGeometry, SinogramGeometry import numpy @@ -89,7 +85,6 @@ elif test_case==2: channels=numchannels) # ASTRA operator using volume and sinogram geometries -print("create Projector") Aop = AstraProjectorMC(vg, pg, 'cpu') -- cgit v1.2.3 From af182a56f321d2d03c23d1627991f539a8e2d4bc Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 8 Mar 2018 14:48:28 +0000 Subject: fix of sinogram axis label --- Wrappers/Python/ccpi/framework.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py index 7324a5c..5a507d9 100644 --- a/Wrappers/Python/ccpi/framework.py +++ b/Wrappers/Python/ccpi/framework.py @@ -543,7 +543,7 @@ class SinogramData(DataSet): 'vertical' , 'horizontal'] else: shape = (channels , num_of_angles, horiz) - dim_labels = ['angle' , 'angle' , + dim_labels = ['channel' , 'angle' , 'horizontal'] else: if vert > 1: @@ -891,7 +891,7 @@ if __name__ == '__main__': vgeometry = geoms.VolumeGeometry(voxel_num_x=2, voxel_num_y=3, channels=2) vol = VolumeData(geometry=vgeometry) - sgeometry = geoms.SinogramGeometry(dimension=2, angles=10, + sgeometry = geoms.SinogramGeometry(dimension=2, angles=numpy.linspace(0, 180, num=20), geom_type='parallel', pixel_num_v=3, pixel_num_h=5 , channels=2) sino = SinogramData(geometry=sgeometry) -- cgit v1.2.3