4D reconstruction LavaLamp

from kblt import recon
import tomopy
#import tomopy.prep.alignment as alignment
import tomopy.prep.stripe
from matplotlib.pyplot import imshow
import matplotlib.pyplot as plt
#import matplotlib.pyplot
import numpy as np
from numpy import inf
from qimtools import visualization, inspection, segment
import os
import shutil, os
from skimage.io.collection import ImageCollection
import time
from PIL import Image
import PIL
from lth import image_processing as imp
import scipy.ndimage as ndi
from skimage.exposure import rescale_intensity
from skimage.util import img_as_uint, img_as_ubyte, img_as_float32

astropy module not found

#sample name
sample = 'LavaLamp/'

#input folder
fname = sample + 'rep_16/'
fname_flat = sample
print(fname)
print(fname_flat)

LavaLamp/rep_16/
LavaLamp/

tomo_start = 1
tomo_end = 100
flat_start = 1
flat_end = 20

ind_flat = range(flat_start, flat_end)
ind_tomo = range(tomo_start, tomo_end)

#Read RGB data
#to do, check min and max values, and let user define them!
tomo = recon.read_rgb_convert(fname, 'tomo', save_out = False)
flat = recon.read_rgb_convert(fname_flat, 'flat', save_out = False)

#Normalized the projection images
data = recon.flat_fielding(fname, tomo, flat, save_out = False)
data[np.isnan(data)] = 0
data[data == -inf] = 0
data[data == inf] = 1

Starting flat-fielding of data.
Starting removal of NaN, Inf and -Inf values.
Removal of NaN, Inf and -Inf took: 0.13 seconds
Flat-fielding took: 0.35 seconds

#jupyterlab
#%matplotlib inline
#jupyter notebook
%matplotlib notebook

#plot figures
prj = plt.figure()
prj.subplots_adjust(hspace=.5)
ax = prj.add_subplot(1, 3, 1)
ax.set_title('projection')
ax.imshow(tomo[1])
ax2 = prj.add_subplot(1, 3, 2, sharex=ax, sharey=ax)
ax2.set_title('flat-field')
ax2.imshow(flat[0])
ax3 = prj.add_subplot(1, 3, 3, sharex=ax, sharey=ax)
ax3.set_title('normalized projection')
ax3.imshow(data[0])
plt.show()

#OPTIONAL
#select cropped region
y1 = 0
y2 = 480
x1 = 75
x2 = 575
%matplotlib notebook

#plot figures
prj = plt.figure()
prj.subplots_adjust(hspace=.5)

ax = prj.add_subplot(1, 3, 1)
ax.set_title('angle 0')
ax.imshow(data[0, y1:y2, x1:x2])

ax2 = prj.add_subplot(1, 3, 2, sharex=ax, sharey=ax)
ax2.set_title('angle 90')
ax2.imshow(data[50, y1:y2, x1:x2])

ax3 = prj.add_subplot(1, 3, 3, sharex=ax, sharey=ax)
ax3.set_title('angle 180')
ax3.imshow(data[99, y1:y2, x1:x2])

plt.show()

#perform the cropping on all projections
data=recon.crop(fname, data, x1, y1, x2, y2, save_out = False)

Starting cropping of data.
cropping of tomo done!
Croping took: 0.0 seconds

%matplotlib notebook

#plot sinograms
prj = plt.figure()
prj.subplots_adjust(hspace=.5)

ax = prj.add_subplot(3, 1, 1)
ax.set_title('sinogram 150')
ax.imshow(data[:, 150, :])

ax2 = prj.add_subplot(3, 1, 2, sharex=ax, sharey=ax)
ax2.set_title('sinogram 220')
ax2.imshow(data[:, 220, :])

ax3 = prj.add_subplot(3, 1, 3, sharex=ax, sharey=ax)
ax3.set_title('sinogram 450')
ax3.imshow(data[:, 450, :])

plt.show()

