From 02d7cc538274bbafca48faa57562d95ec04cbb94 Mon Sep 17 00:00:00 2001
From: Philipp Schlegel <observing@web.de>
Date: Sat, 28 Sep 2024 10:35:40 +0100
Subject: [PATCH] docs: update dotprops tutorial

---
 docs/examples/0_io/tutorial_io_02_dotprops.py | 86 ++++++++++++++++---
 requirements.txt                              |  1 +
 2 files changed, 75 insertions(+), 12 deletions(-)

diff --git a/docs/examples/0_io/tutorial_io_02_dotprops.py b/docs/examples/0_io/tutorial_io_02_dotprops.py
index cab0166d..ba2ba102 100644
--- a/docs/examples/0_io/tutorial_io_02_dotprops.py
+++ b/docs/examples/0_io/tutorial_io_02_dotprops.py
@@ -6,30 +6,92 @@
 
 [`navis.Dotprops`][] are point clouds with associated principal vectors which are mostly used for
 NBLASTing. They are typically derivatives of skeletons or meshes but you can load them straight from
-confocal data using [`navis.read_nrrd`][]:
+e.g. confocal image stacks using [`navis.read_nrrd`][] or [`navis.read_tiff`][].
 """
 
 # %%
 import navis
+import matplotlib.pyplot as plt
 
 # %%
 # ## From image data
 #
-# For this example I downloaded one of Janelia's Fly Light confocal stacks ([link](https://splitgal4.janelia.org/))
-# and converted it to NRRD format using [ImageJ](https://imagej.net/ij/).
+# For this example we will use a stack from [Janelia's split Gal4 collection](https://splitgal4.janelia.org/).
+# This `LH2094` line is also available from [Virtual Fly Brain](https://v2.virtualflybrain.org/org.geppetto.frontend/geppetto?id=VFB_00102926&i=VFB_00101567,VFB_00102926)
+# where, conveniently, they can be downloaded in NRRD format which we can directly read into {{ navis }}.
 #
-# Load NRRD file into Dotprops instead of VoxelNeuron:
-# ```python
-# dp = navis.read_nrrd(
-#    "~/Downloads/JRC_SS86025_JRC_SS86025-20211112_49_B6.nrrd",
-#    output="dotprops",
-#    threshold=3000,
-# )
-# ```
+# Let's do this step-by-step first:
+
+# Load raw NRRD image
+im, header = navis.read_nrrd(
+    "https://v2.virtualflybrain.org/data/VFB/i/0010/2926/VFB_00101567/volume.nrrd",
+    output="raw"
+)
+
+# Plot a maximum projection
+max_proj = im.max(axis=2)
+plt.imshow(
+    max_proj.T,
+    extent=(0, int(0.5189 * 1210), (0.5189 * 566), 0),  # extent is calculated from the spacing (see `header`) times the no of x/y pixels
+    cmap='Greys_r',
+    vmax=10  # make it really bright so we can see neurons + outline of the brain
+    )
+
+# %%
+# At this point we could threshold the image, extract above-threshold voxels and convert them to a Dotprops object.
+# However, the easier option is to use [`navis.read_nrrd`][] with the `output="dotprops"` parameter:
+
+dp = navis.read_nrrd(
+    "https://v2.virtualflybrain.org/data/VFB/i/0010/2926/VFB_00101567/volume.nrrd",
+    output="dotprops",
+    threshold=5,  # threshold to determine which voxels are used for the dotprops
+    thin=True,   # see note below on this parameter!
+    k=10  # number of neighbours to consider when calculating the tangent vector
+)
+
+# %%
+# !!! note "Thinning"
+#     In the above [`read_nrrd`][navis.read_nrrd] call we used `thin=True`. This is a post-processing step that
+#     thins the image to a single pixel width. This will produce "cleaner" dotprops but can also remove denser
+#     neurites thus emphasizing the backbone of the neuron. This option requires the `scikit-image` package:
+#
+#     ```bash
+#     pip install scikit-image
+#     ```
+#
+#  Let's overlay the dotprops on the maximum projection:
+
+fig, ax = plt.subplots()
+ax.imshow(
+    max_proj.T,
+    extent=(0, int(0.5189 * 1210), (0.5189 * 566), 0),
+    cmap='Greys_r',
+    vmax=10
+    )
+navis.plot2d(dp, ax=ax, view=("x", "-y"), method="2d", color="r", linewidth=1.5)
+
+# %%
+# This looks pretty good but we have a bit of little fluff around the brain which we may want to get rid off:
+
+# Drop everything but the two largest connected components
+dp = navis.drop_fluff(dp, n_largest=2)
+
+# Plot again
+fig, ax = plt.subplots()
+ax.imshow(
+    max_proj.T,
+    extent=(0, int(0.5189 * 1210), (0.5189 * 566), 0),
+    cmap='Greys_r',
+    vmax=10
+    )
+navis.plot2d(dp, ax=ax, view=("x", "-y"), method="2d", color="r", linewidth=1.5)
 
 # %%
 # !!! note
-#     Note the threshold parameter? It determines which voxels (by brightness) are used and which are ignored!
+#     To extract the connected components, [`navis.drop_fluff`][] treats all pairs of points within a certain distance
+#     as connected. The distance is determined by the `dp_dist` parameter which defaults to 5 x the average distance
+#     between points. This is a good value ehre but you may need adjust it for your data.
+#
 #
 # ## From other neurons
 #
diff --git a/requirements.txt b/requirements.txt
index 83c0779f..495360f8 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -64,6 +64,7 @@ neuprint-python
 caveclient
 cloud-volume
 flybrains
+scikit-image
 Shapely>=1.6.0
 
 #extra: dev