SNT: how to make a 2D dendrite density map?

Hi All,

I am trying to make some simple analysis but struggle with the implementation.

I have 3D reconstructions (mouse L5 pyramidal cells) in SWC format.
I would like to project them down to 2D and then measure the sum dendrite length in a grid (say at 10x10 µm resolution).
The idea is to get a “2D dendrite density map" for different types of dendrites.

Any thoughts?

Thanks,
Ede

Welcome to the forum, @ede.rancz !

Here is one way of approaching such an analysis using scripting (Jython). We can get a 2D projection of the rasterized skeleton of the reconstruction as an 8-bit binary image (using a nifty one-liner SNT method). Then generate a grid over this image and count the pixels of the projected skeleton that occur in each bin, calculate their probabilities and display the result in a new image.

If you run the below script in the Script Editor, leaving the “Reconstruction file” field empty, it will generate the following images, the first being the projected reconstruction skeleton and the second the density map.

#@ File (label="Reconstruction file (Leave empty for demo)", style="file", required=false) recFile
#@ Float (label="Area per point (spatially calibrated units)^2") boxArea
#@ boolean (label="Center grid on image") centerGrid
#@ SNTService snt
#@ RoiManager roiManager
#@ ResultsTable resultsTable

# Documentation Resources: https://imagej.net/SNT:_Scripting
# Latest SNT API: https://morphonets.github.io/SNT/

import math
from ij import IJ
from ij.gui import NewImage
from ij.plugin import LutLoader
from sc.fiji.snt import Tree


boxW = boxH = math.sqrt(boxArea)

# Retrieve a 2D projection of the rasterized skeleton of the reconstruction
# as an 8-bit binary image (skeleton: 255, background: 0)
if recFile is not None:
	imp = Tree(recFile.getAbsolutePath()).getSkeleton2D()
else :
	# Use dendrites of a mouse SSp-m5 neuron
	# AA0001 in the MouseLight database 
	# https://ml-neuronbrowser.janelia.org/
	imp = snt.demoTrees().get(0).getSkeleton2D()

imp.show()

# Get total sum of pixel intensities over the image
totalIntDen = IJ.getValue(imp, "RawIntDen")

# Create grid overlay using point ROIs, with each grid box having area == boxArea
if centerGrid:
	IJ.run(imp, "Grid...", "grid=Points area=" + str(boxArea) + " center")
else:
	IJ.run(imp, "Grid...", "grid=Points area=" + str(boxArea))

# Get the composite grid ROI from the image overlay and
# add it to the RoiManager
overlay = imp.getOverlay()
roiManager.reset()  # clear any existing rois
roiManager.addRoi(overlay.get(0)) # add the single composite grid ROI to RoiManager
roiManager.select(0)  # select this ROI
roiManager.runCommand("Split")  # split into individual point ROIs
roiManager.runCommand("Delete") # delete selected composite ROI, leaving the point ROIs
imp.setOverlay(None)  # clear overlay from source image

# Collect xy coordinates of the grid box corner points
xs = []
ys = []
for roi in roiManager:
	centroid = roi.getContourCentroid()
	xs.append(int(round(centroid[0])))
	ys.append(int(round(centroid[1])))
xMax = max(xs)
yMax = max(ys)

# Now make the rectangular ROIs using the corner points
roiManager.reset()  # clear existing point ROIs
for x, y in zip(xs, ys):
	if x >= xMax or y >= yMax:
		# Skip boundary points
		continue
	IJ.makeRectangle(x, y, boxW, boxH)
	roiManager.runCommand("Add")
	
# Measure integrated density for each ROI
IJ.run("Set Measurements...", "integrated")
resultsTable.reset()  # clear any existing data
roiManager.runCommand("Select All")
roiManager.runCommand("Measure")
colIdx = resultsTable.getColumnIndex("RawIntDen")
# Base the max pixel intensity (255) on the maximum encountered bin probability
maxProb = max(resultsTable.getColumn(colIdx)) / totalIntDen

# Create image with same dimensions as source image
newImp = NewImage.createImage("Density Map", imp.getWidth(), imp.getHeight(), 0, 8, NewImage.FILL_BLACK)
# Use viridis color palette to paint the density map
lutPath = IJ.getDirectory("luts") + "mpl-viridis.lut"
lut = LutLoader.openLut(lutPath)
ip = newImp.getProcessor()
ip.setLut(lut)

# Fill in ROIs on the new image using bin probabilities
row = 0
for roi in roiManager:
	localProb = resultsTable.getValue("RawIntDen", row) / totalIntDen
	# Map the bin probability onto the interval [0, 255]
	pixelValue = (255 / maxProb) * localProb
	ip.setColor(pixelValue)
	ip.fill(roi)
	row += 1
	
newImp.setOverlay(overlay)  # show grid on new image
newImp.show()

2 Likes

@arshadic, super awesome SNT-fu!
@ede.rancz, a couple of notes on Cameron’s script:
I’m assuming you will need to extend this to a collection of cells (perhaps followed by a Gaussian blur)? In that case you may need to register your cells so that their arbors align. The Script Editor’s Templates>Neuroanatomy>Render Cell Collection (Multipanel Figure) provides a basic example on how to do that.

Also, for sake of completeness: One could also approach the problem without rasterizing the dendritic arbor, and work directly on the reconstruction: Each reconstruction (a Tree in SNT’s API) is associated with a BoundingBox class: One could programmatically tessellate such bounding box, and retrieve the length on each tile. However, the script would likely be more verbose.

2 Likes

Big thank you @arshadic, works like a charm. And thanks for the tips @tferr for averaging across cells, needed indeed.

1 Like