In [ ]:
%matplotlib inline
In [ ]:
import os
import glob
import numpy as np
import pylab as plt
import astra
import tomopy
import cv2
from pprint import pprint
import h5py
import astra
In [ ]:
def log_progress(sequence, every=None, size=None):
from ipywidgets import IntProgress, HTML, VBox
from IPython.display import display
is_iterator = False
if size is None:
try:
size = len(sequence)
except TypeError:
is_iterator = True
if size is not None:
if every is None:
if size <= 200:
every = 1
else:
every = size / 200 # every 0.5%
else:
assert every is not None, 'sequence is iterator, set every'
if is_iterator:
progress = IntProgress(min=0, max=1, value=1)
progress.bar_style = 'info'
else:
progress = IntProgress(min=0, max=size, value=0)
label = HTML()
box = VBox(children=[label, progress])
display(box)
index = 0
try:
for index, record in enumerate(sequence, 1):
if index == 1 or index % every == 0:
if is_iterator:
label.value = '{index} / ?'.format(index=index)
else:
progress.value = index
label.value = u'{index} / {size}'.format(
index=index,
size=size
)
yield record
except:
progress.bar_style = 'danger'
raise
else:
progress.bar_style = 'success'
progress.value = index
label.value = str(index or '?')
In [ ]:
working_dir = '/diskmnt/a/makov/tomo_data/bukreeva/ESRF_osso_2016/56A_bone_/'
In [ ]:
config_file = glob.glob(os.path.join(working_dir,'*.cfg'))[0]
In [ ]:
# %load /diskmnt/a/makov/tomo_data/bukreeva/ESRF_osso_2016/56A_bone_/56A_bone_.cfg
local_tmp_dir /tmp
default_data_dir /tmp
acc_max_expo_time 1
disable_motor_poll 0
read_srcur 0
images_at_end_as_quali 0
save_separate_dark_image 0
accel_disp -1
ref_at_end_after_rotating 0
ref_power 0
rounding_correction 0
secured 0
TOMO_DB 1
XSHUTTER_TIME 7
Z_STEP 0
Y_STEP -25
DARK_N 10
REF_N 10
REF_ON 2000
TOMO_EXPTIME 1
TOMO_N 2000
DISTANCE 11000
ENERGY 30
srcur_stop 92.0568
scan_duration 2148.74
acc_expo_time 1
acc_nb_frames 1
display_start_frame 0
real_final_angles 180.001016
real_start_angles
srcur_start 90.8094
ccd_acq_mode SINGLE
tomo_scan_range 180
latency_time 0.005
readout_time 2e-05
ccd_mode FFM
ct_angle 90
no_ignore_mi_errors 0
use_soft_shutter 0
ftm_noshutter 1
no_dir_moved_check 0
beam_check_thres 20
beam_check 0
no_return_at_end 1
auto_update_ref 0
source_sample_distance 145000
start_pos 0
no_xml_file 0
ref_mot_settle_time 2
open_slits_on_quali 0
quali_images_on_ref 0
mono_tune_on_ref_freq 1
mono_tune_on_ref 0
mono_tune_on_start 0
no_check_safshut 0
beeps_at_end 5
no_ref_at_end 0
shifted_images_at_end 0
no_images_at_end 1
minimum_accel_disp 0.001
interlaced 0
return_btw_turns 1
noref_btw_turns 0
shift_turns 1
shift_angle 0
half_orig_pos 0
half_fov_factor 0
half_acq 0
safe_time 0.005
rot_accel_disp -1
speed_corr_factor -1
icat_on 0
use_safshut 1
transfer_to_nice 2
data_dir /data/visitor/md876/id17/Uccelli
focus_is_modifiable 1
magnification 20
ccd_pixel_size 6.5
inited 1
optics_name UserOptics
optics_type User
scintillator GGG_10
focus_scan_steps 10
focus_scan_range 20
focus_lim_neg -1000
focus_lim_pos 1000
focus_prev_step_size 0
focus_step_size 133.333
focus_type Rotation
focus_mot focus
ccd_flip_horz 0
ccd_flip_vert 1
optics_eyepiece 0
optics_objective 20
pixel_size_correct_ask 1
pixel_size_correct 3.04
unbinned_pixel_size 3.04
manual_pixel_size 0
images_per_step 1
no_accel_corr 0
scan_type 2
live_correction 0
distance_macro ftomo_id17_get_distance
mono_macro ftomo_id17_mono_control
info_macro ftomo_id17_get_info
optics_macro ftomo_id17_optics_control
opiom_restore_macro ftomo_id17_opiom_restore
opiom_setup_macro ftomo_id17_opiom_setup
fs_macro ftomo_id17_fastshutter
saf_shut exp2
def_speed_corr_factor -1
cam_x_mot xc
lat_align_mot yrot
ref_mot yrot
rot_mot srot
image_pixel_size 3.04
In [ ]:
# %load /diskmnt/a/makov/tomo_data/bukreeva/ESRF_osso_2016/56A_bone_/56A_bone_.info
Energy= 30
Distance= 11000
Prefix= 56A_bone_
Directory= /data/visitor/md876/id17/Uccelli/56A_bone_
ScanRange= 180
TOMO_N= 2000
REF_ON= 2000
REF_N= 10
DARK_N= 10
Y_STEP= -25
Z_STEP= 0
Dim_1= 2560
Dim_2= 2160
Count_time= 1
Latency_time (s)= 0.005
Shutter_time= 7
Col_end= 2559
Col_beg= 0
Row_end= 2159
Row_beg= 0
Optic_used= 3.04
PixelSize= 2.83
Date= Sat Apr 30 14:24:43 2016
Scan_Type= continuous
CCD_Mode= FFM
SrCurrent= 90.8094
CTAngle= 90
Comment=
In [ ]:
# %load /diskmnt/a/makov/tomo_data/bukreeva/ESRF_osso_2016/56A_bone_/56A_bone_.rec
Scan_Name radix = 56A_bone_
WF_name = ref0005_0000.edf
Number_of_WF in the file = 10
Number_of_projections = 2000
CT_Range (180/360) = 180
Sample_Positioning (left/right) = right
Half_Acquisition (yes/no)? = no
Reconstruction_starts in line = 1
Reconstruction_ends in line = 2160
Rotation_Axis (H/V) = V
###DEFINE ROI####
y_start = 1
y_end = 2160
Width = 2560
Height = 2160
Pixel_Size = 3.04
BinningXY = 1
BinningZ = 1
###DEFINE PBI PHASE RETRIEVAL###
Do_Paganin (yes/no)? = yes
Paganin_Length = 300
Ring_Correction (yes/no) = no
###MISCELLANEOUS###
Only_CT (yes/no)? = yes
Make_Sinos = no
In [ ]:
ref_files = glob.glob(os.path.join(working_dir,'refHST*.edf'))
ref_files = sorted(ref_files)
print len(ref_files)
pprint(ref_files)
dark_files = glob.glob(os.path.join(working_dir,'dark*.edf'))
dark_files = sorted(dark_files)
print(len(dark_files))
pprint(dark_files)
data_files = glob.glob(os.path.join(working_dir,'56A_bone*.edf'))
data_files = sorted(data_files)
print(len(data_files))
# pprint(data_files)
all_files = glob.glob(os.path.join(working_dir,'*.edf'))
other_files = set(all_files)-set(data_files)-set(dark_files)-set(ref_files)
# pprint(other_files)
In [ ]:
def read_edf(x):
return tomopy.io.reader.EdfFile.EdfFile(x).GetData(0)
In [ ]:
ref = read_edf(ref_files[0])
plt.figure(figsize=(15,10))
plt.imshow(ref, cmap=plt.cm.viridis)
plt.colorbar()
plt.title('reference')
plt.show()
In [ ]:
dark = read_edf(dark_files[0])
plt.figure(figsize=(15,10))
plt.imshow(dark, cmap=plt.cm.viridis)
plt.colorbar()
plt.title('dark')
plt.show()
In [ ]:
ref = ref-dark
ref[ref<0]=0
plt.figure(figsize=(15,10))
plt.imshow(ref, cmap=plt.cm.viridis)
plt.colorbar()
plt.title('reference')
plt.show()
In [ ]:
data = read_edf(data_files[1000])
data = data-dark
data = data/ref
plt.figure(figsize=(15,10))
plt.imshow(data, cmap=plt.cm.viridis)
plt.colorbar()
plt.title('reference')
plt.show()
In [ ]:
xt = np.expand_dims(data, axis=0)
xp = tomopy.prep.phase.retrieve_phase(xt, pixel_size=3.04e-4, dist=1100, energy=30)
In [ ]:
plt.figure(figsize=(15,10))
plt.imshow(np.squeeze(xp), cmap=plt.cm.viridis)
plt.colorbar()
plt.title('reference')
plt.show()
In [ ]:
small_ref = ref
small_dark = dark
tmp_sino_file = 'sino.h5'
with h5py.File(tmp_sino_file,'w') as h5f:
for idf,f in enumerate(log_progress(data_files)):
d = read_edf(f)
d = (d-small_dark)/small_ref
dt = np.expand_dims(d, axis=0)
xp = tomopy.prep.phase.retrieve_phase(dt, pixel_size=3.04e-4, dist=1100, energy=30)
h5f.create_dataset(str(idf), data=np.squeeze(xp), compression="gzip")
In [ ]:
tmp_sino_file = 'sino.h5'
sinogram = []
raw = 1000
with h5py.File(tmp_sino_file) as h5f:
for idk, k in enumerate(log_progress(range(len(h5f.keys())))):
sinogram.append(h5f[str(idk)][raw])
sinogram = np.array(sinogram)
In [ ]:
plt.figure(figsize=(15,10))
plt.imshow(sinogram, cmap=plt.cm.viridis)
plt.colorbar()
plt.title('reference')
plt.show()
In [ ]:
sino_pp = sinogram.T/sinogram.sum(axis=-1)*sinogram.sum(axis=-1).mean()
sino_pp = sino_pp.T
In [ ]:
plt.figure(figsize=(15,10))
plt.imshow(sino_pp, cmap=plt.cm.viridis)
plt.colorbar()
plt.title('normalized')
plt.show()
In [ ]:
sino_pp_log = -np.log(sino_pp)
sino_pp_log[sino_pp_log<0] = 0
plt.figure(figsize=(15,10))
plt.imshow(sino_pp_log, cmap=plt.cm.viridis)
plt.colorbar()
plt.title('log')
plt.show()
In [ ]:
plt.figure(figsize=(15,10))
plt.plot(sino_pp_log[0])
plt.plot(np.flipud(sino_pp_log[-1]))
plt.grid()
plt.show()
In [ ]:
delta = 9
s0 = sino_pp_log[0][delta:][1500:2500]
s1 = np.flipud(sino_pp_log[-1][delta:])[1500:2500]
plt.figure(figsize=(15,10))
plt.plot(s0)
plt.plot(s1)
plt.plot((s0-s1)*10)
plt.grid()
plt.show()
In [ ]:
def build_reconstruction_geomety(detector_size, angles):
proj_geom = astra.create_proj_geom('parallel', 1.0, detector_size, angles)
return proj_geom
In [ ]:
def astra_tomo2d(sinogram, angles):
angles = angles.astype('float64') # hack for astra stability, may be removed in future releases
detector_size = sinogram.shape[1]
rec_size = detector_size # size of reconstruction region
vol_geom = astra.create_vol_geom(rec_size, rec_size)
proj_geom = build_reconstruction_geomety(detector_size, angles)
sinogram_id = astra.data2d.create('-sino', proj_geom, data=sinogram)
# Create a data object for the reconstruction
rec_id = astra.data2d.create('-vol', vol_geom)
# Set up the parameters for a reconstruction algorithm using the GPU
cfg = astra.astra_dict('FBP_CUDA')
cfg['ReconstructionDataId'] = rec_id
cfg['ProjectionDataId'] = sinogram_id
cfg['option'] = {}
# cfg['option']['ShortScan'] = True
cfg['option']['MinConstraint'] = 0
# cfg['option']['MaxConstraint'] = 0.02
# Available algorithms:
# SIRT_CUDA, SART_CUDA, EM_CUDA, FBP_CUDA (see the FBP sample)
# Create the algorithm object from the configuration structure
alg_id = astra.algorithm.create(cfg)
# Run 150 iterations of the algorithm
astra.algorithm.run(alg_id, 1000)
# Get the result
rec = astra.data2d.get(rec_id)
# Clean up. Note that GPU memory is tied up in the algorithm object,
# and main RAM in the data objects.
astra.algorithm.delete(alg_id)
astra.data2d.delete(rec_id)
astra.data2d.delete(sinogram_id)
astra.clear()
return rec, proj_geom, cfg
In [ ]:
slices = np.arange(sinogram.shape[0])
In [ ]:
angles = slices*180./2000/180.*np.pi
In [ ]:
rec, proj_geom, cfg = astra_tomo2d(sino_pp_log[:,delta:], angles)
rec= rec/3e-6
In [ ]:
plt.figure(figsize=(15,10))
plt.imshow(rec[1000:2000,1000:2000], cmap=plt.cm.viridis,vmin=30,vmax=50)
plt.colorbar()
plt.title('rec')
plt.show()
In [ ]:
rc = tomopy.misc.corr.remove_ring(np.expand_dims(rec, axis=0),
thresh=100, theta_min=30, rwidth=50)
rc = np.squeeze(rc)
plt.figure(figsize=(15,10))
# plt.imshow((rc-rec)[1200:1800,1200:1800], cmap=plt.cm.viridis)
plt.imshow(rc[1000:2000,1000:2000], cmap=plt.cm.viridis,vmin=0,vmax=200)
# plt.imshow(rc[1000:2000,1000:2000])
plt.colorbar()
plt.title('rec')
plt.show()
In [ ]:
plt.figure(figsize=(10,10))
plt.plot(rc[1500,1800:2000])
plt.grid(True)
In [ ]:
r = read_edf(os.path.join(working_dir,'Slices','56A_bone_rec__1000.edf'))
plt.figure(figsize=(15,10))
plt.imshow(r[1000:2000,1000:2000], cmap=plt.cm.viridis,vmax=0.4)
plt.colorbar()
plt.title('rec')
plt.show()
In [ ]:
rc = tomopy.misc.corr.remove_ring(np.expand_dims(r, axis=0),
thresh=1, theta_min=30, rwidth=30)
rc = np.squeeze(rc)
plt.figure(figsize=(15,10))
plt.imshow(rc[1000:2000,1000:2000], cmap=plt.cm.viridis, vmax=0.4)
plt.colorbar()
plt.title('rec')
plt.show()
In [ ]: