In [ ]:
%matplotlib inline
In [ ]:
import astra
import numpy as np
import pylab as plt
import os
import glob
import scipy.ndimage
import skimage.segmentation
import matplotlib
font = {'size' : 18}
matplotlib.rc('font', **font)
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 [ ]:
def images_diff(im1, im2):
assert(im1.shape==im2.shape)
rec_diff = np.zeros(shape=(im1.shape[0],im1.shape[1],3), dtype='float32')
im1_t = im1.copy()
im1_t = (im1_t-im1_t.min())/(im1_t.max()-im1_t.min())
im2_t = im2.copy()
im2_t = (im2_t-im2_t.min())/(im2_t.max()-im2_t.min())
# nrecon_rec_t[nrecon_rec_t<0] = 0
diff_rec = im1_t-im2_t
rec_diff[...,0] = diff_rec*(diff_rec>0)
rec_diff[...,1] = -diff_rec*(diff_rec<0)
rec_diff[...,2] = rec_diff[...,1]
return rec_diff
In [ ]:
data_root = '/diskmnt/a/makov/yaivan/MMC_1/'
nrecon_root_folder = os.path.join(data_root,'_tmp','nrecon')
nrecon_folders = glob.glob(os.path.join(nrecon_root_folder, '*'))
nrecon_folders = [nf for nf in nrecon_folders if os.path.isdir(nf)]
print len(nrecon_folders)
In [ ]:
!ls /diskmnt/a/makov/yaivan/Sand/_tmp/nrecon/bh_0_rc_0
In [ ]:
!ls /diskmnt/a/makov/yaivan/MMC_1/_tmp/nrecon/bh_*_rc_0
In [ ]:
!md5sum /diskmnt/a/makov/yaivan/MMC_1/_tmp/nrecon/bh_0_rc_0/*
In [ ]:
bh_images = {}
sr_images = {}
sl_images = {}
for bh in log_progress(np.arange(0,101,10)):
nrecon_folder='/diskmnt/a/makov/yaivan/MMC_1/_tmp/nrecon/bh_{}_rc_0/'.format(bh)
sino_raw = os.path.join(nrecon_folder, 'MMC1_2.82um__sinoraw_0960.tif')
sino_log = os.path.join(nrecon_folder, 'MMC1_2.82um__sino0960.tif')
rec_file = os.path.join(nrecon_folder, 'MMC1_2.82um__rec0960.png')
v_max = 0.52
v_min = -0.18
# nrecon_folder='/diskmnt/a/makov/yaivan/Sand/_tmp/nrecon/bh_{}_rc_0/'.format(bh)
# sino_raw = os.path.join(nrecon_folder, 'Chieftain_Unc_2.8__sinoraw_0980.tif')
# sino_log = os.path.join(nrecon_folder, 'Chieftain_Unc_2.8__sino0980.tif')
# rec_file = os.path.join(nrecon_folder, 'Chieftain_Unc_2.8__rec0980.png')
# v_max = 0.0680
# v_min = -0.0250
if os.path.isfile(rec_file):
bh_images[bh] = np.squeeze(plt.imread(rec_file).astype('float32')[...,0])
sr_images[bh] = plt.imread(sino_raw).astype('float32')
sl_images[bh] = plt.imread(sino_log).astype('float32')
In [ ]:
sr = plt.imread('/diskmnt/a/makov/yaivan/Sand/_tmp/nrecon/noTS/Chieftain_Unc_2.8__sinoraw_0980.tif').astype('float32')
sl = plt.imread('/diskmnt/a/makov/yaivan/Sand/_tmp/nrecon/noTS/Chieftain_Unc_2.8__sino0980.tif').astype('float32')
images = sorted(glob.glob(r'/diskmnt/a/makov/yaivan/Sand/Raw/Chieftain_Unc_2.8_????.tif'))
In [ ]:
sr=sr_images[0]
sl=sl_images[0]
In [ ]:
mx = []
mn = []
me = []
for im in log_progress(images[::10]):
i = plt.imread(im).astype('float32')
mx.append(np.max(i))
mn.append(np.min(i))
me.append(np.mean(i))
In [ ]:
plt.figure(figsize=(10,12))
plt.imshow(sr, cmap=plt.cm.gray)
plt.colorbar(orientation='horizontal')
In [ ]:
plt.figure(figsize=(10,5))
plt.plot(mx, label = 'max')
# plt.plot(mn, label = 'min')
# plt.plot(me, label = 'mean')
plt.grid()
# plt.legend(loc=0)
In [ ]:
srf = sr_images[0]
In [ ]:
slf = sl_images[0]
In [ ]:
plt.figure(figsize=(10,12))
plt.imshow(sr-srf, cmap=plt.cm.gray)
plt.colorbar(orientation='horizontal')
In [ ]:
plt.figure(figsize=(10,12))
plt.imshow(sl-slf, cmap=plt.cm.gray)
plt.colorbar(orientation='horizontal')
In [ ]:
plt.figure(figsize=(10,12))
plt.imshow(sr, cmap=plt.cm.gray)
plt.colorbar(orientation='horizontal')
In [ ]:
print(sr.min(), sr.max() )
In [ ]:
print(np.log(sr.min()), np.log(sr.max()))
In [ ]:
plt.figure(figsize=(10,12))
plt.imshow(sl, cmap=plt.cm.gray)
plt.colorbar(orientation='horizontal')
In [ ]:
print(sl.min(), sl.max() )
In [ ]:
plt.figure(figsize=(10,12))
plt.imshow(-np.log(sr/65535)-(sl/65535), cmap=plt.cm.gray)
plt.colorbar(orientation='horizontal')
In [ ]:
def images_diff(im1, im2):
assert(im1.shape==im2.shape)
rec_diff = np.zeros(shape=(im1.shape[0],im1.shape[1],3), dtype='float32')
im1_t = im1.copy()
im1_t = (im1_t-im1_t.min())/(im1_t.max()-im1_t.min())
im2_t = im2.copy()
im2_t = (im2_t-im2_t.min())/(im2_t.max()-im2_t.min())
# nrecon_rec_t[nrecon_rec_t<0] = 0
diff_rec = im1_t-im2_t
rec_diff[...,0] = diff_rec*(diff_rec>0)
rec_diff[...,1] = -diff_rec*(diff_rec<0)
rec_diff[...,2] = rec_diff[...,1]
return rec_diff
In [ ]:
r=-np.log(sr/65535)
l=sl/65535
plt.figure(figsize=(10,12))
plt.imshow(100*images_diff(r,l))
plt.colorbar(orientation='horizontal')
In [ ]:
r=-np.log(sr/65535)
rr = (r.T-r.min(axis=-1)).T
l=sl/65535
plt.figure(figsize=(10,12))
plt.imshow(100*images_diff(rr,l))
plt.colorbar(orientation='horizontal')
In [ ]:
r=-np.log(sr/65535)
rr = (r-r.min(axis=-1).max())
rr[rr<0]=0
l=sl/65535
plt.figure(figsize=(10,12))
plt.imshow(100*images_diff(rr,l))
plt.colorbar(orientation='horizontal')
In [ ]:
plt.figure(figsize=(10,10))
plt.plot(sr.max(axis=0))
plt.plot(sr.min(axis=0))
plt.grid(True)
In [ ]:
plt.figure(figsize=(10,10))
plt.plot(sr.max(axis=-1))
plt.plot(sr.min(axis=-1))
plt.plot(sr.mean(axis=-1))
plt.grid(True)
In [ ]:
2**16*0.9
In [ ]:
x = sl[100::100]
y=sr.copy()[100::100]
trh = y[:,:10].mean(axis=-1)
for i in range(y.shape[0]):
t=y[i]
t[t>trh[i]]=trh[i]
t/=trh[i]
y = -np.log(y)
# y = (y.T-y.min(axis=-1)).T
plt.figure(figsize=(10,10))
plt.plot(x, y)
# plt.plot(x,y)
plt.grid(True)
In [ ]:
x = sl.copy()
y=sr.copy()
trh = y[:,:11].mean(axis=-1)
for i in range(y.shape[0]):
t=y[i]
t[t>trh[i]]=trh[i]
t/=trh[i]
y = -np.log(y)
In [ ]:
plt.figure(figsize=(10,12))
plt.imshow(1000*images_diff(x,y))
plt.colorbar(orientation='horizontal')
In [ ]:
plt.figure(figsize=(10,12))
plt.imshow(x-y/y.max()*65535)
plt.colorbar(orientation='horizontal')
In [ ]:
np.max(x-y/y.max()*65535)
In [ ]:
trh[np.newaxis,:].shape
In [ ]:
from scipy.optimize import curve_fit
In [ ]:
p = curve_fit(lambda x,a,b: a*x+b, x.ravel(), y.ravel())
print p
In [ ]:
-p[0][1]/p[0][0]
In [ ]:
sr.max()
In [ ]:
sr.max(axis=-1).mean()
In [ ]:
r=sr.copy()
r[sl>0]=65535
plt.figure(figsize=(10,10))
plt.plot(sr[:,:1].mean(axis=-1),r.min(axis=-1),'o')
In [ ]:
r=sr.copy()
r[sl>0]=65535
plt.figure(figsize=(10,10))
for ss in range(8,13):
# plt.plot(r.min(axis=-1))
# plt.plot(sr.max(axis=-1))
plt.plot(sr[:,:ss].mean(axis=-1),r.min(axis=-1), 'o', label=ss)
plt.grid()
plt.legend(loc=0)
# plt.plot(x,y)
# plt.colorbar()
In [ ]:
r=-np.log(sr)
rr = ((r.T-r.min(axis=-1))/r.mean(axis=-1)).T*r.mean()
rr[rr<0]=0
l=sl/60000
plt.figure(figsize=(10,12))
# plt.imshow(100*images_diff(rr,l))
plt.imshow(rr/(1.2*l),cmap=plt.cm.gray,vmin=0.9, vmax=1.1)
plt.colorbar(orientation='horizontal')
In [ ]:
x = sl
y = -np.log(sr)
plt.figure(figsize=(10,10))
plt.plot(x[100], y[100],'x')
plt.grid(True)
In [ ]:
np.log(55000)
In [ ]:
np.exp(0.05)
In [ ]:
plt.figure(figsize=(10,10))
plt.plot(sl.min(axis=-1))
plt.grid(True)
In [ ]:
plt.figure(figsize=(10,10))
plt.plot(sl.max(axis=-1))
plt.grid(True)
In [ ]:
plt.figure(figsize=(10,10))
plt.plot(sl.mean(axis=-1))
plt.grid(True)
In [ ]:
plt.figure(figsize=(10,10))
plt.plot(sr.min(axis=-1))
plt.grid(True)
In [ ]:
In [ ]:
plt.figure(figsize=(10,10))
plt.plot(sr.max(axis=-1))
plt.grid(True)
In [ ]:
plt.figure(figsize=(10,10))
plt.plot(sr.mean(axis=-1))
plt.grid(True)
In [ ]:
np.exp(10.9)
In [ ]:
np.min(sl)
In [ ]:
from skimage.filter import gaussian
In [ ]:
In [ ]:
print(sorted(bh_images.keys()))
In [ ]:
s = []
v_prev = None
for k in log_progress(sorted(bh_images.keys())):
v = bh_images[k]
v = v*(v_max-v_min)/(v.max()-v.min())+v_min
v = medianBlur(v,3)
# if v_prev is None:
# v_prev = v
# continue
r = np.mean(v-gaussian(v,10))/np.mean(v)
v_prev = v
s.append((k,r))
s = np.array(s)
plt.plot(s[:,0],s[:,1],'*')
plt.grid(True)
In [ ]:
s = []
v_prev = None
for k in log_progress(sorted(bh_images.keys())):
v = bh_images[k]
v = v*(v_max-v_min)/(v.max()-v.min())+v_min
# v = medianBlur(v,3)
if v_prev is None:
v_prev = v
continue
t = np.linalg.norm(v-v_prev)
v_prev = v
s.append((k,t))
s = np.array(s)
plt.plot(s[:,0],s[:,1],'*')
plt.grid(True)
In [ ]:
# %load /diskmnt/a/makov/yaivan/Sand/Raw/Chieftain_Unc_2.8_.log
[System]
Scanner=Skyscan1172
Instrument S/N=08G01121
Hardware version=A
Software=Version 1. 5 (build 23)
Home directory=C:\SkyScan
Source Type=Hamamatsu 100/250
Camera=Hamamatsu 10Mp camera
Camera Pixel Size (um)= 11.40
CameraXYRatio=1.0010
Incl.in lifting (um/mm)=-0.4750
[User]
User Name=IYakimchuk
Computer Name=SLB-6BBX74J
[Acquisition]
Data directory=E:\Results\Yakimchuk\2016_Digital Fracture Conductivity\00_Pre-Study\01. Sands\01. Chieftain Sand 20-40 Unconfined\Raw
Filename Prefix=Chieftain_Unc_2.8_
Configuration=C:\Skyscan1172A_10MP_Hamamatsu\std
Number of Files= 2031
Source Voltage (kV)= 100
Source Current (uA)= 100
Number of Rows= 2096
Number of Columns= 4000
Image crop origin X= 0
Image crop origin Y=0
Camera binning=1x1
Image Rotation=0.5200
Gantry direction=CC
Image Pixel Size (um)= 2.83
Object to Source (mm)=56.000
Camera to Source (mm)=225.315
Vertical Object Position (mm)=37.797
Optical Axis (line)= 980
Filter=Al 0.5 mm
Image Format=TIFF
Depth (bits)=16
Screen LUT=0
Exposure (ms)= 1767
Rotation Step (deg)=0.100
Frame Averaging=ON (9)
Random Movement=OFF (10)
Use 360 Rotation=NO
Geometrical Correction=ON
Camera Offset=OFF
Median Filtering=ON
Flat Field Correction=ON
Rotation Direction=CC
Scanning Trajectory=ROUND
Type Of Motion=STEP AND SHOOT
Study Date and Time=Jan 23, 2016 05:45:50
Scan duration=10:06:00
[Reconstruction]
Reconstruction Program=NRecon
Program Version=Version: 1.6.5.8
Program Home Directory=C:\SkyScan\NRecon_GPU
Reconstruction engine=NReconServer
Engine version=Version: 1.6.5
Reconstruction from batch=No
Reconstruction servers= slb-8hlv74j
Option for additional F4F float format=OFF
Dataset Origin=Skyscan1172
Dataset Prefix=Chieftain_Unc_2.8_
Dataset Directory=E:\Results\Yakimchuk\2016_Digital Fracture Conductivity\00_Pre-Study\01. Sands\01. Chieftain Sand 20-40 Unconfined\Raw
Output Directory=E:\Results\Yakimchuk\2016_Digital Fracture Conductivity\00_Pre-Study\01. Sands\01. Chieftain Sand 20-40 Unconfined\Reconstructed
Time and Date=Jan 26, 2016 16:58:16
First Section=98
Last Section=1982
Reconstruction duration per slice (seconds)=5.300265
Total reconstruction time (1885 slices) in seconds=9991.000000
Postalignment=-6.00
Section to Section Step=1
Sections Count=1885
Result File Type=PNG
Result File Header Length (bytes)=Unknown: compressed JPG format (100%)
Result Image Width (pixels)=4000
Result Image Height (pixels)=4000
Pixel Size (um)=2.83356
Reconstruction Angular Range (deg)=203.00
Use 180+=OFF
Angular Step (deg)=0.1000
Smoothing=0
Ring Artifact Correction=20
Draw Scales=OFF
Object Bigger than FOV=OFF
Reconstruction from ROI=OFF
Filter cutoff relative to Nyquisit frequency=100
Filter type=0
Filter type meaning(1)=0: Hamming (Ramp in case of optical scanner); 1: Hann; 2: Ramp; 3: Almost Ramp;
Filter type meaning(2)=11: Cosine; 12: Shepp-Logan; [100,200]: Generalized Hamming, alpha=(iFilter-100)/100
Undersampling factor=1
Threshold for defect pixel mask (%)=0
Beam Hardening Correction (%)=60
CS Static Rotation (deg)=0.0
Minimum for CS to Image Conversion=-0.0250
Maximum for CS to Image Conversion=0.0680
HU Calibration=OFF
BMP LUT=0
Cone-beam Angle Horiz.(deg)=11.557156
Cone-beam Angle Vert.(deg)=6.070880
In [ ]:
# %load /diskmnt/a/makov/yaivan/Sand/_tmp/nrecon/bh_0_rc_0/tomo_config.log
[System]
Scanner = Skyscan1172
Instrument S/N = 08G01121
Hardware version = A
Software = Version 1. 5 (build 23)
Home directory = C:\SkyScan
Source Type = Hamamatsu 100/250
Camera = Hamamatsu 10Mp camera
Camera Pixel Size (um) = 11.40
CameraXYRatio = 1.0010
Incl.in lifting (um/mm) = -0.4750
[User]
User Name = IYakimchuk
Computer Name = SLB-6BBX74J
[Acquisition]
Data directory = E:\Results\Yakimchuk\2016_Digital Fracture Conductivity\00_Pre-Study\01. Sands\01. Chieftain Sand 20-40 Unconfined\Raw
Filename Prefix = Chieftain_Unc_2.8_
Configuration = C:\Skyscan1172A_10MP_Hamamatsu\std
Number of Files = 2031
Source Voltage (kV) = 100
Source Current (uA) = 100
Number of Rows = 2096
Number of Columns = 4000
Image crop origin X = 0
Image crop origin Y = 0
Camera binning = 1x1
Image Rotation = 0.5200
Gantry direction = CC
Image Pixel Size (um) = 2.83
Object to Source (mm) = 56.000
Camera to Source (mm) = 225.315
Vertical Object Position (mm) = 37.797
Optical Axis (line) = 980
Filter = Al 0.5 mm
Image Format = TIFF
Depth (bits) = 16
Screen LUT = 0
Exposure (ms) = 1767
Rotation Step (deg) = 0.100
Frame Averaging = ON (9)
Random Movement = OFF (10)
Use 360 Rotation = NO
Geometrical Correction = ON
Camera Offset = OFF
Median Filtering = ON
Flat Field Correction = ON
Rotation Direction = CC
Scanning Trajectory = ROUND
Type Of Motion = STEP AND SHOOT
Study Date and Time = Jan 23, 2016 05:45:50
Scan duration = 10:06:00
[Reconstruction]
Reconstruction Program = NRecon
Program Version = Version: 1.6.5.8
Program Home Directory = C:\SkyScan\NRecon_GPU
Reconstruction engine = NReconServer
Engine version = Version: 1.6.5
Reconstruction from batch = No
Reconstruction servers = slb-8hlv74j
Option for additional F4F float format = OFF
Dataset Origin = Skyscan1172
Dataset Prefix = Chieftain_Unc_2.8_
Dataset Directory = f:\big\yaivan\Sand\Raw
Output Directory = C:\Users\makov\Desktop\NRecon_out\Chieftain_Unc_2.8_\bh_0_rc_0
Time and Date = Jan 26, 2016 16:58:16
First Section = 980
Last Section = 980
Reconstruction duration per slice (seconds) = 5.300265
Total reconstruction time (1885 slices) in seconds = 9991.000000
Postalignment = -6.00
Section to Section Step = 1
Sections Count = 1885
Result File Type = PNG
Result File Header Length (bytes) = Unknown: compressed JPG format (100%)
Result Image Width (pixels) = 4000
Result Image Height (pixels) = 4000
Pixel Size (um) = 2.83356
Reconstruction Angular Range (deg) = 203.00
Use 180+ = OFF
Angular Step (deg) = 0.1000
Smoothing = 0
Ring Artifact Correction = 0
Draw Scales = OFF
Object Bigger than FOV = OFF
Reconstruction from ROI = OFF
Filter cutoff relative to Nyquisit frequency = 100
Filter type = 0
Filter type meaning(1) = 0: Hamming (Ramp in case of optical scanner); 1: Hann; 2: Ramp; 3: Almost Ramp;
Filter type meaning(2) = 11: Cosine; 12: Shepp-Logan; [100,200]: Generalized Hamming, alpha=(iFilter-100)/100
Undersampling factor = 1
Threshold for defect pixel mask (%) = 0
Beam Hardening Correction (%) = 0
CS Static Rotation (deg) = 0.0
Minimum for CS to Image Conversion = -0.0250
Maximum for CS to Image Conversion = 0.0680
HU Calibration = OFF
BMP LUT = 0
Cone-beam Angle Horiz.(deg) = 11.557156
Cone-beam Angle Vert.(deg) = 6.070880
In [ ]:
tmp_bh =bh_images[70]
In [ ]:
plt.figure(figsize=(10,12))
plt.imshow(tmp_bh, cmap=plt.cm.gray)
plt.colorbar(orientation='horizontal')
In [ ]:
plt.figure(figsize=(10,12))
plt.imshow(bh_images[60]-bh_images[0], cmap=plt.cm.gray)
plt.colorbar(orientation='horizontal')
In [ ]:
roi = tmp_bh > (np.percentile(tmp_bh.flatten(),95))
In [ ]:
roi = tmp_bh > (np.mean(tmp_bh)*1.5)
In [ ]:
from skimage.morphology import erosion
In [ ]:
roi = erosion(roi)
In [ ]:
plt.figure(figsize=(10,12))
plt.imshow(roi, cmap=plt.cm.gray)
plt.colorbar(orientation='horizontal')
In [ ]:
filt_im = roi*tmp_bh
plt.figure(figsize=(10,12))
plt.imshow(filt_im, cmap=plt.cm.gray)
plt.colorbar(orientation='horizontal')
In [ ]:
grad = np.gradient(np.ntmp_bh)
In [ ]:
In [ ]:
grad_abs = np.sqrt(grad[0]**2+grad[1]**2)
In [ ]:
plt.figure(figsize=(10,12))
plt.imshow(grad_abs*tmp_bh, cmap=plt.cm.gray,vmax=0.)
plt.colorbar(orientation='horizontal')
In [ ]:
grad_median = np.median(grad_abs.flatten())
print(grad_median)
In [ ]:
plt.figure(figsize=(10,12))
plt.imshow(grad_abs>grad_median, cmap=plt.cm.gray,vmax=0.01)
plt.colorbar(orientation='horizontal')
In [ ]:
plt.hist(grad_abs[grad_abs>grad_median].flatten(),bins=100);
In [ ]:
plt.figure(figsize=(10,12))
plt.imshow((bh_images[0]-bh_images[92]), cmap=plt.cm.gray)
plt.colorbar(orientation='horizontal')
In [ ]:
from cv2 import medianBlur
In [ ]:
In [ ]: