Updated docs and added one more example.

doc/source/examples/ex_vessel_graph.txt

@@ -0,0 +1,232 @@
+.. _vessel_graph_example:
+
+Vessel graph
+============
+
+This example demonstrates generation of a vessel graph structure from a tomographic image of vasculature in mouse brain. The analysis process is similar to that used in [1], where it was used to trace vessels in tomographic images of corrosion cast models of cerebral vasculature.
+
+The analysis process consists of thresholding the original image in order to convert it to binary form, where the vessels are white and everything else is black. The binary image is skeletonized into single pixel wide lines, and the lines are traced into a graph structure. Spurious branches are removed using a custom condition based on the ratio of length and diameter of the branch. The average diameter of the branches is calculated and all the information is saved into a .vtk file for further analysis. The figures below represent the same slice through the tomographic image in various processing phases. Also, 3D visualization of the original image and the traced vessels are provided below.
+
+For running this example code, a demonstration image package available at the `GitHub releases page <https://github.com/arttumiettinen/pi2/releases>`__ is required.
+
+::
+
+	def vessel_tracing():
+
+		# Make binary image where vessels are foreground
+		img = pi.read(input_file('orig_basal_ganglia_217'))
+		pi.threshold(img, 18000) # The threshold value should be chosen using e.g. Otsu method
+		pi.closingfilter(img, 1, True, NeighbourhoodType.BOX) # This closes small holes in the vessels
+		pi.convert(img, ImageDataType.UINT8)
+		pi.multiply(img, 255)
+
+		# In order to fill cavities in the vessels, fill the background by a temporary color.
+		# Here we perform the background fill by a flood fill starting from non-vessel points
+		# in the image edge. Here we try to fill from all corners of the image.
+		# Note that it is possible that all corners are occupied by vessels and in this case
+		# this filling method does not work, but in practice that situation is very rare.
+		fill(img, 0, 0, 0)
+		fill(img, img.get_width() - 1, 0, 0)
+		fill(img, 0, img.get_height() - 1, 0)
+		fill(img, img.get_width() - 1, img.get_height() - 1, 0)
+		fill(img, 0, 0, img.get_depth() - 1)
+		fill(img, img.get_width() - 1, 0, img.get_depth() - 1)
+		fill(img, 0, img.get_height() - 1, img.get_depth() - 1)
+		fill(img, img.get_width() - 1, img.get_height() - 1, img.get_depth() - 1)
+		# Here the vessels, cavities and background have colors 255, 0, and 128, respectively.
+		# Set cavities and vessels to 255 and background to 0.
+		pi.replace(img, 255, 0)
+		pi.replace(img, 0, 255)
+		pi.replace(img, 128, 0)
+
+		pi.writeraw(img, output_file('vessel_bin'))
+
+		# Now calculate skeleton. We use the surfaceskeleton function and specify False as second
+		# argument to indicate that we want a line skeleton. In many cases this method makes
+		# more stable skeletons than the lineskeleton function.
+		pi.surfaceskeleton(img, False)
+		pi.writeraw(img, output_file('vessel_skele'))
+
+		# Below we will need the distance map of vessel phase, so we calculate it here.
+		img = pi.read(output_file('vessel_bin'))
+		dmap = pi.newimage(ImageDataType.FLOAT32)
+		pi.dmap(img, dmap)
+		pi.writeraw(dmap, output_file('vessel_dmap'))
+		dmap_data = dmap.get_data()
+
+
+		# Trace skeleton
+		skeleton = pi.read(output_file('vessel_skele'))
+		smoothing_sigma = 2
+		max_displacement = 2
+		vertices = pi.newimage(ImageDataType.FLOAT32)
+		edges = pi.newimage(ImageDataType.UINT64)
+		measurements = pi.newimage(ImageDataType.FLOAT32)
+		points = pi.newimage(ImageDataType.INT32)
+		pi.tracelineskeleton(skeleton, vertices, edges, measurements, points, True, 1, smoothing_sigma, max_displacement)
+
+
+		# Next, we will remove all edges that has at least one free and and whose
+		# L/r < 2.
+		# First, get edges, vertices, and branch length as NumPy arrays.
+		old_edges = edges.get_data()
+		vert_coords = vertices.get_data()
+
+		# The tracelineskeleton measures branch length by anchored convolution and returns it in the
+		# measurements image.
+		meas_data = measurements.get_data()
+		length_data = meas_data[:, 1]
+
+
+		# Calculate degree of each vertex
+		deg = {}
+		for i in range(0, vert_coords.shape[0]):
+			deg[i] = 0
+
+		for i in range(0, old_edges.shape[0]):
+			deg[old_edges[i, 0]] += 1
+			deg[old_edges[i, 1]] += 1
+
+		# Determine which edges should be removed
+		remove_flags = []
+		for i in range(0, old_edges.shape[0]):
+			n1 = old_edges[i, 0]
+			n2 = old_edges[i, 1]
+
+			# Remove edge if it has at least one free end point, and if L/r < 2, where
+			# r = max(r_1, r_2) and r_1 and r_2 are radii at the end points or the edge.
+			should_remove = False
+			if deg[n1] == 1 or deg[n2] == 1:
+
+				p1 = vert_coords[n1, :]
+				p2 = vert_coords[n2, :]
+
+				r1 = dmap_data[int(p1[1]), int(p1[0]), int(p1[2])]
+				r2 = dmap_data[int(p2[1]), int(p2[0]), int(p2[2])]
+
+				r = max(r1, r2)
+				L = length_data[i]
+
+				if L < 2 * r:
+					should_remove = True
+
+			# Remove very short isolated branches, too.
+			if deg[n1] == 1 and deg[n2] == 1:
+				L = length_data[i]
+				if L < 5 / 0.75: # (5 um) / (0.75 um/pixel)
+					should_remove = True
+
+			remove_flags.append(should_remove)
+
+		remove_flags = np.array(remove_flags).astype(np.uint8)
+		print(f"Before dynamic pruning: {old_edges.shape[0]} edges")
+		print(f"Removing {np.sum(remove_flags)} edges")
+
+		# This call adjusts the vertices, edges, and measurements images such that
+		# the edges for which remove_flags entry is True are removed from the graph.
+		pi.removeedges(vertices, edges, measurements, points, remove_flags, True, True)
+
+
+
+		# Convert to vtk format in order to get radius for each point and line
+		vtkpoints = pi.newimage()
+		vtklines = pi.newimage()
+		pi.getpointsandlines(vertices, edges, measurements, points, vtkpoints, vtklines);
+
+		# Get radius for each point
+		points_data = vtkpoints.get_data()
+		radius_points = np.zeros([points_data.shape[0]])
+		for i in range(0, points_data.shape[0]):
+			p = points_data[i, :]
+			r = dmap_data[int(p[1]), int(p[0]), int(p[2])]
+			radius_points[i] = r
+
+
+		# Get average radius for each branch
+		# Notice that the vtklines image has a special format that is detailed in
+		# the documentation of getpointsandlines function.
+		lines_data = vtklines.get_data()
+		radius_lines = []
+		i = 0
+		edge_count = lines_data[i]
+		i += 1
+		for k in range(0, edge_count):
+			count = lines_data[i]
+			i += 1
+
+			R = 0
+			for n in range(0, count):
+				index = lines_data[i]
+				i += 1
+				p = points_data[index, :]
+				R += dmap_data[int(p[1]), int(p[0]), int(p[2])]
+			R /= count
+
+			radius_lines.append(R)
+
+		radius_lines = np.array(radius_lines)
+
+
+		# Convert to vtk format again, now with smoothing the point coordinates to get non-jagged branches.
+		vtkpoints = pi.newimage()
+		vtklines = pi.newimage()
+		pi.getpointsandlines(vertices, edges, measurements, points, vtkpoints, vtklines, smoothing_sigma, max_displacement);
+
+
+		# Write to file
+		pi.writevtk(vtkpoints, vtklines, output_file('vessels'), "radius", radius_points, "radius", radius_lines)
+
+
+
+.. figure:: figures/vessels_orig.png
+   :scale: 100 %
+   :alt: Input image.
+
+   Input image. The bright regions are corrosion-cast blood vessels.
+
+
+.. figure:: figures/vessels_bin.png
+   :scale: 100 %
+   :alt: Segmented vessels.
+
+   Segmented vessels.
+
+
+.. figure:: figures/vessels_dmap.png
+   :scale: 100 %
+   :alt: Distance map of the segmented vessels.
+
+   Distance map of the segmented vessels.
+
+
+.. figure:: figures/vessels_skele.png
+   :scale: 100 %
+   :alt: Skeleton of the segmented vessels.
+
+   Skeleton of the segmented vessels.
+
+
+.. figure:: figures/vessels_orig_3d.png
+   :scale: 50 %
+   :alt: Volume visualization.
+
+   Volume visualization of the original vessel image. This visualization has been generated with the `Volume Viewer <https://imagej.net/plugins/volume-viewer.html>`__ plugin of `Fiji <https://fiji.sc>`__.
+
+
+.. figure:: figures/vessels_traced_3d.png
+   :scale: 50 %
+   :alt: Traced vessels.
+
+   Visualization of the traced vessels. Each vessel branch has been colored based on the average diameter of that branch. Bright colors correspond to large diameters. This visualization has been generated with `ParaView <https://www.paraview.org/>`__.
+
+
+
+
+
+
+References
+----------
+
+[1] Wälchli T., Bisschop J., Miettinen A. et al. Hierarchical imaging and computational analysis of three-dimensional vascular network architecture in the entire postnatal and adult mouse brain. Preprint available at https://doi.org/10.1101/2020.10.19.344903
+
+

doc/source/index.txt

@@ -7,7 +7,7 @@ Documentation of pi2
 ====================
 
 This is documentation for pi2/itl2, an image analysis program/library designed for processing and analysis of
-terapixel-scale volume images :ref:`locally or on a computer cluster <cluster_config>`.
+terapixel-scale volume images :ref:`locally or on a computer cluster <cluster_config>`. Some typical uses of this program are, e.g., :ref:`tracing blood vessels in tomographic images <vessel_graph_example>`, :ref:`analysis of fibre orientation in composite materials <orientation_example>`, or :ref:`stitching of volumetric mosaic images <nr_stitcher>`.
 
 The pi2 environment can be used from e.g.
 :ref:`Python <usage_python>`,
@@ -21,7 +21,7 @@ Some of the capabilities of pi2 include
 * :ref:`Fast non-rigid stitching <nr_stitcher>`
 * :ref:`Filtering <filtering_example>` for noise reduction: :ref:`bilateral <bilateralfilter>`, :ref:`approximate bilateral <bilateralfilterapprox>`, :ref:`variance weighted mean <vawefilter>`, :ref:`median <medianfilter>`, :ref:`Gaussian <gaussfilter>`, :ref:`high-pass <highpassfilter>`, ...
 * Post-processing of segmented images: :ref:`small region removal <regionremoval>` (volume opening), :ref:`minimum <minfilter>` and :ref:`maximum <maxfilter>` filtering (optionally with fast approximate spherical structuring element), :ref:`opening <openingfilter>` and :ref:`closing <closingfilter>`, ...
-* :ref:`Skeletonization <skeleton_example>` (:ref:`line <lineskeleton>` and :ref:`plate+line <surfaceskeleton>` skeletons), skeleton :ref:`tracing <tracelineskeleton>` with measurements (branch length, cross-sectional area...), skeleton :ref:`visualization <fillskeleton>`, :ref:`pruning <pruneskeleton>`,...
+* :ref:`Skeletonization <skeleton_example>` (:ref:`line <lineskeleton>` and :ref:`plate+line <surfaceskeleton>` skeletons), skeleton :ref:`tracing <vessel_graph_example>` with measurements (branch length, cross-sectional area...).
 * :ref:`Detection of interfaces <carpet_example>` using an Edwards-Wilkinson surface model
 * :ref:`Orientation analysis <orientation_example>` using the structure tensor approach
 * Fast separable :ref:`distance map <dmap>` and :ref:`local thickness map <local_thickness_example>`

doc/source/license

@@ -676,6 +676,34 @@ Public License instead of this License.  But first, please read
 <https://www.gnu.org/licenses/why-not-lgpl.html>.
 
 
+-----------------------------------------------------------------------------------------
+
+This program contains code from the C++ Mathematical Expression Toolkit Library, licensed
+under the MIT license available at http://www.opensource.org/licenses/MIT.
+See also http ://www.partow.net/programming/exprtk/index.html
+
+Copyright 2021 Arash Partow
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of
+this software and associated documentation files(the "Software"), to deal in
+the Software without restriction, including without limitation the rights to
+use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
+the Software, and to permit persons to whom the Software is furnished to do so,
+subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
+FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
+COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
+IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
+CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+
+
+
 -----------------------------------------------------------------------------------------
 
 This program contains code from PStreams library, licensed under the following license: