In [ ]:
%matplotlib inline
In [ ]:
from __future__ import division
import os
import glob
try:
import ConfigParser
except ImportError:
import configparser as ConfigParser
import numpy as np
import h5py
import pylab as plt
import cv2
import astra
from ipywidgets import interact
from pprint import pprint
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 = unicode(index or '?')
In [ ]:
# Data directory
# data_root = '/diskmnt/a/makov/yaivan/MMC_1/'
data_root = '/media/makov/buext4/yaivan/MMC_1'
raw_root = os.path.join(data_root,'Raw')
out_file = os.path.join(data_root, 'raw.h5')
In [ ]:
# files = !ls {raw_root}
# pprint(sorted(files))
In [ ]:
# SkyScan config output
def print_config(config):
for k,v in config._sections.items():
print('[{}]'.format(k))
for kk, vv in v.items():
print('\t{} = {}'.format(kk,vv))
config = ConfigParser.RawConfigParser()
config.optionxform = str
config.read(os.path.join(raw_root,'MMC1_2.82um_.log'))
print_config(config)
In [ ]:
object_name = config._sections['Acquisition']['Filename Prefix']
print('Object name:', object_name)
In [ ]:
data_files = glob.glob(os.path.join(raw_root,object_name+'[0-9]'*4+'.tif'))
data_files = sorted(data_files)
print('Data files found:', len(data_files))
print('First file:', data_files[0])
print('Last file:', data_files[-1])
In [ ]:
test_data = plt.imread(data_files[0])
print('Size of image:', test_data.shape)
print('Image data type:', test_data.dtype)
In [ ]:
# if not packed raw HDF5 file found, create it
#TODO: add adittional reference iif files
if not os.path.exists(out_file):
with h5py.File(out_file,'w-') as h5f:
h5f.create_dataset('data', (len(data_files), test_data.shape[0], test_data.shape[1]),
chunks=True, dtype=test_data.dtype)
for idf, df in log_progress(enumerate(data_files)):
df = plt.imread(data_files[idf])
h5f['data'][idf]=df
In [ ]:
def build_reconstruction_geomety(detector_size, angles):
# proj_geom = astra.create_proj_geom('parallel', 1.0, detector_size, angles)
#Object to Source (mm) = 56.135
#Camera to Source (mm) = 225.082
# All distances in [pixels]
pixel_size = 2.82473e-3
os_distance = 56.135/pixel_size
ds_distance = 225.082/pixel_size
proj_geom = astra.create_proj_geom('fanflat', ds_distance/os_distance, detector_size, angles,
os_distance, (ds_distance-os_distance))
return proj_geom
In [ ]:
def astra_tomo2d_fanflat_fbp(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'] = 5
# 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
def astra_tomo2d_fanflat_sirt(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('SIRT_CUDA')
cfg['ReconstructionDataId'] = rec_id
cfg['ProjectionDataId'] = sinogram_id
cfg['option'] = {}
# cfg['option']['MinConstraint'] = 0
# cfg['option']['MaxConstraint'] = 5
# 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, 2000)
# 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
# Define the plugin class (has to subclass astra.plugin.base)
# Note that usually, these will be defined in a separate package/module
class SIRTPlugin(astra.plugin.base):
"""Example of an ASTRA plugin class, implementing a simple 2D SIRT algorithm.
Options:
'rel_factor': relaxation factor (optional)
"""
# The astra_name variable defines the name to use to
# call the plugin from ASTRA
astra_name = "SIRT-PLUGIN"
def initialize(self,cfg, rel_factor = 1):
self.W = astra.OpTomo(cfg['ProjectorId'])
self.vid = cfg['ReconstructionDataId']
self.sid = cfg['ProjectionDataId']
self.rel = rel_factor
def run(self, its):
v = astra.data2d.get_shared(self.vid)
s = astra.data2d.get_shared(self.sid)
print(s.shape)
W = self.W
for i in range(its):
v[:] += self.rel*(W.T*(s - (W*v).reshape(s.shape))).reshape(v.shape)/s.size
# from plugin import SIRTPlugin
def astra_tomo2d_fanflat_plugin(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)
proj_id = astra.create_projector('cuda',proj_geom,vol_geom)
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)
astra.plugin.register(SIRTPlugin)
print(astra.plugin.get_registered())
# Set up the parameters for a reconstruction algorithm using the GPU
cfg = astra.astra_dict('SIRT-PLUGIN')
cfg['ProjectorId'] = proj_id
cfg['ReconstructionDataId'] = rec_id
cfg['ProjectionDataId'] = sinogram_id
cfg['option'] = {}
cfg['option']['rel_factor'] = 1.5
# cfg['option']['MinConstraint'] = 0
# cfg['option']['MaxConstraint'] = 5
# 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, 10)
# 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
def create_sinogram(data, angles):
angles = angles.astype('float64') # hack for astra stability, may be removed in future releases
detector_size = data.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)
W = astra.OpTomo(proj_id)
P = data
sinogram = W * P
sinogram = sinogram.reshape([len(angles), detector_size])
return np.rot90(sinogram,3)
In [ ]:
try:
from functools32 import lru_cache
except ImportError:
from functools import lru_cache
def cv_rotate(x,angle):
"""
Rotate square array using OpenCV2 around center of the array
:param x: numpy array
:param angle: angle in degrees
:return: rotated array
"""
x_center=tuple(np.array((x.shape[1],x.shape[0]),dtype='float32')/2.0-0.5)
rot_mat=cv2.getRotationMatrix2D(x_center,angle,1.0)
xro=cv2.warpAffine(x,rot_mat,(x.shape[1],x.shape[0]),flags=cv2.INTER_LINEAR)
return xro
@lru_cache(maxsize=32)
def load_sinogram(sinogram_number):
# Without thermal correction
with h5py.File(out_file,'r') as h5f:
sinogram = h5f['data'][:,-sinogram_number,:]
# With thermal correction
# Loading geometrical correction map generated by NRecon
# corrections_map = np.loadtxt(os.path.join(raw_root,'MMC1_2.82um__TS.crv'),skiprows=2)
# with h5py.File(out_file,'r') as h5f:
# sinogram = np.zeros_like(h5f['data'][:,-sinogram_number,:])
# for ii in tqdm.tqdm(range(sinogram.shape[0])):
# sinogram[ii] = np.roll(
# h5f['data'][ii,-sinogram_number-int(corrections_map[ii,2]),:],
# -int(corrections_map[ii,1])
# )
# select only central 180 deg for reconstruction
# sinogram = sinogram[int(sinogram.shape[0]/2-900):1800-int(sinogram.shape[0]/2-900)]
# sinogram = sinogram[:1800]
sinogram = np.array(sinogram/2.**16)
return sinogram
In [ ]:
def get_sinogram(sinogram_number):
sinogram = load_sinogram(sinogram_number)
# sinogram_max = np.mean(sinogram, axis=1)
# for ii in range(sinogram.shape[0]):
# sinogram[ii] *= sinogram_max[ii]/sinogram_max.mean()
# plt.plot(np.max(sinogram, axis=1))
# plt.show()
# np.savetxt('sinogram.txt', sinogram)
# angles_reduce = 1
# size_reduce = 1
# sinogram = cv2.resize(sinogram,(sinogram.shape[1]//size_reduce, sinogram.shape[0]//angles_reduce))
# print sinogram.min(),sinogram.max()
# sinogram[sinogram<1e-3] = 1
sinogram = np.rot90(sinogram)
sinogram = -np.log(sinogram)
# Corrections here !!!
sinoram_max = sinogram.max()
# for i in range(sinogram.shape[0]):
# sinogram[i] = sinogram[i] - sinogram[i].min()
sinogram = sinogram/sinoram_max
return sinogram
def get_reconstruction(sinogram, reconstruction_function):
angles = np.arange(sinogram.shape[1])*0.1-11.493867*2
angles = angles.astype('float64')/180.*np.pi
# angles = angles + (angles.max()+angles.min())/2.
# angles = angles - np.pi
astra_data = np.rot90(sinogram)
astra_rec = reconstruction_function(astra_data, angles)
astra_rec = np.rot90(astra_rec)
return astra_rec
def get_reconstruction_fbp(sinogram):
return get_reconstruction(sinogram, astra_tomo2d_fanflat_fbp)
def get_reconstruction_sirt(sinogram):
return get_reconstruction(sinogram, astra_tomo2d_fanflat_sirt)
def get_reconstruction_plugin(sinogram):
return get_reconstruction(sinogram, astra_tomo2d_fanflat_plugin)
astra_tomo2d_fanflat_plugin
def show_sinogram(sinogram_number):
sinogram = get_sinogram(sinogram_number)
plt.figure(figsize=(7,7))
plt.imshow(sinogram, cmap=plt.cm.gray_r)
plt.colorbar()
plt.axis('tight')
plt.title('Original')
plt.show()
def show_reconstruction(rec, roi, title=None):
# plt.figure(figsize=(12,10))
plt.imshow(rec[roi], cmap=plt.cm.gray, interpolation='nearest')
# plt.colorbar(orientation='vertical')
if title is None:
plt.title('Astra')
else:
plt.title('Astra {}'.format(title))
# plt.show()
def show_original(sinogram_number, roi):
# plt.figure(figsize=(12,10))
source_rec = plt.imread(os.path.join(data_root,
'Reconstructed','{}_rec{:04d}.png'.format(object_name, sinogram_number)))[...,0]
source_rec = np.rot90(source_rec)
source_rec = np.array(source_rec)*(0.52+0.18)-0.18
plt.imshow(source_rec[roi], cmap=plt.cm.gray, interpolation='nearest')
# plt.colorbar(orientation='vertical')
plt.title('Skyscan')
# plt.show()
In [ ]:
N=960
sinogram_orig = get_sinogram(N)
rec_fbp = get_reconstruction_fbp(sinogram_orig)
angles = np.arange(sinogram_orig.shape[1])*0.1+90-11.493867*2
sinogram_fbp = create_sinogram(rec_fbp, angles*np.pi/180)
# # rec_plugin = get_reconstruction_plugin(sinogram)
# rec_sirt = get_reconstruction_sirt(sinogram)
rois = []
# rois.append(np.ix_(np.r_[1900:2400],np.r_[2000:2500]))
# rois.append(np.ix_(np.r_[1100:1700],np.r_[2500:3200]))
# rois.append(np.ix_(np.r_[2660:2800],np.r_[2200:2400]))
# rois.append(np.ix_(np.r_[2750:2900],np.r_[2600:2800]))
# rois.append(np.ix_(np.r_[3000:3500],np.r_[2000:3000]))
# rois.append(np.ix_(np.r_[2700:3500],np.r_[1300:2100]))
# rois.append(np.ix_(np.r_[1000:1600],np.r_[1600:2200]))
rois.append(np.ix_(np.r_[0:4000],np.r_[0:4000]))
# for roi in rois:
# plt.figure(figsize = (20,20))
# plt.subplot(121)
# show_original(N ,roi)
# plt.subplot(122)
# show_reconstruction(rec_fbp, roi, 'FBP')
# # plt.subplot(224)
# # show_reconstruction(rec_sirt, roi, 'SIRT')
# # plt.show()
# # show_reconstruction(rec_plugin, roi, 'SIRT-PLUGIN')
# # show_reconstruction(rec_sirt, roi, 'SIRT')
for roi in rois:
plt.figure(figsize = (10,10))
show_reconstruction(rec_fbp, roi, 'FBP')
plt.savefig('fbp_no_corrections.png')
# show_sinogram(N)
# plt.figure(figsize=(7,7))
# plt.imshow(sinogram_fbp, cmap=plt.cm.gray_r)
# plt.colorbar()
# plt.axis('tight')
# plt.title('FBP')
# plt.show()
# save_marina_txt(sinogram_orig,'sinogram_original_no_corrections')
# save_marina_txt(sinogram_fbp, 'sinogram_fbp_no_corrections')
# save_marina_txt(rec_fbp, 'rec_fbp_no_corrections')
In [ ]:
# gamma corractions
from skimage.exposure import adjust_gamma
N=960
for gamma in np.arange(0.5,2,0.1):
sinogram_orig = get_sinogram(N)
sinogram_orig = adjust_gamma(sinogram_orig, gamma)
rec_fbp = get_reconstruction_fbp(sinogram_orig)
angles = np.arange(sinogram_orig.shape[1])*0.1+90-11.493867*2
sinogram_fbp = create_sinogram(rec_fbp, angles*np.pi/180)
# # rec_plugin = get_reconstruction_plugin(sinogram)
# rec_sirt = get_reconstruction_sirt(sinogram)
rois = []
# rois.append(np.ix_(np.r_[1900:2400],np.r_[2000:2500]))
# rois.append(np.ix_(np.r_[1100:1700],np.r_[2500:3200]))
# rois.append(np.ix_(np.r_[2660:2800],np.r_[2200:2400]))
# rois.append(np.ix_(np.r_[2750:2900],np.r_[2600:2800]))
# rois.append(np.ix_(np.r_[3000:3500],np.r_[2000:3000]))
# rois.append(np.ix_(np.r_[2700:3500],np.r_[1300:2100]))
# rois.append(np.ix_(np.r_[1000:1600],np.r_[1600:2200]))
rois.append(np.ix_(np.r_[0:4000],np.r_[0:4000]))
for roi in rois:
plt.figure(figsize = (20,20))
plt.subplot(121)
show_original(N ,roi)
plt.subplot(122)
show_reconstruction(rec_fbp, roi, 'FBP')
# # plt.subplot(224)
# # show_reconstruction(rec_sirt, roi, 'SIRT')
# # plt.show()
# # show_reconstruction(rec_plugin, roi, 'SIRT-PLUGIN')
# # show_reconstruction(rec_sirt, roi, 'SIRT')
show_sinogram(N)
plt.figure(figsize=(7,7))
plt.imshow(sinogram_fbp, cmap=plt.cm.gray_r)
plt.colorbar()
plt.axis('tight')
plt.title('FBP')
plt.show()
save_marina_txt(sinogram_orig,'sinogram_original_gamma_{}'.format(gamma))
save_marina_txt(sinogram_fbp, 'sinogram_fbp_gamma_{}'.format(gamma))
save_marina_txt(rec_fbp, 'rec_fbp_gamma_{}'.format(gamma))
In [ ]:
def save_marina_txt(data, name):
file_name = '{}_{}_{}.txt'.format(name, data.shape[0], data.shape[1])
print name, data.shape, data.dtype, file_name
np.savetxt(file_name, data.reshape(np.prod(data.shape),1), fmt='%0.3e')
In [ ]:
np.save('rec_sirt',rec_sirt)
np.save('rec_fbp',rec_fbp)
source_rec = plt.imread(os.path.join(data_root,
'Reconstructed','{}_rec{:04d}.png'.format(object_name, N)))[...,0]
np.save('rec_orig',source_rec)
In [ ]:
rec_sirt = np.load('rec_sirt.npy')
rec_fbp = np.load('rec_fbp.npy')
In [ ]:
for name in ['rec_sirt', 'rec_fbp', 'rec_orig']:
x = np.load(name+'.npy')
print name, x.shape, x.dtype
np.savetxt('{}_{}_{}.txt'.format(name, x.shape[0], x.shape[1]), x.reshape(np.prod(x.shape),1), fmt='%0.3e')
In [ ]:
!head rec_sirt_4000_4000.txt
In [ ]:
!gzip *.txt rec.zip
In [ ]:
save_marina_txt(sinogram_orig,'sinogram_original_normalized')
In [ ]:
import tomopy
In [ ]:
data_orig = np.load('rec_sirt.npy')
print data_orig.shape
In [ ]:
data = np.zeros(shape=(1, data_orig.shape[0], data_orig.shape[1]))
data[0] = data_orig*1e3
In [ ]:
data_corr = tomopy.misc.corr.remove_ring(data)
In [ ]:
plt.imshow(data_corr[0,1900:2100,1900:2100])
plt.colorbar()
In [ ]:
plt.imshow(data[0,1900:2100,1900:2100])
plt.colorbar()
In [ ]:
save_marina_txt(data[0],'rec_sirt_norings')
In [ ]:
!gzip rec_sirt_norings_4000_4000.txt
In [ ]:
N=980
sinogram_orig = get_sinogram(N)
rec_fbp = get_reconstruction_fbp(sinogram_orig)
In [ ]:
save_marina_txt(rec_fbp,'rec_fbp_980')
In [ ]:
!gzip rec_fbp_980_4000_4000.txt
In [ ]: