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 = '/home/makov/diskmnt/big/tomo_data/bukreeva/ESRF_osso_2016/spine_mf_10x_mf'
In [ ]:
config_file = '/diskmnt/a/makov/tomo_data/bukreeva/ESRF_osso_2016/spine_mf_10x/tif/spine_mf_10x_.log'
In [ ]:
# %load /diskmnt/a/makov/tomo_data/bukreeva/ESRF_osso_2016/spine_mf_10x/tif/spine_mf_10x_.log
User ID : e14251
FAST-TOMO scan of sample spine_mf_10x_ started on Fri Jun 21 11:34:56 2013
--------------------Beamline Settings-------------------------
Ring current [mA] : 401.846
Beam energy [keV] : 23.999
Monostripe : W/Si
FE-Filter : Filter 50%
OP-Filter 1 : 100um Al
OP-Filter 2 : 10um Cu
OP-Filter 3 : 10um Fe
--------------------Detector Settings-------------------------
Camera : PCO.Edge
Microscope : Opt.Peter MB op
Magnification : 10.00
Scintillator : LuAG:Ce 20um
Exposure time [ms] : 500
Delay time [ms] : 0
Millisecond shutter [ms] : not used
X-ROI : 1 - 2560
Y-ROI : 1 - 2160
Actual pixel size [um] : 0.65
------------------------Scan Settings-------------------------
Sample folder : /sls/X02DA/data/e14251/Data10/disk1/spine_mf_10x_
File Prefix : spine_mf_10x_
Number of scans : 1
Number of projections : 1601
Number of darks : 32
Number of flats : 160
Number of inter-flats : 0
Flat frequency : 0
Readout frequency : 0
Rot Y min position [deg] : 0.000
Rot Y max position [deg] : 180.000
Angular step [deg] : 0.113
Sample In [um] : 0
Sample Out [um] : 10000
-----------------------Sample coordinates---------------------
X-coordinate : 9633.50
Y-coordinate : 10700.00
Z-coordinate : 26390.00
XX-coordinate : -119.40
ZZ-coordinate : 2265.36
--------------------------------------------------------------
TOMOGRAPHIC SCAN FINISHED at Fri Jun 21 11:52:18 2013
Lines used for estimate: 1, 51, 101, 151, 201, 251, 301, 351, 401, 451, 501, 551, 601, 651, 701, 751, 801, 851, 901, 951, 1001, 1051, 1101, 1151, 1201, 1251, 1301, 1351, 1401, 1451, 1501, 1551, 1601, 1651, 1701, 1751, 1801, 1851, 1901, 1951, 2001, 2051, 2101, 2151
Rotation center: 1277.99
Rotation center (non padded): 1277.99
------------------------------------------------------------
-------------- Projection information -------------------
Original tif projections 2560x2160 pixels
Total size of Tiff projections 21599012 kB
------------- Reconstruction parameters -----------------
Used algorithm GridRec
Reconstruction directory /sls/X02DA/data/e14251/Data10/disk1/spine_mf_10x_/rec_8bit/
Center 1280.50
Filter Parzen
Rotation 0.00
Geometry Homogeneous angular projection distribution between 0 and pi
Binning for reconstruction 1
Roi for reconstruction 0,0,0,0
Ring removal Off
Zinger removal Off
Output format TIFF8 Scaling parameters -0.0003325,0.0003269
Zero Padding 0.50
------------- Reconstruction Information -----------------
Size of reconstructed slice 6.55 MB
Size of reconstructed dataset 14.16 GB
In [ ]:
flat_files = glob.glob(os.path.join(working_dir,'flat1','*.tif'))
flat_files = sorted(flat_files)
print len(flat_files)
# pprint(flat_files)
dark_files = glob.glob(os.path.join(working_dir,'dark','*.tif'))
dark_files = sorted(dark_files)
print(len(dark_files))
# pprint(dark_files)
data_files = glob.glob(os.path.join(working_dir,'proj', '*.tif'))
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 [ ]:
mean_flat = None
for flat_file in flat_files:
tf = plt.imread(flat_files[0]).astype('float32')
if mean_flat is None:
mean_flat = np.zeros_like(tf)
mean_flat = mean_flat+tf
mean_flat = mean_flat/len(flat_files)
plt.figure(figsize=(15,10))
plt.imshow(mean_flat, cmap=plt.cm.viridis)
plt.colorbar()
plt.title('reference')
plt.show()
In [ ]:
mean_dark = None
for dark_file in dark_files:
td = plt.imread(dark_files[0]).astype('float32')
if mean_dark is None:
mean_dark = np.zeros_like(td)
mean_dark = mean_dark+td
mean_dark = mean_dark/len(dark_files)
plt.figure(figsize=(15,10))
plt.imshow(mean_dark, cmap=plt.cm.viridis)
plt.colorbar()
plt.title('dark')
plt.show()
In [ ]:
mean_flat = mean_flat-mean_dark
mean_flat[mean_flat<0]=0
plt.figure(figsize=(15,10))
plt.imshow(mean_flat, cmap=plt.cm.viridis)
plt.colorbar()
plt.title('reference')
plt.show()
In [ ]:
data = plt.imread(data_files[500]).astype('float32')
data = data-mean_dark
data = data/mean_flat
plt.figure(figsize=(15,10))
plt.imshow(data, cmap=plt.cm.viridis)
plt.colorbar()
plt.title('reference')
plt.show()
In [ ]:
import numba
fast_retrieve_phase=numba.jit()(tomopy.prep.phase.retrieve_phase)
In [ ]:
xt = np.expand_dims(data[::], axis=0)
xp = fast_retrieve_phase(xt[:,:], pixel_size=0.65e-4,
dist=21,
energy=24) # distances in [cm], energy [kev]
In [ ]:
plt.figure(figsize=(15,10))
plt.imshow(np.squeeze(xp), cmap=plt.cm.viridis)
plt.colorbar()
plt.title('reference')
plt.show()
In [ ]:
distance = 21
tmp_sino_file = '/home/makov/diskmnt/fast/inna/sino_spine_{}.h5'.format(distance)
In [ ]:
for distance in log_progress([21,]):
tmp_sino_file = '/home/makov/diskmnt/fast/inna/sino_spine_{}.h5'.format(distance)
if not os.path.exists(tmp_sino_file):
with h5py.File(tmp_sino_file,'w') as h5f:
for idf,f in enumerate(log_progress(data_files)):
d = plt.imread(f)
d = (d-mean_dark)/mean_flat
dt = np.expand_dims(d, axis=0)
# xp = fast_retrieve_phase(dt, pixel_size=0.65e-4, dist=distance, energy=24)
xp = dt
h5f.create_dataset(str(idf), data=np.squeeze(xp), compression="lzf")
In [ ]:
# if not os.path.exists(tmp_sino_file):
# with h5py.File(tmp_sino_file,'w') as h5f:
# for idf,f in enumerate(log_progress(data_files)):
# d = plt.imread(f)
# d = (d-mean_dark)/mean_flat
# dt = np.expand_dims(d[1800:1810], axis=0)
# xp = fast_retrieve_phase(dt, pixel_size=0.65e-4, dist=10, energy=24)
# # xp = dt
# h5f.create_dataset(str(idf), data=np.squeeze(xp), compression="lzf")
In [ ]:
sinogram_list = []
raw = 9
with h5py.File(tmp_sino_file) as h5f:
for idk, k in enumerate(log_progress(range(len(h5f.keys())))):
sinogram_list.append(h5f[str(idk)][raw])
sinogram = np.array(sinogram_list)
%xdel sinogram_list
In [ ]:
np.save('sinogram.pkl', sinogram)
In [ ]:
plt.figure(figsize=(15,10))
plt.imshow(sinogram, cmap=plt.cm.gray, vmin=0)
plt.axis('tight')
plt.colorbar()
plt.title('reference')
plt.show()
In [ ]:
plt.figure(figsize=(15,10))
plt.imshow(sinogram, cmap=plt.cm.gray)
plt.axis('tight')
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.gray,vmin=0)
plt.axis('tight')
plt.colorbar()
plt.title('reference')
plt.show()
In [ ]:
plt.figure(figsize=(15,10))
plt.subplot(2,2,1)
plt.imshow(sinogram, cmap=plt.cm.gray,vmin=0, interpolation='nearest')
plt.axis('tight')
plt.colorbar(orientation='horizontal')
plt.title('Without normalizaion')
plt.ylabel('Angle number')
plt.xlabel('Pixel number')
plt.subplot(2,2,2)
plt.imshow(sino_pp, cmap=plt.cm.gray,vmin=0, interpolation='nearest')
plt.axis('tight')
plt.colorbar(orientation='horizontal')
plt.title('With normalizaion')
plt.ylabel('Angle number')
plt.xlabel('Pixel number')
plt.subplot(212)
plt.plot(sinogram.sum(axis=-1)*sinogram.sum(axis=-1).mean(), 'k', label='Without normalizaion')
plt.plot(sino_pp.sum(axis=-1)*sinogram.sum(axis=-1).mean(), 'k--', label='With normalizaion')
plt.grid()
plt.ylabel('Radon invariant')
plt.xlabel('Angle number')
plt.legend(loc=0)
plt.show()
In [ ]:
plt.figure(figsize=(15,10))
plt.subplot(2,3,1)
plt.imshow(sinogram, cmap=plt.cm.gray,vmin=0, interpolation='nearest')
plt.axis('tight')
plt.colorbar(orientation='horizontal')
plt.title('Without normalizaion')
plt.ylabel('Angle number')
plt.xlabel('Pixel number')
plt.subplot(2,3,2)
plt.imshow(sino_pp, cmap=plt.cm.gray,vmin=0, interpolation='nearest')
plt.axis('tight')
plt.colorbar(orientation='horizontal')
plt.title('With normalizaion')
plt.ylabel('Angle number')
plt.xlabel('Pixel number')
plt.subplot(2,3,3)
plt.imshow(sino_pp-sinogram, cmap=plt.cm.gray,vmin=0, interpolation='nearest')
plt.axis('tight')
plt.colorbar(orientation='horizontal')
plt.title('Sinograms difference')
plt.ylabel('Angle number')
plt.xlabel('Pixel number')
plt.subplot(212)
plt.plot(sinogram.sum(axis=-1)*sinogram.sum(axis=-1).mean(), 'k', label='Without normalizaion')
plt.plot(sino_pp.sum(axis=-1)*sinogram.sum(axis=-1).mean(), 'k--', label='With normalizaion')
plt.grid()
plt.ylabel('Radon invariant')
plt.xlabel('Angle number')
plt.legend(loc=0)
plt.show()
In [ ]:
sino_pp = sinogram.T/sinogram.sum(axis=-1)*sinogram.sum(axis=-1).mean()
sino_pp = sino_pp.T
# sino_pp = sinogram.T-sinogram.min(axis=-1)
# sino_pp = sino_pp.T+0.01
# sino_pp = sino_pp - sino_pp.min()+0.01
sino_pp_log = -np.log(sino_pp)
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 = 1
s0 = sino_pp_log[0][delta:]
s1 = np.flipud(sino_pp_log[-1][delta:])
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(sino_pp_log.shape[0],dtype='float32')
angles = slices/len(slices-1)*180*np.pi/180.
In [ ]:
tomo_source = sino_pp_log.copy()
# tomo_source_t = tomo_source
tomo_source_t = tomo_source-tomo_source.min()
In [ ]:
plt.figure(figsize=(15,10))
plt.imshow(tomo_source_t, cmap=plt.cm.viridis)
plt.colorbar(orientation='horizontal')
plt.title('normalized')
plt.show()
In [ ]:
# trh =16
# tomo_source_t[tomo_source_t>trh]=trh
tomo_source_t = tomo_source_t/tomo_source_t.sum(axis=0)*np.sqrt(np.sum(tomo_source_t))
In [ ]:
plt.figure(figsize=(15,10))
plt.imshow(tomo_source_t, cmap=plt.cm.viridis)
plt.colorbar(orientation='horizontal')
plt.title('normalized')
plt.show()
In [ ]:
rec, proj_geom, cfg = astra_tomo2d(tomo_source_t, angles)
In [ ]:
rec=rec/np.mean(rec)
print rec.shape
In [ ]:
from scipy.ndimage import gaussian_filter
rs = np.array(rec.shape)//2
radius=1150
rec_t = rec[rs[0]-radius:rs[0]+radius,rs[1]-radius:rs[1]+radius]
mask = gaussian_filter(rec_t,50)
plt.figure(figsize=(15,15))
# plt.imshow(rec,vmin=0.1, vmax=0.2)
plt.imshow(rec_t-mask, cmap=plt.cm.gray,vmin=-0.15,vmax=0.25)
plt.colorbar(orientation='horizontal')
plt.title('rec')
plt.imsave('1809.tiff', rec_t-mask, cmap=plt.cm.gray,vmin=-0.15,vmax=0.25)
plt.show()
In [ ]:
# rs = rec.shape[0]//2
# radius=900
rc = tomopy.misc.corr.remove_ring(np.expand_dims(rec_t-mask, axis=0),
thresh=0.3, theta_min=90, rwidth=30)
rc = np.squeeze(rc)
plt.figure(figsize=(15,15))
# plt.imshow(rec,vmin=0.1, vmax=0.2)
plt.imshow(rc, cmap=plt.cm.gray,vmin=-0.15,vmax=0.25)
plt.colorbar(orientation='horizontal')
plt.imsave('1809_ring_postproc.tiff', rc, cmap=plt.cm.gray,vmin=-0.15,vmax=0.25)
plt.title('rec')
plt.show()
In [ ]:
from skimage.io import imread
In [ ]:
ref = imread(
'/home/makov/diskmnt/big/tomo_data/bukreeva/ESRF_osso_2016/spine_mf_10x_mf/tom_pulito_girato/Reslice of Reslice0181.tif')
In [ ]:
plt.figure(figsize=(15,10))
# plt.imshow((rc-rec)[1200:1800,1200:1800], cmap=plt.cm.viridis)
plt.imshow(ref, cmap=plt.cm.viridis)
# plt.imshow(rc[1000:2000,1000:2000])
plt.colorbar()
plt.title('ref')
plt.show()
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 [ ]:
sinogram = np.load('sinogram.pkl.npy')
In [ ]:
plt.figure(figsize=(10,10))
plt.imshow(sinogram, cmap=plt.cm.gray, interpolation='nearest')
plt.show()
In [ ]:
sino_log = -np.log(sinogram)
slices = np.arange(sino_log.shape[0],dtype='float32')
angles = slices/len(slices-1)*180*np.pi/180.
# TODO: fix rings
In [ ]:
plt.figure(figsize=(10,10))
plt.imshow(sino_log, cmap=plt.cm.gray, interpolation='nearest')
plt.colorbar(orientation='horizontal')
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
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, 1)
# 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 [ ]:
def astra_build_sinogram(rec, angles):
angles = angles.astype('float64') # hack for astra stability, may be removed in future releases
detector_size = rec.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)
proj_id = astra.create_projector('cuda',proj_geom,vol_geom)
sinogram_id, sinogram = astra.create_sino(rec, proj_id)
astra.data2d.delete(sinogram_id)
astra.clear()
return sinogram
In [ ]:
In [ ]:
padsize = 2000
sinogram_padded = np.zeros((sino_log.shape[0],sino_log.shape[1]+padsize*2), dtype='float32')
sinogram_padded[:,padsize:-padsize] = sino_log
rec, proj_geom, cfg = astra_tomo2d(sinogram_padded, angles)
sino = astra_build_sinogram(rec, angles)
sino[:,padsize:-padsize] = sino_log
MU = rec.sum()*1.5
X,Y = np.meshgrid(np.arange(rec.shape[0]),np.arange(rec.shape[1]))
X-=rec.shape[0]//2
Y-=rec.shape[1]//2
mask = (X**2+Y**2)<(rec.shape[0]//2)**2-10
for i in log_progress(range(10)):
rec, proj_geom, cfg = astra_tomo2d(sino, angles)
rec*=rec>0
rec*=mask
rec[rec>1] = 1
if rec.sum()>MU:
rec = rec/rec.sum()*MU
sino = astra_build_sinogram(rec, angles)
sino[:,padsize:-padsize] = sino_log
rec, proj_geom, cfg = astra_tomo2d(sino, angles)
plt.figure(figsize=(10,10))
plt.imshow(sino, cmap=plt.cm.gray, interpolation='nearest', vmin=0, vmax=1)
plt.colorbar(orientation='horizontal')
plt.show()
plt.figure(figsize=(10,10))
# plt.imshow(rec,vmin=0.1, vmax=0.2)
plt.imshow(rec, interpolation='nearest', cmap=plt.cm.gray)
plt.colorbar(orientation='horizontal')
plt.show()
In [ ]:
plt.figure(figsize=(10,10))
plt.imshow(rec[padsize:-padsize,padsize:-padsize], cmap=plt.cm.gray, interpolation='nearest')
plt.colorbar(orientation='horizontal')
plt.show()
In [ ]:
sino[:,padsize:-padsize] -= sino_log
In [ ]:
plt.figure(figsize=(10,10))
plt.imshow(sino, cmap=plt.cm.gray, interpolation='nearest')
plt.colorbar(orientation='horizontal')
plt.show()
In [ ]:
rec, proj_geom, cfg = astra_tomo2d(sino_log, angles)
plt.figure(figsize=(15,15))
# plt.imshow(rec,vmin=0.1, vmax=0.2)
plt.imshow(rec, interpolation='nearest', cmap=plt.cm.gray)
plt.colorbar(orientation='horizontal')
plt.show()
In [ ]:
rec.sum()
In [ ]:
sino_log.sum(axis=-1).mean()
In [ ]:
In [ ]: