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(-) (limited to 'Wrappers/Python') 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.3