/*
Cross Profiles macro
N. Vischer, 14.12.18, 14:21

Creates "Two Peaks" plot of all 2-peak cross profiles.
Creates "One Peak" plot of all other cross profiles.
Creates a "OneP+TwoP" plot which is a weighted sum of the above.
Center/border ratio in 1-peak cells appear in columns "CB1-Ratio" and "CB-Ratio"
Center/border ratio in 2-peak cells appear in columns "CB2-Ratio" and "CB-Ratio"
Cross profiles are sampled perpendicularly at two positions of the cell:
at 20..40 % and at 60..80% of the cell axis.
Profiles are averaged with their mirror images, i.e. symmetrical
Collective results go to the Log window
While running, keep Caps Lock key down to observe perpendicular line rois
Works in combination with Coli-Inspector project files containing marked objects 
*/

macro "Cross-Profiles"{
	print("\\Clear");
	close("* Peak*");
	fluorChannel = 2;
	lWidthUM = 2; //cross profile evaluation length microns
	cellDiaUM = 0.95;//expected cell diameter in microns, used for 1 peak plots
	relTolerance = 1/25;//of the profiles max value for find Maxima
	ojHideResults();
	ojInitColumn("CB1-Ratio, CB2-Ratio, CB-Ratio");//center-border
	profile2P = newArray(100);//two peaks
	profile1P = newArray(100);//one peak
	profileT = newArray(100);//total
	nTwo = 0;
	nOne = 0;
	setBatchMode(true);
	setBatchMode("hide");
	for(img = 1; img <= ojNImages(); img++){
		pixSize = ojGetVoxelSize(img, "x");
		lWidth = lWidthUM/pixSize;
		cellDia = cellDiaUM/pixSize;
		

		ojShowImage(img);
		for(obj = ojFirstObject(img); obj <= ojLastObject(img); obj++){
			showProgress(obj/ojNObjects());
			ojShowObject(obj);
			ojItemToRoi();
			Stack.setChannel(fluorChannel);
			run("Interpolate", "interval=1 adjust");
			getSelectionCoordinates(xx, yy);
			id = getImageID;
			thisValley = 0;
			thisPeak = 0;
			thisBorder = 0;
			thisCenter = 0;
			for(loop = 1; loop <=2; loop++){
				selectImage(id);
				pos1 = round(xx.length * 0.2);
				pos2 = maxOf(round(xx.length * 0.4), pos1+2);
				run("Line Width...", "line=&lWidth");
				setLineWidth(lWidth);
				makeLine(xx[pos1], yy[pos1], xx[pos2], yy[pos2]);
	
				run("Straighten...", "title=straight line=&lWidth");
				run("Select All");
				setKeyDown("alt");
				profile = getProfile;
				close;
				Array.getStatistics(profile, min, max, mean, stdDev);
				tolerance = max * relTolerance;
				profile = Array.resample(profile, 100);
				indexesP = Array.findMaxima(profile,  tolerance, 1);
				indexesV = Array.findMinima(profile,  tolerance, 1);
				offset = 0;
				twoPeaks = (indexesP.length == 2);
				if(twoPeaks){//two p
					nTwo++;
					centerPos = round((indexesP[0] + indexesP[1])/2);
					offset = 100/2 - centerPos;//align center
					if(offset > 10) offset = 10;
					if(offset < -10) offset = -10;
					if(indexesV.length == 1){
						thisValley += profile[indexesV[0]];
						thisPeak += (profile[indexesP[0]] + profile[indexesP[1]])/2;
					}
				}
				else{
					thisCenter += profile[100/2];
					thisBorder += (profile[100/2 - cellDia/2] + profile[100/2  + cellDia/2])/2;
					nOne++;
				}
				for(jj = 0; jj < 100; jj++){
					kk = maxOf(0, jj - offset);
					kk = minOf(100-1, kk);
					if(twoPeaks)
						profile2P[jj] += profile[kk];
					else
						profile1P[jj] += profile[kk];

				}
		
				if(is("Caps Lock Set")){
					selectImage(id);
					setBatchMode("show");			
					waitForUser;
					setBatchMode("hide");
				}			
		
				Array.reverse(xx);//other end of cell 
				Array.reverse(yy);
			}
			
			if(twoPeaks){
				thisRatio2P  = thisValley/thisPeak;//can be from one or both cell halves
				ojSetResult("CB2-Ratio", obj, thisRatio2P);
				ojSetResult("CB-Ratio", obj, thisRatio2P);
			}
			else{
				thisRatio1P  = thisCenter/thisBorder;//can be from one or both cell halves
				ojSetResult("CB1-Ratio", obj, thisRatio1P);
				ojSetResult("CB-Ratio", obj, thisRatio1P);
			}
		
		}
	}
// debug;
 	setBatchMode(false);
	for(jj = 0; jj < 100; jj++)
		profileT[jj] = profile1P[jj] + profile2P[jj];//weighted
	Array.getStatistics(profile1P, min, max1, mean, stdDev);
	Array.getStatistics(profile2P, min, max2, mean, stdDev);
	Array.getStatistics(profileT, min, maxT, mean, stdDev);


	for(jj = 0; jj < 100/2; jj++){//make profile1P and profile2P symmetrical
		profile2P[jj] = (profile2P[jj] + profile2P[100 - jj -1])/2/max2;
		profile2P[100 - jj -1] = profile2P[jj];
		profile1P[jj] = (profile1P[jj] + profile1P[100 - jj -1])/2/max1;
		profile1P[100 - jj -1] = profile1P[jj];
		profileT[jj] = (profileT[jj] + profileT[100 - jj -1])/2/maxT;
		profileT[100 - jj -1] = profileT[jj];
	}
	print("\n---\npath=", ojGetProjectPath()+ojGetProjectName());
	xx = newArray(2);//create x-axis
	xx[0] = -lWidthUM/2;
	xx[1] = lWidthUM/2;
	xx = Array.resample(xx, 100);
	score2 = nTwo/(nTwo + nOne);
	score1 = d2s((1-score2) * 100, 1) + "%";
	score2 = d2s(score2 * 100, 1) + "%";
	setFont("SanSerif", 12, "antialiased");
  	setColor("magenta");
	if(nOne > 0){
		ratio1 = 1/profile1P[25];
		Plot.create("One Peak", "Off-Axis [um]", "Rel. Fluorescence", xx, profile1P);
		print("\One-Peak score: " + score1);
		print("Center/border ratio R1 = " + d2s(ratio1, 2));
		Plot.show;
		Overlay.drawString("Score1 = " + score1  + "  CB1-ratio = " +  d2s(ratio1, 2), 64, 14);
		Overlay.show;
	}
	if(nTwo > 0){
		ratio2 = profile2P[100/2];//the center
		Plot.create("Two Peaks", "Off-Axis [um]", "Rel. Fluorescence", xx, profile2P);
		print("Two-peak score: ", score2);
		print("Center/border ratio: R2= ", d2s(ratio2, 2));
		Plot.show;
		Overlay.drawString("Score2 = " + score2 + "  CB2-ratio = " +  d2s(ratio2, 2), 64, 14);
		Overlay.show;
	}
	if(nOne > 0 && nTwo > 0){
		Plot.create("OneP+TwoP", "Off-Axis [um]", "Rel. Fluorescence", xx, profileT);
		Plot.show;
		Overlay.show;
	}
	selectWindow("Log");
	close("Plot Values");
}