Skip to content

pwra / NumericalExperiments

Go to a project
Commit ce100ab0 by pwra
Options

Add reconstruction_analysis.py.

parent f6238281
Showing with 186 additions and 24 deletions
dynamic_reconstruction_FBP.py
......@@ -35,10 +35,10 @@ def main(args):
print_step = 0.0
for iFrame in range(num_frame):
if np.isclose(iFrame/num_frame, print_step):
print('Reconstructing time step {} out of {} steps'.format(iFrame, num_frame))
print_step += 0.1
for iSlice in range(num_slice):
if np.isclose(iFrame/num_frame, print_step):
print('Reconstructing time step {} out of {} steps'.format(iFrame, num_frame))
print_step += 0.1
sinograms_id = astra.data2d.create('-sino', proj_geom, projections[iFrame][iSlice, proj_used_idx, :])
rec_id = astra.data2d.create('-vol', vol_geom)
......@@ -56,10 +56,10 @@ def main(args):
astra.algorithm.delete(alg_id)
astra.data2d.delete(rec_id)
astra.data2d.delete(sinograms_id)
file_path = '/home/pwra/simulation/output_files/noise{}/num_proj_used{}/reco_type{}/'.format(noise, args['num_proj_used'], 0)
file_name = 'dynamic_reconstruction_Noise{}_NumProjUsed{}_RecoType{}_NumItr{}.npy'.format(noise, args['num_proj_used'], 0, num_itr)
np.save(file_path+file_name, reconstruction)
# %%
file_path = '/home/pwra/simulation/output_files/noise{}/num_proj_used{}/reco_type{}/'.format(noise, args['num_proj_used'], 0)
file_name = 'dynamic_reconstruction_Noise{}_NumProjUsed{}_RecoType{}_NumItr{}.npy'.format(noise, args['num_proj_used'], 0, num_itr)
np.save(file_path+file_name, reconstruction)
# %%
if __name__ == '__main__':
# Construct the argument parser
......
dynamic_reconstruction_SIRT.py
......@@ -148,7 +148,7 @@ def main(args):
astra.algorithm.delete(alg_id)
astra.data3d.delete(rec_id)
astra.data3d.delete(sinograms_id)
# %%
file_path = '/home/pwra/simulation/output_files/noise{}/num_proj_used{}/reco_type{}/'.format(noise, args['num_proj_used'], reco_type)
file_name = 'dynamic_reconstruction_Noise{}_NumProjUsed{}_RecoType{}_NumItr{}.npy'.format(noise, args['num_proj_used'], reco_type, num_itr)
np.save(file_path+file_name, reconstruction)
......
reconstruction_analysis.py 0 → 100755
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Apr 8 14:22:55 2020
@author: Peter Winkel Rasmussen
"""
# import os
# import sys
import numpy as np
import cupy as cp
import argparse
from skimage.draw import circle
mempool = cp.get_default_memory_pool()
file_path = '/home/pwra/simulation/output_files/'
file_path_plots = '/home/pwra/simulation/plots/'
ground_truth = np.load(file_path+'simulation_attenuation.npy')
num_frame, num_slice, num_row, num_col = ground_truth.shape
# %%
mask_outer = np.ones((num_row, num_col), dtype='bool')
mask_inner = np.zeros((num_row, num_col), dtype='bool')
rr, cc = circle(num_row//2-1, num_col//2-1, 127)
mask_outer[rr, cc] = False
mask_inner[rr, cc] = True
mask_inner_rock = np.tile(mask_inner, (100, 256, 1, 1))
# %%
def main(args):
# %%
noise = args['noise']
num_proj_used = args['num_proj_used']
reco_type = args['reco_type']
if reco_type == 0:
num_itr_list = [None]
elif any(reco_type == np.array([1,2,3])):
num_itr_list = [1, 5, 10, 15, 20, 50, 100, 500]
else:
raise ValueError('Reco type does not exist!\n'
'Select 0,1,2,3 instead of {}'.format(reco_type))
mode = args['mode']
# data = []
res_norm = np.zeros((len(num_itr_list), num_frame), dtype=np.float32)
vol_elements = np.count_nonzero(mask_inner_rock)//num_frame
num_bins = 400
bin_values = np.zeros((len(num_itr_list), num_bins), dtype=np.float32)
min_edge = -2.5
max_edge = 2.5
bin_edges = np.linspace(min_edge, max_edge, num_bins+1, endpoint=True, dtype=np.float32)
tmp_edge = cp.array(bin_edges)
# bin_centers = (bin_edges[:-1]+bin_edges[1:])/2
# %%
if mode == 'cpu':
ground_truth_inner = ground_truth[mask_inner_rock]
for iData, num_itr in enumerate(num_itr_list):
# if iData == 0:
# file_path = '/home/pwra/simulation/output_files/noise{}/num_proj_used{}/reco_type{}/'.format(noise, num_proj_used, 0)
# file_name = 'dynamic_reconstruction_Noise{}_NumProjUsed{}_RecoType{}_NumItr{}.npy'.format(noise, num_proj_used, 0, num_itr)
# else:
file_path = '/home/pwra/simulation/output_files/noise{}/num_proj_used{}/reco_type{}/'.format(noise, num_proj_used, reco_type)
file_name = 'dynamic_reconstruction_Noise{}_NumProjUsed{}_RecoType{}_NumItr{}.npy'.format(noise, num_proj_used, reco_type, num_itr)
data = np.load(file_path+file_name)[mask_inner_rock]
res = data-ground_truth_inner
tmp_bin, tmp_edge = np.histogram(res.ravel(), bins=num_bins, range=(min_edge, max_edge))
bin_values[iData] = tmp_bin
res_norm[iData] = np.linalg.norm(res.reshape(num_frame, vol_elements), axis=1)
ground_truth_inner = None
data = None
res = None
tmp_bin = None
tmp_edge = None
elif mode == 'gpu':
with cp.cuda.Device(args['gpu']):
ground_truth_inner = cp.array(ground_truth[mask_inner_rock])
for iData, num_itr in enumerate(num_itr_list):
if iData == 0:
file_path = '/home/pwra/simulation/output_files/noise{}/num_proj_used{}/reco_type{}/'.format(noise, num_proj_used, 0)
file_name = 'dynamic_reconstruction_Noise{}_NumProjUsed{}_RecoType{}_NumItr{}.npy'.format(noise, num_proj_used, 0, num_itr)
else:
file_path = '/home/pwra/simulation/output_files/noise{}/num_proj_used{}/reco_type{}/'.format(noise, num_proj_used, reco_type)
file_name = 'dynamic_reconstruction_Noise{}_NumProjUsed{}_RecoType{}_NumItr{}.npy'.format(noise, num_proj_used, reco_type, num_itr)
data = cp.array(np.load(file_path+file_name)[mask_inner_rock])
res = data-ground_truth_inner
data = None
mempool.free_all_blocks()
tmp_bin, tmp_edge = cp.histogram(res.ravel(), bins=tmp_edge)
bin_values[iData] = cp.asnumpy(tmp_bin)
res_norm[iData] = cp.asnumpy(cp.linalg.norm(res.reshape(num_frame, vol_elements), axis=1))
res = None
mempool.free_all_blocks()
ground_truth_inner = None
data = None
res = None
tmp_bin = None
tmp_edge = None
mempool.free_all_blocks()
else:
raise ValueError("Chosen mode of calculation does not exist!\n"
"Select 'gpu' or 'cpu' instead of {}".format(mode))
# %%
file_name = 'res_norm_Noise{}_NumProjUsed{}_RecoType{}.npy'.format(noise, num_proj_used, reco_type)
np.save(file_path+file_name, res_norm)
file_name = 'res_binned_Noise{}_NumProjUsed{}_RecoType{}.npy'.format(noise, num_proj_used, reco_type)
np.save(file_path+file_name, bin_values)
# %%
if __name__ == '__main__':
# Construct the argument parser
parser = argparse.ArgumentParser()
# Add the arguments to the parser
parser.add_argument("-num_proj_used", "--num_proj_used", required=False,
nargs='?', const=1, default=360, type=int, help="first operand")
parser.add_argument("-noise", "--noise", required=False,
nargs='?', const=1, default=1.1, type=float, help="second operand")
parser.add_argument("-reco_type", "--reco_type", required=False,
nargs='?', const=1, default=0, type=int, help="third operand")
parser.add_argument("-mode", "--mode", required=False,
nargs='?', const=1, default='cpu', type=str, help="fourth operand")
parser.add_argument("-gpu", "--gpu", required=False,
nargs='?', const=1, default=0, type=int, help="fifth operand")
args = vars(parser.parse_args())
# %%
main(args)
\ No newline at end of file
reconstructions.sh
nohup python dynamic_reconstruction_SIRT.py -reco_type 1 -num_itr 1 -num_proj_used 360 -gpu 1 -noise 1.1 &
nohup python dynamic_reconstruction_SIRT.py -reco_type 1 -num_itr 5 -num_proj_used 360 -gpu 1 -noise 1.1 &
nohup python dynamic_reconstruction_SIRT.py -reco_type 1 -num_itr 10 -num_proj_used 360 -gpu 1 -noise 1.1 &
nohup python dynamic_reconstruction_SIRT.py -reco_type 1 -num_itr 15 -num_proj_used 360 -gpu 1 -noise 1.1 &
nohup python dynamic_reconstruction_SIRT.py -reco_type 1 -num_itr 20 -num_proj_used 360 -gpu 1 -noise 1.1 &
nohup python dynamic_reconstruction_SIRT.py -reco_type 1 -num_itr 50 -num_proj_used 360 -gpu 1 -noise 1.1 &
nohup python dynamic_reconstruction_SIRT.py -reco_type 1 -num_itr 100 -num_proj_used 360 -gpu 0 -noise 1.1 &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 1 -num_itr 500 -num_proj_used 90 -gpu 0 -noise 1.1 &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 1 -num_itr 1 -num_proj_used 180 -gpu 0 -noise 1.1 &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 1 -num_itr 5 -num_proj_used 180 -gpu 0 -noise 1.1 &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 1 -num_itr 10 -num_proj_used 180 -gpu 0 -noise 1.1 &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 1 -num_itr 15 -num_proj_used 180 -gpu 0 -noise 1.1 &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 1 -num_itr 20 -num_proj_used 180 -gpu 0 -noise 1.1 &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 1 -num_itr 50 -num_proj_used 180 -gpu 0 -noise 1.1 &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 1 -num_itr 100 -num_proj_used 180 -gpu 0 -noise 1.1 &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 1 -num_itr 500 -num_proj_used 180 -gpu 1 -noise 1.1 &
num_proj_used_arr=(20 45 90 180 360)
reco_type_arr=(0 1 2 3)
for num_proj_used in "${num_proj_used_arr[@]}"; do
# echo "num_proj_used: $num_proj_used"
for reco_type in "${reco_type_arr[@]}"; do
nohup python reconstruction_analysis.py -num_proj_used $num_proj_used -reco_type $reco_type &
done
done
# nohup python dynamic_reconstruction_SIRT.py -reco_type 0 -num_proj_used 20 -gpu 1 -noise 1.1 &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 0 -num_proj_used 45 -gpu 1 -noise 1.1 &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 0 -num_proj_used 90 -gpu 0 -noise 1.1 &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 0 -num_proj_used 180 -gpu 0 -noise 1.1 &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 1 -num_proj_used 360 -gpu 1 -noise 1.1 > reco_type_3_num_itr_1_num_proj_used_360_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 5 -num_proj_used 360 -gpu 1 -noise 1.1 > reco_type_3_num_itr_5_num_proj_used_360_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 10 -num_proj_used 360 -gpu 1 -noise 1.1 > reco_type_3_num_itr_10_num_proj_used_360_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 15 -num_proj_used 360 -gpu 1 -noise 1.1 > reco_type_3_num_itr_15_num_proj_used_360_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 20 -num_proj_used 360 -gpu 1 -noise 1.1 > reco_type_3_num_itr_20_num_proj_used_360_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 50 -num_proj_used 360 -gpu 1 -noise 1.1 > reco_type_3_num_itr_50_num_proj_used_360_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 100 -num_proj_used 360 -gpu 1 -noise 1.1 > reco_type_3_num_itr_100_num_proj_used_360_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 500 -num_proj_used 360 -gpu 0 -noise 1.1 > reco_type_3_num_itr_500_num_proj_used_360_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 1 -num_proj_used 45 -gpu 0 -noise 1.1 > reco_type_3_num_itr_1_num_proj_used_45_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 5 -num_proj_used 45 -gpu 0 -noise 1.1 > reco_type_3_num_itr_5_num_proj_used_45_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 10 -num_proj_used 45 -gpu 0 -noise 1.1 > reco_type_3_num_itr_10_num_proj_used_45_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 15 -num_proj_used 45 -gpu 0 -noise 1.1 > reco_type_3_num_itr_15_num_proj_used_45_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 20 -num_proj_used 45 -gpu 0 -noise 1.1 > reco_type_3_num_itr_20_num_proj_used_45_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 50 -num_proj_used 45 -gpu 0 -noise 1.1 > reco_type_3_num_itr_50_num_proj_used_45_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 100 -num_proj_used 45 -gpu 0 -noise 1.1 > reco_type_3_num_itr_100_num_proj_used_45_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 500 -num_proj_used 45 -gpu 1 -noise 1.1 > reco_type_3_num_itr_500_num_proj_used_45_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 1 -num_proj_used 90 -gpu 1 -noise 1.1 > reco_type_3_num_itr_1_num_proj_used_90_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 5 -num_proj_used 90 -gpu 1 -noise 1.1 > reco_type_3_num_itr_5_num_proj_used_90_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 10 -num_proj_used 90 -gpu 1 -noise 1.1 > reco_type_3_num_itr_10_num_proj_used_90_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 15 -num_proj_used 90 -gpu 1 -noise 1.1 > reco_type_3_num_itr_15_num_proj_used_90_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 20 -num_proj_used 90 -gpu 1 -noise 1.1 > reco_type_3_num_itr_20_num_proj_used_90_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 50 -num_proj_used 90 -gpu 1 -noise 1.1 > reco_type_3_num_itr_50_num_proj_used_90_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 100 -num_proj_used 90 -gpu 0 -noise 1.1 > reco_type_3_num_itr_100_num_proj_used_90_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 500 -num_proj_used 90 -gpu 0 -noise 1.1 > reco_type_3_num_itr_500_num_proj_used_90_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 1 -num_proj_used 180 -gpu 0 -noise 1.1 > reco_type_3_num_itr_1_num_proj_used_180_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 5 -num_proj_used 180 -gpu 0 -noise 1.1 > reco_type_3_num_itr_5_num_proj_used_180_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 10 -num_proj_used 180 -gpu 0 -noise 1.1 > reco_type_3_num_itr_10_num_proj_used_180_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 15 -num_proj_used 180 -gpu 0 -noise 1.1 > reco_type_3_num_itr_15_num_proj_used_180_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 20 -num_proj_used 180 -gpu 0 -noise 1.1 > reco_type_3_num_itr_20_num_proj_used_180_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 50 -num_proj_used 180 -gpu 0 -noise 1.1 > reco_type_3_num_itr_50_num_proj_used_180_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 100 -num_proj_used 180 -gpu 0 -noise 1.1 > reco_type_3_num_itr_100_num_proj_used_180_noise_1.1.out &
# nohup python dynamic_reconstruction_SIRT.py -reco_type 3 -num_itr 500 -num_proj_used 180 -gpu 1 -noise 1.1 > reco_type_3_num_itr_500_num_proj_used_180_noise_1.1.out &
# for f in output_files/dynamic_reconstruction*.npy
# do old_name=$f
# echo $old_name
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or sign in to comment