// Varve counter by Janne Jääskeläinen
// Department of geology and mineralogy
// University of Turku 2006

showStatus("Starting varveCounter");
// ImageJ requirements
requires("1.34m");
	
// If "Log" window is open, we close it now...
if (isOpen("Log")) {
selectWindow("Log");
  	run("Close");
}

// First we open dialog for user input.
// User defines the smallest distance between different varves
// (in pixels).
width=512; height=512;
Dialog.create("Varve counter");
Dialog.addNumber("Smallest distance between varves (in pixels):", 50);
Dialog.addCheckbox("Show full result data", false);
Dialog.addCheckbox("Show plot", false); 	
Dialog.addCheckbox("Draw varves", true);
Dialog.addNumber("Line thickness (in pixels):", 4);  	 	
Dialog.show();
scope = Dialog.getNumber();
showFullData = Dialog.getCheckbox();
showPlot = Dialog.getCheckbox();
showVarves = Dialog.getCheckbox();
thick = Dialog.getNumber(); 
 	 		
// Get profile and calculate the slope data in "Log" window
profile = getProfile();
// Average values of x in the defined scope area
xAvg = newArray(profile.length);
// Average values of y in the defined scope area
yAvg = newArray(profile.length);
// y[i] - yAvg(in scope area)
ySumma = newArray(profile.length);
// x[i] - xAvg(in scope area) 
xSumma = newArray(profile.length);
// (x[i] - xAvg(in scope area)) * (y[i] - yAvg(in scope area))
xy = newArray(profile.length);
// (x[i] - xAvg(in scope area)) * (x[i] - xAvg(in scope area))
xx = newArray(profile.length);
// [sigma(x[i] - xAvg(in scope area)) * (y[i] - yAvg(in scope area))] / 
// [sigma(x[i] - xAvg(in scope area)) * (x[i] - xAvg(in scope area))]
slope = newArray(profile.length);
// 0 or 1, depending on the sign of the slope ("+" or "-")
slopeChange = newArray(profile.length);
// if 0 followed 0 or 1 followed by 1, then 0, else 1.
varveEstim = newArray(profile.length);
showProgress(0.0);

// Calculates average of x and y in the user defined scope area
showStatus("Calculating slope data...");
for (i=0; i<profile.length; i++){
xTotal = 0;
 	yTotal = 0;	
 	amis = scope+i;
 	if(amis > profile.length){amis = profile.length;}
  	for (j=i; j< amis; j++){
 		xTotal += j;
 		yTotal += profile[j]; 
 	}
 	xAvg[i] = xTotal/scope;
 	yAvg[i] = yTotal/scope; 
}
showProgress(0.14);
// Calculates "x[i] - xAvg(in scope area)" and "y[i] ? 
// yAvg(in scope area)".
for (i=0; i<profile.length; i++){
	xSumma[i] = i-xAvg[i];
	ySumma[i] = profile[i]-yAvg[i];
}
	
// Calculates "(x[i] - xAvg(in scope area)) * (y[i] - yAvg(in 
// scope area))" and "(x[i] - xAvg(in scope area)) * (x[i] ? xAvg
// (in scope area))".
for (i=0; i<profile.length; i++){
	xy[i] = xSumma[i]*ySumma[i];
	xx[i] = xSumma[i]*xSumma[i];
}
	
showProgress(0.28); 	
// Calculates the slope parameter in the user defined scope area	
for (i=0; i<profile.length; i++){ 	
	xySummis = 0;
	xxSummis = 0;
	damis = scope+i;
// Here we make sure that the method doesn't return error.
	if(damis > profile.length){damis = profile.length;} 
  	for (j=i; j< damis; j++){
  		xySummis += xy[j];
  		xxSummis += xx[j];
 	}
 	slope[i] = xySummis/xxSummis;
} 	
showProgress(0.42);
	
// Now we can calculate from each slope parameter 0 or 1 (is it bigger
// or smaler than zero).
for (i=0; i<profile.length; i++){
	if(slope[i] <= 0){
	slopeChange[i] = 0;
 	}else{
		slopeChange[i] = 1;
 	}
}
showProgress(0.56);
showStatus("Counting varves...");	
// If the slope changes from negative to positive (or otherways), we
// add 1 to array.
varveCount = 0;
for (i=0; i<profile.length; i++){
	adeis = i-1;
// Here we make sure that the method doesn't return error.
	if(adeis < 0){adeis = 0;}
	if(slopeChange[adeis] != slopeChange[i]){
		varveEstim[i] = 1;
		varveCount += 1;
	}else{
		varveEstim[i] = 0;
	}
}
	
showProgress(0.70);
// And final varve count is changes in plot directions /2.
finalVarveCount = floor(varveCount/2);
 	
// Output window
showStatus("Printing output...");
print("Varve counter by Janne Jääskeläinen");
print("Department of geology and mineralogy");
print("University of Turku 2006");
 	
showProgress(0.89);
print("-----------------------------------------------------");	
if(showFullData==true){
print("X: | Y: | xSumma: | ySumma: | XY: | XX: | Slope: | 
SlopeChange:|");
print("-----------------------------------------------------");
for (i=0; i<profile.length; i++){
print(i+"|"+profile[i]+"|"+xSumma[i]+"|"+ySumma[i]+"|"+
xy[i]+"|"+xx[i]+"|"+slope[i]+"|"+varveEstim[i]);
	} 
print("-----------------------------------------------------");
} 
print("Estimated varves: "+finalVarveCount); 	 	
showProgress(1.0);
 	
// Plot profile
if(showPlot==true){
	showStatus("Plotting profile...");
 	Plot.create("Profile", "X", "Value", profile);
}
	
// Drawing varves
if(showVarves==true){ 	
setColor(255,0,0);
	setLineWidth(thick);
	for (i=0; i<profile.length; i++){
		if(varveEstim[i] == 1){
		if(slopeChange[i]<=0){
				abo = getHeight();
				moveTo(i,0);
				lineTo(i,abo);				
			} 
		}
	}		
}	
showStatus("Done.");