# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# %load_ext autoreload
# %autoreload 2

from dtoolbioimage import ImageDataSet, Image3D, Image
import numpy as np

imageds = ImageDataSet('azure://jicimagedatasets/cc9d757a-b149-4626-83be-abc5a3765cb8')

wall_stack = imageds.get_stack('fca-3_FLC-Venus_root01', 'fca-3_FLC-Venus_root01 #1', 1)

venus_stack = imageds.get_stack('fca-3_FLC-Venus_root01', 'fca-3_FLC-Venus_root01 #1', 0)

from dtoolbioimage.segment import sitk_watershed_segmentation, filter_segmentation_by_size

from dtoolbioimage import zoom_to_match_scales
zoomed_stack = zoom_to_match_scales(wall_stack)

zoomed_stack

raw_segmentation = sitk_watershed_segmentation(zoomed_stack)

raw_segmentation.view(Image3D)

filtered = filter_segmentation_by_size(raw_segmentation)

filtered.view(Image3D)

from dtoolbioimage.segment import Segmentation3D

uci = filtered.view(Segmentation3D).unique_color_image

uci[:,:,10].view(Image)

filtered.view(Segmentation3D).save('root01segmentation.tif')

import os; os.getcwd()

pretty = filtered.view(Segmentation3D).pretty_color_image

pretty[:,:,10].view(Image)

# +
from dtoolbioimage.util.array import pretty_color_array
from ipywidgets import interactive, IntSlider, Output, HBox
from IPython.display import clear_output, display, display_png

def simple_segmentation_viewer(stack):

    _, _, max_z = stack.shape

    slider = IntSlider(min=0, max=max_z, step=1, description='Z plane:')

    def show_z_plane(z):
        display_png(pretty_color_array(stack[:,:,z]).view(Image))
    
    return interactive(show_z_plane, z=slider)



# -

simple_segmentation_viewer(filtered)

from io import BytesIO
from imageio import imsave

png_byte_arrays = []
for z in range(filtered.shape[2]):
    b = BytesIO()
    imsave(b, pretty_color_array(filtered[:,:,z]), 'PNG', compress_level=0)
    png_byte_arrays.append(b)


class CachedImageView(object):
    
    def _repr_png_(self):
        return png_byte_arrays[0].getvalue()


CachedImageView()

# +
from dtoolbioimage.util.array import pretty_color_array
from ipywidgets import interactive, IntSlider, Output, HBox
from IPython.display import clear_output, display, display_png

def cached_segmentation_viewer(stack):

    _, _, max_z = stack.shape
    
    slider = IntSlider(min=0, max=max_z, step=1, description='Z plane:')

    png_byte_arrays = {}
    def show_z_plane(z):
        if z in png_byte_arrays:
            raw_png_rep = png_byte_arrays[z]
        else:
            b = BytesIO()
            imsave(b, pretty_color_array(stack[:,:,z]), 'PNG', compress_level=0)
            raw_png_rep = b.getvalue()
            png_byte_arrays[z] = raw_png_rep

        display({'image/png': raw_png_rep}, raw=True)
    
    return interactive(show_z_plane, z=slider)


# -

cached_segmentation_viewer(filtered)

filtered[np.where(filtered > 50)] = 0

filtered.view(Segmentation3D)

fpath = 'root01segmentation.tif'

from imageio import volread

unique_color_image = volread(fpath)

unique_color_image.shape

print(unique_color_image[18,256,512])

unique_color = unique_color_image[18,256,512]

from dtoolbioimage.util.color import identifier_from_unique_color

identifier_from_unique_color(unique_color)

planes = []
for z in range(unique_color_image.shape[0]):
    pass

zdim, xdim, ydim, _ = unique_color_image.shape
planes = []
for z in range(zdim):
    segmentation = np.zeros((xdim, ydim), dtype=np.uint32)
    segmentation += unique_color_image[z,:,:,2]
    segmentation += unique_color_image[z,:,:,1] * 256
    segmentation += unique_color_image[z,:,:,0] * 256 * 256
    planes.append(segmentation)

np.dstack(planes).view(Segmentation3D)

Segmentation3D.from_file('root01segmentation.tif')

fpath = 'root01segmentation.tif'

unique_color_image = volread(fpath)

unique_color_image.shape

from skimage.measure import regionprops

filtered.view(Segmentation3D)

zoomed_venus_stack = zoom_to_match_scales(venus_stack)

zoomed_venus_stack

from skimage.measure import regionprops

by_label = {r.label: r for r in regionprops(filtered)}

by_label[5].bbox

rmin, cmin, zmin, rmax, cmax, zmax = by_label[5].bbox

selected = filtered[rmin:rmax, cmin:cmax, zmin:zmax]

selected.view(Segmentation3D)

region = selected == 5

from scipy.ndimage.morphology import binary_erosion

eroded = binary_erosion(region)

eroded.view(Image3D)

surface = region ^ eroded

surface.view(Image3D)

np.sum(surface)

np.sum(region)


def spherality(region):
    eroded = binary_erosion(region)
    surface = region ^ eroded
    
    S = np.sum(surface)
    V = np.sum(region)
    
    mult = 4.5 * np.sqrt(np.pi)
    
    return mult * V / np.power(S, 1.5)
    


spherality(region)

np.power(2, 1.5)

from skimage.morphology import ball

ball(20).view(Image3D)

spherality(ball == 1)

ballbool = ball(100) == 1

ballbool.view(Image3D)

eroded = binary_erosion(ballbool)

surface = ballbool ^ eroded

surface.view(Image3D)

spherality(ballbool)

ballbool.shape

(4 / 3) * np.pi * np.power(100, 3)

np.sum(ballbool)

np.sum(surface)

4 * np.pi * 100 * 100

V = (4 / 3) * np.pi * np.power(100, 3)
S = 4 * np.pi * 100 * 100
mult = 6 * np.sqrt(np.pi)
mult * V / np.power(S, 1.5)

from dtoolbioimage.segment import select_region, spherality

r = select_region(filtered, 50)

r.view(Image3D)

wonky = select_region(filtered, 50)

good = select_region(filtered, 8)

spherality(good)

spherality(wonky)

print(np.unique(filtered))

bad = [l for l in range(3, 51) if spherality(select_region(filtered, l)) < 0.5]

bad

select_region(filtered, 3).view(Image3D)

for l in bad:
    filtered[np.where(filtered == l)] = 0

filtered.view(Segmentation3D)

larger_seg = Segmentation3D.from_file('root01segmentation.tif')

larger_seg


def filter_by_sphericity(segmentation):
    
    filtered = segmentation.copy()
    
    ids = set(np.unique(segmentation)) - set([0])

    bad = [l for l in ids if spherality(select_region(segmentation, l)) < 0.5]
    
    for l in bad:
        filtered[np.where(segmentation == l)] = 0
        
    return filtered


new_filtered = filter_by_sphericity(larger_seg)

new_filtered

labels = set(np.unique(new_filtered)) - set([0])

spheralities = [spherality(select_region(new_filtered, l)) for l in labels]

spheralities

eroded_sps = [spherality(binary_erosion(select_region(new_filtered, l))) for l in labels]

for s1, s2, s3 in zip(dilated_sps, spheralities, eroded_sps):
    print(s1, s2, s3)

from scipy.ndimage.morphology import binary_dilation

dilated_sps = [spherality(binary_dilation(select_region(new_filtered, l))) for l in labels]



StarcoderdataJupyterScriptsDedupFiltered