// Olivier Seynnes, June 19, 2019

// _Thickness_BB v3.0
// This macro finds the thickness of the brachialis and biceps brachii muscles in a cross-sectional image.

// OPTION PARAMETERS ------------------------------------------------------------
/*									
// String(value = "--------------- Arm Flexor Muscles Thickness ---------------", visibility="MESSAGE") title
// Boolean (label="Scaled measurements", value=false, persist=false) scaling
// String (label = "Scaling", choices = {"Philips HD11", "no scaling"}) app	
// Integer(label="Scan depth (cm)", min=3, max=7, style="scroll bar", value = 5) depth 
*/
//--------------------------------------------------------------------------------
macro "Arm Flexor Muscles Thickness Action Tool - C0a0L18f8L818f" {
		/*
	if (app=="Philips HD11") {		
		if (depth=3) {
			SF = 13.3333;
		}
		if (depth=4) {
			SF = 11.3;
		}
		if (depth=5) {
			SF = 9.04;
		}
		if (depth=6) {
			SF = 7.5333;
		}
		if (depth=7) {
			SF = 6.4571;
		}
	
	} else if (app=="Telemed LS128") {
		if (depth=="3cm") {
			SF = ?;
		}
		if (depth=="4cm") {
			SF = 15.7;
		}
		if (depth=="5cm") {
			SF = 12.62;
		}
		if (depth=="6cm") {
			SF = 10.57;
		}
		if (depth=="7cm") {
			SF = ?;
		}
	*/	
	SF = 1; // no calibration

	run("Set Measurements...", "centroid redirect=None decimal=1");
	IDraw = getImageID();
	width = getWidth();
	height = getHeight();
	
	//Close results window (obtained from previous analysis, should be kept out of loop for stacks)
	if (isOpen("Results")) {
		selectWindow("Results");
	    run("Close");
	}
	if (isOpen("ROI Manager")) {
		selectWindow("ROI Manager");
		roiManager("Show None");
	 	run("Close");
	}
	run("Select None");
	
	// prompt user to draw lines over features of interest
	selectImage(IDraw);
	setTool("polyline"); //in case of a single image
	waitForUser("Draw a line over the bone edge. Click OK when done");
	roiManager("Add");
	run("Select None");
	waitForUser("Draw a line over the biceps brachii upper aponeurosis. Click OK when done");
	roiManager("Add");
	run("Select None");
	waitForUser("Draw a line over the brachialis upper aponeurosis. Click OK when done");
	roiManager("Add");
	roiManager("Deselect");
	if(roiManager("count") < 2)
		exit("The macro needs at least 2 entries in the ROI Manager");
	
	n = roiManager("count");
	roiManager("select", 0); // Bone ROI
	run("Fit Spline");
	run("Interpolate");
	roiManager("update");
	run("Fit Circle"); // fit circle to find center
	List.setMeasurements; //to pull measurements of centroid coordinates from the "Measure" function
	x1 = List.getValue("X");
	y1 = List.getValue("Y");
	
	roiManager("select", 1); // BB upper aponeurosis ROI
	run("Fit Spline");
	run("Interpolate");
	roiManager("update");
	run("Fit Circle"); // fit circle to find center
	List.setMeasurements; //to pull measurements of centroid coordinates from the "Measure" function
	u1 = List.getValue("X");
	v1 = List.getValue("Y");

	// Define measurement line (y = mx + b) between the centre of each fitted circles
	m = (y1-v1)/(x1-u1); //slope 
	b = v1 - m * u1; 
	u2 = -b/m; //coordinates of intercept are (-b/m,0)
	makeLine(x1, y1, u2, 0);
	roiManager("Add");
	
	alpha = 180-(atan(m)*(180/PI)); // get angle of main measurement line
	beta1 = alpha + 10; // set secondary measurement line 10 deg away from the main one
	n1 = tan((180-beta1)*(PI/180)); // slope
	c1 = y1 - n1 * x1;
	u3 = -c1/n1;
	makeLine(x1, y1, u3, 0);
	roiManager("Add");
	beta2 = alpha - 10;
	n2 = tan((180-beta2)*(PI/180));
	c2 = y1 - n2 * x1;
	u4 = -c2/n2;
	makeLine(x1, y1, u4, 0);
	roiManager("Add");
	//roiManager("Show All");
 
	// intersection between measurement lines and other ROIs. Method from https://forum.image.sc/t/intersection-of-a-spline-fitted-curve-and-a-straight-line/3268
	strokeWd = 1;
	n = roiManager("count");
	for (i=0; i<n; i++){ 
		lineToArea(i);	//convert line ROIs to area
	}	

	intersXb1 = newArray(2); 
	intersXb2 = newArray(2);
	intersXb3 = newArray(2);
	intersXb1 = intCoordinates (0, 3, intersXb1); 
	xbmean = intersXb1[0]; // ROIs intersection coordinates (bone edge and line of measurement 1)
	ybmean = intersXb1[1];
	intersXb2 = intCoordinates (0, 4, intersXb2); 
	xb2mean = intersXb2[0]; // ROIs intersection coordinates (bone edge and line of measurement 2)
	yb2mean = intersXb2[1];
	intersXb3 = intCoordinates (0, 5, intersXb3); 
	xb3mean = intersXb3[0]; // ROIs intersection coordinates (bone edge and line of measurement 3)
	yb3mean = intersXb3[1];

	intersXbb1 = newArray(2); 
	intersXbb2 = newArray(2);
	intersXbb3 = newArray(2);
	intersXbb1 = intCoordinates (1, 3, intersXbb1); 
	xbbmean = intersXbb1[0]; // ROIs intersection coordinates (BB and line of measurement 1)
	ybbmean = intersXbb1[1];
	intersXbb2 = intCoordinates (1, 4, intersXbb2); 
	xbb2mean = intersXbb2[0]; // ROIs intersection coordinates (BB and line of measurement 2)
	ybb2mean = intersXbb2[1];
	intersXbb3 = intCoordinates (1, 5, intersXbb3); 
	xbb3mean = intersXbb3[0]; // ROIs intersection coordinates (BB edge and line of measurement 3)
	ybb3mean = intersXbb3[1];

	intersXbr1 = newArray(2); 
	intersXbr2 = newArray(2);
	intersXbr3 = newArray(2);
	intersXbr1 = intCoordinates (2, 3, intersXbr1); 
	xbrmean = intersXbr1[0]; // ROIs intersection coordinates (BB and line of measurement 1)
	ybrmean = intersXbr1[1];
	intersXbr2 = intCoordinates (2, 4, intersXbr2); 
	xbr2mean = intersXbr2[0]; // ROIs intersection coordinates (BB and line of measurement 2)
	ybr2mean = intersXbr2[1];
	intersXbr3 = intCoordinates (2, 5, intersXbr3); 
	xbr3mean = intersXbr3[0]; // ROIs intersection coordinates (BB edge and line of measurement 3)
	ybr3mean = intersXbr3[1];
		
	// Draw and measure lines
	makeLine(xbmean, ybmean, xbbmean, ybbmean);
	roiManager("Add");
	makeLine(xb2mean, yb2mean, xbb2mean, ybb2mean);
	roiManager("Add");
	makeLine(xb3mean, yb3mean, xbb3mean, ybb3mean);
	roiManager("Add");

	ThTot1 = calcLength (xbmean, ybmean, xbbmean, ybbmean);
	ThTot2 = calcLength (xb2mean, yb2mean, xbb2mean, ybb2mean);
	ThTot3 = calcLength (xb3mean, yb3mean, xbb3mean, ybb3mean);
	ThBr1 = calcLength (xbmean, ybmean, xbrmean, ybrmean);
	ThBr2 = calcLength (xb2mean, yb2mean, xbr2mean, ybr2mean);
	ThBr3 = calcLength (xb3mean, yb3mean, xbr3mean, ybr3mean);

	ThTot = (ThTot1+ThTot2+ThTot3)/3;
	ThBr = (ThBr1+ThBr2+ThBr3)/3;	
	ThBB = ThTot - ThBr;		
	
	print( "Total thickness: " + ThTot);
	print( "Brachialis thickness: " + ThBr);
	print( "Biceps brachii thickness: " + ThBB);
	
	roiManager("select", newArray(3,4,5));
	roiManager("Delete");
	roiManager("Show All");
		
} // end of macro

// FUNCTIONS ---------------------------------------------------------------

function lineToArea(ROI) {
	roiManager("select", ROI);
	run("Line to Area");
	roiManager("update");
}

function intCoordinates(ROI1, ROI2, intersX) {
	roiManager("select", newArray(ROI1,ROI2)); // select ROIs that intersect
	roiManager("and");
	getSelectionCoordinates(xI, yI); // coordinates of polygon at intersection 
	Array.getStatistics(xI, xImin, xImax, xImean, xIstdDev);
	Array.getStatistics(yI, yImin, yImax, yImean, yIstdDev);
	intersX = newArray(xImean,yImean);
	return intersX;
}

function calcLength(x1, y1, x2, y2) {
	Length = round((sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)))/SF);
	return Length;
}