7. dev with dask#

from collections import defaultdict

import dask.array as da
import holoviews as hv
import hvplot
import hvplot.pandas
import hvplot.xarray
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import skimage
import tifffile as tff
from scipy import ndimage

from nima import io, nima, utils

%load_ext autoreload
%autoreload 2

hvplot.pandas.Interactive
hvplot.xarray.Widget

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

daimg = da.from_zarr(tff.imread(fp, aszarr=True))
daimg
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
4 1 512 512 3 
utils.bg(daimg[0, 0].compute())

(np.float64(457.8380994040984), np.float64(48.50287340564226))

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 1 2
0 457.838099 257.010244 289.378226
1 457.295254 259.072941 289.627118
2 457.760167 260.182049 290.268666
3 453.995203 257.189940 285.613624 
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 *= 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)

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/46be0613d488fe974570edf813f72dd5ba293df975e91e99c2790c6d2f4921f0.png

NEXT:

  • make utils.bg and utils.prob work with dask arrays

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 *= 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

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

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

np.float64(195.49311541527658)

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

47

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

array([[122, 329],
       [122, 510],
       [475, 125],
       [341, 116],
       [421,   1]])
../_images/19039d36e0b71dc50ba8b08fa5fe4154654519ce0e72640b745ff702a16362a8.png 
masks = [dmask(daimg[t]) for t in range(4)]

masks = np.stack(masks)
masks.shape

(4, 512, 512)

tff.imshow(masks)

(<Figure size 988.8x604.8 with 3 Axes>,
 <Axes: >,
 <matplotlib.image.AxesImage at 0x7f81f86d9590>)
../_images/aa531aaa6fb2de1d7a2d757e1cbe0faed77ea756ac39a7ff42ddd9829fecd0ec.png 
distance = ndimage.distance_transform_edt(masks)

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

import impy

impy.array(distance).imshow()
../_images/426dad1849827a239f854abdb0667bc0efd466fe8f3d82566ab8e7ce607312cb.png
nameNo name
shape4(t), 512(y), 512(x)
label shapeNo label
dtypefloat64
sourceNone
scaleScaleView(t=1.0000, y=1.0000, x=1.0000)
for t in range(4):
    coords = skimage.feature.peak_local_max(distance[t], footprint=np.ones((130, 130)))
    print(coords)

[[114 346]
 [473 128]
 [344 110]]
[[114 346]
 [473 128]
 [344 110]]
[[114 346]
 [473 128]
 [344 110]]
[[114 346]
 [473 128]
 [344 110]]

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

coords.T

array([[114, 473, 344],
       [346, 128, 110]])

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)

<matplotlib.image.AxesImage at 0x7f81f8204410>
../_images/bcdff561e3271a7bbe2d0947bc0a4217b001fda24638379489565bb07aee91fd.png 
img = tff.imread(fp)

dim = io.read_image(fp, channels=["R", "G", "C"])

/home/runner/work/nima/nima/.venv/lib/python3.14/site-packages/requests/__init__.py:113: RequestsDependencyWarning: urllib3 (2.6.3) or chardet (7.0.1)/charset_normalizer (3.4.5) doesn't match a supported version!
  warnings.warn(

/home/runner/work/nima/nima/src/nima/io.py:53: UserWarning: Channel wavelength validation failed: Expected λ_C < λ_G < λ_R. Got C=458.0nm, G=563.0nm, R=482.0nm.
  return nio.read_image(str(fp), channels=channels)

bg_params = nima.BgParams()
res = nima.bg(dim, bg_params)
bgs = res[1]

def ratio(t, roi):
    if not np.any(labels[t] == roi):
        return np.nan, np.nan
    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)

(nan, nan)

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])

[<matplotlib.lines.Line2D at 0x7f81a7c01550>]
../_images/5bdb3e96ed11490cb3e71fae46b8a8ca0c70b824a09b89123bad02ec8a47ee0b.png 
plt.plot(rcl[1])
plt.plot(rcl[2])
plt.plot(rcl[3])
plt.plot(rcl[4])

