BPsCentroid_2.groovy

# @ImageJ ij
# @LUTService lut

import java.text.DecimalFormat
import net.imglib2.display.ColorTable
import org.scijava.util.Colors;
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("Thalamus")
/* Criteria 3. Only cells w/ at least this no. of BPs are considered */
bpsCutoff = 10

class Nuclei {
	SNTPoint centerL	//mesh barycente L hemisphere
	SNTPoint centerR	//mesh barycente R hemisphere
	ColorRGB color		// annotation color
	AllenCompartment compartment
	OBJMesh mesh
}

class Centroid {
	SNTPoint soma, bps	// barycenters for soma and axonal BPs 
	String label		// identifier
	ColorRGB color		// annotation color
	double mValue		// mapping value
	AllenCompartment bpsCompartment // the closest compartment to Axonal BPs
}

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 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
}


children = thalamus.getChildren(3)
colors = SNTColor.getDistinctColors(children.size())
inputMap = [:]
children.eachWithIndex { child, index ->
	if (child.getMesh())
		inputMap << ["${child.acronym()}":colors[index]]
}

nucleiMap = [:]
inputMap.each { label, color ->
	compartment =  AllenUtils.getCompartment(label)
	mesh = compartment.getMesh()
	if (mesh) {
		mesh.setColor(color, 90)
		centerL = mesh.getBarycentre("left")
		centerR = mesh.getBarycentre("right")
		nuclei = new Nuclei()
		nuclei.centerL = centerL
		nuclei.centerR = centerR
		nuclei.color = color
		nuclei.compartment = compartment
		nuclei.mesh = mesh
		nucleiMap << ["${label}":nuclei]
		println "${label}:${compartment.name()}   color:$color"
	}
}

// Parse cells
centroidMap = [:].withDefault {key -> return []}
getIDs().each() { id ->
	print("Parsing $id...")
	loader = new MouseLightLoader(id)
	bpsCentroid = getCentroids(loader)
	if (!bpsCentroid) {
			println(" Not enough branch points associated with $thalamus")
	} else {

		isLeft = AllenUtils.isLeftHemisphere(bpsCentroid.soma);
		Nuclei closestNuclei
		dx = Double.MAX_VALUE
		for (n in nucleiMap) {
			nuclei = n.getValue()			
			d = bpsCentroid.bps.distanceTo((isLeft) ? nuclei.centerL : nuclei.centerR)
			if (d < dx) {
				closestNuclei = nuclei
				dx = d
			}
		}
		bpsCentroid.color = closestNuclei.color
		bpsCentroid.bpsCompartment = closestNuclei.compartment
		println("{$id}: Closest nuclei: ${closestNuclei.compartment.name()}, color: ${closestNuclei.color}")
		if (AllenUtils.isLeftHemisphere(bpsCentroid.soma))
			centroidMap."left" << bpsCentroid
		else
			centroidMap."right" << bpsCentroid
	}
	
}
println("Filtered cells: ${centroidMap."left".size()} + ${centroidMap."right".size()}")



// 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 ->
	id = centroid.bpsCompartment.id()
	annot = viewer.annotatePoint(centroid.bps, centroid.label + " bps")
	annot.setSize(20)
	annot.setColor(centroid.color, 10)
	annotMap."${id}" << annot
	annot = viewer.annotatePoint(centroid.soma, centroid.label + " soma")
	annot.setSize(30)
	annot.setColor(centroid.color, 10)
	annotMap."${id}" << annot
	annot = viewer.annotateLine([centroid.soma, centroid.bps], centroid.label + " vector")
	annot.setSize(10)
	annot.setColor(centroid.color, 50)
	annotMap."${id}vec" << annot
}

// Add meshes and annotations to viewer
viewer.add(AllenUtils.getCompartment("Whole Brain").getMesh())
for (centroid in centroids.unique {c1, c2 -> c1.bpsCompartment.id() <=> c2.bpsCompartment.id()}) {

	// Add annotations
	id = centroid.bpsCompartment.id()
	acronym = centroid.bpsCompartment.acronym()
	viewer.mergeAnnotations(annotMap."${id}", "Centroids - ${acronym}")
	viewer.mergeAnnotations(annotMap."${id}vec", "Vectors - ${acronym}")
	// Add meshes
	mesh = centroid.bpsCompartment.getMesh()
	mesh.setColor(centroid.color, 90)
	viewer.add(mesh)
}
viewer.show()
println("Done.")