BPsCentroid_1.groovy

# @ImageJ ij
# @LUTService lut

import java.text.DecimalFormat
import net.imglib2.display.ColorTable
import sc.fiji.snt.*
import sc.fiji.snt.io.*
import sc.fiji.snt.analysis.*
import sc.fiji.snt.annotation.*
import sc.fiji.snt.viewer.*
import sc.fiji.snt.util.*
import org.scijava.util.*

/* Criteria 1. Only cells with somata in this area are considered */
motorAreas = AllenUtils.getCompartment("Somatomotor areas")
/* Criteria 2. Only cells w/ axonal BPs in this area are considered */
thalamus = AllenUtils.getCompartment("DORsm")
/* Criteria 3. Only cells w/ at least this no. of BPs are considered */
bpsCutoff = 10


class Centroid {
	SNTPoint soma, bps	// barycenters for soma and axonal BPs 
	String label		// identifier
	ColorRGB color		// annotation color
	double mValue		// mapping value
}

def getIDs() {
	// Read pre-filtered ids from file, otherwise retrieve all ids in Database
	ids = []
	file = new File("/groups/mousebrainmicro/home/ferreirat/code/analysis/TH_targeting_cells")
	if (file.exists()) {
		file.eachLine { ids << it }
	} else {
		df = new DecimalFormat("0000")
		for (id in 1..MouseLightLoader.getNeuronCount()) {
			id = "AA" + df.format(id)
			ids << id
		}
	}
	println("Will be parsing ${ids.size()} cells")
	ids
}

def assignUniqueColors(centroidList) {
	// Assigns (random) high constrast colors to a Centroid list
	colors = SNTColor.getDistinctColors(centroidList.size())
	centroidList.eachWithIndex { centroid, index ->
		centroid.color = colors[index]
	}
}

def assignMappingValues(centroidList) {
	// Assigns the mapping value for color coding
	meshCenterL = thalamus.getMesh().getBarycentre("left")
	meshCenterR = thalamus.getMesh().getBarycentre("right")
	for (centroid in centroidList) {
		meshCenter = (AllenUtils.isLeftHemisphere(centroid.bps)) ? meshCenterL : meshCenterR
		centroid.mValue = ((PointInImage)centroid.bps).distanceTo((PointInImage)meshCenter)
	}
}

def assignLUT(centroidList) {
	// Assigns a Centroid's mapping value to a LUT entry
	if (centroidList.isEmpty()) return
	colorTable = lut.loadLUT(lut.findLUTs().get("Ice.lut"))
	Comparator mc = { a, b -> a.mValue == b.mValue ? 0 : a.mValue < b.mValue ? -1 : 1 }
	min = centroidList.min(mc).mValue
	max = centroidList.max(mc).mValue
	int idx
	for (centroid in centroidList) {
		if (centroid.mValue <= min) idx = 0
		else if (centroid.mValue > max) idx = colorTable.getLength() - 1
		else idx = (int) Math.round((colorTable.getLength() - 1) * (centroid.mValue - min) / (max - min))
		color = new ColorRGB(colorTable.get(ColorTable.RED, idx), colorTable.get(ColorTable.GREEN, idx), 
				colorTable.get(ColorTable.BLUE, idx))
		centroid.color = color
	}
}

def getCentroids(loader) {

	if (!loader.idExists()) {
		println(" id not found. Skipping...")
		return null
	}

	// Retrieve the 1st node of the soma and its annotated compartment
	soma = loader.getNodes("soma")[0]
	compartment = (AllenCompartment) soma.getAnnotation()
	if (compartment && !motorAreas.contains(compartment)) {
		println(" Soma not associated with " + motorAreas + ". Skipping")
		return null
	}
	println(" Id matches soma location requirements!")

	// Retrieve the axonal arbor as a Tree object. Instantiate a TreeAnalyzer
	// so that we can conveniently access all of the axonal branch points
	axonalTree = loader.getTree("axon")
	analyzer = new TreeAnalyzer(axonalTree)
	branchPoints = analyzer.getBranchPoints()

	// Iterate through all the branch points (PointInImage objects)
	// and extract those in thalamus
	println(" Assessing location of ${branchPoints.size()} branch points...")
	filteredBranchPoints = []
	for (branchPoint in branchPoints) {
		compartment = (AllenCompartment) branchPoint.getAnnotation()
		if (compartment && thalamus.contains(compartment))
			filteredBranchPoints.add(branchPoint)
	}
	println(" Found ${filteredBranchPoints.size()} match(es)")

	// Ignore cells with less than the specified cutoff
	if (filteredBranchPoints.size() < bpsCutoff)
		return null

	centroid = new Centroid()
	centroid.soma = soma
	centroid.bps = SNTPoint.average(filteredBranchPoints)
	centroid.label = loader.id

	centroid
}

// Parse cells
centroidMap = [:].withDefault {key -> return []}
getIDs().stream().forEach {
	print("Parsing $it...")
	loader = new MouseLightLoader(it)
	centroid = getCentroids(loader)
	if (!centroid) {
		println(" Not enough branch points associated with $thalamus")
	} else {
		if (AllenUtils.isLeftHemisphere(centroid.soma))
			centroidMap."left" << centroid
		else
			centroidMap."right" << centroid
	}
}
println("Filtered cells: ${centroidMap."left".size()} + ${centroidMap."right".size()}")

// Perform the color mapping on each hemisphere
for (e in centroidMap) {
	println("Mapping ${e.key} hemisphere...")
	assignMappingValues(e.value)
	assignLUT(e.value)
}

// Merge centroid from both hemispheres
centroids = centroidMap."left" + centroidMap."right"
if (centroids.isEmpty()) return

// Render annotations
println("Rendering...")
viewer = new Viewer3D(ij.context())

// Organize annotations by category
annotMap = [:].withDefault {key -> return []}
centroids.each { centroid ->
	annot = viewer.annotatePoint(centroid.bps, centroid.label + " bps")
	annot.setSize(20)
	annot.setColor(centroid.color, 10)
	annotMap."bps" << annot
	annot = viewer.annotatePoint(centroid.soma, centroid.label + " soma")
	annot.setSize(30)
	annot.setColor(centroid.color, 10)
	annotMap."soma" << annot
	annot = viewer.annotateLine([centroid.soma, centroid.bps], centroid.label + " vector")
	annot.setSize(10)
	annot.setColor(centroid.color, 50)
	annotMap."vector" << annot
}
viewer.mergeAnnotations(annotMap."soma", "Somas")
viewer.mergeAnnotations(annotMap."bps", "BPs centroids")
viewer.mergeAnnotations(annotMap."vector", "vectors")

// Add meshes
["Whole Brain", "Somatomotor areas", "Thalamus", "VENT", "VAL", "VM", "VP"].each {
	viewer.add(AllenUtils.getCompartment(it).getMesh())
}

viewer.show()
println("Done.")