[<matplotlib.lines.Line2D at 0x7f81a5bee900>]
../_images/aa2c202fb7b1a30b99e21d7f6dc889d16fc5a3b0881c51b4c0c882ee51111f14.png 
t = 2
mask = dmask(daimg[t])
# skimage.io.imshow(mask)
lab, nlab = ndimage.label(mask)
lab[~mask] = -1
# lab[lab==1] = -1
if np.any(lab == 0):
    labels_ws = skimage.segmentation.random_walker(
        daimg[t, 1].compute(), lab, beta=1e10, mode="bf"
    )
else:
    labels_ws = lab
# 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)

<matplotlib.image.AxesImage at 0x7f81a67865d0>
../_images/fd086c0413a81076e15c0d477fd2c974780638a6eff5bf47910b94d4bec79bba.png 
imar = impy.imread(fp)

imar.label_threshold()
name1b_c16_15.tif
shape4(t), 3(c), 512(y), 512(x)
dtypeuint16
source../../tests/data/1b_c16_15.tif
scaleScaleView(t=1.0000, c=1.0000, y=0.2000, x=0.2000)
imar[:, 2].imshow(label=1)
../_images/cc214f3e784447358b189ee05a9e7e9d473c16a56432966016c84da668850a92.png
name1b_c16_15.tif
shape4(t), 512(y), 512(x)
label shape4(t), 512(y), 512(x)
dtypeuint16
source../../tests/data/1b_c16_15.tif
scaleScaleView(t=1.0000, y=0.2000, x=0.2000)
def dmask0(im, erf_pvalue=1e-100, size=10):
    p = utils.prob(im[0], *utils.bg(im[0]))
    for img in im[1:]:
        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)

dmask0(imar[1])

array([[False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       ...,
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False]], shape=(512, 512))

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

<matplotlib.image.AxesImage at 0x7f81a6b68050>
../_images/cc608c1412ca98c56f1e6e6c9babef46e5b3ef6f59b82a58d9e954c681f91a2b.png 
distance = skimage.filters.gaussian(distance, sigma=30)
tff.imshow(distance)

(<Figure size 988.8x604.8 with 1 Axes>,
 <Axes: >,
 <matplotlib.image.AxesImage at 0x7f81a6b21d10>)
../_images/a8eb61cf9b76ee6bc881d60989f052f561ab367938aa9d0feae664474271c444.png 
np.transpose(np.nonzero(skimage.morphology.local_maxima(distance)))

array([[  3, 110, 353],
       [  3, 346, 107],
       [  3, 456,  18],
       [  3, 486, 128]])

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

(<Figure size 988.8x604.8 with 2 Axes>,
 <Axes: >,
 <matplotlib.image.AxesImage at 0x7f81a48bf890>)
../_images/7701447e18567644e1e645c090f5b5a0bbef1af5cd9ae84cd72ffabb2d446714.png 
res[1]
R G C
0 461.0 257.0 289.0
1 461.0 260.0 290.0
2 461.0 261.0 291.0
3 457.0 258.0 286.0 
res[2]["G"][2][0]
../_images/e6a22dff7c3cb284cea0f9057a3e34d90ed49fe8ad5ba385538beb1402b10355.png 
res[1].plot()

<Axes: >
../_images/ce3b7004bd4c2de32b10279578c724f82f2a9467db3e87eaaee7c091f16af635.png 
res[1].hvplot()

xim = dim
xim.data
Array Chunk
Bytes 6.00 MiB 512.00 kiB
Shape (4, 3, 1, 512, 512) (1, 1, 1, 512, 512)
Dask graph 12 chunks in 44 graph layers
Data type uint16 numpy.ndarray
3 4 512 512 1 
xim.coords["Y"]

<xarray.DataArray 'Y' (Y: 512)> Size: 4kB
array([  0. ,   0.2,   0.4, ..., 101.8, 102. , 102.2], shape=(512,))
Coordinates:
  * Y        (Y) float64 4kB 0.0 0.2 0.4 0.6 0.8 ... 101.6 101.8 102.0 102.2
xim.sel(C="G")[1, 0].hvplot(width=400, height=300)

xim[0, :, 0].hvplot(
    width=300,
    subplots=True,
    by="C",
)

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

img = xim.sel(C="R")[0][0]

hvimg = hv.Image(img)

hvimg

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

f = xim.sel(C="R")[:, 0].hvplot(
    frame_width=300,
    frame_height=200,
    subplots=True,
    col="T",
    yaxis=False,
    colorbar=False,
    xaxis=False,
    cmap="Reds",
) + xim.sel(C="C")[:, 0].hvplot(
    subplots=True, col="T", yaxis=False, colorbar=False, xaxis=False, cmap="Greens"
)
f

import bioio

bioio.__version__

'3.2.0'

hv.opts.defaults(
    hv.opts.Image(aspect=1388 / 1038),
    hv.opts.Image("Cyan", cmap=plt.cm.Blues),
    hv.opts.Image("Green", cmap=plt.cm.Greens),
    hv.opts.Image("Red", cmap=plt.cm.Reds),
)

dim.sel(C="C", Z=0)

<xarray.DataArray 'transpose-c0830e1a0a4ea379f5b5a7d5a76a504e' (T: 4, Y: 512,
                                                                X: 512)> Size: 2MB
dask.array<getitem, shape=(4, 512, 512), dtype=uint16, chunksize=(1, 512, 512), chunktype=numpy.ndarray>
Coordinates:
  * T        (T) float64 32B 13.38 73.38 133.4 193.4
  * Y        (Y) float64 4kB 0.0 0.2 0.4 0.6 0.8 ... 101.6 101.8 102.0 102.2
  * X        (X) float64 4kB 0.0 0.2 0.4 0.6 0.8 ... 101.6 101.8 102.0 102.2
    Z        float64 8B 0.0
    C        <U1 4B 'C'
Attributes:
    unprocessed:   {254: <FILETYPE.PAGE: 2>, 256: 512, 257: 512, 258: 16, 259...
    processed:     experiments=[{'description': '', 'experimenter_ref': {'id'...
    metadata:      Metadata(S=1, T=[4], C=[3], Z=[1], Y=[512], X=[512]\n  bit...
    ome_metadata:  experiments=[{'description': '', 'experimenter_ref': {'id'...
chans = (
    hv.Image(dim.sel(C="C", Z=0)[0], group="cyan")
    + hv.Image(dim.sel(C="G", Z=0)[2], group="green")
    + hv.Image(dim.sel(C="R", Z=0)[1], group="red")
)

chans

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

8. Holoviews#

hv.extension("bokeh")
cm = plt.cm.inferno_r
channels = ["G", "R", "C"]
dim = io.read_image(fp, channels)

dimm = nima.median(dim)
f = nima.plot_img(dimm, cmap=cm)

f
../_images/622403422675338a2ade99399eb330da56ba9077b5d2834ac1be822dc439b5c3.png 
hv.opts.defaults(
    hv.opts.Image(aspect=512 / 512),
    hv.opts.Image("Cyan", cmap=plt.cm.Blues),
    hv.opts.Image("Green", cmap=plt.cm.Greens),
    hv.opts.Image("Red", cmap=plt.cm.Reds),
)

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

chans

nima.plot_img(dimm.sel(Z=0) - dim.sel(Z=0))
../_images/60b1bc0d0b6fe78a65020ea0e0893333da8e94648e59ec93ed25adb9b024f361.png 
c = [(i, hv.Image(im)) for i, im in enumerate(dim.sel(C="C", Z=0))]
c = hv.HoloMap(c, kdims=["Frame"])
g = [(i, hv.Image(im)) for i, im in enumerate(dim.sel(C="G", Z=0))]
g = hv.HoloMap(g, kdims=["Frame"])
r = [(i, hv.Image(im)) for i, im in enumerate(dim.sel(C="R", Z=0))]
r = hv.HoloMap(r, kdims=["Frame"])

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

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

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)

hv.opts.defaults(hv.opts.Image(cmap="viridis", frame_width=300, aspect="equal"))
channels[::5].grid(["Frame", "excitation wavelength"])

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

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

[str(k) for k in dim.C.data]

['G', 'R', 'C']

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

hv.opts.defaults(hv.opts.Image(cmap="viridis"), hv.opts.Image("A", aspect=2))
im = hv.Image(dim.sel(Z=0, C="G")[1])
im2 = hv.Image(dim.sel(Z=0, C="C")[1])
im3 = hv.Image(dimm.sel(Z=0, C="C")[1])
(
    (im * hv.HLine(y=350 * 0.2))
    + im.sample(Y=350 * 0.2)
    + (im2 * hv.VLine(x=30))
    + im.sample(X=76) * im3.sample(X=30)
).cols(3)
On this page

This Page