In [ ]:
%matplotlib inline

Измерения кристалла кремния с одиночными дислокациями

  1. Трубка - Mo,режим 50 Кв, 30 мА
  2. Геометрия эксперимента:
    • Угол Брэгга - 10.66
    • Угол детектора - 21.32 (выстовлено приблизительно)
  3. Данные SI*.fit
    • 181 проекция (разбиение на 2 дня)
    • экспозиция - 10х20 сек. (складывались в матрице)
  4. Темновой ток FREE*.fit
    • 9 проекций (разбиение на 2 дня),
    • экспозиция - 10х20 сек. (складывались в матрице)

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 [ ]: