In [ ]:
%pylab inline
In [ ]:
import numpy as np
import EdfFile
import glob
import os
import h5py
import cv2
import cPickle as pickle
import pylab as plt
from pprint import pprint
In [ ]:
#data_fbp_root = './prepFBP_proj_topolamino_Cr7_220_transmission'
data_source_root = './proj_topolamino_Cr7_220_transmission/'
data_root = data_source_root
edf_files = glob.glob(data_root+'/*.edf')
out_hdf5 = data_root+'/hdf5/all.h5'
print len(edf_files)
print edf_files[0]
In [ ]:
def print_edf_header(edf_header):
print 'Number of images in EDF file {}'.format(edf_structure.GetNumImages())
print 'EDF header:'
for k,v in edf_header.items():
print '\t{:20}:\t{}'.format(k,v)
print edf_structure.GetData(0).dtype
return
def get_array_from_edf(edf_file_name):
edf_structure = EdfFile.EdfFile(edf_file_name,'r')
edf_header = edf_structure.GetHeader(0)
#print_edf_header(edf_header)
image_data = edf_structure.GetData(0).astype('float32')
return edf_header, image_data
In [ ]:
!mkdir -p {data_root+'/png'}
!mkdir -p {data_root+'/hdf5'}
In [ ]:
def get_motors_values(header):
motor_names = header['motor_mne'].split(' ')
motor_names = ['motor_'+ motor_name for motor_name in motor_names]
motor_values = header['motor_pos'].split(' ')
motor_values = [np.float32(f) for f in motor_values]
motor_header = dict(zip(motor_names, motor_values))
return motor_header
In [ ]:
#Generate HDF5 and PNGs
if not os.path.isfile(out_hdf5):
with h5py.File(out_hdf5,'w') as h5f:
for id_edf_file, edf_file in enumerate(edf_files):
header, res = get_array_from_edf(edf_file)
file_prefix = os.path.splitext(os.path.basename(edf_file))[0]
attributes = {}
attributes.update(header)
attributes.update(get_motors_values(header))
del attributes['motor_pos']
del attributes['motor_mne']
ds0 = h5f.create_dataset('/original/'+file_prefix,
data=res, compression="gzip",compression_opts=4)
for k,v in attributes.iteritems():
ds0.attrs[k] = v
small_image = cv2.resize(res,(res.shape[0]//3,res.shape[1]//3))
ds1 = h5f.create_dataset('/compressed_3/'+file_prefix,
data=small_image, compression="gzip",compression_opts=4)
plt.imsave(data_root+'/png/'+file_prefix+'.png', small_image, cmap = plt.cm.Greys_r)
for k,v in attributes.iteritems():
ds1.attrs[k] = v
if id_edf_file%100 == 0:
print
if id_edf_file%10 == 0:
print id_edf_file,
else:
print 'File {} already exists. Skip.'.format(out_hdf5)
In [ ]:
scan_angles = []
with h5py.File(out_hdf5 ,'r') as h5f:
h5g = h5f['compressed_3']
for k in h5g:
scan_angles.append(h5g[k].attrs['motor_pmo'])
scan_angles = np.array(scan_angles)
In [ ]:
plt.plot(scan_angles)
In [ ]:
from ipywidgets import interact
In [ ]:
def find_angles(h5g):
scan_angles = []
for k in h5g:
scan_angles.append((k.encode(), h5g[k].attrs['motor_pmo']))
scan_angles = dict(scan_angles)
return scan_angles
def get_key(d, angle):
for k,v in d.iteritems():
if np.isclose(v, angle):
return k
def browse_hdf5():
with h5py.File(out_hdf5,'r') as h5f:
h5g = h5f['compressed_3']
angles_dict = find_angles(h5g)
angles_list=sorted(np.array(angles_dict.values()))
def show_key(i):
angle = angles_list[i]
# print angle, get_key(angles_dict, angle)
plt.figure(figsize=(15,7))
plt.subplot(121)
plt.plot(sorted(angles_list))
plt.hold(True)
plt.plot(i, angle, 'r*')
plt.subplot(122)
image_key = filter(lambda item: np.isclose(item[1], angle), angles_dict.iteritems())[0][0]
with h5py.File(out_hdf5,'r') as h5f:
h5g = h5f['compressed_3']
image = h5g[image_key].value
#print image_key
image = image[100:-60, 100:-150] # Remove overlap regions
plt.imshow(image, cmap=plt.cm.gray)
plt.show()
interact(show_key, i=(0, len(angles_list)-1))
browse_hdf5()
In [ ]:
data = None
with h5py.File(out_hdf5,'r') as h5f:
name_ang = []
for k,v in h5f['compressed_3'].iteritems():
name_ang.append([k, v.attrs['motor_pmo']])
angles = [x[1] for x in name_ang]
keys = [x[0] for x in name_ang]
angles_indexes = np.argsort(angles)
for ai in angles_indexes:
frame = h5f['compressed_3'][keys[ai]].value
if data is None:
data = np.zeros(shape=(len(angles),frame.shape[0],frame.shape[1]), dtype='float32')
data[ai]=frame
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))
plt.colorbar(orientation='horizontal')
plt.show()
local_data = data3d
interact(show_slice, i=(0,data3d.shape[axis_numb]-1))
# data_filtered = np.array(data[240:488:5,120:-100, int(1006.5/3)-230:int(1006.5/3)+230])
data_filtered = np.array(data[::5,120:-100, int(1006.5/3)-230:int(1006.5/3)+230])
data_filtered[data_filtered<0] = 0
dfs = data_filtered.shape
zp = 200
data_filtered_zeros = np.zeros(shape=(dfs[0], dfs[1],dfs[2]+zp))
data_filtered_zeros[:,:,zp/2+30:dfs[2]+zp/2]= np.log(1+data_filtered[:,:,30:])-4.0
data_filtered_zeros[data_filtered_zeros<0] = 0
# angles_filtered = np.array(angles[240:488:5])
angles_filtered = np.array(angles[::5])
show_array(data_filtered_zeros, 0)
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)
# print proj_geom
# print sinogram.shape
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 = 4.637*np.pi/180
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)
# print proj_geom
# print sinogram.shape
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
In [ ]:
s1 = np.require(data_filtered_zeros[:,:,:].swapaxes(0,1),dtype=np.float32, requirements=['C'])
In [ ]:
rec_3d = astra_topotomo3d(s1, angles_filtered*np.pi/180)
In [ ]:
show_array(rec_3d,0)
In [ ]:
data_file_parts = os.path.join(data_root,'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, 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 [ ]:
image_key = 'TL_CR7_pos1_int1_0466'
with h5py.File(out_hdf5,'r') as h5f:
h5g = h5f['compressed_3']
image = h5g[image_key].value
print image_key
image = image[100:-60, 100:-150]
plt.figure(figsize=(10,10))
plt.imshow(image, cmap=plt.cm.gray);
plt.colorbar()
In [ ]:
def find_angles(h5g):
scan_angles = []
for k in h5g:
scan_angles.append((k.encode(), h5g[k].attrs['motor_pmo']))
scan_angles = dict(scan_angles)
return scan_angles
def get_data_by_angle(h5g, angle):
for k,v in h5g.iteritems():
if np.isclose(angle,v.attrs['motor_pmo']):
return v.value[100:-60, 100:-150]
with h5py.File(out_hdf5,'r') as h5f:
h5g = h5f['compressed_3']
angles_dict = find_angles(h5g)
angles_list=sorted(np.array(angles_dict.values()))
angles_list = filter(lambda x: x >-20 and x<40, angles_list)
angles_list = np.array(angles_list)
In [ ]:
with h5py.File(out_hdf5,'r') as h5f:
h5g = h5f['compressed_3']
d = get_data_by_angle(h5g, angles_list[0])
d.shape
In [ ]:
import astra
In [ ]:
# vol_geom = astra.create_vol_geom(d.shape[1], d.shape[0], d.shape[0])
In [ ]:
tomo_angles = angles_list[::10]
len(tomo_angles)
sinogram_volume = np.zeros(shape=(d.shape[0],len(tomo_angles),d.shape[1]), dtype='float32')
print sinogram_volume.shape
with h5py.File(out_hdf5,'r') as h5f:
h5g = h5f['compressed_3']
for iang,ang in enumerate(tomo_angles):
raw_sinogram = get_data_by_angle(h5g, ang)
raw_sinogram = (raw_sinogram - raw_sinogram.min())/(raw_sinogram.max() - raw_sinogram.min())
sinogram_volume[:,iang,:] = raw_sinogram
In [ ]:
def show_sinogram(i):
plt.imshow(sinogram_volume[i,:,:], cmap=plt.cm.gray)
plt.axis('tight')
plt.colorbar()
interact(show_sinogram, i=(0, sinogram_volume.shape[0]-1))
In [ ]:
# We generate the same geometry as the circular one above.
astra_angles = tomo_angles/180*np.pi
vectors = np.zeros((len(astra_angles), 12))
for i in range(len(astra_angles)):
# ray direction
vectors[i,0] = np.sin(astra_angles[i])
vectors[i,1] = -np.cos(astra_angles[i])
vectors[i,2] = 0
# 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', sinogram_volume.shape[0], sinogram_volume.shape[-1], vectors)
In [ ]:
# Create a data object for the reconstruction
sinogram_id = astra.data3d.create('-sino', proj_geom, data=sinogram_volume)
vol_geom = astra.create_vol_geom(d.shape[1], d.shape[1], d.shape[0])
rec_id = astra.data3d.create('-vol', vol_geom)
In [ ]:
# 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
# Create the algorithm object from the configuration structure
alg_id = astra.algorithm.create(cfg)
# Run 150 iterations of the algorithm
# Note that this requires about 750MB of GPU memory, and has a runtime
# in the order of 10 seconds.
astra.algorithm.run(alg_id, 150)
# Get the result
rec = astra.data3d.get(rec_id)
astra.algorithm.info()
astra.data3d.info()
astra.projector.info()
astra.matrix.info()
astra.data3d.clear()
In [ ]:
pylab.figure(figsize=(5,5))
def show_rec(i):
plt.imshow(rec[i,:,:], cmap=plt.cm.gray)
plt.colorbar()
interact(show_rec, i=(0, rec.shape[0]-1))
In [ ]:
# print tmp_sino.shape
print vol_geom
print proj_geom
print sinogram_volume.shape
In [ ]:
vol_geom = astra.create_vol_geom(128, 128, 128)
angles = np.linspace(0, np.pi, 180,False)
proj_geom = astra.create_proj_geom('parallel3d', 1.0, 1.0, 128, 192, angles)
# Create a simple hollow cube phantom
cube = np.zeros((128,128,128))
cube[17:113,17:113,17:113] = 1
cube[33:97,33:97,33:97] = 0
# Create projection data from this
proj_id, proj_data = astra.create_sino3d_gpu(cube, proj_geom, vol_geom)
sinogram_id = astra.data3d.create('-sino', proj_geom, data=proj_data)
In [ ]:
plt.imshow(proj_data[64,:,:], cmap=plt.cm.gray)
plt.colorbar()
In [ ]:
proj_data.shape
In [ ]:
astra.algorithm.info()
astra.data3d.info()
astra.projector.info()
astra.matrix.info()
# astra.data3d.clear()
In [ ]:
1006.5/3
In [ ]: