Imat Demo

2 files changed, 350 insertions, 0 deletions

diff --git a/Wrappers/Python/wip/Demos/IMAT_Reconstruction/TV_WhiteBeam_reconstruction.py b/Wrappers/Python/wip/Demos/IMAT_Reconstruction/TV_WhiteBeam_reconstruction.py
new file mode 100644
index 0000000..ef7a9ec
--- /dev/null
+++ b/Wrappers/Python/wip/Demos/IMAT_Reconstruction/TV_WhiteBeam_reconstruction.py
@@ -0,0 +1,164 @@
+# -*- 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-2019 Evangelos Papoutsellis and 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 ccpi.framework import ImageGeometry, AcquisitionGeometry, AcquisitionData
+from astropy.io import fits
+import numpy as np 
+import numpy                          
+import matplotlib.pyplot as plt
+
+from ccpi.optimisation.algorithms import PDHG
+
+from ccpi.optimisation.operators import BlockOperator, Gradient
+from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, L2NormSquared,\
+                      MixedL21Norm, BlockFunction
+
+from ccpi.astra.ops import AstraProjectorSimple
+
+
+# load IMAT file
+#filename_sino = '/media/newhd/shared/DataProcessed/IMAT_beamtime_Feb_2019/preprocessed_test_flat/sino/rebin_slice_350/sino_log_rebin_141.fits' 
+filename_sino = '/media/newhd/shared/DataProcessed/IMAT_beamtime_Feb_2019/preprocessed_test_flat/sino/rebin_slice_350/sino_log_rebin_564.fits' 
+
+sino_handler = fits.open(filename_sino)
+sino_tmp = numpy.array(sino_handler[0].data, dtype=float)
+# reorder sino coordinate: channels, angles, detectors
+sinogram = numpy.rollaxis(sino_tmp, 2)
+sino_handler.close()
+#%%
+# white beam data
+sinogram_wb = sinogram.sum(axis=0)
+
+pixh = sinogram_wb.shape[1] # detectors
+pixv = sinogram_wb.shape[1] # detectors
+
+# WhiteBeam Geometry
+igWB = ImageGeometry(voxel_num_x = pixh, voxel_num_y = pixv)
+
+# Load Golden angles
+with open("/media/newhd/vaggelis/CCPi/CCPi-Block/CCPi-Framework/Wrappers/Python/wip/Demos/IMAT_Reconstruction/golden_angles.txt") as f:
+    angles_string = [line.rstrip() for line in f]
+    angles = numpy.array(angles_string).astype(float)
+agWB = AcquisitionGeometry('parallel', '2D',  angles * numpy.pi / 180, pixh)
+op_WB = AstraProjectorSimple(igWB, agWB, 'gpu')
+sinogram_aqdata = AcquisitionData(sinogram_wb, agWB)
+
+# BackProjection
+result_bp = op_WB.adjoint(sinogram_aqdata)
+
+plt.imshow(result_bp.subset(channel=50).array)
+plt.title('BackProjection')
+plt.show()
+
+
+
+#%%
+
+# Regularisation Parameter
+alpha = 10
+
+# Create operators
+op1 = Gradient(igWB)
+op2 = op_WB
+
+# Create BlockOperator
+operator = BlockOperator(op1, op2, shape=(2,1) ) 
+
+# Create functions
+      
+f1 = alpha * MixedL21Norm()
+f2 = KullbackLeibler(sinogram_aqdata)    
+#f2 = L2NormSquared(b = sinogram_aqdata)  
+f = BlockFunction(f1, f2)  
+                                      
+g = ZeroFunction()
+
+diag_precon =  True
+
+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:
+    # Compute operator Norm
+    normK = operator.norm()
+    print ("normK", normK)
+    # Primal & dual stepsizes
+    sigma = 0.1
+    tau = 1/(sigma*normK**2)
+
+#%%
+
+
+## Primal & dual stepsizes
+sigma = 0.1
+tau = 1/(sigma*normK**2)
+#
+#
+## Setup and run the PDHG algorithm
+pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True)
+pdhg.max_iteration = 10000
+pdhg.update_objective_interval = 500
+
+def circ_mask(h, w, center=None, radius=None):
+
+    if center is None: # use the middle of the image
+        center = [int(w/2), int(h/2)]
+    if radius is None: # use the smallest distance between the center and image walls
+           radius = min(center[0], center[1], w-center[0], h-center[1])
+
+    Y, X = numpy.ogrid[:h, :w]
+    dist_from_center = numpy.sqrt((X - center[0])**2 + (Y-center[1])**2)
+
+    mask = dist_from_center <= radius
+    return mask
+
+def show_result(niter, objective, solution):
+    
+    mask = circ_mask(pixh, pixv, center=None, radius = 220)  # 55 with 141,     
+    plt.imshow(solution.as_array() *  mask)
+    plt.colorbar()
+    plt.title("Iter: {}".format(niter))
+    plt.show()
+    
+
+    print( "{:04}/{:04} {:<5} {:.4f} {:<5} {:.4f} {:<5} {:.4f}".\
+                      format(niter, pdhg.max_iteration,'', \
+                             objective[0],'',\
+                             objective[1],'',\
+                             objective[2]))
+
+pdhg.run(10000, callback = show_result)
+
+#%%
+
+mask = circ_mask(pixh, pixv, center=None, radius = 210) #  55 with 141,
+plt.figure(figsize=(15,15))
+plt.subplot(3,1,1)
+plt.imshow(pdhg.get_output().as_array() *  mask)
+plt.title('Ground Truth')
+plt.colorbar()
+plt.show()
diff --git a/Wrappers/Python/wip/Demos/IMAT_Reconstruction/golden_angles.txt b/Wrappers/Python/wip/Demos/IMAT_Reconstruction/golden_angles.txt
new file mode 100644
index 0000000..95ce73a
--- /dev/null
+++ b/Wrappers/Python/wip/Demos/IMAT_Reconstruction/golden_angles.txt
@@ -0,0 +1,186 @@
+0
+0.9045
+1.809
+2.368
+3.2725
+4.736
+5.6405
+6.1995
+7.104
+8.5675
+9.472
+10.9356
+11.8401
+12.3991
+13.3036
+14.7671
+15.6716
+16.2306
+17.1351
+18.0396
+18.5986
+19.5031
+20.9666
+21.8711
+22.4301
+23.3346
+24.7981
+25.7026
+27.1661
+28.0706
+28.6297
+29.5342
+30.9977
+31.9022
+32.4612
+33.3657
+34.8292
+35.7337
+37.1972
+38.1017
+38.6607
+39.5652
+41.0287
+41.9332
+42.4922
+43.3967
+44.3012
+44.8602
+45.7647
+47.2283
+48.1328
+48.6918
+49.5963
+51.0598
+51.9643
+53.4278
+54.3323
+54.8913
+55.7958
+57.2593
+58.1638
+58.7228
+59.6273
+60.5318
+61.0908
+61.9953
+63.4588
+64.3633
+64.9224
+65.8269
+67.2904
+68.1949
+69.6584
+70.5629
+71.1219
+72.0264
+73.4899
+74.3944
+74.9534
+75.8579
+77.3214
+78.2259
+79.6894
+80.5939
+81.1529
+82.0574
+83.521
+84.4255
+84.9845
+85.889
+86.7935
+87.3525
+88.257
+89.7205
+90.625
+91.184
+92.0885
+93.552
+94.4565
+95.92
+96.8245
+97.3835
+98.288
+99.7516
+100.656
+101.215
+102.12
+103.583
+104.488
+105.951
+106.856
+107.415
+108.319
+109.783
+110.687
+111.246
+112.151
+113.055
+113.614
+114.519
+115.982
+116.887
+117.446
+118.35
+119.814
+120.718
+122.182
+123.086
+123.645
+124.55
+126.013
+126.918
+127.477
+128.381
+129.286
+129.845
+130.749
+132.213
+133.117
+133.676
+134.581
+136.044
+136.949
+138.412
+139.317
+139.876
+140.78
+142.244
+143.148
+143.707
+144.612
+146.075
+146.98
+148.443
+149.348
+149.907
+150.811
+152.275
+153.179
+153.738
+154.643
+155.547
+156.106
+157.011
+158.474
+159.379
+159.938
+160.842
+162.306
+163.21
+164.674
+165.578
+166.137
+167.042
+168.505
+169.41
+169.969
+170.873
+172.337
+173.242
+174.705
+175.609
+176.168
+177.073
+178.536
+179.441