<jupyter_start><jupyter_text>Adaptive Neuron Thresholding  Load data from AWSData is downloaded via cloud-volume. Here, we select the 400th segment from the second brain segment. We scale the default chunk size of 50x50x50 by a factor of 2 in each dimension to download a 100x100x100 voxel volume.<jupyter_code>import brainlit
from brainlit.utils.ngl_pipeline import NeuroglancerSession

url = "s3://open-neurodata/brainlit/brain1"
ngl = NeuroglancerSession(url, mip=6)
img, bounds, voxel = ngl.pull_chunk(2, 100)<jupyter_output>Downloading: 100%|██████████| 1/1 [00:00<00:00,  2.44it/s]
Downloading:   0%|          | 0/1 [00:00<?, ?it/s]<jupyter_text>Prepare plotting functionsWe add some useful plotting functions to view the downloaded volume in matplotlib (2D) and Napari (3D).<jupyter_code># Reference: https://github.com/NeuroDataDesign/mouselit/blob/master/bijan/mouse_test/final%20notebook.ipynb

import matplotlib.pyplot as plt
import SimpleITK as sitk
import numpy as np
import napari
%gui qt


def plot_2d(img, title=None, margin=0.05, dpi=80):
    nda = sitk.GetArrayFromImage(img)
    spacing = img.GetSpacing()

    if nda.ndim == 3:
        c = nda.shape[-1]

        if c not in (3, 4):
            nda = nda[nda.shape[0] // 2, :, :]

    elif nda.ndim == 4:
        c = nda.shape[-1]

        if c not in (3, 4):
            raise RuntimeError("Unable to show 3D-vector Image")

        nda = nda[nda.shape[0] // 2, :, :, :]

    xsize = nda.shape[1] * 2
    ysize = nda.shape[0] * 2

    figsize = (1 + margin) * xsize / dpi, (1 + margin) * ysize / dpi

    plt.figure(figsize=figsize, dpi=dpi, tight_layout=True)
    ax = plt.gca()

    extent = (0, xsize * spacing[0], ysize * spacing[1], 0)

    t = ax.imshow(nda, extent=extent, interpolation=None)

    if nda.ndim == 2:
        t.set_cmap("gray")

    if title:
        plt.title(title)

    plt.show()


def plot_3d(img, xslices=[], yslices=[], zslices=[], title=None, margin=0.05, dpi=80):
    img_xslices = [img[s, :, :] for s in xslices]
    img_yslices = [img[:, s, :] for s in yslices]
    img_zslices = [img[:, :, s] for s in zslices]

    maxlen = max(len(img_xslices), len(img_yslices), len(img_zslices))

    img_null = sitk.Image([0, 0], img.GetPixelID(), img.GetNumberOfComponentsPerPixel())

    img_slices = []
    d = 0

    if len(img_xslices):
        img_slices += img_xslices + [img_null] * (maxlen - len(img_xslices))
        d += 1

    if len(img_yslices):
        img_slices += img_yslices + [img_null] * (maxlen - len(img_yslices))
        d += 1

    if len(img_zslices):
        img_slices += img_zslices + [img_null] * (maxlen - len(img_zslices))
        d += 1

    if maxlen != 0:
        if img.GetNumberOfComponentsPerPixel() == 1:
            img = sitk.Tile(img_slices, [maxlen, d])
        else:
            img_comps = []
            for i in range(0, img.GetNumberOfComponentsPerPixel()):
                img_slices_c = [sitk.VectorIndexSelectionCast(s, i) for s in img_slices]
                img_comps.append(sitk.Tile(img_slices_c, [maxlen, d]))
            img = sitk.Compose(img_comps)

    plot_2d(img, title, margin, dpi)


def napari_viewer(img, labels=None, label_name="Segmentation"):
    viewer = napari.view_image(np.squeeze(np.array(img)))
    if labels is not None:
        viewer.add_labels(labels, name=label_name)
    return viewer<jupyter_output><empty_output><jupyter_text>Let's take a look at the downloaded volume. Napari will open in a new window.<jupyter_code>print(np.sum(img))
print(img.shape)

plot_3d(sitk.GetImageFromArray(np.squeeze(img), isVector=False), zslices=range(48,53), title = "Downloaded Mouselight Volume")
napari_viewer(img)
# viewer = napari.view_image(np.squeeze(np.array(img)))
from skimage.data import astronaut

# create the viewer and display the image
viewer = napari.view_image(astronaut(), rgb=True)<jupyter_output><empty_output><jupyter_text>Generate SegmentationsNext, we will generate segmentations using the following algorithms: \Gaussian Mixture Model Thresholding \Otsu\Connected Threshold\Confidence-Connected Threshold\Neighborhood-Connected Threshold\Level-Set Segmentation\Fast-Marching SegmentationThe segmentations from the functions are output as numpy arrays, so segmentation results can easily be displayed locally and are compatible with Neuroglancer.To start, we'll import adaptive_thresh to access the segmentation algorithms and convert the cloudvolume volume of the brain to a SimpleITK image format for segmentation.<jupyter_code>from brainlit.algorithms.generate_fragments import adaptive_thresh
img_T1, img_T1_255 = adaptive_thresh.get_img_T1(img)<jupyter_output><empty_output><jupyter_text>Gaussian Mixture Model ThresholdingIn this algorithm, a 2-component Gaussian mixture model is trained on the volume. Since neurons have higher voxel intensities in the Mouselight data than the background intensity, we threshold the entire volume by the *least* intensity in the component containing the *highest* intensity.<jupyter_code>labels = adaptive_thresh.gmm_seg(img)
print(labels.shape)
label_img = sitk.GetImageFromArray(labels, isVector=False)
plot_3d(sitk.LabelOverlay(img_T1_255, label_img), zslices=range(48,53), title="2-Component Gaussian Mixture Model")<jupyter_output>(66, 50, 52)<jupyter_text>Otsu's MethodOtsu's method thresholds the volume into two classes (segments) by maximizing inter-class variance.<jupyter_code>labels = adaptive_thresh.otsu(img)
label_img = sitk.GetImageFromArray(labels, isVector=False)
plot_3d(sitk.LabelOverlay(img_T1_255, label_img), zslices=range(48,53), title="Otsu")<jupyter_output><empty_output><jupyter_text>Region Growing MethodsRegion growing methods iteratively examine pixels neighboring seed points, and determine if these neighboring pixels should be added to a region.Three such methods are available: \In *closed-connected thresholding*, pixels in a region are connected to a seed and lie within a range of values. Here, the lower threshold is automatically calculated based on the Gaussian Mixture model thresholding above. \In *neighborhood-connected thrsholding*, pixels in a region are connected to a seed and lie within a neighborhood. Again, the lower threshold is automatically calculated based on the Gaussian Mixture model thresholding above.\In *confidence-connected thresholding*, pixels in a region have intensities consistent with a seed point; the pixel intensities lie within a confidence interval such that they are within some "multiplier" number of standard deviations from the mean of the neighborhood of the seed point. In this example, the multiplier is set to 1. If multiple iterations are used, the mean and standard deviation are calculated again at each iteration using pixels in the previous segmentation.In each of these methods, binary morphological closing is used to connect separate segments of the neuron / fill holes after the algorithm has run.We base  the algorithms on a seed set to the center of the downloaded volume, since the volume is centered on a point known to be within a neuron.<jupyter_code>_, seed = adaptive_thresh.get_seed(voxel)
labels = adaptive_thresh.connected_threshold(img, seed)
print(labels.shape)
label_img = sitk.GetImageFromArray(labels, isVector=False)
plot_3d(sitk.LabelOverlay(img_T1_255, label_img), zslices=range(48,53), title="Closed-Connected Threshold")
labels = adaptive_thresh.neighborhood_connected_threshold(img, seed)
label_img = sitk.GetImageFromArray(labels, isVector=False)
plot_3d(sitk.LabelOverlay(img_T1_255, label_img), zslices=range(48,53), title="Neighborhood-Connected Threshold")
labels = adaptive_thresh.confidence_connected_threshold(img, seed, num_iter=1, multiplier=1)
label_img = sitk.GetImageFromArray(labels, isVector=False)
plot_3d(sitk.LabelOverlay(img_T1_255, label_img), zslices=range(48,53), title="Confidence-Connected Threshold")<jupyter_output><empty_output><jupyter_text>Since this thresholding doesn't appear to segment the image, we'll try again with a smaller multiplier, 80% of the value previously used.<jupyter_code>labels = adaptive_thresh.confidence_connected_threshold(img, seed, num_iter=1, multiplier=0.8)
label_img = sitk.GetImageFromArray(labels, isVector=False)
plot_3d(sitk.LabelOverlay(img_T1_255, label_img), zslices=range(48,53), title="Confidence-Connected Threshold")<jupyter_output><empty_output><jupyter_text>Level-Set SegmentationThe level set method segments an image based on its contours. The main idea is to begin with a region based on an initial contour, and expand the region outward until edges are reached. When root mean squared change in the level set function for an iteration is below a threshold, or the maximum number of iteration have elapsed, the algorithm is said to have converged. Lower and upper thresholds may be set to the mean of the image minus or plus its standard deviation times some factor (1.5 here), if the thresholds are not set explicitly.<jupyter_code>labels = adaptive_thresh.level_set_seg(img, seed, factor = 1.5)
label_img = sitk.GetImageFromArray(labels, isVector=False)
plot_3d(sitk.LabelOverlay(img_T1_255, label_img), zslices=range(48,53), title="Level-Set Segmentation")<jupyter_output><empty_output><jupyter_text>Unfortunately, this method did effectively segment this neuron.  Fast-Marching SegmentationThe fast marching method is another form of level-set segmentation. Only a positive speed term is used to govern the differential equation in the method.<jupyter_code>labels = adaptive_thresh.fast_marching_seg(img, seed)
label_img = sitk.GetImageFromArray(labels, isVector=False)
plot_3d(sitk.LabelOverlay(img_T1_255, label_img), zslices=range(48,53), title="Fast-Marching Segmentation")<jupyter_output><empty_output><jupyter_text>Unfortunately, this method did effectively segment this neuron.  Napari to View SegmentationsWe can also view the segmentations in 3D using Napari.<jupyter_code>%gui qt
labels = adaptive_thresh.confidence_connected_threshold(img, seed, num_iter=1, multiplier=0.8)
napari_viewer(img, labels, label_name="Confidence-Connected Threshold")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09295189071470995, 'fraction_numerical': 0.020392001583844785, 'mean_word_length': 5.201964395334561, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 43, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '60232', 'license': 'Apache-2.0', 'n_tokens_mistral': 3424, 'n_tokens_neox': 3180, 'n_words': 960, 'path': 'docs/notebooks/algorithms/tutorial_notebook_adaptive_thresh.ipynb', 'repo_name': 'Jwright707/brainl', 'size': '10102'}