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