diff options
| -rw-r--r-- | Wrappers/Python/demos/PDHG_examples/ColorbayDemo.py | 266 |
1 files changed, 266 insertions, 0 deletions
diff --git a/Wrappers/Python/demos/PDHG_examples/ColorbayDemo.py b/Wrappers/Python/demos/PDHG_examples/ColorbayDemo.py new file mode 100644 index 0000000..0a2126e --- /dev/null +++ b/Wrappers/Python/demos/PDHG_examples/ColorbayDemo.py @@ -0,0 +1,266 @@ +#======================================================================== +# Copyright 2019 Science Technology Facilities Council +# Copyright 2019 University of Manchester +# +# This work is part of the Core Imaging Library developed by Science Technology +# Facilities Council and University of Manchester +# +# 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.txt +# +# 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 ccpi.framework import ImageGeometry, ImageData, AcquisitionGeometry, AcquisitionData, BlockDataContainer + +import numpy as numpy +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG, CGLS +from ccpi.optimisation.algs import CGLS as CGLS_old + +from ccpi.optimisation.operators import BlockOperator, Gradient +from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \ + MixedL21Norm, BlockFunction + +from ccpi.astra.operators import AstraProjectorMC +from scipy.io import loadmat +import h5py + +#%% + +phantom = 'powder' + +if phantom == 'carbon': + pathname = '/media/newhd/shared/Data/ColourBay/spectral_data_sets/CarbonPd/' + filename = 'carbonPd_full_sinogram_stripes_removed.mat' + X = loadmat(pathname + filename) + X = numpy.transpose(X['SS'],(3,1,2,0)) + +elif phantom == 'powder': + pathname = '/media/newhd/shared/DataProcessed/' + filename = 'S_180.mat' + path = pathname + filename + arrays = {} + f = h5py.File(path) + for k, v in f.items(): + arrays[k] = numpy.array(v) + XX = arrays['S'] + X = numpy.transpose(XX,(0,2,1,3)) + X = X[0:250] + + + +#%% Setup Geometry of Colorbay + +num_channels = X.shape[0] +num_pixels_h = X.shape[3] +num_pixels_v = X.shape[2] +num_angles = X.shape[1] + +# Display a single projection in a single channel +plt.imshow(X[100,5,:,:]) +plt.title('Example of a projection image in one channel' ) +plt.show() + +# Set angles to use +angles = numpy.linspace(-numpy.pi,numpy.pi,num_angles,endpoint=False) + +# Define full 3D acquisition geometry and data container. +# Geometric info is taken from the txt-file in the same dir as the mat-file +ag = AcquisitionGeometry('cone', + '3D', + angles, + pixel_num_h=num_pixels_h, + pixel_size_h=0.25, + pixel_num_v=num_pixels_v, + pixel_size_v=0.25, + dist_source_center=233.0, + dist_center_detector=245.0, + channels=num_channels) +data = AcquisitionData(X, geometry=ag) + +# Reduce to central slice by extracting relevant parameters from data and its +# geometry. Perhaps create function to extract central slice automatically? +data2d = data.subset(vertical=40) +ag2d = AcquisitionGeometry('cone', + '2D', + ag.angles, + pixel_num_h=ag.pixel_num_h, + pixel_size_h=ag.pixel_size_h, + pixel_num_v=1, + pixel_size_v=ag.pixel_size_h, + dist_source_center=ag.dist_source_center, + dist_center_detector=ag.dist_center_detector, + channels=ag.channels) +data2d.geometry = ag2d + +# Set up 2D Image Geometry. +# First need the geometric magnification to scale the voxel size relative +# to the detector pixel size. +mag = (ag.dist_source_center + ag.dist_center_detector)/ag.dist_source_center +ig2d = ImageGeometry(voxel_num_x=ag2d.pixel_num_h, + voxel_num_y=ag2d.pixel_num_h, + voxel_size_x=ag2d.pixel_size_h/mag, + voxel_size_y=ag2d.pixel_size_h/mag, + channels=X.shape[0]) + +# Create GPU multichannel projector/backprojector operator with ASTRA. +Aall = AstraProjectorMC(ig2d,ag2d,'gpu') + +# Compute and simple backprojction and display one channel as image. +Xbp = Aall.adjoint(data2d) +plt.imshow(Xbp.subset(channel=100).array) +plt.show() + +#%% CGLS + +x_init = ig2d.allocate() +cgls1 = CGLS(x_init=x_init, operator=Aall, data=data2d) +cgls1.max_iteration = 100 +cgls1.update_objective_interval = 1 +cgls1.run(5,verbose=True) + +plt.imshow(cgls1.get_output().subset(channel=100).array) +plt.title('CGLS') +plt.show() + +#%% Tikhonov + +alpha = 2.5 +Grad = Gradient(ig2d, correlation=Gradient.CORRELATION_SPACE) # use also CORRELATION_SPACECHANNEL + +# Form Tikhonov as a Block CGLS structure +op_CGLS = BlockOperator( Aall, alpha * Grad, shape=(2,1)) +block_data = BlockDataContainer(data2d, Grad.range_geometry().allocate()) + +cgls2 = CGLS(x_init=x_init, operator=op_CGLS, data=block_data) +cgls2.max_iteration = 100 +cgls2.update_objective_interval = 1 + +cgls2.run(10,verbose=True) + +plt.imshow(cgls2.get_output().subset(channel=100).array) +plt.title('Tikhonov') +plt.show() + +#%% Total Variation + +# Regularisation Parameter +alpha_TV = 50 + +# Create operators +op1 = Gradient(ig2d, correlation=Gradient.CORRELATION_SPACE) +op2 = Aall + +# Create BlockOperator +operator = BlockOperator(op1, op2, shape=(2,1) ) + +# Create functions +f1 = alpha * MixedL21Norm() +f2 = 0.5 * L2NormSquared(b=data2d) +f = BlockFunction(f1, f2) +g = ZeroFunction() + +# Compute operator Norm +normK = 8.70320267279591 # Run one time no need to compute again takes time + +# Primal & dual stepsizes +sigma = 10 +tau = 1/(sigma*normK**2) + +# Setup and run the PDHG algorithm +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) +pdhg.max_iteration = 2000 +pdhg.update_objective_interval = 100 +pdhg.run(1000, verbose =True) + + +#%% Show sinograms +channel_ind = [25,75,125] + +plt.figure(figsize=(15,15)) + +plt.subplot(4,3,1) +plt.imshow(data2d.subset(channel = channel_ind[0]).as_array()) +plt.title('Channel {}'.format(channel_ind[0])) +plt.colorbar() + +plt.subplot(4,3,2) +plt.imshow(data2d.subset(channel = channel_ind[1]).as_array()) +plt.title('Channel {}'.format(channel_ind[1])) +plt.colorbar() + +plt.subplot(4,3,3) +plt.imshow(data2d.subset(channel = channel_ind[2]).as_array()) +plt.title('Channel {}'.format(channel_ind[2])) +plt.colorbar() + +############################################################################### +# Show CGLS +plt.subplot(4,3,4) +plt.imshow(cgls1.get_output().subset(channel = channel_ind[0]).as_array()) +plt.colorbar() + +plt.subplot(4,3,5) +plt.imshow(cgls1.get_output().subset(channel = channel_ind[1]).as_array()) +plt.colorbar() + +plt.subplot(4,3,6) +plt.imshow(cgls1.get_output().subset(channel = channel_ind[2]).as_array()) +plt.colorbar() + +############################################################################### +# Show Tikhonov + +plt.subplot(4,3,7) +plt.imshow(cgls2.get_output().subset(channel = channel_ind[0]).as_array()) +plt.colorbar() + +plt.subplot(4,3,8) +plt.imshow(cgls2.get_output().subset(channel = channel_ind[1]).as_array()) +plt.colorbar() + +plt.subplot(4,3,9) +plt.imshow(cgls2.get_output().subset(channel = channel_ind[2]).as_array()) +plt.colorbar() + + +############################################################################### +# Show Total variation + +plt.subplot(4,3,10) +plt.imshow(pdhg.get_output().subset(channel = channel_ind[0]).as_array()) +plt.colorbar() + +plt.subplot(4,3,11) +plt.imshow(pdhg.get_output().subset(channel = channel_ind[1]).as_array()) +plt.colorbar() + +plt.subplot(4,3,12) +plt.imshow(pdhg.get_output().subset(channel = channel_ind[2]).as_array()) +plt.colorbar() + + +############################################################################### + + + + + + + + + + + + |