4. dev with dask#

[1]:

import os
from collections import defaultdict
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import tifffile as tff
import holoviews as hv
import hvplot
import dask.array as da
import skimage
import scipy.stats
from scipy import ndimage
from dask_image import ndfilters, ndmeasure, ndmorph
from nima import nima
from nima import utils

%load_ext autoreload
%autoreload 2

fp = "../../tests/data/1b_c16_15.tif"

[2]:

daimg = da.from_zarr(tff.imread(fp, aszarr=True))
daimg

[2]:

	Array	Chunk
Bytes	6.00 MiB	512.00 kiB
Shape	(4, 3, 512, 512)	(1, 1, 512, 512)
Dask graph	12 chunks in 2 graph layers
Data type	uint16 numpy.ndarray

[3]:

utils.bg(daimg[0, 0].compute())

[ ] | 0% Completed | 198.63 us: </pre>
[ ] | 0% Completed | 198.63 us: end{sphinxVerbatim}

[ ] | 0% Completed | 198.63 us

[########################################] | 100% Completed | 102.08 ms: </pre>
[########################################] | 100% Completed | 102.08 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 102.08 ms

[3]:

(457.8380944892146, 48.50286684025243)

[4]:

def dabg(daimg):
    r = defaultdict(list)
    n_t, n_c = daimg.shape[:2]
    for t in range(n_t):
        for c in range(n_c):
            r[c].append(utils.bg(daimg[t, c].compute())[0])
    return pd.DataFrame(r)


dabg(daimg)

[ ] | 0% Completed | 157.26 us: </pre>
[ ] | 0% Completed | 157.26 us: end{sphinxVerbatim}

[ ] | 0% Completed | 157.26 us

[########################################] | 100% Completed | 101.31 ms: </pre>
[########################################] | 100% Completed | 101.31 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.31 ms

[ ] | 0% Completed | 163.37 us: </pre>
[ ] | 0% Completed | 163.37 us: end{sphinxVerbatim}

[ ] | 0% Completed | 163.37 us

[########################################] | 100% Completed | 103.10 ms: </pre>
[########################################] | 100% Completed | 103.10 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 103.10 ms

[ ] | 0% Completed | 273.96 us: </pre>
[ ] | 0% Completed | 273.96 us: end{sphinxVerbatim}

[ ] | 0% Completed | 273.96 us

[########################################] | 100% Completed | 101.21 ms: </pre>
[########################################] | 100% Completed | 101.21 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.21 ms

[ ] | 0% Completed | 119.56 us: </pre>
[ ] | 0% Completed | 119.56 us: end{sphinxVerbatim}

[ ] | 0% Completed | 119.56 us

[########################################] | 100% Completed | 101.04 ms: </pre>
[########################################] | 100% Completed | 101.04 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.04 ms

[ ] | 0% Completed | 146.70 us: </pre>
[ ] | 0% Completed | 146.70 us: end{sphinxVerbatim}

[ ] | 0% Completed | 146.70 us

[########################################] | 100% Completed | 101.33 ms: </pre>
[########################################] | 100% Completed | 101.33 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.33 ms

[ ] | 0% Completed | 127.39 us: </pre>
[ ] | 0% Completed | 127.39 us: end{sphinxVerbatim}

[ ] | 0% Completed | 127.39 us

[########################################] | 100% Completed | 101.09 ms: </pre>
[########################################] | 100% Completed | 101.09 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.09 ms

[ ] | 0% Completed | 122.75 us: </pre>
[ ] | 0% Completed | 122.75 us: end{sphinxVerbatim}

[ ] | 0% Completed | 122.75 us

[########################################] | 100% Completed | 101.02 ms: </pre>
[########################################] | 100% Completed | 101.02 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.02 ms

[ ] | 0% Completed | 155.11 us: </pre>
[ ] | 0% Completed | 155.11 us: end{sphinxVerbatim}

[ ] | 0% Completed | 155.11 us

[########################################] | 100% Completed | 102.03 ms: </pre>
[########################################] | 100% Completed | 102.03 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 102.03 ms

[ ] | 0% Completed | 135.41 us: </pre>
[ ] | 0% Completed | 135.41 us: end{sphinxVerbatim}

[ ] | 0% Completed | 135.41 us

[########################################] | 100% Completed | 101.13 ms: </pre>
[########################################] | 100% Completed | 101.13 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.13 ms

[ ] | 0% Completed | 129.10 us: </pre>
[ ] | 0% Completed | 129.10 us: end{sphinxVerbatim}

[ ] | 0% Completed | 129.10 us

[########################################] | 100% Completed | 114.50 ms: </pre>
[########################################] | 100% Completed | 114.50 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 114.50 ms

[ ] | 0% Completed | 151.49 us: </pre>
[ ] | 0% Completed | 151.49 us: end{sphinxVerbatim}

[ ] | 0% Completed | 151.49 us

[########################################] | 100% Completed | 102.67 ms: </pre>
[########################################] | 100% Completed | 102.67 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 102.67 ms

[ ] | 0% Completed | 148.72 us: </pre>
[ ] | 0% Completed | 148.72 us: end{sphinxVerbatim}

[ ] | 0% Completed | 148.72 us

[########################################] | 100% Completed | 101.31 ms: </pre>
[########################################] | 100% Completed | 101.31 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.31 ms

[4]:

	0	1	2
0	457.838094	257.010233	289.378226
1	457.295251	259.072941	289.627119
2	457.760165	260.182107	290.268666
3	453.995204	257.189941	285.613624

[5]:

def dabg_fg(daimg, erf_pvalue=1e-100, size=10):
    n_t, n_c = daimg.shape[:2]
    bgs = defaultdict(list)
    fgs = defaultdict(list)
    for t in range(n_t):
        p = np.ones(daimg.shape[-2:])
        multichannel = daimg[t].compute()
        for c in range(n_c):
            av, sd = utils.bg(multichannel[c])
            p = p * utils.prob(multichannel[c], av, sd)
            bgs[c].append(av)
        mask = ndimage.median_filter((p) ** (1 / n_c), size=size) < erf_pvalue
        for c in range(n_c):
            fgs[c].append(np.ma.mean(np.ma.masked_array(multichannel[c], mask=~mask)))
    return pd.DataFrame(bgs), pd.DataFrame(fgs)


dfb, dff = dabg_fg(daimg)

[ ] | 0% Completed | 135.07 us: </pre>
[ ] | 0% Completed | 135.07 us: end{sphinxVerbatim}

[ ] | 0% Completed | 135.07 us

[########################################] | 100% Completed | 101.23 ms: </pre>
[########################################] | 100% Completed | 101.23 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.23 ms

[ ] | 0% Completed | 168.33 us: </pre>
[ ] | 0% Completed | 168.33 us: end{sphinxVerbatim}

[ ] | 0% Completed | 168.33 us

[########################################] | 100% Completed | 102.13 ms: </pre>
[########################################] | 100% Completed | 102.13 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 102.13 ms

[ ] | 0% Completed | 156.53 us: </pre>
[ ] | 0% Completed | 156.53 us: end{sphinxVerbatim}

[ ] | 0% Completed | 156.53 us

[########################################] | 100% Completed | 102.23 ms: </pre>
[########################################] | 100% Completed | 102.23 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 102.23 ms

[ ] | 0% Completed | 161.92 us: </pre>
[ ] | 0% Completed | 161.92 us: end{sphinxVerbatim}

[ ] | 0% Completed | 161.92 us

[########################################] | 100% Completed | 101.62 ms: </pre>
[########################################] | 100% Completed | 101.62 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.62 ms

[6]:

plt.subplot(121)
((dff - dfb)[0] / (dff - dfb)[2]).plot(marker="s")
plt.grid()
plt.subplot(122)
((dff - dfb)[2] / (dff - dfb)[1]).plot(marker="o")
plt.grid()

../_images/tutorials_dev_xr_hvplot_holoviews_6_0.png

NEXT: - make utils.bg and utils.prob work with dask arrays

[7]:

def dmask(daim, erf_pvalue=1e-100, size=10):
    n_c = daim.shape[0]
    im = daim[0].compute()
    p = utils.prob(im, *utils.bg(im))
    for c in range(1, n_c):
        im = daim[c].compute()
        p = p * utils.prob(im, *utils.bg(im))
    p = ndimage.median_filter((p) ** (1 / n_c), size=size)
    mask = p < erf_pvalue
    return skimage.morphology.remove_small_objects(mask)
    # mask = skimage.morphology.remove_small_holes(mask)
    # return np.ma.masked_array(plane, mask=~mask), np.ma.masked_array(plane, mask=mask)


mask = dmask(daimg[2])

lab, nlab = ndimage.label(mask)
lab, nlab

[ ] | 0% Completed | 167.06 us: </pre>
[ ] | 0% Completed | 167.06 us: end{sphinxVerbatim}

[ ] | 0% Completed | 167.06 us

[########################################] | 100% Completed | 101.25 ms: </pre>
[########################################] | 100% Completed | 101.25 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.25 ms

[ ] | 0% Completed | 154.34 us: </pre>
[ ] | 0% Completed | 154.34 us: end{sphinxVerbatim}

[ ] | 0% Completed | 154.34 us

[########################################] | 100% Completed | 103.12 ms: </pre>
[########################################] | 100% Completed | 103.12 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 103.12 ms

[ ] | 0% Completed | 159.02 us: </pre>
[ ] | 0% Completed | 159.02 us: end{sphinxVerbatim}

[ ] | 0% Completed | 159.02 us

[########################################] | 100% Completed | 101.45 ms: </pre>
[########################################] | 100% Completed | 101.45 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.45 ms

[7]:

(array([[0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        ...,
        [2, 2, 2, ..., 0, 0, 0],
        [2, 2, 2, ..., 0, 0, 0],
        [2, 2, 2, ..., 0, 0, 0]], dtype=int32),
 2)

[8]:

pr = skimage.measure.regionprops(lab, intensity_image=daimg[0][0])
pr[1].equivalent_diameter

[8]:

205.6159504584483

[9]:

max_diameter = pr[0].equivalent_diameter
size = int(max_diameter * 0.3)
size

[9]:

[10]:

t = 0
mask = dmask(daimg[t])
# skimage.io.imshow(mask)
lab, nlab = ndimage.label(mask)

distance = ndimage.distance_transform_edt(mask)
# distance = skimage.filters.gaussian(distance, sigma=0)   min_distance=size,
coords = skimage.feature.peak_local_max(
    distance, footprint=np.ones((size, size)), labels=lab
)
mm = np.zeros(distance.shape, dtype=bool)
mm[tuple(coords.T)] = True
# markers, _ = ndimage.label(mm)
markers = skimage.measure.label(mm)

labels = skimage.segmentation.watershed(-distance, markers, mask=mask)

_, (ax0, ax1, ax2) = plt.subplots(1, 3)
ax0.imshow(distance)
ax1.imshow(labels)
ax2.imshow(labels == 3)
coords

[ ] | 0% Completed | 139.94 us: </pre>
[ ] | 0% Completed | 139.94 us: end{sphinxVerbatim}

[ ] | 0% Completed | 139.94 us

[########################################] | 100% Completed | 101.19 ms: </pre>
[########################################] | 100% Completed | 101.19 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.19 ms

[ ] | 0% Completed | 143.34 us: </pre>
[ ] | 0% Completed | 143.34 us: end{sphinxVerbatim}

[ ] | 0% Completed | 143.34 us

[########################################] | 100% Completed | 102.95 ms: </pre>
[########################################] | 100% Completed | 102.95 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 102.95 ms

[ ] | 0% Completed | 127.66 us: </pre>
[ ] | 0% Completed | 127.66 us: end{sphinxVerbatim}

[ ] | 0% Completed | 127.66 us

[########################################] | 100% Completed | 101.10 ms: </pre>
[########################################] | 100% Completed | 101.10 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.10 ms

[10]:

array([[119, 332],
       [121, 329],
       [127, 321],
       [123, 510],
       [487, 108],
       [440,   1],
       [341, 117]])

../_images/tutorials_dev_xr_hvplot_holoviews_11_10.png

[11]:

masks = [dmask(daimg[t]) for t in range(4)]

[ ] | 0% Completed | 144.58 us: </pre>
[ ] | 0% Completed | 144.58 us: end{sphinxVerbatim}

[ ] | 0% Completed | 144.58 us

[########################################] | 100% Completed | 101.14 ms: </pre>
[########################################] | 100% Completed | 101.14 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.14 ms

[ ] | 0% Completed | 130.75 us: </pre>
[ ] | 0% Completed | 130.75 us: end{sphinxVerbatim}

[ ] | 0% Completed | 130.75 us

[########################################] | 100% Completed | 102.54 ms: </pre>
[########################################] | 100% Completed | 102.54 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 102.54 ms

[ ] | 0% Completed | 115.42 us: </pre>
[ ] | 0% Completed | 115.42 us: end{sphinxVerbatim}

[ ] | 0% Completed | 115.42 us

[########################################] | 100% Completed | 101.30 ms: </pre>
[########################################] | 100% Completed | 101.30 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.30 ms

[ ] | 0% Completed | 155.00 us: </pre>
[ ] | 0% Completed | 155.00 us: end{sphinxVerbatim}

[ ] | 0% Completed | 155.00 us

[########################################] | 100% Completed | 101.82 ms: </pre>
[########################################] | 100% Completed | 101.82 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.82 ms

[ ] | 0% Completed | 146.71 us: </pre>
[ ] | 0% Completed | 146.71 us: end{sphinxVerbatim}

[ ] | 0% Completed | 146.71 us

[########################################] | 100% Completed | 102.54 ms: </pre>
[########################################] | 100% Completed | 102.54 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 102.54 ms

[ ] | 0% Completed | 136.44 us: </pre>
[ ] | 0% Completed | 136.44 us: end{sphinxVerbatim}

[ ] | 0% Completed | 136.44 us

[########################################] | 100% Completed | 101.10 ms: </pre>
[########################################] | 100% Completed | 101.10 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.10 ms

[ ] | 0% Completed | 139.83 us: </pre>
[ ] | 0% Completed | 139.83 us: end{sphinxVerbatim}

[ ] | 0% Completed | 139.83 us

[########################################] | 100% Completed | 101.94 ms: </pre>
[########################################] | 100% Completed | 101.94 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.94 ms

[ ] | 0% Completed | 138.37 us: </pre>
[ ] | 0% Completed | 138.37 us: end{sphinxVerbatim}

[ ] | 0% Completed | 138.37 us

[########################################] | 100% Completed | 101.72 ms: </pre>
[########################################] | 100% Completed | 101.72 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.72 ms

[ ] | 0% Completed | 142.30 us: </pre>
[ ] | 0% Completed | 142.30 us: end{sphinxVerbatim}

[ ] | 0% Completed | 142.30 us

[########################################] | 100% Completed | 101.55 ms: </pre>
[########################################] | 100% Completed | 101.55 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.55 ms

[ ] | 0% Completed | 150.93 us: </pre>
[ ] | 0% Completed | 150.93 us: end{sphinxVerbatim}

[ ] | 0% Completed | 150.93 us

[########################################] | 100% Completed | 102.85 ms: </pre>
[########################################] | 100% Completed | 102.85 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 102.85 ms

[ ] | 0% Completed | 156.92 us: </pre>
[ ] | 0% Completed | 156.92 us: end{sphinxVerbatim}

[ ] | 0% Completed | 156.92 us

[########################################] | 100% Completed | 102.84 ms: </pre>
[########################################] | 100% Completed | 102.84 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 102.84 ms

[ ] | 0% Completed | 138.11 us: </pre>
[ ] | 0% Completed | 138.11 us: end{sphinxVerbatim}

[ ] | 0% Completed | 138.11 us

[########################################] | 100% Completed | 101.10 ms: </pre>
[########################################] | 100% Completed | 101.10 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.10 ms

[12]:

masks = np.stack(masks)
masks.shape

[12]:

(4, 512, 512)

[13]:

tff.imshow(masks)

[13]:

(<Figure size 988.8x604.8 with 3 Axes>,
 <Axes: >,
 <matplotlib.image.AxesImage at 0x7fb01ca608d0>)

../_images/tutorials_dev_xr_hvplot_holoviews_14_1.png

[14]:

distance = ndimage.distance_transform_edt(masks)

distance = skimage.filters.gaussian(distance, sigma=5)

[15]:

import impy

impy.array(distance).imshow()

../_images/tutorials_dev_xr_hvplot_holoviews_16_0.png

[15]:

name	No name
shape	4(t), 512(y), 512(x)
label shape	No label
dtype	float32
source	None
scale	ScaleView(t=1.0000px, y=1.0000px, x=1.0000px)

[16]:

for t in range(4):
    coords = skimage.feature.peak_local_max(distance[t], footprint=np.ones((130, 130)))
    print(coords)

[[497 108]
 [445   2]
 [113 345]
 [344 110]]
[[497 108]
 [445   2]
 [113 345]
 [344 110]]
[[497 108]
 [445   2]
 [113 345]
 [344 110]]
[[497 108]
 [445   2]
 [113 345]
 [344 110]]

[17]:

co = np.stack([coords, coords, coords, coords])

[18]:

coords.T

[18]:

array([[497, 445, 113, 344],
       [108,   2, 345, 110]])

[19]:

mm = np.zeros(masks[0].shape, dtype=bool)
mm[tuple(co.T)] = True
# markers, _ = ndimage.label(mm)
markers = skimage.measure.label(np.stack([mm, mm, mm, mm]))

labels = skimage.segmentation.watershed(-distance, markers, mask=masks)

_, (ax1, ax2) = plt.subplots(1, 2)
ax1.imshow(labels[3])
ax2.imshow(labels[3] == 4)

[19]:

<matplotlib.image.AxesImage at 0x7fb01cc6c950>

../_images/tutorials_dev_xr_hvplot_holoviews_20_1.png

[20]:

img = tff.imread(fp)

[21]:

dim, _, _ = nima.read_tiff(fp, channels=["R", "G", "C"])

(4, 512, 512)

[22]:

res = nima.d_bg(dim)
bgs = res[1]

[23]:

def ratio(t, roi):
    g = img[t, 0][labels[t] == roi].mean() - bgs["G"][t]
    r = img[t, 1][labels[t] == roi].mean() - bgs["R"][t]
    c = img[t, 2][labels[t] == roi].mean() - bgs["C"][t]
    return g / c, c / r


ratio(1, 4)

[23]:

(6.347092058582823, 1.312316699431818)

[24]:

rph = defaultdict(list)
rcl = defaultdict(list)
for roi in range(1, 5):
    for t in range(4):
        ph, cl = ratio(t, roi)
        rph[roi].append(ph)
        rcl[roi].append(cl)

plt.plot(rph[1])
plt.plot(rph[2])
plt.plot(rph[3])
plt.plot(rph[4])

[24]:

[<matplotlib.lines.Line2D at 0x7fb00f033790>]

../_images/tutorials_dev_xr_hvplot_holoviews_25_1.png

[25]:

plt.plot(rcl[1])
plt.plot(rcl[2])
plt.plot(rcl[3])
plt.plot(rcl[4])

[25]:

[<matplotlib.lines.Line2D at 0x7fb00f0fdd90>]

../_images/tutorials_dev_xr_hvplot_holoviews_26_1.png

[26]:

t = 2
mask = dmask(daimg[t])
# skimage.io.imshow(mask)
lab, nlab = ndimage.label(mask)
lab[~mask] = -1
# lab[lab==1] = -1
labels_ws = skimage.segmentation.random_walker(
    daimg[t, 1].compute(), lab, beta=1e10, mode="bf"
)
# labels_ws = skimage.segmentation.random_walker(-distance, lab, beta=10000, mode="bf")

_, (ax1, ax2) = plt.subplots(1, 2)
ax1.imshow(labels_ws)
ax2.imshow(labels_ws == 2)

[ ] | 0% Completed | 152.02 us: </pre>
[ ] | 0% Completed | 152.02 us: end{sphinxVerbatim}

[ ] | 0% Completed | 152.02 us

[########################################] | 100% Completed | 101.24 ms: </pre>
[########################################] | 100% Completed | 101.24 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.24 ms

[ ] | 0% Completed | 147.40 us: </pre>
[ ] | 0% Completed | 147.40 us: end{sphinxVerbatim}

[ ] | 0% Completed | 147.40 us

[########################################] | 100% Completed | 101.20 ms: </pre>
[########################################] | 100% Completed | 101.20 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.20 ms

[ ] | 0% Completed | 142.19 us: </pre>
[ ] | 0% Completed | 142.19 us: end{sphinxVerbatim}

[ ] | 0% Completed | 142.19 us

[########################################] | 100% Completed | 101.44 ms: </pre>
[########################################] | 100% Completed | 101.44 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.44 ms

[ ] | 0% Completed | 180.37 us: </pre>
[ ] | 0% Completed | 180.37 us: end{sphinxVerbatim}

[ ] | 0% Completed | 180.37 us

[########################################] | 100% Completed | 101.79 ms: </pre>
[########################################] | 100% Completed | 101.79 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 101.79 ms

/home/docs/checkouts/readthedocs.org/user_builds/nima/envs/tut2/lib/python3.11/site-packages/skimage/segmentation/random_walker_segmentation.py:475: UserWarning: Random walker only segments unlabeled areas, where labels == 0. No zero valued areas in labels were found. Returning provided labels.
  inds_isolated_seeds, isolated_values) = _preprocess(labels)

[26]:

<matplotlib.image.AxesImage at 0x7fb00f3233d0>

../_images/tutorials_dev_xr_hvplot_holoviews_27_14.png

[27]:

imar = impy.imread(fp)

imar.label_threshold()

[ ] | 0% Completed | 149.79 us: </pre>
[ ] | 0% Completed | 149.79 us: end{sphinxVerbatim}

[ ] | 0% Completed | 149.79 us

[########################################] | 100% Completed | 102.34 ms: </pre>
[########################################] | 100% Completed | 102.34 ms: end{sphinxVerbatim}

[########################################] | 100% Completed | 102.34 ms

[27]:

name	1b_c16_15.tif
shape	4(t), 3(c), 512(y), 512(x)
dtype	uint16
source	../../tests/data/1b_c16_15.tif
scale	ScaleView(t=1.0000px, c=1.0000px, y=0.2000px, x=0.2000px)

[28]:

imar[:, 2].imshow(label=1)

../_images/tutorials_dev_xr_hvplot_holoviews_29_0.png

[28]:

name	1b_c16_15.tif
shape	4(t), 512(y), 512(x)
label shape	4(t), 512(y), 512(x)
dtype	uint16
source	../../tests/data/1b_c16_15.tif
scale	ScaleView(t=1.0000px, y=0.2000px, x=0.2000px)

[29]:

def dmask0(im, erf_pvalue=1e-100, size=10):
    p = utils.prob(im[0], *utils.bg(im[0]))
    for img in im[1:]:
        p = p * utils.prob(img, *utils.bg(img))
    p = ndimage.median_filter((p) ** (1 / len(im)), size=size)
    mask = p < erf_pvalue
    return skimage.morphology.remove_small_objects(mask)

[30]:

dmask0(imar[1])

[30]:

array([[False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       ...,
       [ True,  True,  True, ..., False, False, False],
       [ True,  True,  True, ..., False, False, False],
       [ True,  True,  True, ..., False, False, False]])

[31]:

plt.imshow(skimage.measure.label(mask))

[31]:

<matplotlib.image.AxesImage at 0x7fb00ed28890>

../_images/tutorials_dev_xr_hvplot_holoviews_32_1.png

[32]:

distance = skimage.filters.gaussian(distance, sigma=30)
tff.imshow(distance)

[32]:

(<Figure size 988.8x604.8 with 1 Axes>,
 <Axes: >,
 <matplotlib.image.AxesImage at 0x7fb00f0547d0>)

../_images/tutorials_dev_xr_hvplot_holoviews_33_1.png

[33]:

np.transpose(np.nonzero(skimage.morphology.local_maxima(distance)))

[33]:

array([[  0, 446,   0],
       [  3, 110, 352],
       [  3, 347, 107],
       [  3, 511, 108]])

[34]:

tff.imshow(ndimage.label(mask)[0])

[34]:

(<Figure size 988.8x604.8 with 2 Axes>,
 <Axes: >,
 <matplotlib.image.AxesImage at 0x7fb00db7b910>)

../_images/tutorials_dev_xr_hvplot_holoviews_35_1.png

[35]:

res[1]

[35]:

	R	G	C
0	459.0	245.0	272.0
1	462.0	249.0	275.0
2	463.0	250.0	277.0
3	460.0	248.0	273.0

[36]:

res[2]["G"][2][0]

[36]:

../_images/tutorials_dev_xr_hvplot_holoviews_37_0.png

[37]:

res[1].plot()

[37]:

<Axes: >

../_images/tutorials_dev_xr_hvplot_holoviews_38_1.png

[38]:

import hvplot.pandas

[39]:

res[1].hvplot()

[39]:

[40]:

import xarray as xr

[41]:

xim = xr.DataArray(
    data=[dim["G"], dim["R"], dim["C"]],
    dims=["channel", "time", "y", "x"],
    coords=dict(
        channel=["Green", "Red", "Cyan"], time=[0, 1, 2, 3], y=range(512), x=range(512)
    ),
)

[42]:

import hvplot.xarray

[43]:

xim.sel(time=0, channel="Green").hvplot(width=400, height=300)

[43]:

[44]:

xim.sel(time=0).hvplot(
    width=300,
    subplots=True,
    by="channel",
)

[44]:

[45]:

hvplot.extension(
    "bokeh",
    "matplotlib",
)

[46]:

img = xim.sel(time=0).sel(channel="Red")

[47]:

hvimg = hv.Image(img)

[48]:

%%opts Image style(cmap='viridis')
#%%opts Image [aspect=1388/1038]
#%%output size=300
f = xim.hvplot(frame_width=170,
               subplots=True, row="channel", col="time", yaxis=False, colorbar=False, xaxis=False)
f

[48]:

[49]:

hv.save(f, "b.png", dpi=200)

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[49], line 1
----> 1 hv.save(f, "b.png", dpi=200)

File ~/checkouts/readthedocs.org/user_builds/nima/envs/tut2/lib/python3.11/site-packages/holoviews/util/__init__.py:807, in save(obj, filename, fmt, backend, resources, toolbar, title, **kwargs)
    805     if formats[-1] in supported:
    806         filename = '.'.join(formats[:-1])
--> 807 return renderer_obj.save(obj, filename, fmt=fmt, resources=resources,
    808                          title=title)

File ~/checkouts/readthedocs.org/user_builds/nima/envs/tut2/lib/python3.11/site-packages/holoviews/plotting/renderer.py:599, in Renderer.save(self_or_cls, obj, basename, fmt, key, info, options, resources, title, **kwargs)
    596     plot.layout.save(basename, embed=True, resources=resources, title=title)
    597     return
--> 599 rendered = self_or_cls(plot, fmt)
    600 if rendered is None: return
    601 (data, info) = rendered

File ~/checkouts/readthedocs.org/user_builds/nima/envs/tut2/lib/python3.11/site-packages/holoviews/plotting/renderer.py:199, in Renderer.__call__(self, obj, fmt, **kwargs)
    197     return self.static_html(plot), info
    198 else:
--> 199     data = self._figure_data(plot, fmt, **kwargs)
    200     data = self._apply_post_render_hooks(data, obj, fmt)
    201     return data, info

File ~/checkouts/readthedocs.org/user_builds/nima/envs/tut2/lib/python3.11/site-packages/holoviews/plotting/bokeh/renderer.py:125, in BokehRenderer._figure_data(self, plot, fmt, doc, as_script, **kwargs)
    123 elif fmt == 'png':
    124     from bokeh.io.export import get_screenshot_as_png
--> 125     img = get_screenshot_as_png(plot.state, driver=state.webdriver)
    126     imgByteArr = BytesIO()
    127     img.save(imgByteArr, format='PNG')

File ~/checkouts/readthedocs.org/user_builds/nima/envs/tut2/lib/python3.11/site-packages/bokeh/io/export.py:274, in get_screenshot_as_png(obj, driver, timeout, resources, width, height, scale_factor, state)
    271         raise ValueError(f'Expected the web driver to have a device pixel ratio greater than {scale_factor}. '
    272                          f'Was given a web driver with a device pixel ratio of {device_pixel_ratio}.')
    273 else:
--> 274     web_driver = webdriver_control.get(scale_factor=scale_factor)
    275 web_driver.maximize_window()
    276 web_driver.get(f"file://{tmp.path}")

File ~/checkouts/readthedocs.org/user_builds/nima/envs/tut2/lib/python3.11/site-packages/bokeh/io/webdriver.py:176, in _WebdriverState.get(self, scale_factor)
    173 if not self.reuse or self.current is None or not scale_factor_less_than_web_driver_device_pixel_ratio(
    174         scale_factor, self.current):
    175     self.reset()
--> 176     self.current = self.create(scale_factor=scale_factor)
    177 return self.current

File ~/checkouts/readthedocs.org/user_builds/nima/envs/tut2/lib/python3.11/site-packages/bokeh/io/webdriver.py:180, in _WebdriverState.create(self, kind, scale_factor)
    179 def create(self, kind: DriverKind | None = None, scale_factor: float = 1) -> WebDriver:
--> 180     driver = self._create(kind, scale_factor=scale_factor)
    181     self._drivers.add(driver)
    182     return driver

File ~/checkouts/readthedocs.org/user_builds/nima/envs/tut2/lib/python3.11/site-packages/bokeh/io/webdriver.py:198, in _WebdriverState._create(self, kind, scale_factor)
    195         self.kind = "firefox"
    196         return driver
--> 198     raise RuntimeError("Neither firefox and geckodriver nor a variant of chromium browser and " \
    199                        "chromedriver are available on system PATH. You can install the former " \
    200                        "with 'conda install -c conda-forge firefox geckodriver'.")
    201 elif driver_kind == "chromium":
    202     return create_chromium_webdriver(scale_factor=scale_factor)

RuntimeError: Neither firefox and geckodriver nor a variant of chromium browser and chromedriver are available on system PATH. You can install the former with 'conda install -c conda-forge firefox geckodriver'.

[50]:

# %%opts Image [aspect=1388/1038]

f = xim.sel(channel="Red").hvplot(
    frame_width=300,
    frame_height=200,
    subplots=True,
    col="time",
    yaxis=False,
    colorbar=False,
    xaxis=False,
    cmap="Reds",
) + xim.sel(channel="Cyan").hvplot(
    subplots=True, col="time", yaxis=False, colorbar=False, xaxis=False, cmap="Greens"
)
f

[50]:

[51]:

import aicsimageio

aicsimageio.__version__

[51]:

'4.14.0'

[52]:

reader = aicsimageio.readers.tiff_reader.TiffReader
aim1 = aicsimageio.AICSImage(
    "/home/dati/dt-evolv/data/2022-06-17/images/Vero-Hek/2022-06-14/13080/TimePoint_1/6w-20Xph1-SpikeTest3_A02_s570_w14510D534-71A3-4EB5-B48F-F4331FE96517.tif",
    reader=reader,
)
aim2 = aicsimageio.AICSImage(
    "/home/dati/dt-evolv/data/2022-06-17/images/Vero-Hek/2022-06-14/13080/TimePoint_1/6w-20Xph1-SpikeTest3_A02_s570_w25049D5AC-5888-492F-891D-8BECC1AB67DF.tif",
    reader=reader,
)

---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
Cell In[52], line 2
      1 reader = aicsimageio.readers.tiff_reader.TiffReader
----> 2 aim1 = aicsimageio.AICSImage(
      3     "/home/dati/dt-evolv/data/2022-06-17/images/Vero-Hek/2022-06-14/13080/TimePoint_1/6w-20Xph1-SpikeTest3_A02_s570_w14510D534-71A3-4EB5-B48F-F4331FE96517.tif",
      4     reader=reader,
      5 )
      6 aim2 = aicsimageio.AICSImage(
      7     "/home/dati/dt-evolv/data/2022-06-17/images/Vero-Hek/2022-06-14/13080/TimePoint_1/6w-20Xph1-SpikeTest3_A02_s570_w25049D5AC-5888-492F-891D-8BECC1AB67DF.tif",
      8     reader=reader,
      9 )

File ~/checkouts/readthedocs.org/user_builds/nima/envs/tut2/lib/python3.11/site-packages/aicsimageio/aics_image.py:283, in AICSImage.__init__(self, image, reader, reconstruct_mosaic, fs_kwargs, **kwargs)
    280     ReaderClass = reader
    282 # Init and store reader
--> 283 self._reader = ReaderClass(image, fs_kwargs=fs_kwargs, **kwargs)
    285 # Store delayed modifiers
    286 self._reconstruct_mosaic = reconstruct_mosaic

File ~/checkouts/readthedocs.org/user_builds/nima/envs/tut2/lib/python3.11/site-packages/aicsimageio/readers/tiff_reader.py:83, in TiffReader.__init__(self, image, chunk_dims, dim_order, channel_names, fs_kwargs, **kwargs)
     73 def __init__(
     74     self,
     75     image: types.PathLike,
   (...)
     81 ):
     82     # Expand details of provided image
---> 83     self._fs, self._path = io_utils.pathlike_to_fs(
     84         image,
     85         enforce_exists=True,
     86         fs_kwargs=fs_kwargs,
     87     )
     89     # Store params
     90     if isinstance(chunk_dims, str):

File ~/checkouts/readthedocs.org/user_builds/nima/envs/tut2/lib/python3.11/site-packages/aicsimageio/utils/io_utils.py:56, in pathlike_to_fs(uri, enforce_exists, fs_kwargs)
     54 if enforce_exists:
     55     if not fs.exists(path):
---> 56         raise FileNotFoundError(f"{fs.protocol}://{path}")
     58 # Get and store details
     59 # We do not return an AbstractBufferedFile (i.e. fs.open) as we do not want to have
     60 # any open file buffers _after_ any API call. API calls must themselves call
     61 # fs.open and complete their function during the context of the opened buffer.
     62 return fs, path

FileNotFoundError: file:///home/dati/dt-evolv/data/2022-06-17/images/Vero-Hek/2022-06-14/13080/TimePoint_1/6w-20Xph1-SpikeTest3_A02_s570_w14510D534-71A3-4EB5-B48F-F4331FE96517.tif

[53]:

x1 = aim1.xarray_data
x2 = aim2.xarray_data

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[53], line 1
----> 1 x1 = aim1.xarray_data
      2 x2 = aim2.xarray_data

NameError: name 'aim1' is not defined

[54]:

# Create a new Dataset with new coordinates
ds = xr.Dataset({"c1": x1, "c2": x2})

# Assuming ds is your Dataset
new_coords = {"Frame": [1, 2], "excitation_wavelength": [400, 500]}

# Use assign_coords to set new coordinates
ds_assigned = ds.assign_coords(**new_coords)

ds_assigned

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[54], line 2
      1 # Create a new Dataset with new coordinates
----> 2 ds = xr.Dataset({"c1": x1, "c2": x2})
      4 # Assuming ds is your Dataset
      5 new_coords = {"Frame": [1, 2], "excitation_wavelength": [400, 500]}

NameError: name 'x1' is not defined

[55]:

aim2.metadata[220:230] == aim1.metadata[220:230]

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[55], line 1
----> 1 aim2.metadata[220:230] == aim1.metadata[220:230]

NameError: name 'aim2' is not defined

[56]:

im = x1.to_numpy()[0, 0, 0]

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[56], line 1
----> 1 im = x1.to_numpy()[0, 0, 0]

NameError: name 'x1' is not defined

[57]:

im1 = tff.imread("/home/dati/dt-evolv/data/2022-06-17/flat_w1.tif")
im2 = tff.imread("/home/dati/dt-evolv/data/2022-06-17/flat_w2.tif")

---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
/tmp/ipykernel_715/2721790190.py in ?()
----> 1 im1 = tff.imread("/home/dati/dt-evolv/data/2022-06-17/flat_w1.tif")
      2 im2 = tff.imread("/home/dati/dt-evolv/data/2022-06-17/flat_w2.tif")

~/checkouts/readthedocs.org/user_builds/nima/envs/tut2/lib/python3.11/site-packages/tifffile/tifffile.py in ?(files, aszarr, key, series, level, squeeze, maxworkers, mode, name, offset, size, pattern, axesorder, categories, imread, sort, container, chunkshape, dtype, axestiled, ioworkers, chunkmode, fillvalue, zattrs, multiscales, omexml, out, out_inplace, _multifile, _useframes, **kwargs)
   1036
   1037         if isinstance(files, str) or not isinstance(
   1038             files, collections.abc.Sequence
   1039         ):
-> 1040             with TiffFile(
   1041                 files,
   1042                 mode=mode,
   1043                 name=name,

~/checkouts/readthedocs.org/user_builds/nima/envs/tut2/lib/python3.11/site-packages/tifffile/tifffile.py in ?(self, file, mode, name, offset, size, omexml, _multifile, _useframes, _parent, **is_flags)
   3929
   3930         if mode not in (None, 'r', 'r+', 'rb', 'r+b'):
   3931             raise ValueError(f'invalid mode {mode!r}')
   3932
-> 3933         fh = FileHandle(file, mode=mode, name=name, offset=offset, size=size)
   3934         self._fh = fh
   3935         self._multifile = True if _multifile is None else bool(_multifile)
   3936         self._files = {fh.name: self}

~/checkouts/readthedocs.org/user_builds/nima/envs/tut2/lib/python3.11/site-packages/tifffile/tifffile.py in ?(self, file, mode, name, offset, size)
  13631         self._offset = -1 if offset is None else offset
  13632         self._size = -1 if size is None else size
  13633         self._close = True
  13634         self._lock = NullContext()
> 13635         self.open()
  13636         assert self._fh is not None

~/checkouts/readthedocs.org/user_builds/nima/envs/tut2/lib/python3.11/site-packages/tifffile/tifffile.py in ?(self)
  13646         if isinstance(self._file, str):
  13647             # file name
  13648             self._file = os.path.realpath(self._file)
  13649             self._dir, self._name = os.path.split(self._file)
> 13650             self._fh = open(self._file, self._mode)  # type: ignore
  13651             self._close = True
  13652             if self._offset < 0:
  13653                 self._offset = 0

FileNotFoundError: [Errno 2] No such file or directory: '/home/dati/dt-evolv/data/2022-06-17/flat_w1.tif'

[58]:

%%opts Image [aspect=1388/1038]

%%opts Image.Cyan style(cmap=plt.cm.Blues)
%%opts Image.Green style(cmap=plt.cm.Greens)
%%opts Image.Red style(cmap=plt.cm.Reds)

[59]:

chans = (
    hv.Image(dim["C"][0], group="cyan")
    + hv.Image(dim["G"][2], group="green")
    + hv.Image(dim["R"][1], group="red")
)

chans

[59]:

[60]:

hv.save(chans, "a.png")

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[60], line 1
----> 1 hv.save(chans, "a.png")

File ~/checkouts/readthedocs.org/user_builds/nima/envs/tut2/lib/python3.11/site-packages/holoviews/util/__init__.py:807, in save(obj, filename, fmt, backend, resources, toolbar, title, **kwargs)
    805     if formats[-1] in supported:
    806         filename = '.'.join(formats[:-1])
--> 807 return renderer_obj.save(obj, filename, fmt=fmt, resources=resources,
    808                          title=title)

File ~/checkouts/readthedocs.org/user_builds/nima/envs/tut2/lib/python3.11/site-packages/holoviews/plotting/renderer.py:599, in Renderer.save(self_or_cls, obj, basename, fmt, key, info, options, resources, title, **kwargs)
    596     plot.layout.save(basename, embed=True, resources=resources, title=title)
    597     return
--> 599 rendered = self_or_cls(plot, fmt)
    600 if rendered is None: return
    601 (data, info) = rendered

File ~/checkouts/readthedocs.org/user_builds/nima/envs/tut2/lib/python3.11/site-packages/holoviews/plotting/renderer.py:199, in Renderer.__call__(self, obj, fmt, **kwargs)
    197     return self.static_html(plot), info
    198 else:
--> 199     data = self._figure_data(plot, fmt, **kwargs)
    200     data = self._apply_post_render_hooks(data, obj, fmt)
    201     return data, info

File ~/checkouts/readthedocs.org/user_builds/nima/envs/tut2/lib/python3.11/site-packages/holoviews/plotting/bokeh/renderer.py:125, in BokehRenderer._figure_data(self, plot, fmt, doc, as_script, **kwargs)
    123 elif fmt == 'png':
    124     from bokeh.io.export import get_screenshot_as_png
--> 125     img = get_screenshot_as_png(plot.state, driver=state.webdriver)
    126     imgByteArr = BytesIO()
    127     img.save(imgByteArr, format='PNG')

File ~/checkouts/readthedocs.org/user_builds/nima/envs/tut2/lib/python3.11/site-packages/bokeh/io/export.py:274, in get_screenshot_as_png(obj, driver, timeout, resources, width, height, scale_factor, state)
    271         raise ValueError(f'Expected the web driver to have a device pixel ratio greater than {scale_factor}. '
    272                          f'Was given a web driver with a device pixel ratio of {device_pixel_ratio}.')
    273 else:
--> 274     web_driver = webdriver_control.get(scale_factor=scale_factor)
    275 web_driver.maximize_window()
    276 web_driver.get(f"file://{tmp.path}")

File ~/checkouts/readthedocs.org/user_builds/nima/envs/tut2/lib/python3.11/site-packages/bokeh/io/webdriver.py:176, in _WebdriverState.get(self, scale_factor)
    173 if not self.reuse or self.current is None or not scale_factor_less_than_web_driver_device_pixel_ratio(
    174         scale_factor, self.current):
    175     self.reset()
--> 176     self.current = self.create(scale_factor=scale_factor)
    177 return self.current

File ~/checkouts/readthedocs.org/user_builds/nima/envs/tut2/lib/python3.11/site-packages/bokeh/io/webdriver.py:180, in _WebdriverState.create(self, kind, scale_factor)
    179 def create(self, kind: DriverKind | None = None, scale_factor: float = 1) -> WebDriver:
--> 180     driver = self._create(kind, scale_factor=scale_factor)
    181     self._drivers.add(driver)
    182     return driver

File ~/checkouts/readthedocs.org/user_builds/nima/envs/tut2/lib/python3.11/site-packages/bokeh/io/webdriver.py:198, in _WebdriverState._create(self, kind, scale_factor)
    195         self.kind = "firefox"
    196         return driver
--> 198     raise RuntimeError("Neither firefox and geckodriver nor a variant of chromium browser and " \
    199                        "chromedriver are available on system PATH. You can install the former " \
    200                        "with 'conda install -c conda-forge firefox geckodriver'.")
    201 elif driver_kind == "chromium":
    202     return create_chromium_webdriver(scale_factor=scale_factor)

RuntimeError: Neither firefox and geckodriver nor a variant of chromium browser and chromedriver are available on system PATH. You can install the former with 'conda install -c conda-forge firefox geckodriver'.

5. Holoviews#

[61]:

hv.notebook_extension()
cm = plt.cm.inferno_r
channels = ["G", "R", "C"]
dim, n_ch, times = nima.read_tiff(fp, channels)

dimm = nima.d_median(dim)
f = nima.d_show(dimm, cmap=cm)

(4, 512, 512)

[62]:

%%opts Image [aspect=512/512]

%%opts Image.Cyan style(cmap=plt.cm.Blues)
%%opts Image.Green style(cmap=plt.cm.Greens)
%%opts Image.Red style(cmap=plt.cm.Reds)

chans = hv.Image(dim['C'][0], group='cyan') \
    + hv.Image(dim['G'][0], group='green') \
    + hv.Image(dim['R'][0], group='red')

chans

[62]:

[63]:

c = [(i, hv.Image(im)) for i, im in enumerate(dim["C"])]
c = hv.HoloMap(c, kdims=["Frame"])
g = [(i, hv.Image(im)) for i, im in enumerate(dim["G"])]
g = hv.HoloMap(g, kdims=["Frame"])
r = [(i, hv.Image(im)) for i, im in enumerate(dim["R"])]
r = hv.HoloMap(r, kdims=["Frame"])

[64]:

%%output holomap='auto'
%%opts Image style(cmap='viridis')
(c + g).select(Frame={0,5,6,7,10,30}).cols(2)

[64]:

[65]:

c[::20].overlay("Frame")

[65]:

[66]:

wl = hv.Dimension("excitation wavelength", unit="nm")
c = c.add_dimension(wl, 1, 458)
g = g.add_dimension(wl, 1, 488)
r = r.add_dimension(wl, 1, 561)

channels = c.clone()
channels.update(g)
channels.update(r)

[67]:

%%opts Image style(cmap='viridis')
%%output size=300
channels[::5].grid(['Frame', 'excitation wavelength'])

[67]:

[68]:

t = [(i, hv.Image(im)) for i, im in enumerate(dim["C"])]

[69]:

hv.HoloMap([(i, hv.Image(im)) for i, im in enumerate(dim["C"])], kdims=["frame"])

[69]:

[70]:

hv.NdLayout(
    {
        k: hv.HoloMap(
            [(i, hv.Image(im)) for i, im in enumerate(dim[k])], kdims=["frame"]
        )
        for k in dim
    },
    kdims=["channels"],
)[::4]

[70]:

[71]:

%%opts Image (cmap='viridis')
%%opts Image.A [aspect=2]
im = hv.Image(dim["G"][1], bounds=(0, 0, 512, 512))
im2 = hv.Image(dim['C'][1], bounds=(0, 0, 512, 512))
im3 = hv.Image(dimm['C'][1], bounds=(0, 0, 512, 512))
((im * hv.HLine(y=350)) + im.sample(y=350) + (im2 * hv.HLine(y=150)) + im2.sample(y=150) * im3.sample(y=150)).cols(3)

[71]: