From 885160057882b61862f46b177efbc7d8101f890e Mon Sep 17 00:00:00 2001 From: "Jakob Jorgensen, WS at HMXIF" Date: Thu, 10 May 2018 09:26:49 +0100 Subject: Fix old imports to allow running --- Wrappers/Python/wip/demo_nexus.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Wrappers/Python/wip/demo_nexus.py b/Wrappers/Python/wip/demo_nexus.py index 4dcc9f8..1b40e2b 100755 --- a/Wrappers/Python/wip/demo_nexus.py +++ b/Wrappers/Python/wip/demo_nexus.py @@ -8,10 +8,10 @@ Created on Wed Mar 21 14:26:21 2018 from ccpi.framework import ImageData , AcquisitionData, ImageGeometry, AcquisitionGeometry from ccpi.optimisation.algs import FISTA, FBPD, CGLS from ccpi.optimisation.funcs import Norm2sq, Norm1 -from ccpi.reconstruction.ccpiops import CCPiProjectorSimple +from ccpi.plugins.ops import CCPiProjectorSimple from ccpi.reconstruction.parallelbeam import alg as pbalg -from ccpi.reconstruction.processors import CCPiForwardProjector, CCPiBackwardProjector , \ -Normalizer , CenterOfRotationFinder , AcquisitionDataPadder +from ccpi.plugins.processors import CCPiForwardProjector, CCPiBackwardProjector +from ccpi.processors import Normalizer , CenterOfRotationFinder , AcquisitionDataPadder from ccpi.io.reader import NexusReader @@ -31,7 +31,7 @@ def avg_img(image): return avg -reader = NexusReader(os.path.join(".." ,".." ,".." , "data" , "24737_fd.nxs" )) +reader = NexusReader(os.path.join(".." ,".." ,".." , "..", "CCPi-ReconstructionFramework","data" , "24737_fd.nxs" )) dims = reader.get_projection_dimensions() print (dims) -- cgit v1.2.1 From e9905059de627978bdcb8a6eda118d10f445d8a7 Mon Sep 17 00:00:00 2001 From: "Jakob Jorgensen, WS at HMXIF" Date: Thu, 10 May 2018 10:11:01 +0100 Subject: Add comment etc. --- Wrappers/Python/wip/demo_nexus.py | 110 ++++++++++++++++++-------------------- 1 file changed, 53 insertions(+), 57 deletions(-) diff --git a/Wrappers/Python/wip/demo_nexus.py b/Wrappers/Python/wip/demo_nexus.py index 1b40e2b..a51ec2d 100755 --- a/Wrappers/Python/wip/demo_nexus.py +++ b/Wrappers/Python/wip/demo_nexus.py @@ -1,27 +1,28 @@ -# -*- coding: utf-8 -*- -""" -Created on Wed Mar 21 14:26:21 2018 -@author: ofn77899 -""" +# This script demonstrates how to load a parallel beam data set in Nexus +# format, apply dark and flat field correction and reconstruct using the +# modular optimisation framework. +# +# The data set is available from +# https://github.com/DiamondLightSource/Savu/blob/master/test_data/data/24737_fd.nxs +# and should be downloaded to a local directory to be specified below. -from ccpi.framework import ImageData , AcquisitionData, ImageGeometry, AcquisitionGeometry +# All own imports +from ccpi.framework import ImageData, AcquisitionData, ImageGeometry, AcquisitionGeometry from ccpi.optimisation.algs import FISTA, FBPD, CGLS from ccpi.optimisation.funcs import Norm2sq, Norm1 from ccpi.plugins.ops import CCPiProjectorSimple from ccpi.reconstruction.parallelbeam import alg as pbalg from ccpi.plugins.processors import CCPiForwardProjector, CCPiBackwardProjector from ccpi.processors import Normalizer , CenterOfRotationFinder , AcquisitionDataPadder - from ccpi.io.reader import NexusReader +# All external imports import numpy import matplotlib.pyplot as plt - import os -import pickle - +# Define utility function to average over flat and dark images. def avg_img(image): shape = list(numpy.shape(image)) l = shape.pop(0) @@ -30,116 +31,115 @@ def avg_img(image): avg += image[i] / l return avg - +# Set up a reader object pointing to the Nexus data set. Revise path as needed. reader = NexusReader(os.path.join(".." ,".." ,".." , "..", "CCPi-ReconstructionFramework","data" , "24737_fd.nxs" )) +# Read and print the dimensions of the raw projections dims = reader.get_projection_dimensions() print (dims) +# Load and average all flat and dark images in preparation for normalising data. flat = avg_img(reader.load_flat()) dark = avg_img(reader.load_dark()) +# Set up normaliser object for normalising data by flat and dark images. norm = Normalizer(flat_field=flat, dark_field=dark) +# Load the raw projections and pass as input to the normaliser. norm.set_input(reader.get_acquisition_data()) +# Set up CenterOfRotationFinder object to center data. cor = CenterOfRotationFinder() + +# Set the output of the normaliser as the input and execute to determine center. cor.set_input(norm.get_output()) center_of_rotation = cor.get_output() -voxel_per_pixel = 1 +# Set up AcquisitionDataPadder to pad data for centering using the computed +# center, set the output of the normaliser as input and execute to produce +# padded/centered data. padder = AcquisitionDataPadder(center_of_rotation=center_of_rotation) padder.set_input(norm.get_output()) padded_data = padder.get_output() -pg = padded_data.geometry +# Create Acquisition and Image Geometries for setting up projector. +voxel_per_pixel = 1 +ag = padded_data.geometry geoms = pbalg.pb_setup_geometry_from_acquisition(padded_data.as_array(), - pg.angles, + ag.angles, center_of_rotation, voxel_per_pixel ) -vg = ImageGeometry(voxel_num_x=geoms['output_volume_x'], +ig = ImageGeometry(voxel_num_x=geoms['output_volume_x'], voxel_num_y=geoms['output_volume_y'], voxel_num_z=geoms['output_volume_z']) -#data = numpy.reshape(reader.getAcquisitionData()) -print ("define projector") -Cop = CCPiProjectorSimple(vg, pg) + +# Define the projector object +print ("Define projector") +Cop = CCPiProjectorSimple(ig, ag) + # Create least squares object instance with projector and data. print ("Create least squares object instance with projector and data.") f = Norm2sq(Cop,padded_data,c=0.5) + +# Set initial guess print ("Initial guess") -# Initial guess -x_init = ImageData(geometry=vg, dimension_labels=['horizontal_x','horizontal_y','vertical']) +x_init = ImageData(geometry=ig, dimension_labels=['horizontal_x','horizontal_y','vertical']) -#%% -print ("run FISTA") -# Run FISTA for least squares without regularization +# Run FISTA reconstruction for least squares without regularization +print ("run FISTA for least squares") opt = {'tol': 1e-4, 'iter': 10} x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt=opt) -pickle.dump(x_fista0, open("fista0.pkl", "wb")) - plt.imshow(x_fista0.subset(horizontal_x=80).array) plt.title('FISTA0') -#plt.show() +plt.show() -# Now least squares plus 1-norm regularization +# Set up 1-norm function for FISTA least squares plus 1-norm regularisation +print ("Run FISTA for least squares plus 1-norm regularisation") 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=opt) -pickle.dump(x_fista1, open("fista1.pkl", "wb")) plt.imshow(x_fista0.subset(horizontal_x=80).array) plt.title('FISTA1') -#plt.show() - -plt.semilogy(criter1) -#plt.show() +plt.show() # Run FBPD=Forward Backward Primal Dual method on least squares plus 1-norm +print ("Run FBPD for least squares plus 1-norm regularisation") x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt=opt) -pickle.dump(x_fbpd1, open("fbpd1.pkl", "wb")) plt.imshow(x_fbpd1.subset(horizontal_x=80).array) plt.title('FBPD1') -#plt.show() - -plt.semilogy(criter_fbpd1) -#plt.show() +plt.show() -# Run CGLS, which should agree with the FISTA0 +# Run CGLS, which should agree with the FISTA least squares +print ("Run CGLS for least squares") x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Cop, padded_data, opt=opt) -pickle.dump(x_CGLS, open("cgls.pkl", "wb")) plt.imshow(x_CGLS.subset(horizontal_x=80).array) plt.title('CGLS') -plt.title('CGLS recon, compare FISTA0') -#plt.show() - -plt.semilogy(criter_CGLS) -plt.title('CGLS criterion') -#plt.show() - +plt.show() +# Display all reconstructions and decay of objective function cols = 4 rows = 1 current = 1 fig = plt.figure() -# projections row current = current a=fig.add_subplot(rows,cols,current) -a.set_title('FISTA0') +a.set_title('FISTA LS') imgplot = plt.imshow(x_fista0.subset(horizontal_x=80).as_array()) current = current + 1 a=fig.add_subplot(rows,cols,current) -a.set_title('FISTA1') +a.set_title('FISTA LS+1') imgplot = plt.imshow(x_fista1.subset(horizontal_x=80).as_array()) current = current + 1 a=fig.add_subplot(rows,cols,current) -a.set_title('FBPD1') +a.set_title('FBPD LS+1') imgplot = plt.imshow(x_fbpd1.subset(horizontal_x=80).as_array()) current = current + 1 @@ -149,16 +149,12 @@ imgplot = plt.imshow(x_CGLS.subset(horizontal_x=80).as_array()) plt.show() - -#%% fig = plt.figure() -# projections row b=fig.add_subplot(1,1,1) b.set_title('criteria') -imgplot = plt.loglog(criter0 , label='FISTA0') -imgplot = plt.loglog(criter1 , label='FISTA1') -imgplot = plt.loglog(criter_fbpd1, label='FBPD1') +imgplot = plt.loglog(criter0 , label='FISTA LS') +imgplot = plt.loglog(criter1 , label='FISTA LS+1') +imgplot = plt.loglog(criter_fbpd1, label='FBPD LS+1') imgplot = plt.loglog(criter_CGLS, label='CGLS') -#imgplot = plt.loglog(criter_fbpdtv, label='FBPD TV') b.legend(loc='right') plt.show() \ No newline at end of file -- cgit v1.2.1 From c2062f3713fce3dc414a5a214a80c1c29f757dd4 Mon Sep 17 00:00:00 2001 From: "Jakob Jorgensen, WS at HMXIF" Date: Thu, 10 May 2018 10:17:03 +0100 Subject: Exclude unnec import --- Wrappers/Python/wip/demo_nexus.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/Wrappers/Python/wip/demo_nexus.py b/Wrappers/Python/wip/demo_nexus.py index a51ec2d..ffa5741 100755 --- a/Wrappers/Python/wip/demo_nexus.py +++ b/Wrappers/Python/wip/demo_nexus.py @@ -13,8 +13,7 @@ from ccpi.optimisation.algs import FISTA, FBPD, CGLS from ccpi.optimisation.funcs import Norm2sq, Norm1 from ccpi.plugins.ops import CCPiProjectorSimple from ccpi.reconstruction.parallelbeam import alg as pbalg -from ccpi.plugins.processors import CCPiForwardProjector, CCPiBackwardProjector -from ccpi.processors import Normalizer , CenterOfRotationFinder , AcquisitionDataPadder +from ccpi.processors import Normalizer, CenterOfRotationFinder, AcquisitionDataPadder from ccpi.io.reader import NexusReader # All external imports @@ -86,12 +85,12 @@ print ("Initial guess") x_init = ImageData(geometry=ig, dimension_labels=['horizontal_x','horizontal_y','vertical']) # Run FISTA reconstruction for least squares without regularization -print ("run FISTA for least squares") +print ("Run FISTA for least squares") opt = {'tol': 1e-4, 'iter': 10} x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt=opt) plt.imshow(x_fista0.subset(horizontal_x=80).array) -plt.title('FISTA0') +plt.title('FISTA LS') plt.show() # Set up 1-norm function for FISTA least squares plus 1-norm regularisation @@ -103,7 +102,7 @@ g0 = Norm1(lam) x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0,opt=opt) plt.imshow(x_fista0.subset(horizontal_x=80).array) -plt.title('FISTA1') +plt.title('FISTA LS+1') plt.show() # Run FBPD=Forward Backward Primal Dual method on least squares plus 1-norm @@ -111,7 +110,7 @@ print ("Run FBPD for least squares plus 1-norm regularisation") x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt=opt) plt.imshow(x_fbpd1.subset(horizontal_x=80).array) -plt.title('FBPD1') +plt.title('FBPD LS+1') plt.show() # Run CGLS, which should agree with the FISTA least squares -- cgit v1.2.1 From a3dd21598a1d61fec1bc8b7e17abe74d5dbd6a51 Mon Sep 17 00:00:00 2001 From: "Jakob Jorgensen, WS at HMXIF" Date: Thu, 10 May 2018 11:29:41 +0100 Subject: Remove geoms and pbalg dependence. --- Wrappers/Python/wip/demo_nexus.py | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/Wrappers/Python/wip/demo_nexus.py b/Wrappers/Python/wip/demo_nexus.py index ffa5741..03739b1 100755 --- a/Wrappers/Python/wip/demo_nexus.py +++ b/Wrappers/Python/wip/demo_nexus.py @@ -12,7 +12,6 @@ from ccpi.framework import ImageData, AcquisitionData, ImageGeometry, Acquisitio from ccpi.optimisation.algs import FISTA, FBPD, CGLS from ccpi.optimisation.funcs import Norm2sq, Norm1 from ccpi.plugins.ops import CCPiProjectorSimple -from ccpi.reconstruction.parallelbeam import alg as pbalg from ccpi.processors import Normalizer, CenterOfRotationFinder, AcquisitionDataPadder from ccpi.io.reader import NexusReader @@ -62,15 +61,10 @@ padder.set_input(norm.get_output()) padded_data = padder.get_output() # Create Acquisition and Image Geometries for setting up projector. -voxel_per_pixel = 1 ag = padded_data.geometry -geoms = pbalg.pb_setup_geometry_from_acquisition(padded_data.as_array(), - ag.angles, - center_of_rotation, - voxel_per_pixel ) -ig = ImageGeometry(voxel_num_x=geoms['output_volume_x'], - voxel_num_y=geoms['output_volume_y'], - voxel_num_z=geoms['output_volume_z']) +ig = ImageGeometry(voxel_num_x=ag.pixel_num_h, + voxel_num_y=ag.pixel_num_h, + voxel_num_z=ag.pixel_num_v) # Define the projector object print ("Define projector") -- cgit v1.2.1