%matplotlib inline
import matplotlib.pyplot as plt
from glob import glob
import pandas as pd
import os
from skimage.io import imread
import numpy as np
all_images = glob(os.path.join('..', 'data', '*.tif'))
img_df = pd.DataFrame([dict(path = c_file) for c_file in all_images])
img_df['basename'] = img_df['path'].map(lambda x: os.path.basename(x).split('.')[0])
img_df['frame_id'] = img_df['basename'].map(lambda x: '_'.join(x.split('_')[0:-1]))
img_df['sample'] = img_df['basename'].map(lambda x: x.split('_')[0])
img_df['slice'] = img_df['basename'].map(lambda x: int(x.split('_')[3]))
img_df['frame'] = img_df['basename'].map(lambda x: int(x.split('_')[1]))
# group by prefix
img_3d_df = img_df.sort_values(['sample', 'frame_id', 'slice']).groupby(['sample', 'frame', 'frame_id'])['path'].apply(list).reset_index()
img_3d_df
| sample | frame | frame_id | path | |
|---|---|---|---|---|
| 0 | MI5 | 0 | MI5_0_00 | [../data/MI5_0_00_1500.rec.8bit.tif, ../data/M... |
| 1 | MI5 | 1 | MI5_1_00 | [../data/MI5_1_00_1500.rec.8bit.tif, ../data/M... |
| 2 | MI5 | 1 | MI5_1_01 | [../data/MI5_1_01_1500.rec.8bit.tif, ../data/M... |
| 3 | MI5 | 1 | MI5_1_02 | [../data/MI5_1_02_1500.rec.8bit.tif, ../data/M... |
| 4 | MI7 | 0 | MI7_0_00 | [../data/MI7_0_00_1300.rec.8bit.tif, ../data/M... |
| 5 | MI7 | 1 | MI7_1_00 | [../data/MI7_1_00_1300.rec.8bit.tif, ../data/M... |
| 6 | MI7 | 1 | MI7_1_01 | [../data/MI7_1_01_1300.rec.8bit.tif, ../data/M... |
| 7 | MI7 | 1 | MI7_1_02 | [../data/MI7_1_02_1300.rec.8bit.tif, ../data/M... |
from skimage.filters import threshold_minimum
from skimage.morphology.convex_hull import convex_hull_image
from skimage.morphology import binary_closing, binary_opening, disk, label
from skimage.measure import regionprops
from skimage.color import label2rgb
from scipy import ndimage as ndi
def threshold_slice_2d(raw_slice):
in_slice = ndi.median_filter(raw_slice, [7, 7])
# basic threshold
thresh_min = threshold_minimum(in_slice)
# make solid component
solid_part = binary_closing(in_slice>thresh_min, disk(4))
solid_part = ndi.binary_fill_holes(solid_part)
# calculate porous component inside solid
porous_part = (in_slice < thresh_min)*solid_part
#porous_part = binary_opening(porous_part, disk(2))
porous_part = binary_closing(porous_part, disk(2))
# create label by combining
label_slice = solid_part.astype(np.uint8).clip(0,1)
label_slice += (porous_part).astype(np.uint8).clip(0,1)
return label_slice, thresh_min
def paths_to_vol(in_paths):
return np.stack([threshold_slice_2d(imread(x))[0] for x in in_paths],0)
def paths_to_raw_vol(in_paths):
return np.stack([imread(x) for x in in_paths], 0)
import ipywidgets as ipw
import ipyvolume as p3
_, c_row = next(img_3d_df.query('sample == "MI5"').iterrows())
n_vol = paths_to_raw_vol(c_row['path'])
fig = p3.figure()
# create a custom LUT
temp_tf = plt.cm.gray(np.linspace(0, 1, 256))
temp_tf[:,3] = np.linspace(-.3, 0.5, 256).clip(0, 1) # make transparency more aggressive
tf = p3.transferfunction.TransferFunction(rgba = temp_tf)
p3.volshow(n_vol.astype(np.float32),
lighting = True,
max_opacity=0.5,
tf = tf,
controls = True)
p3.show()
p3.save('mi5_raw.html')
/Users/mader/anaconda/envs/qbi2018/lib/python3.6/site-packages/ipyvolume/serialize.py:66: RuntimeWarning: invalid value encountered in true_divide gradient = gradient / np.sqrt(gradient[0]**2 + gradient[1]**2 + gradient[2]**2)
Failed to display Jupyter Widget of type VBox.
If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean that the widgets JavaScript is still loading. If this message persists, it likely means that the widgets JavaScript library is either not installed or not enabled. See the Jupyter Widgets Documentation for setup instructions.
If you're reading this message in another frontend (for example, a static rendering on GitHub or NBViewer), it may mean that your frontend doesn't currently support widgets.
/Users/mader/anaconda/envs/qbi2018/lib/python3.6/site-packages/ipyvolume/serialize.py:66: RuntimeWarning: invalid value encountered in true_divide gradient = gradient / np.sqrt(gradient[0]**2 + gradient[1]**2 + gradient[2]**2)
n_vol = paths_to_vol(c_row['path'])
p3.quickvolshow(n_vol)
/Users/mader/anaconda/envs/qbi2018/lib/python3.6/site-packages/ipyvolume/serialize.py:66: RuntimeWarning: invalid value encountered in true_divide gradient = gradient / np.sqrt(gradient[0]**2 + gradient[1]**2 + gradient[2]**2)
Failed to display Jupyter Widget of type VBox.
If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean that the widgets JavaScript is still loading. If this message persists, it likely means that the widgets JavaScript library is either not installed or not enabled. See the Jupyter Widgets Documentation for setup instructions.
If you're reading this message in another frontend (for example, a static rendering on GitHub or NBViewer), it may mean that your frontend doesn't currently support widgets.
_, c_row = next(img_3d_df.query('sample == "MI7"').iterrows())
n_vol = paths_to_vol(c_row['path'])
p3.quickvolshow(n_vol, lighting=True)
/Users/mader/anaconda/envs/qbi2018/lib/python3.6/site-packages/ipyvolume/serialize.py:66: RuntimeWarning: invalid value encountered in true_divide gradient = gradient / np.sqrt(gradient[0]**2 + gradient[1]**2 + gradient[2]**2)
Failed to display Jupyter Widget of type VBox.
If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean that the widgets JavaScript is still loading. If this message persists, it likely means that the widgets JavaScript library is either not installed or not enabled. See the Jupyter Widgets Documentation for setup instructions.
If you're reading this message in another frontend (for example, a static rendering on GitHub or NBViewer), it may mean that your frontend doesn't currently support widgets.