def histoplot(data, no_bins, parameter):
    fig = plt.figure()
    plt.hist(data, bins=no_bins)
    plt.ylabel(ylabel='frequency')
    plt.xlabel(xlabel=parameter)

    fig.show()

#OPTIONAL
#cut off is given in decimal values between 0.0 to 1.0 of the mean values between x ranges x3 to x4
cut_off = 0.9

x1_sino = 0
x2_sino = 20

broken_projs_collection = []

for i in range(len(data[0])):
    broken_projs, sino_row, cut_off_value = recon.search_broken_proj(data, height=i, x1=x1_sino, x2=x2_sino, cut_off=cut_off, sign=False)
    broken_projs_collection = np.append(broken_projs_collection, broken_projs)
    #print(i)
    #print(broken_projs)

print(broken_projs_collection)

%matplotlib notebook
histoplot(broken_projs_collection, len(data[:,:,:]), parameter='broken projections')

[ 9. 14. 42. ... 79. 84. 95.]

%matplotlib notebook

fig = plt.figure()
plt.plot(sino_row)
plt.ylabel(ylabel='Grey Level [32 bit]')
plt.xlabel(xlabel='projection')

#if list broken_projs is not empty do...
if broken_projs:
    plt.plot(broken_projs, np.max(sino_row), 'go')

    vector_cut_off_values = []
    for i in range(len(sino_row)):
        vector_cut_off_values = np.append(vector_cut_off_values, cut_off_value)

    plt.plot(vector_cut_off_values, 'r-')
    fig.show()

data_repaired = recon.repair_sino_horiz(fname, data, broken_projs, save_out=False)

%matplotlib notebook

#plot sinograms
prj = plt.figure()
prj.subplots_adjust(hspace=.5)

ax = prj.add_subplot(3, 1, 1)
ax.set_title('sinogram 150')
ax.imshow(data[:, 150, :])

ax2 = prj.add_subplot(3, 1, 2, sharex=ax, sharey=ax)
ax2.set_title('sinogram 220')
ax2.imshow(data[:, 220, :])

ax3 = prj.add_subplot(3, 1, 3, sharex=ax, sharey=ax)
ax3.set_title('sinogram 450')
ax3.imshow(data[:, 450, :])

plt.show()

#OPTIONAL
#remove vertical strips in the sinogram regime
data_ring_rem_fw = tomopy.prep.stripe.remove_stripe_fw(data, level=None, wname='db5', sigma=2, pad=True, ncore=None, nchunk=None)

%matplotlib notebook

#plot figures
sino = plt.figure()
sino.subplots_adjust(hspace=.5)

ax = sino.add_subplot(1, 2, 1)
ax.set_title('Sinogram before')
ax.imshow(data[:,400,:])

ax8 = sino.add_subplot(1, 2, 2, sharex=ax, sharey=ax)
ax8.set_title('Sinogram after fw')
ax8.imshow(data_ring_rem_fw[:,400,:])

plt.show()

data = data_ring_rem_fw
#recon.save_projections(fname, '5_PROJ_ring_rem', data)

data_filt = ndi.gaussian_filter( data, sigma=3)

fig = plt.figure()
ax1 = fig.add_subplot(1, 2, 1)
ax1.set_title('Original projection, normalized')
ax1.imshow(data[1])

ax2 = fig.add_subplot(1, 2, 2, sharex=ax1, sharey=ax1)
ax2.set_title('Gaussian filtered projection')
ax2.imshow(data_filt[1])
plt.show()

data = data_filt
del data_filt
#recon.save_projections(fname, '7_PROJ_gaussian', data)

proj_no, proj_y, proj_x = data.shape
print(proj_no)
print(proj_y)
print(proj_x)

100
480
500

#the CR will be automatically applied below
rot_center_m = recon.test_cr(data)

#for 180 scans
theta = tomopy.angles(data.shape[0], 0, 180)

The central point of the image is: 250.0
The automatically detected CoR is: 247.25

#jupyterlab
#%matplotlib inline
#jupyter notebook
%matplotlib notebook

#slice height (top to bottom), choose a slice between 0 and the height of your image in pixels:
height = 200

#recon_name is a global parameter from kblt package
name = '7_RECON_cr_test_height_' + str(height)

#reconstruct over a range of pixel values around the automatic CR:
recon_range = range(-40, 40, 1)

#manually remove any old recon test data first
recon.delete_saved_out_data(fname+name)

#reconstruct +/- 20 pixels around the automatic CR
for i in recon_range:
    print(i)
    #32 bit is better to find the optimal CoR
    rec_test_32bit = recon.recon_gridrec(fname+name, data, slice_start = height, slice_end = height+1, theta=theta, rot_center=rot_center_m+i, save_out = True, rem_old_rec = False, cr_in_filename = True)

-40
Reconstructed using GRIDREC done!
Reconstruction took: 0.11 seconds
-39
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
-38
Reconstructed using GRIDREC done!
Reconstruction took: 0.05 seconds
-37
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
-36
Reconstructed using GRIDREC done!
Reconstruction took: 0.05 seconds
-35
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
-34
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
-33
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
-32
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
-31
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
-30
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
-29
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
-28
Reconstructed using GRIDREC done!
Reconstruction took: 0.05 seconds
-27
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
-26
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
-25
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
-24
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
-23
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
-22
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
-21
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
-20
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
-19
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
-18
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
-17
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
-16
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
-15
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
-14
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
-13
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
-12
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
-11
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
-10
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
-9
Reconstructed using GRIDREC done!
Reconstruction took: 0.05 seconds
-8
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
-7
Reconstructed using GRIDREC done!
Reconstruction took: 0.06 seconds
-6
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
-5
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
-4
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
-3
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
-2
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
-1
Reconstructed using GRIDREC done!
Reconstruction took: 0.05 seconds
0
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
1
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
2
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
3
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
4
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
5
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
6
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
7
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
8
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
9
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
10
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
11
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
12
Reconstructed using GRIDREC done!
Reconstruction took: 0.05 seconds
13
Reconstructed using GRIDREC done!
Reconstruction took: 0.05 seconds
14
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
15
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
16
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
17
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
18
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
19
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
20
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
21
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
22
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
23
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
24
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
25
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
26
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
27
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
28
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
29
Reconstructed using GRIDREC done!
Reconstruction took: 0.04 seconds
30
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
31
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
32
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
33
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
34
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
35
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
36
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
37
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
38
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds
39
Reconstructed using GRIDREC done!
Reconstruction took: 0.03 seconds

#OPTIONAL
#delete old test reconstructions, only if needed
for k in range(1,2):
    path = fname + str(name)
    print(path)
    shutil.rmtree(path, ignore_errors=True)

LavaLamp/rep_16/7_RECON_cr_test_height_200

#read in cr test data into the memory
ind_slice = '*.tiff'
slice_path_in = os.path.join(fname+name, ind_slice)
print(slice_path_in)

LavaLamp/rep_16/7_RECON_cr_test_height_200/*.tiff

rec_test = ImageCollection(slice_path_in)

#concatenate the cr test data into a virtual 3D-volume stack
rec_test = rec_test.concatenate()

# compute total variation
#The lower a total variation, the more smooth the image is.

p = len(recon_range)
TV = np.zeros((p))

#print('Total variation for slice recon is: ')
for i in range(p):
    TV[i] = segment.totvar( rec_test[i,:,:] )
    #print("Test slice no.: " + str(i) + " "+ str(np.round( segment.totvar( rec_test[i,:,:] ),3 ) ))

#print("Min value is: " + str(round(np.min(TV),3)))

#%matplotlib inline
%matplotlib notebook
plt.plot(TV)
plt.ylabel('Total Variation')
plt.xlabel('Center of Rotation Range')
plt.plot(np.argmax(TV), np.max(TV), 'ro')
plt.plot(np.argmin(TV), np.min(TV), 'go')
plt.show()
print("CoR range number: " + str(np.argmax(TV)))
print("TV automatically detected CoR: " + str(rot_center_m + recon_range[np.argmax(TV)] ))

CoR range number: 79
TV automatically detected CoR: 286.25

#if using jupyterlab
#%matplotlib inline
#if using jupyter notebook
%matplotlib notebook
visualization.show_vol( rec_test, axis=3 )

Instructions:
    Controls
    --------
        'up-arrow'    : Next slice in volume
        'down-arrow'  : Previous slice in volume
        'right-arrow' : +10 slices
        'left-arrow'  : -10 slices
        'PgUp'        : +50 slices
        'PgDown'      : -50 slices

        scroll wheel  : previous/next slice

#double check, rot_center value
rot_center_m = recon.test_cr(data)
#change number [-] to the optimal slice number (actually meaning shift in X in number of pixels) in the stack above
rc_range = 41
rot_center_m = rot_center_m + recon_range[rc_range]
print('The manually set CoR is: ' + str(rot_center_m))

The central point of the image is: 250.0
The automatically detected CoR is: 247.25
The manually set CoR is: 248.25

#investigate conversion values to 8-bit, 2D test slice
l_val = -0.06
u_val = 0.02

slice_orig = rec_test[rc_range]

test_slice = rescale_intensity(slice_orig, in_range=(l_val, u_val) )
test_slice = img_as_ubyte(test_slice)

# histogram of voxel values
visualization.plot_histogram( slice_orig, nbin=255 )
print('Values in 32-bit image stack:')
print('Min: ' + str(np.min(slice_orig)))
print('Max: ' + str(np.max(slice_orig)))

visualization.plot_histogram( test_slice, nbin=255 )
print('Values in 8-bit image stack:')
print('Min: ' + str(np.min(test_slice)))
print('Max: ' + str(np.max(test_slice)))

Values in 32-bit image stack:
Min: -0.012099617
Max: 0.012226308

Values in 8-bit image stack:
Min: 50
Max: 205

#jupyterlab
#%matplotlib inline
#jupyter notebook
%matplotlib notebook

#Before correction
fig = plt.figure()
ax1 = fig.add_subplot(1, 2, 1)
ax1.set_title('32bit slice, before')
ax1.imshow(slice_orig)

#After manual test correction
#fig.subplots_adjust(hspace=.5)
ax2 = fig.add_subplot(1, 2, 2, sharex=ax1, sharey=ax1)
ax2.set_title('8bit slice, after')
ax2.imshow(test_slice)
plt.show()

#reconstruct all stack of slices
name = '8_RECON_FULL_rot_center_' + str(rot_center_m)

#reconstruct all slices
slice_start = 0
slice_end = len(data[0])

#rec = recon.recon_gridrec(fname + name + '_32bit', data, slice_start = slice_start, slice_end = slice_end, theta=theta, rot_center=rot_center_m, save_out = True, rem_old_rec = True, cr_in_filename = False)
rec = recon.recon_gridrec_8bit(fname + name + '_8bit', data, slice_start = slice_start, slice_end = slice_end, theta=theta, rot_center=rot_center_m, l_val=l_val, u_val=u_val, save_out = True, rem_old_rec = True, cr_in_filename = False)

#imshow(rec[100])

Reconstructed using GRIDREC done!
Reconstruction took: 8.11 seconds

#jupyterlab
#%matplotlib inline
#jupyter notebook
%matplotlib notebook

visualization.plot_orthogonal_slices(rec, 75, 250, 250, layout=(1,3))

#batch reconstruction

print('Flat fields')
fname_flat = sample
print(fname_flat)

t0 = time.time()

for p in range(10,11):
    t1 = time.time()

    #input folder
    fname = sample + 'rep_' + str(p) + '/'
    print(fname)

    #Read RGB data
    tomo = recon.read_rgb_convert(fname, 'tomo', save_out = False)
    flat = recon.read_rgb_convert(fname_flat, 'flat', save_out = False)
    #read in of KBLT 16bit or real NCT/XCT data
    #tomo, flat = recon.read_kblt(fname, ind_tomo, ind_flat, proj=None, sino=None)

    #---------------------

    #Normalized the projection images
    data = recon.flat_fielding(fname, tomo, flat, save_out = False)
    data[np.isnan(data)] = 0
    data[data == -inf] = 0
    data[data == inf] = 1

    #perform the cropping on all projections
    data=recon.crop(fname, data, x1, y1, x2, y2, save_out = False)

    #perform tilt correction on all projection images
    #print(tilt_step)
    #data = recon.rotate_tomo(fname, data, tilt=tilt_step, normalized = True, save_out = False)

    #search broken projections
    broken_projs_collection = []
    for i in range(len(data[0])):
        broken_projs, sino_row, cut_off_value = recon.search_broken_proj(data, height=i, x1=x1_sino, x2=x2_sino, cut_off=cut_off, sign=False)
        broken_projs_collection = np.append(broken_projs_collection, broken_projs)

    #repair projections
    data = recon.repair_sino_horiz(fname, data, broken_projs, save_out = True)

    #remove stripes in the sinogram regime
    data = tomopy.prep.stripe.remove_stripe_fw(data, level=None, wname='db5', sigma=2, pad=True, ncore=None, nchunk=None)

    #gaussian filter the projections
    data = ndi.gaussian_filter(data, sigma=3)

    #for 180 scans
    #theta = tomopy.angles(data.shape[0], 0, 180)
    if p % 2 == 0:
        theta = tomopy.angles(data.shape[0], 0, 180)
    else:
        theta = tomopy.angles(data.shape[0], 180, 360)

    #reconstruct all stack of slices
    name = '8_RECON_FULL_rot_center_' + str(rot_center_m)
    #reconstruct all slices
    slice_start = 0
    slice_end = len(data[0])

    #rec = recon.recon_gridrec(fname + name + '_32bit', data, slice_start = slice_start, slice_end = slice_end, theta=theta, rot_center=rot_center_m, save_out = True, rem_old_rec = True, cr_in_filename = False)
    rec = recon.recon_gridrec_8bit(fname + name + '_8bit', data, slice_start = slice_start, slice_end = slice_end, theta=theta, rot_center=rot_center_m, l_val=l_val, u_val=u_val, save_out = True, rem_old_rec = True, cr_in_filename = False)

    print('Recon per sample took: ' + str(round(time.time()-t1,2)) + ' seconds')

print('Entire reconstruction script for all samples took: ' + str(round(time.time()-t0,2)) + ' seconds')

Flat fields
LavaLamp/
LavaLamp/rep_10/
Starting flat-fielding of data.
Starting removal of NaN, Inf and -Inf values.
Removal of NaN, Inf and -Inf took: 0.1 seconds
Flat-fielding took: 0.31 seconds
Starting cropping of data.
cropping of tomo done!
Croping took: 0.0 seconds
LavaLamp/rep_10/4_PROJ_repaired/tomo
Saved down repaired projections tomo
Reconstructed using GRIDREC done!
Reconstruction took: 7.69 seconds
Recon per sample took: 60.48 seconds
Entire reconstruction script for all samples took: 60.48 seconds

visualization.plot_orthogonal_slices(rec, 75, 250, 250, layout=(1,3))

#OPTIONAL
#delete old reconstructions, only if needed

cr = 244.5

for k in range(16,17):
    path = sample + 'rep_' + str(k) + "/" + '8_RECON_FULL_rot_center_' + str(cr) + '_8bit/'
    print(path)
    shutil.rmtree(path, ignore_errors=True)

LavaLamp/rep_16/8_RECON_FULL_rot_center_244.5_8bit/