In [ ]:
%matplotlib inline
In [ ]:
import os
import glob
import numpy as np
import pylab as plt
import h5py
import astra
import cv2
from astropy.io import fits
from pprint import pprint
import warnings
import ipywidgets
In [ ]:
import re
def natural_key(string_):
"""See http://www.codinghorror.com/blog/archives/001018.html"""
return [int(s) if s.isdigit() else s for s in re.split(r'(\d+)', string_)]
In [ ]:
#lodad files and sort it ina natural way
data_dir = os.path.expanduser('/home/makov/diskmnt/big/topo_tomo/2015_11_09/Si/exp')
data_files = glob.glob(os.path.join(data_dir,'SI*.fit'))
data_files = sorted(data_files, key=natural_key)
print len(data_files)
empty_files = glob.glob(os.path.join(data_dir,'FREE*.fit'))
empty_files = sorted(empty_files, key=natural_key)
print len(empty_files)
In [ ]:
#loading data files
data = None
for idf, df in enumerate(data_files):
with warnings.catch_warnings():
warnings.simplefilter("ignore")
with fits.open(df) as fit:
tmp_data = fit[0].data.astype('uint16')
if data is None:
data = np.ndarray(shape = (len(data_files),tmp_data.shape[0],tmp_data.shape[1]),
dtype = 'float32')
data[idf] = tmp_data
In [ ]:
#loading empty files
empty = None
for idf, df in enumerate(empty_files):
with warnings.catch_warnings():
warnings.simplefilter("ignore")
with fits.open(df) as fit:
tmp_data = fit[0].data.astype('uint16')
if empty is None:
empty = np.ndarray(shape = (len(empty_files),tmp_data.shape[0],tmp_data.shape[1]),
dtype = 'float32')
empty[idf] = tmp_data
In [ ]:
angles = np.arange(len(data_files))*2.0 # Angles staep = 2 deg
In [ ]:
# remove unused regions (crop ROI)
data = data[:,320:620,103:403]
empty = empty[:,320:620,103:403]
In [ ]:
empty_frame = empty.sum(axis=0)/empty.shape[0] # avearaging emplty frames
In [ ]:
data = data-empty_frame
In [ ]:
data[data<0]=0
In [ ]:
# #take only images with good contrast
# good_slices=np.concatenate((np.arange(0,25),np.arange(40,90),np.arange(130,180)))
# data = data[good_slices]
# angles = angles[good_slices]
In [ ]:
def show_array(data3d, axis_numb=0):
def show_slice(i):
plt.figure(figsize=(10,10))
plt.imshow(local_data.take(i,axis=axis_numb),vmin=0)
plt.colorbar(orientation='horizontal')
plt.show()
local_data = data3d
ipywidgets.interact(show_slice, i=(0,data3d.shape[axis_numb]-1))
show_array(data, 1)
In [ ]:
import astra
def astra_tomo2d_parallel(sinogram, angles):
angles = angles.astype('float64')
detector_size = sinogram.shape[1]
rec_size = detector_size
vol_geom = astra.create_vol_geom(rec_size, rec_size)
proj_geom = astra.create_proj_geom('parallel', 1.0, 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)
# proj_id = astra.create_projector('strip', proj_geom, vol_geom) # for CPU reconstruction only
# 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['ProjectorId'] = proj_id # for CPU reconstruction only
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, 150)
# 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_tomo3d_parallel(sinogram, angles):
angles = angles.astype('float64')
detector_size = sinogram.shape[-1]
slices_number = sinogram.shape[0]
rec_size = detector_size
vol_geom = astra.create_vol_geom(rec_size, rec_size, slices_number)
proj_geom = astra.create_proj_geom('parallel3d', 1.0, 1.0, slices_number, detector_size, angles)
sinogram_id = astra.data3d.create('-sino', proj_geom, data=sinogram)
# Create a data object for the reconstruction
rec_id = astra.data3d.create('-vol', vol_geom)
# proj_id = astra.create_projector('strip', proj_geom, vol_geom) # for CPU reconstruction only
# Set up the parameters for a reconstruction algorithm using the GPU
cfg = astra.astra_dict('SIRT3D_CUDA')
cfg['ReconstructionDataId'] = rec_id
cfg['ProjectionDataId'] = sinogram_id
# cfg['ProjectorId'] = proj_id # for CPU reconstruction only
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, 150)
# Get the result
rec = astra.data3d.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.data3d.delete(rec_id)
astra.data3d.delete(sinogram_id)
astra.clear()
return rec
def astra_topotomo3d(sinogram, angles):
astra_angles = angles.astype('float64')
detector_size = sinogram.shape[-1]
slices_number = sinogram.shape[0]
rec_size = detector_size
vol_geom = astra.create_vol_geom(rec_size, rec_size, slices_number)
# We generate the same geometry as the circular one above.
vectors = np.zeros((len(astra_angles), 12))
alpha = -10.66*np.pi/180 # define bragg angle
for i in range(len(astra_angles)):
# ray direction
vectors[i,0] = np.sin(astra_angles[i])*np.cos(alpha)
vectors[i,1] = -np.cos(astra_angles[i])*np.cos(alpha)
vectors[i,2] = np.sin(alpha)
# center of detector
vectors[i,3:6] = 0
# vector from detector pixel (0,0) to (0,1)
vectors[i,6] = np.cos(astra_angles[i])
vectors[i,7] = np.sin(astra_angles[i])
vectors[i,8] = 0;
# vector from detector pixel (0,0) to (1,0)
vectors[i,9] = 0
vectors[i,10] = 0
vectors[i,11] = 1
# Parameters: #rows, #columns, vectors
proj_geom = astra.create_proj_geom('parallel3d_vec', slices_number, detector_size, vectors)
# proj_geom = astra.create_proj_geom('parallel3d', 1.0, 1.0, slices_number, detector_size, angles)
sinogram_id = astra.data3d.create('-sino', proj_geom, data=sinogram)
# Create a data object for the reconstruction
rec_id = astra.data3d.create('-vol', vol_geom)
# proj_id = astra.create_projector('strip', proj_geom, vol_geom) # for CPU reconstruction only
# Set up the parameters for a reconstruction algorithm using the GPU
cfg = astra.astra_dict('SIRT3D_CUDA')
cfg['ReconstructionDataId'] = rec_id
cfg['ProjectionDataId'] = sinogram_id
# cfg['ProjectorId'] = proj_id # for CPU reconstruction only
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, 100)
# Get the result
rec = astra.data3d.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.data3d.delete(rec_id)
astra.data3d.delete(sinogram_id)
astra.clear()
return rec
In [ ]:
# rec_slice = astra_tomo2d_parallel(np.log(data[:,150,:]+1), angles*np.pi/180)
# rec_slice = astra_tomo2d_parallel(data[:,150,:], angles*np.pi/180)
In [ ]:
# plt.imshow(rec_slice)
# plt.colorbar()
In [ ]:
s1 = np.require(data[:,:,:].swapaxes(0,1),dtype=np.float32, requirements=['C'])
In [ ]:
# rec_slice = astra_tomo3d_parallel(np.log(s1+1), angles*np.pi/180)
In [ ]:
# show_array(rec_slice)
In [ ]:
rec_3d = astra_topotomo3d(s1, angles*np.pi/180)
In [ ]:
rec_3d_masked = rec_3d
rec_3d_masked[rec_3d_masked<100]=0
show_array(rec_3d_masked,2)
In [ ]:
sl = rec_3d[rec_3d.shape[0]/3]
mask = np.zeros_like(sl)
X,Y = meshgrid(np.arange(mask.shape[0]),np.arange(mask.shape[1]))
X = X-mask.shape[0]/2
Y = Y-mask.shape[1]/2
mask = (X*X+Y*Y)**0.5<(mask.shape[0]/2)-3
rec_3d_masked = rec_3d*mask
In [ ]:
show_array(rec_3d_masked)
In [ ]:
data_file_parts = os.path.join(data_dir,'result.h5')
out_file = '{}_reconstruction.h5'.format(''.join(data_file_parts[:-1]))
print out_file
with h5py.File(out_file,'w') as h5f:
h5f.create_dataset('Result', data=rec_3d_masked, dtype=np.float32)
In [ ]:
def save_amira(result_file):
"""
Функция сохраняет реконструированные слои в формате Amira raw file
Inputs:
data_path - путь к директории, где находиться файл res_tomo.hdf5 в формате HDF5
в этом файде должен быть раздел (node) /Results в котором в виде 3D массива
записаны реконструированный объём
Outputs:
Файлы amira.raw и tomo.hx. Если файлы уже существуют, то они перезаписываются.
Тип данных: float32 little endian
"""
data_path = os.path.dirname(result_file)
with open(os.path.join(data_path, 'amira.raw'), 'wb') as amira_file:
with h5py.File(result_file, 'r') as h5f:
x = h5f['Result']
for i in range(x.shape[0]):
numpy.array(x[i, :, :]).tofile(amira_file)
file_shape = x.shape
with open(os.path.join(data_path, 'tomo.hx'), 'w') as af:
af.write('# Amira Script\n')
af.write('remove -all\n')
af.write(r'[ load -raw ${SCRIPTDIR}/amira.raw little xfastest float 1 '+
str(file_shape[1])+' '+str(file_shape[2])+' '+str(file_shape[0])+
' 0 '+str(file_shape[1]-1)+' 0 '+str(file_shape[2]-1)+' 0 '+str(file_shape[0]-1)+
' ] setLabel tomo.raw\n')
In [ ]:
save_amira(out_file)
In [ ]:
%matplotlib?
In [ ]:
ls /home/makov/diskmnt/big/topo_tomo/2015_11_09/Si/exp
In [ ]: