<jupyter_start><jupyter_text>Fiber segmentation<jupyter_code>import matplotlib.pylab as plt
%matplotlib inline
from itkwidgets import view, compare

from pathlib import Path
import numpy as np
import time

from skimage.morphology import skeletonize
from skimage.filters import threshold_otsu
import itk
import skan

import sys
import os
import warnings

import zarr
import xarray as xr

from dask.distributed import Client, LocalCluster
from fsspec.implementations.http import HTTPFileSystem
import dask_image.ndmeasure
import dask.array as da
import dask
# Use a Coiled.io Dask cloud cluster
# See https://github.com/dani-lbnl/SC20_pyHPC/blob/master/coiled/README.md
use_coiled = False
# use_coiled = True

# If we are on the NERSC cluster and using the NERSC Scratch filesystem
def on_nersc():
    return 'SCRATCH' in os.environ
if on_nersc():
    os.chdir(os.environ['SCRATCH'])
if use_coiled:
    # Option 1: Coiled AWS cluster
    #
    # See https://github.com/dani-lbnl/SC20_pyHPC/blob/master/coiled/README.md
    #
    # You must first log into Coiled.
    import coiled
    # Set to re-use a running cluster when re-running the notebook. Listed at https://cloud.coiled.io/<username>/clusters.
    name = None
    cluster = coiled.Cluster(n_workers=16,
                             worker_cpu=4,
                             worker_memory='12G',
                             name=name,
                             configuration='thewtex/sc20-pyhpc',
                             software='thewtex/sc20-pyhpc')
    
    client = Client(cluster)
    
elif on_nersc():
    # Option 2: NERSC Cori Slurm dask-mpi cluster
    #
    # See SC20_pyHPC/nersc/README.md
    #
    # You must first start the dask cluster (Step 1).
    scheduler_file = os.path.join(os.environ["SCRATCH"], "scheduler.json")
    dask.config.config["distributed"]["dashboard"]["link"] = "{JUPYTERHUB_SERVICE_PREFIX}proxy/{host}:{port}/status"

    client = Client(scheduler_file=scheduler_file)
else:
    # Option 3: Local "cluster"
    #
    local_cluster = LocalCluster(n_workers=2, processes=False, memory_limit='6G')
    client = Client(local_cluster)
    client
if on_nersc():
    # Generated by fibers_to_xarray_zarr.ipynb
    store = './rec20160318_191511_232p3_2cm_cont__4097im_1500ms_ML17keV_6.zarr'
else:
    # We can access our API using fsspec's HTTPFileSystem
    fs = HTTPFileSystem()
    # The http mapper gives us a dict-like interface to the API
    store = fs.get_mapper("https://fiber-bed-zarr.netlify.com/rec20160318_191511_232p3_2cm_cont__4097im_1500ms_ML17keV_6.zarr")
ds = xr.open_zarr(store, consolidated=True)
ds
fibers = ds.rec20160318_191511_232p3_2cm_cont__4097im_1500ms_ML17keV_6
fibers
view(fibers[1000,:,:].values)
view(fibers[950:1050,1000:1200,1000:1200].values)
image = itk.image_view_from_array(fibers[950:1050,1000:1200,1000:1200].values.astype(np.float64))
sigma = 5.25
hessian = itk.hessian_recursive_gaussian_image_filter(image, sigma=sigma)
vesselness_filter = itk.Hessian3DToVesselnessMeasureImageFilter[itk.ctype('float')].New(hessian)
vesselness_filter.Update()
vesselness = vesselness_filter.GetOutput()
compare(image, vesselness, gradient_opacity=0.2, shadow=False, mode='z')
thresh = itk.otsu_threshold_image_filter(vesselness)
inverted = itk.invert_intensity_image_filter(thresh)
rescaled = itk.rescale_intensity_image_filter(inverted, ttype=(type(inverted), itk.Image[itk.UC, 3]))
view(image=rescaled)
label, nlabels = dask_image.ndmeasure.label(np.asarray(rescaled)[50,:,:])
dask_image.ndmeasure.center_of_mass(label, label, da.unique(label).compute()).compute()
view(label_image=label.compute().astype(np.uint16))
threshold_value = threshold_otsu(np.asarray(vesselness))
threshold_value
# To run on a smaller slab
fiber_slab = fibers[950:1150,512:1024,512:1024]
#fiber_slab = fibers
fiber_slab
fiber_slab = fiber_slab.chunk((100,128,128))
fiber_slab
def compute_vesselness(image_chunk, threshold=0.15):
    import itkConfig
    itkConfig.LazyLoading = False
    import itk
    
    image = itk.image_view_from_array(image_chunk.astype(np.float64))

    sigma = 5.25
    hessian = itk.hessian_recursive_gaussian_image_filter(image, sigma=sigma)
    vesselness_filter = itk.Hessian3DToVesselnessMeasureImageFilter[itk.ctype('float')].New(hessian)
    vesselness_filter.Update()
    vesselness = vesselness_filter.GetOutput()
    binary = itk.binary_threshold_image_filter(vesselness, upper_threshold=float(threshold))
    inverted = itk.invert_intensity_image_filter(binary)
    rescaled = itk.rescale_intensity_image_filter(inverted, ttype=(type(inverted), itk.Image[itk.UC, 3]))
    return np.asarray(rescaled)
seg = fiber_slab.data.map_overlap(compute_vesselness,
                                depth=32,
                                trim=True,
                                dtype=np.uint8,
                                threshold=threshold_value)
start = time.time()

seg = seg.compute()

elapsed = time.time() - start
print(elapsed, 'seconds')
label, nlabels = dask_image.ndmeasure.label(seg)
nlabels.compute()
label_slice = label[50,:,:]
label_slice_result = label_slice.compute()
label_slice_result
view(label_image=label_slice_result.astype(np.uint64))
centers = dask_image.ndmeasure.center_of_mass(label_slice_result, label_slice_result, np.unique(label_slice_result)).compute()
centers.shape
import pandas as pd
df = pd.DataFrame(centers[1:,:], columns=('X', 'Y'))
df
df.to_csv('./Results-slice1000-fiber_segmentation.csv')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10589519650655022, 'fraction_numerical': 0.04112081513828239, 'mean_word_length': 4.300867888138862, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 6, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '661970', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2089, 'n_tokens_neox': 1915, 'n_words': 341, 'path': 'code/03_fiber_segmentation.ipynb', 'repo_name': 'dani-lbnl/SC20_pyHPC', 'size': '5496'}