diff options
Diffstat (limited to 'Wrappers')
| -rw-r--r-- | Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py | 273 |
1 files changed, 273 insertions, 0 deletions
diff --git a/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py b/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py new file mode 100644 index 0000000..ea1ad08 --- /dev/null +++ b/Wrappers/Python/demos/pdhg_TV_tomography2Dccpi.py @@ -0,0 +1,273 @@ +# -*- coding: utf-8 -*- + +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri Feb 22 14:53:03 2019 + +@author: evangelos +""" + +from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer, \ + AcquisitionGeometry, AcquisitionData + +import numpy as np +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG, PDHG_old + +from ccpi.optimisation.operators import BlockOperator, Identity, Gradient +from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ + MixedL21Norm, BlockFunction, ScaledFunction + +#from ccpi.astra.ops import AstraProjectorSimple +from ccpi.plugins.ops import CCPiProjectorSimple +#from skimage.util import random_noise +from timeit import default_timer as timer + +#%% + +#%%############################################################################### +# Create phantom for TV tomography + +#import os +#import tomophantom +#from tomophantom import TomoP2D +#from tomophantom.supp.qualitymetrics import QualityTools + +#model = 1 # select a model number from the library +#N = 150 # set dimension of the phantom +## one can specify an exact path to the parameters file +## path_library2D = '../../../PhantomLibrary/models/Phantom2DLibrary.dat' +#path = os.path.dirname(tomophantom.__file__) +#path_library2D = os.path.join(path, "Phantom2DLibrary.dat") +##This will generate a N_size x N_size phantom (2D) +#phantom_2D = TomoP2D.Model(model, N, path_library2D) +#ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) +#data = ImageData(phantom_2D, geometry=ig) + +N = 150 +#x = np.zeros((N,N)) + +vert = 4 +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N, voxel_num_z=vert) +data = ig.allocate() +Phantom = data +# Populate image data by looping over and filling slices +i = 0 +while i < vert: + if vert > 1: + x = Phantom.subset(vertical=i).array + else: + x = Phantom.array + x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 + x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 0.98 + if vert > 1 : + Phantom.fill(x, vertical=i) + i += 1 + + +#%% +#detectors = N +#angles = np.linspace(0,np.pi,100) +angles_num = 100 +det_w = 1.0 +det_num = N + +angles = np.linspace(0,np.pi,angles_num,endpoint=False,dtype=np.float32)*\ + 180/np.pi + +# Inputs: Geometry, 2D or 3D, angles, horz detector pixel count, +# horz detector pixel size, vert detector pixel count, +# vert detector pixel size. +ag = AcquisitionGeometry('parallel', + '3D', + angles, + N, + det_w, + vert, + det_w) + +from ccpi.reconstruction.parallelbeam import alg as pbalg +#from ccpi.plugins.processors import setupCCPiGeometries +def setupCCPiGeometries(ig, ag, counter): + + #vg = ImageGeometry(voxel_num_x=voxel_num_x,voxel_num_y=voxel_num_y, voxel_num_z=voxel_num_z) + #Phantom_ccpi = ImageData(geometry=vg, + # dimension_labels=['horizontal_x','horizontal_y','vertical']) + ##.subset(['horizontal_x','horizontal_y','vertical']) + ## ask the ccpi code what dimensions it would like + Phantom_ccpi = ig.allocate(dimension_labels=[ImageGeometry.HORIZONTAL_X, + ImageGeometry.HORIZONTAL_Y, + ImageGeometry.VERTICAL]) + + voxel_per_pixel = 1 + angles = ag.angles + geoms = pbalg.pb_setup_geometry_from_image(Phantom_ccpi.as_array(), + angles, + voxel_per_pixel ) + + pg = AcquisitionGeometry('parallel', + '3D', + angles, + geoms['n_h'], 1.0, + geoms['n_v'], 1.0 #2D in 3D is a slice 1 pixel thick + ) + + center_of_rotation = Phantom_ccpi.get_dimension_size('horizontal_x') / 2 + #ad = AcquisitionData(geometry=pg,dimension_labels=['angle','vertical','horizontal']) + ad = pg.allocate(dimension_labels=[AcquisitionGeometry.ANGLE, + AcquisitionGeometry.VERTICAL, + AcquisitionGeometry.HORIZONTAL]) + geoms_i = pbalg.pb_setup_geometry_from_acquisition(ad.as_array(), + angles, + center_of_rotation, + voxel_per_pixel ) + + counter+=1 + + if counter < 4: + print (geoms, geoms_i) + if (not ( geoms_i == geoms )): + print ("not equal and {} {} {}".format(counter, geoms['output_volume_z'], geoms_i['output_volume_z'])) + X = max(geoms['output_volume_x'], geoms_i['output_volume_x']) + Y = max(geoms['output_volume_y'], geoms_i['output_volume_y']) + Z = max(geoms['output_volume_z'], geoms_i['output_volume_z']) + return setupCCPiGeometries(X,Y,Z,angles, counter) + else: + print ("happy now {} {} {}".format(counter, geoms['output_volume_z'], geoms_i['output_volume_z'])) + + return geoms + else: + return geoms_i + + + +#voxel_num_x, voxel_num_y, voxel_num_z, angles, counter +print ("###############################################") +print (ig) +print (ag) +g = setupCCPiGeometries(ig, ag, 0) +print (g) +print ("###############################################") +print ("###############################################") +#ag = AcquisitionGeometry('parallel','2D',angles, detectors) +#Aop = AstraProjectorSimple(ig, ag, 'cpu') +Aop = CCPiProjectorSimple(ig, ag, 'cpu') +sin = Aop.direct(data) + +plt.imshow(sin.subset(vertical=0).as_array()) +plt.title('Sinogram') +plt.colorbar() +plt.show() + + +#%% +# Add Gaussian noise to the sinogram data +np.random.seed(10) +n1 = np.random.random(sin.shape) + +noisy_data = sin + ImageData(5*n1) + +plt.imshow(noisy_data.subset(vertical=0).as_array()) +plt.title('Noisy Sinogram') +plt.colorbar() +plt.show() + + +#%% Works only with Composite Operator Structure of PDHG + +#ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) + +# Create operators +op1 = Gradient(ig) +op2 = Aop + +# Form Composite Operator +operator = BlockOperator(op1, op2, shape=(2,1) ) + +alpha = 50 +f = BlockFunction( alpha * MixedL21Norm(), \ + 0.5 * L2NormSquared(b = noisy_data) ) +g = ZeroFunction() + +normK = Aop.norm() + +# Compute operator Norm +normK = operator.norm() + +## Primal & dual stepsizes +diag_precon = False + +if diag_precon: + + def tau_sigma_precond(operator): + + tau = 1/operator.sum_abs_row() + sigma = 1/ operator.sum_abs_col() + + return tau, sigma + + tau, sigma = tau_sigma_precond(operator) + +else: + sigma = 1 + tau = 1/(sigma*normK**2) + +# Compute operator Norm +normK = operator.norm() + +# Primal & dual stepsizes +sigma = 1 +tau = 1/(sigma*normK**2) + +opt = {'niter':2000} +opt1 = {'niter':2000, 'memopt': True} + + +pdhg1 = PDHG(f=f,g=g, operator=operator, tau=tau, sigma=sigma) +pdhg1.max_iteration = 2000 +pdhg1.update_objective_interval = 200 +pdhg2 = PDHG(f=f,g=g, operator=operator, tau=tau, sigma=sigma, memopt=True) +pdhg2.max_iteration = 2000 +pdhg2.update_objective_interval = 200 +t1 = timer() +#res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) +pdhg1.run(200) +res = pdhg1.get_output().subset(vertical=0) +t2 = timer() + +t3 = timer() +#res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) +pdhg2.run(200) +res1 = pdhg2.get_output().subset(vertical=0) +t4 = timer() +# +print ("No memopt in {}s, memopt in {}s ".format(t2-t1, t4 -t3)) + +#%% +plt.figure(figsize=(15,15)) +plt.subplot(3,1,1) +plt.imshow(res.as_array()) +plt.title('no memopt') +plt.colorbar() +plt.subplot(3,1,2) +plt.imshow(res1.as_array()) +plt.title('memopt') +plt.colorbar() +plt.subplot(3,1,3) +plt.imshow((res1 - res).abs().as_array()) +plt.title('diff') +plt.colorbar() +plt.show() +# +plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'memopt') +plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'no memopt') +plt.legend() +plt.show() +# +print ("Time: No memopt in {}s, \n Time: Memopt in {}s ".format(t2-t1, t4 -t3)) +diff = (res1 - res).abs().as_array().max() +# +print(" Max of abs difference is {}".format(diff)) + |