// Global variables
var n_buttons = 9;
var code = newArray(n_buttons), choice = newArray(n_buttons); 
var D = "       ";
var next = "";
var n_contour=0, add_contour = "";
var leftButton=16, rightButton=4;
var imORIGINAL=0, imMULTITEMP=0, imMULTIMASK=0;
var folder="", filename = "", filename_noextension = "", filename_extension = ""; 
var MaxCoord = 8000, joinedX = newArray(8000), joinedY = newArray(8000);
var step = 50; // step of the color-dashed line of joined contours
var maxpoints = 10000;
var EQUI_ptsX = newArray(8000), EQUI_ptsY = newArray(8000);
var width, height;
var SRFT = 0.25; // semiRange of the first threshold
var SRST = 0.15; // semiRange of the second threshold

k=0 ; code[k] = "OLM"; choice[k] = "Open leaf image"; 
k++ ; code[k] = "OSC"; choice[k] = "Open single contour";
k++ ; code[k] = "STS"; choice[k] = "no OL: single thresholding segmentation & contour acquisition";
k++ ; code[k] = "DTS"; choice[k] = "OL: double thresholding segmentation";
k++ ; code[k] = "TIC"; choice[k] = "OL: take inner contours - the outer last";
k++ ; code[k] = "SAL"; choice[k] = "OL: show/open all contours";
k++ ; code[k] = "JAC"; choice[k] = "OL: join all contours";
k++ ; code[k] = "DAC"; choice[k] = "OL: delete all saved contour files of this leaf";
k++ ; code[k] = "END"; choice[k] = "End";

next = "OLM,END"; // start
opz = create_dialog(next);

//   LOOP   //////////////////////////////////////////////////////////////////////////////////

while(opz != ">END")
{
	if(opz=="OLM") {	run("Close All"); open_image(); next = "OLM,OSC,STS,DTS,SAL,JAC,END"; opz = create_dialog(next);	}
	if(opz=="OSC") {	open_contour(); next = "OLM,OSC,END"; opz = create_dialog(next);	}
	if(opz=="STS") {	single_threshold(); next = "OLM,END"; opz = create_dialog(next); }
	if(opz=="DTS") {	double_threshold(); next = "OLM,TIC,END"; opz = create_dialog(next); }
	if(opz=="TIC") {	do { draw_single_contour(); } while (add_contour == "yes"); next = "OLM,TIC,SAL,JAC,END"; opz = create_dialog(next); }
	if(opz=="SAL") {  show_all_contours(); next = "OLM,JAC,END"; opz = create_dialog(next); }
	if(opz=="JAC") {	join_all_contours(); next = "OLM,DAC,END"; opz = create_dialog(next); }
	if(opz=="DAC") {	delete_all_contours(); next = "OLM,OSC,STS,DTS,END"; opz = create_dialog(next); }
	if(opz=="END") {  exit(); }
}

function clear_Roi_manager()
{
	n = roiManager("count");
	// print(n);
	if(n > 0) 
	{
		for(j = n-1; j >=0; j--) 
		{
			roiManager("select", j);
			roiManager("Delete");
		}
	}
}
			
function open_image()
{
	if (isOpen("Log")) 
	{
	  selectWindow("Log");
		run("Close");
	} 
	clear_Roi_manager();
	run("Close All");
	open("");		
	run("8-bit");
	run("Set Scale...", "distance=1 known=1 pixel=NaN unit=pixel global");		
	imORIGINAL = getImageID();
	getDimensions(width, height, channels, slices, frames);
	folder = getInfo("image.directory");
	filename = getInfo("image.filename");
	filename_noextension = substring(filename, 0, lengthOf(filename)-4);
	filename_extension = substring(filename, lengthOf(filename)-3, lengthOf(filename));
	//print(folder);
	//print(filename);
	//print(filename_noextension);
	//print(filename_extension);
}

function open_contour()
{
	selectImage(imORIGINAL);
	run("RGB Color");
	run("XY Coordinates... ","open"); 
	Roi.getCoordinates(XXsinglecontour, YYsinglecontour);
	// re-load the contour file and store the coordinates in array
	nXXsinglecontour = lengthOf(XXsinglecontour);
	nYYsinglecontour = lengthOf(YYsinglecontour);
	// base 0!!!
	contour_choice = newArray("black line", "color line (red>green>blue)", "smooth", "add contour", "back to the main menu");
	run("Select None");
	do
	{
		Dialog.create(""); Dialog.addRadioButtonGroup("Menu", contour_choice, 5, 1, contour_choice[1]); Dialog.show();
		D = Dialog.getRadioButton();
		color = getValue("color.foreground");
		
		if(D=="black line")
		{
			setColor(0,0,0);
			moveTo(XXsinglecontour[0], YYsinglecontour[0]);
			for(n = 1; n < nXXsinglecontour; n++) lineTo(XXsinglecontour[n], YYsinglecontour[n]);
			lineTo(XXsinglecontour[0], YYsinglecontour[0]);
			setForegroundColor(255,255,0);
			setForegroundColor(color);	
			// salvare?		
		}
		
		if(D=="color line (red>green>blue)")
		{
			run("Line Width...", "line=5");
			color = getValue("color.foreground");
			//print(color);
			setForegroundColor(255, 0, 0);	
			setColor(255,0,0);			
			u = 0;
			moveTo(XXsinglecontour[0], YYsinglecontour[0]);
			for(n = 1; n < nXXsinglecontour; n++) 
			{
				u++;
				if(u == step) setColor(0,255,0);		
				if(u == (step*2)) setColor(0,0,255);		
				if(u == (step*3))
				{
					setColor(255,0,0);
					u = 0;		
				}
				lineTo(XXsinglecontour[n], YYsinglecontour[n]);
			}
			lineTo(XXsinglecontour[0], YYsinglecontour[0]);
			setForegroundColor(color);		
			run("Line Width...", "line=1");
			// to be saved?
		}
	
		if(D=="smooth") // equally spaced coordinates
		{
			stringa  = "actual number of points: " + toString(nXXsinglecontour) + " - how many equally spaced points?";
			NumEquiPoints = getNumber(stringa, 1000);
			NumEquiPoints = equispac(XXsinglecontour, YYsinglecontour, nXXsinglecontour, NumEquiPoints);
			makeSelection("freehand", EQUI_ptsX, EQUI_ptsY, NumEquiPoints);
			new_filename = folder + filename_noextension + " - " + toString(NumEquiPoints) + "-equicontour.txt";
			saveAs("XY Coordinates", new_filename);
			setForegroundColor(color);					
			setTool("hand");		

			/*
			moveTo(EQUI_ptsX[0], EQUI_ptsY[0]);
			for(n = 1; n < NumEquiPoints; n++) lineTo(EQUI_ptsX[n], EQUI_ptsY[n]);
			lineTo(EQUI_ptsX[0], EQUI_ptsY[0]);
			setForegroundColor(color);		
			*/
		}
		if(D=="add contour")
		{
			open_contour(); // recursive!!!
		}
	} while (D != "back to the main menu")	
}
	
function single_threshold()
{
	clear_Roi_manager();
	selectImage(imORIGINAL);
	run("Duplicate...", " ");
	run("Auto Threshold", "method=Default");
	imSOGLIA = getImageID();
	risp = getBoolean("OK?");
	if(risp == 0)
	{
		selectImage(imSOGLIA);	
		close();
	}
	if(risp == 1)
	{
		new_filename = folder + filename_noextension + " -mask.tif";
		// print(new_filename);
		saveAs("Tiff", new_filename);
		run("Clear Results");	
		run("Set Measurements...", "centroid redirect=None decimal=3");
		run("Measure");
		centrX = getResult("X", 0);
		centrY = getResult("Y", 0);
		// run("Invert");
		setTool("wand");
		doWand(centrX, centrY);
		print("centrX=" + centrX + "   centrY=" + centrY);
		waitForUser("if the contour didn't appear automatically, make it now using the 'magic wand' and then click in this window");
		new_filename = folder + filename_noextension + " -contour.txt";
		saveAs("XY Coordinates", new_filename);
		setTool("hand");
		

	}
}
		
function double_threshold()
{
	clear_Roi_manager();
	selectImage(imORIGINAL);
	// HISTOGRAM
	// =============================================================
	Classi = newArray(256);
	Freq = newArray(256);
	nBins = 256;
	MaxGray = 255;
	getHistogram(Classi, Freq, nBins);
	
	// SMOOTHING DELLE FREQUENZE
	// =============================================================
	FreqSmooth = newArray(255);
	nSmooth = 17; // must be odd
	// the range will be from Smooth_from to Smooth_to
	Smooth_from = floor(nSmooth/2);
	Smooth_to = MaxGray - floor(nSmooth/2) +1;
	Smooth_wing = Smooth_from - 1;
	for (j = Smooth_from; j <= Smooth_to; j+=1)
	{	
		sum = 0;
		for(q = j - Smooth_wing; q <= j + Smooth_wing; q++)
		{
			sum = sum + Freq[q];
		}
		FreqSmooth[j] = sum/nSmooth;  // average
	}
	
	// save the file of the original histogram and the smoothed histogram
	new_filename = folder + filename_noextension + " -original hist.txt";
	f = File.open(new_filename);
	for (j = 0; j <= 255; j+=1)
	{
		print(f, j + "\t" + Freq[j]);
	}
	File.close(f); // output files must be closed as soon as possible
	
	new_filename = folder + filename_noextension + " -smoothed hist.txt";
	f = File.open(new_filename);
	for (j = Smooth_from; j <= Smooth_to; j+=1)
	{
		print(f, j + "\t" + FreqSmooth[j]);
	}
	File.close(f);

/*                              
 17terms smoothed histogram  
the darkest first peak is due to lobes overlapping
the mid-dark second peak is due to the normal leaf density
the clear third peak is due to the empty background

			1 peak                    2 peak                                 3 peak	

                                 #
                                ***
       #                       ******                                   *
      **                     **********                               *****
     ****                 ***************                           **********
    *********            ********************                     **************
  ************  ****** ******************************   ****   ***********************   
**************#*******#******************************#*******#****************************
_____max1___min1a___min1b______max2________________min2a___min2b______________________________
*/

	// FROM THE START, LOOK FOR max1
	// =============================================================
	max1 = 0;
	for (j = Smooth_from+3; j <= Smooth_to; j+=1)
	{	
		if (FreqSmooth[j] < FreqSmooth[j-1])
		{
			j = Smooth_to + 1; // to exit the for loop 
			// having found FreqSmooth[j] < FreqSmooth[j-1], 
			// this means that the maximum has been already found
		}	
		if (j <= Smooth_to)
		{	
			// print(j + "\t" + "FreqSmooth[j]" + "\t" + FreqSmooth[j] + "\t" + "max1" + "\t" + max1);
			
			if (FreqSmooth[j] > max1)
			{
				max1 = FreqSmooth[j];
				max1_n = j;  // max1_n is the index
			}
		}
	}	
	if (max1 == 0)
	{
		// error	
		exit("explored up to the gray level " + j +" - max1 not found");
	}
	print ("max1: " + max1 + "with index: "+ max1_n);
	
	// FORM max1_n TO max1_n x 2, LOOK FOR min1a
	// =============================================================
	min1a = 1000000;
	min1a_n = 1000000;
	for (j = max1_n +1; j <= max1_n * 2; j+=1)
	{	
		// print("value =" + FreqSmooth[j]);
		if (FreqSmooth[j] > FreqSmooth[j-1])
		{
			j = max1_n * 2 + 1; // to exit the for loop 
			// having found FreqSmooth[j] > FreqSmooth[j-1]), 
			// this means that the min1a has been already found
		}	
		if (j <= max1_n * 2)
		{
			if (FreqSmooth[j] < min1a)
			{
				min1a = FreqSmooth[j];
				min1a_n = j;  // min1a_n is the index of the first minimum
			}
		}
	}
	print ("min1a: " + "\t" + min1a + "\t" + "with index: " + "\t"+ min1a_n);
	
	if (min1a == 1000000)
	{
		// errore
		exit("explored up to the gray level " + j +" - min1a not found");
	}
	 
	// FORM min1_n TO max_n x 2, LOOK FOR min2
	// =============================================================
	min1b = 1000000;
	min1b_n = 1000000;
	for (j = min1a_n+1; j <= max1_n * 2; j+=1)
	{	
		if(FreqSmooth[j] < min1b)
		{
			min1b = FreqSmooth[j];
			min1b_n = j;  // min1b_n is the index of min1b
		}
	}
	
	if (min1b == 1000000)
	{
		print ("min1b not found - use only min1a");	
		first_threshold = (min1_n);
	}
	else
	{
		print ("min1b: " + "\t" + min1b + "\t" + "with index: " + "\t" + min1b_n );
		first_threshold = (min1a_n + min1b_n)/2;
	}
	
	// FIRST THRESHOLD
	print("first threshold: " + "\t" + first_threshold);
	
	// FROM max1_n * 2, LOOK FOR max2
	// =============================================================
	max2 = 0;
	for (j = max1_n * 2; j <= 150; j+=1)
	{	
		if (FreqSmooth[j] < FreqSmooth[j-1])
		{
			j = Smooth_to + 1; // to exit the for loop 
			// having found FreqSmooth[j] < FreqSmooth[j-1]), 
			// this means that max2 has been already found
		}	
		if (j <= Smooth_to)
		{
			if (FreqSmooth[j] > max2)
			{
				max2 = FreqSmooth[j];
				max2_n = j;  // max_n is the index
			}
		}
	}
	print ("max2: " + "\t" + max2 + "\t" + "with index: " + "\t" + max2_n);
	
	if (max2 == 0)
	{
		// errore	
		exit("explored up to the gray level " + j +" - max2 not found");		
	}
	
	// FROM max2_n TO max2_n x 2, LOOK FOR min2a
	// =============================================================
	min2a = 1000000;
	min2a_n = 1000000;
	for (j = max2_n +1; j <= max2_n * 2; j+=1)
	{	
		// print("value =" + FreqSmooth[j]);
		if (FreqSmooth[j] > FreqSmooth[j-1])
		{
			j = max2_n * 2 + 1; // to exit the for loop 
			// having found FreqSmooth[j] > FreqSmooth[j-1]), 
			// this means that min2a has been already found
		}	
		if (j <= max2_n * 2)
		{
			if (FreqSmooth[j] < min2a)
			{
				min2a = FreqSmooth[j];
				min2a_n = j;  // min2a_n is the index
			}
		}
	}
	print ("min2a: " + "\t" + min2a + "\t" + "with index: " + "\t" + min2a_n);
	
	if (min2a == 1000000)
	{
		// errore
		exit("explored up to the gray level " + j +" - min2a not found");		
	}
	
	// waitForUser("  ");
	
	// FROM min2a_n TO max2_n x 2, LOOK FOR min2b
	// =============================================================
	min2b = 1000000;
	min2b_n = 1000000;
	for (j = min2a_n+1; j <= max2_n * 2; j+=1)
	{	
		// DEVE TROVARE IL MINIMO GENERALE (min2), SENZA ALTRE CONDIZIONI
		if(FreqSmooth[j] < min2b)
		{
			min2b = FreqSmooth[j];
			min2b_n = j;  // min2_n è l'indice del secondo minimo
		}
	}
	
	if (min2b == 1000000)
	{
		print ("min2b not found - use only min2a");	
		second_threshold = (min2a_n);
	}
	else
	{
		print ("min2b: " + "\t" + min2b + "\t"  + "with index: " + "\t" + min2b_n);
		second_threshold = (min2a_n + min2b_n)/2;
	}	
	// SECOND THRESHOLD
	print ("second_threshold " + "\t" + second_threshold);
	
    /*
     nine images combining the ranges of the first and second thrashold
     
	                                     first threshold
	                                  x       x         x
	                               -SRFT      1      +SRFT  
	                            ____________________________  
	                   x -SRST  |    1        4        7
	                            |
	 seconda soglia    x     1  |    2        5        8  
	                            |      
	                   x +SRST  |    3        6        9 

  SRFT is assigned in the header; // semiRange of the first threshold
  SRST is assigned in the header; // semiRange of the second threshold
	                         
	 */
	
	deltas1 = first_threshold * SRFT;
	deltas2 = second_threshold * SRST;
	imageOpz = newArray(10);
	for(q = 1; q <= 3; q++) 
	{
		s1 = first_threshold - 2 * deltas1 + q * deltas1;
		for(h = 1; h<= 3; h++)
		{
			s2 = second_threshold - 2 * deltas2 + h * deltas2;
			selectImage(imORIGINAL);
			run("Duplicate...", " ");
			imMULTITEMP = getImageID();
			setThreshold(0, s1);
			setOption("BlackBackground", false);
			run("Convert to Mask");
			//waitForUser("1");
			// run("Image Calculator...");
			imageCalculator("OR create", imORIGINAL, imMULTITEMP);
			//waitForUser("OR");
			if(q == 1 && h == 1) imageOpz[1] = getImageID(); 
			if(q == 1 && h == 2) imageOpz[2] = getImageID();
			if(q == 1 && h == 3) imageOpz[3] = getImageID();						
			if(q == 2 && h == 1) imageOpz[4] = getImageID();						
			if(q == 2 && h == 2) imageOpz[5] = getImageID();						
			if(q == 2 && h == 3) imageOpz[6] = getImageID();					
			if(q == 3 && h == 1) imageOpz[7] = getImageID();						
			if(q == 3 && h == 2) imageOpz[8] = getImageID();						
			if(q == 3 && h == 3) imageOpz[9] = getImageID();
			s1 = round(s1); s2 = round(s2); 	
			stringa = "n=" + toString((q-1)*3+h) + ": "+ toString(s1) + "-" + toString(s2);
			rename(stringa);					
			setThreshold(0, s2);
			run("Convert to Mask");
			//waitForUser("2");			
			run("Set... ", "zoom=75 x=0 y=0");
			selectImage(imMULTITEMP);  // closes the first threshold mask
			close();
		}
	}
	run("Tile");
	showMessage("click OK and then select the best mask");		
	esce_dal_loop = 0;
	while (esce_dal_loop == 0) 
	{
		getCursorLoc(x, y, z, flags);
		if (flags&leftButton!=0)
		{
			im = getImageID();
			for(p = 1; p <= 9; p++)
			{ 
				if(im == imageOpz[p])
				{
					imMULTIMASK = im;
					for(z = 1; z < p; z++)
					{
						selectImage(imageOpz[z]);
						close();
						imageOpz[z] = 999;
					}
					for(z = p + 1; z <= 9; z++)
					{
						selectImage(imageOpz[z]);
						close();
						imageOpz[z] = 999;
					}
					esce_dal_loop = 1; 
				}
			}
		}
	}
	selectImage(imORIGINAL);
	run("View 100%");
	run("Scale to Fit");
	setLocation(30,30);
	selectImage(imMULTIMASK);
	// salva 
	new_filename = folder + filename_noextension + " -multimask.tif";
	saveAs("Tiff", new_filename);	
	run("View 100%");
	run("Scale to Fit");
	setLocation(1,1);
}
			
function draw_single_contour()
{
	selectImage(imMULTIMASK);			
	run("Select None");
	setTool("wand");
	// print("ok wand");
	x = 0; y = 0;
	while(x == 0 && y == 0)
	{
		getSelectionBounds(x, y, w, h);	
		//print("Roi");
		//print(x);
		//print(y);								
	}
	risp = getBoolean("to be saved?");
	if(risp == 0)
	{
		selectImage(imMULTIMASK);	
		run("Select None");
	}
	if(risp == 1)
	{
		selectImage(imMULTIMASK);	
		roiManager("Add");
		n_contour++;
		new_filename = folder + filename_noextension + " -contour " + toString(n_contour) + ".txt";
		saveAs("XY Coordinates", new_filename);
		run("Set...", "value=150");
		run("Select None");
	}
	risp = getBoolean("one more contour?");
	if(risp == 1) add_contour = "yes";
	else add_contour = "no";
}
			
function show_all_contours()
{	
	n_saved_contours = 0;	
	run("RGB Color");	
	for(j = 1; j <= 100; j++)
	{
		new_filename = folder + filename_noextension + " -contour " + toString(j) + ".txt";
		file_esiste = File.exists(new_filename);
		if (file_esiste == 1) n_saved_contours++;
	}
	for(j = 1; j <= n_saved_contours; j++)
	{
		// Read all txt files which contin the x and y coordinates
		// previously obtained with the command SaveAs XY Coordinates
		new_filename = folder + filename_noextension + " -contour " + toString(j) + ".txt";
		string = File.openAsString(new_filename);
		lines = split(string, "\t\n");
		n_lines = lengthOf(lines);
		n_coord = n_lines/2;
		xcoord = newArray(n_coord+1);
		ycoord = newArray(n_coord+1);
		z = 0;
		sum = 0;
		for(n = 0; n < n_lines; n+=2) 
		// the [z] index of contours is base 1
		// the element [0] has the number of coordinates
		{
			z++;
			xcoord[z] = lines[n];
			ycoord[z] = lines[n+1];
		}
		xcoord[0] = z;
		ycoord[0] = z;
		
		// draw the contour
		color = getValue("color.foreground");
		//print(color);
		setColor(255,0,0);
		setForegroundColor(255, 0, 0);
		run("Colors...", "foreground=red background=white selection=red");
		moveTo(xcoord[1], ycoord[1]);
		for (i = 2; i <= xcoord[0]; i++) 
		{
			lineTo(xcoord[i], ycoord[i]);
		}
		// to close the contour
		lineTo(xcoord[1], ycoord[1]);
		setForegroundColor(color);
	}
	risp = getBoolean("image to be saved?");
	if(risp == 1)
	{
		new_filename = folder + filename_noextension + " -contours.tif";
		saveAs("Tiff", new_filename);	
	}
}
	
function join_all_contours()
{
	run("Close All");
	n_saved_contours = 0;
	new_filename = folder + filename_noextension + "." + filename_extension;
	open(new_filename);
	zoom = getZoom();
	run("8-bit");	// to make a gray image that can be eventually colored
	run("RGB Color");	
	imSET_CROSSING_POINTS = getImageID();
	rename("image to set crossing points");	
	run("Duplicate...", " ");
	imJOINED_CONTOURS = getImageID();	
	rename("image to check the correct joining of contours");
		
	selectImage(imSET_CROSSING_POINTS);		
	for(j = 1; j <= 100; j++)
	{
		new_filename = folder + filename_noextension + " - contour " + toString(j) + ".txt";
		file_esiste = File.exists(new_filename);
		if (file_esiste == 1) n_saved_contours++;
	}
	for(j = 1; j <= n_saved_contours; j++)
	{
		// Read all txt files which contin the x and y coordinates
		// previously obtained with the command SaveAs XY Coordinates
		new_filename = folder + filename_noextension + " -contour " + toString(j) + ".txt";
		string = File.openAsString(new_filename);
		lines = split(string, "\t\n");
		n_lines = lengthOf(lines);
		n_coord = n_lines/2;
		xcoord = newArray(n_coord+1);
		ycoord = newArray(n_coord+1);
		z = 0;
		sum = 0;
		for(n = 0; n < n_lines; n+=2) 
		// the [z] index of contours is base 1
		// the element [0] has the number of coordinates
		{
			z++;
			xcoord[z] = lines[n];
			ycoord[z] = lines[n+1];
			if(z > 1) sum = sum + ((xcoord[z] - xcoord[z-1]) * (ycoord[z] + ycoord[z-1]));
		}
		xcoord[0] = z;
		ycoord[0] = z;

		// draw the contour
		color = getValue("color.foreground");
		//print(color);
		setColor(0,255,0);
		setForegroundColor(0, 255, 0);
		run("Colors...", "foreground=green background=white selection=red");
		moveTo(xcoord[1], ycoord[1]);
		for (i = 2; i <= xcoord[0]; i++) 
		{
			//if(i < 15) waitForUser("  "); // to stop the drawing
			lineTo(xcoord[i], ycoord[i]);
			//updateDisplay();
		}
		// to close the contour
		lineTo(xcoord[1], ycoord[1]);
		setForegroundColor(color);
	}
	// ImageJ doesn't support 2D arrays. To obviate, I create two very long arrays XX and YY with 
	// a number of elements given by: MaxCoord x the number of arrays  
	// I suppose that the total contour (produced by joining all single contours) 
	// is not longer than MaxCoord (now set to 8000).
	// The indices of each array will be from 1 to MaxCoord + (MaxCoord x (n° array -1))
	// that is, for example, for 5 contours, having MaxCoord = 8000,
	// the indices of the 5 arrays will be 
	// 1-8000, 8001-16000, 16001-24000, 24001-32000, 32001-40000.
	
	// For es. the x e y incoordinates at the position 64 of the q.th array will be in 
	// XX[64 + MaxCoord * (q-1)]  and  YY[64 + MaxCoord * (q-1)]
	
	// Note that: array = contour
	
	XX = newArray(MaxCoord * n_saved_contours + 1); 
	YY = newArray(MaxCoord * n_saved_contours + 1);

	nXX = newArray(n_saved_contours + 1);
	nYY = newArray(n_saved_contours + 1);
	// If the q.th array contains 300 elements
	// nXX[q] = 300 and nYY[q] = 300;	
	
	// Reads the files and put data in arrays
	for(j = 1; j <= n_saved_contours; j++)
	{
		new_filename = folder + filename_noextension + " -contour " + toString(j) + ".txt";
		string = File.openAsString(new_filename);
		lines = split(string, "\t\n");
		n_lines = lengthOf(lines);
		z = 0;
		for(n = 0; n < n_lines; n+=2) 
		{
			z++;
			XX[z + MaxCoord * (j-1)] = lines[n];
			YY[z + MaxCoord * (j-1)] = lines[n+1];
		}
		nXX[j] = z;
		nYY[j] = z;
	}
	// draw contours
	// ... first as selections, adding labels 
	run("Colors...", "foreground=blue background=yellow selection=red");		
	for(j = 1; j <= n_saved_contours; j++)
	{
		xcoord=newArray(nXX[j]+1);
		ycoord=newArray(nXX[j]+1);
		for(n = 1; n <= nXX[j]; n++)
		{ 
			xcoord[n-1] = XX[n + MaxCoord * (j-1)];
			ycoord[n-1] = YY[n + MaxCoord * (j-1)];
		}
		makeSelection("freehand", xcoord, ycoord, nXX[j]);
		// waitForUser(" ");
		run("Set Measurements...", "centroid redirect=None decimal=3");
		run("Measure");
		centrX = getResult("X", 0);
		centrY = getResult("Y", 0);
		print(centrX);
		print(centrY);
		IJ.deleteRows(0, 0);
		setFont("SansSerif", 22, " antialiased");
		setColor("green");
		stringa = toString(j);									
		Overlay.drawString(stringa, centrX-10, centrY+15, 0.0);
		Overlay.show();
		setForegroundColor(0, 0, 255);
		run("Draw", "slice");
	}	
	run("Select None");
	// waitForUser("    ");
	
	// ... then as lines
	run("Colors...", "foreground=white background=white selection=red");				
	for(j = 1; j <= n_saved_contours; j++)
	{
		moveTo(XX[1 + MaxCoord * (j-1)], YY[1 + MaxCoord * (j-1)]);
		for(n = 2; n <= nXX[j]; n++) 
		{
			lineTo(XX[n + MaxCoord * (j-1)], YY[n + MaxCoord * (j-1)]);
		}
		lineTo(XX[1 + MaxCoord * (j-1)], YY[1 + MaxCoord * (j-1)]);
	}
	setForegroundColor(color);
	// array per i confronti
	contour=newArray(n_saved_contours + 1);
	distance=newArray(n_saved_contours + 1);
	location=newArray(n_saved_contours + 1);	

	// here starts the big loop to join all contours
	contours_to_join = n_saved_contours;
	while(contours_to_join > 1)
	{
		selectImage(imJOINED_CONTOURS);
		setLocation(1,1);	
		selectImage(imSET_CROSSING_POINTS);
		setLocation((zoom*width)+30,1);	
		point_ok = 0;
		setTool("point");

		while(point_ok == 0) 
		{
			waitForUser("clock OK\n then click on the image on the right\n the crossing point between two contours to join\n the order is arbitrary\n the joined contours will appear in the image on the left");
			// indico il punto 
			selectImage(imSET_CROSSING_POINTS);
			do
			{
				getCursorLoc(x, y, z, flags);
				if(flags == 16) 
				{
					xSelez = x; ySelez = y;
					risp = getBoolean("ok?");
					if(risp == 1) 
					{
						point_ok = 1;

					}
					run("Select None");	
				}
			} while(flags != 16)
		}
		print("crossing point selected", xSelez, ySelez);
		
		// esamino tutti i contorni
		for(j = 1; j <= n_saved_contours; j++)
		{
			// to join two contours, 
			// one of them extends updating nXX[j] = ...
			// the other one is deleted by setting nXX[j] = 0 
			if(nXX[j] != 0)
			{
				distmin = 10000000; 
				
				// Find the distance between each contour and the crossing point indicated 
				for(n = 1; n <= nXX[j]; n++) 
				{
					// to calculate the shortest distance it is not necessary the extract the square root
					// but doing so the values will be much more big
					dist = (pow((XX[n + (MaxCoord * (j-1))] - xSelez),2) + pow((YY[n + (MaxCoord * (j-1))] - ySelez),2));
					if(dist < distmin) 
					{
						distmin = dist;
						pdistmin = n;
						// waitForUser("riga 808");
					}
				}
				distance[j] = distmin;
				location[j] = pdistmin;
			}
			print("contour ", j, "  -  distance ", distance[j], "   -   location ", location[j]);
		}
		
		// Find the two contours most close to the crossing point
		// contornA is the first closest
		distA = 1000; contornA = 0;
		for(j = 1; j <= n_saved_contours; j++)
		{			
			if(nXX[j] != 0)
			{
				if(distance[j] < distA)
				{
					distA = distance[j];
					contornA = j;
					locationA = location[j];
				}
			}
		}
		if(contornA == 0) exit("contour A not found"); 
		print("primo contour =", contornA);

		// contornB is the second closest
		distB = 1000; contornB = 0;
		for(j = 1; j <= n_saved_contours; j++)
		{			
			if((nXX[j] != 0) && (j != contornA))
			{
				if(distance[j] < distB)
				{
					distB = distance[j];
					contornB = j;
					locationB = location[j];
				}
			}
		}
		if(contornB == 0) exit("contour B not found");  
		print("secondo contour =", contornB);
		
		// once the contours A and B have been identified as the closests contours to the crossing point,						
		// in A and B find two quasi-adjacent (5 pre- and 5 post-, in the direction of the array) points 
		// to the closests point of A and B to the crossing point 

		// xA1,yA1 e xA2,yA2 
		xA1 = parseFloat(XX[(locationA - 5) + MaxCoord * (contornA - 1)]);
		yA1 = parseFloat(YY[(locationA - 5) + MaxCoord * (contornA - 1)]);
		xA2 = parseFloat(XX[(locationA + 5) + MaxCoord * (contornA - 1)]);
		yA2 = parseFloat(YY[(locationA + 5) + MaxCoord * (contornA - 1)]);
		
		// xB1,yB1 e xB2,yB2 
		xB1 = parseFloat(XX[(locationB - 5) + MaxCoord * (contornB - 1)]);
		yB1 = parseFloat(YY[(locationB - 5) + MaxCoord * (contornB - 1)]);			
		xB2 = parseFloat(XX[(locationB + 5) + MaxCoord * (contornB - 1)]);
		yB2 = parseFloat(YY[(locationB + 5) + MaxCoord * (contornB - 1)]);	

		// to join the two contours...

		// consider the segment S1 between A1 and B2
		// coordinates
		xS1p1 = xA1; 
		yS1p1 = yA1;
		xS1p2 = xB2;
		yS1p2 = yB2;
		// range
		xS1min = minOf(xS1p1, xS1p2);
		yS1min = minOf(yS1p1, yS1p2);
		xS1max = maxOf(xS1p1, xS1p2);
		yS1max = maxOf(yS1p1, yS1p2);
		slopeS1 = (yS1p1 - yS1p2)/(xS1p1 - xS1p2);
		constS1 = (yS1p1 + yS1p2)/2 - slopeS1 * (xS1p1 + xS1p2)/2;
					
		// consider the segment S2 between A2 and B1
		// coordinates
		xS2p1 = xA2; 
		yS2p1 = yA2;
		xS2p2 = xB1;
		yS2p2 = yB1;
		// range
		xS2min = minOf(xS2p1, xS2p2);
		yS2min = minOf(yS2p1, yS2p2);
		xS2max = maxOf(xS2p1, xS2p2);
		yS2max = maxOf(yS2p1, yS2p2);			
		slopeS2 = (yS2p1 - yS2p2)/(xS2p1 - xS2p2);
		constS2 = (yS2p1 + yS2p2)/2 - slopeS2 * (xS2p1 + xS2p2)/2;
					
		// calculate the intersection between the two segments
		// x_inters, y_inters
		x_inters = (constS1 - constS2)/(slopeS2 - slopeS1);
		y_inters = constS1 + slopeS1 * x_inters;

		// if there is intersection 
		// contour A is open in A0 
		// and continues in contour B from punto B0, B+1, B+2,... ut to go back to B0					
		// and then to continue to A0, A+1, A+2,... up to close again in A0
		// contour B will be deleted by setting nXX[B] = 0
		 
		// otherwise, if there is no intersection
		// contour A is open in A0 
		// il contour A si spezza in A0 
		// and continues in contour B from punto B0, B-1, B-2,... ut to go back to B0				
		// per riprendere con A0, A+1, A+2,... sino a tornare su A0
		// il contour B sarà eliminato con nXX[B] = 0 

		// evaluates if there is intersection
		if(x_inters >= xS1min && x_inters <= xS1max &&
		   x_inters >= xS2min && x_inters <= xS2max &&
		   y_inters >= yS1min && y_inters <= yS1max &&
		   y_inters >= yS2min && y_inters <= yS2max)
		{ 
			t = 0;
			for(w = 1; w <= locationA; w++)
			{
				t++; 
				joinedX[t] = XX[w + MaxCoord * (contornA - 1)];
				joinedY[t] = YY[w + MaxCoord * (contornA - 1)];
			}
			for(w = locationB; w <= nXX[contornB]; w++)
			{
				t++;
				joinedX[t] = XX[w + MaxCoord * (contornB - 1)];
				joinedY[t] = YY[w + MaxCoord * (contornB - 1)];
			}
			for(w = 1; w <= locationB; w++)
			{
				t++;
				joinedX[t] = XX[w + MaxCoord * (contornB - 1)];
				joinedY[t] = YY[w + MaxCoord * (contornB - 1)];
			}
			for(w = locationA; w <= nXX[contornA]; w++)
			{
				t++; 
				joinedX[t] = XX[w + MaxCoord * (contornA - 1)];
				joinedY[t] = YY[w + MaxCoord * (contornA - 1)];
			}
			njoinedX = t; njoinedY = t;
		}				

		else
		{
			t = 0;
			for(w = 1; w <= locationA; w++)
			{
				t++; 
				joinedX[t] = XX[w + MaxCoord * (contornA - 1)];
				joinedY[t] = YY[w + MaxCoord * (contornA - 1)];
			}
			for(w = locationB; w >= 1; w--)
			{
				t++;
				joinedX[t] = XX[w + MaxCoord * (contornB - 1)];
				joinedY[t] = YY[w + MaxCoord * (contornB - 1)];
			}
			for(w = nXX[contornB]; w >= locationB; w--)
			{
				t++;
				joinedX[t] = XX[w + MaxCoord * (contornB - 1)];
				joinedY[t] = YY[w + MaxCoord * (contornB - 1)];
			}
			for(w = locationA; w <= nXX[contornA]; w++)
			{
				t++; 
				joinedX[t] = XX[w + MaxCoord * (contornA - 1)];
				joinedY[t] = YY[w + MaxCoord * (contornA - 1)];
			}
			njoinedX = t; njoinedY = t;
		}
		for(h = 1; h <= njoinedX; h++)
		{
			XX[h + MaxCoord * (contornA - 1)] = joinedX[h];
			YY[h + MaxCoord * (contornA - 1)] = joinedY[h];
		}
		nXX[contornA] = njoinedX;
		nYY[contornA] = njoinedY;
		nXX[contornB] = 0;
		nYY[contornB] = 0;
		contours_to_join--;

		// draw the joined contours
		selectImage(imJOINED_CONTOURS);
		run("Select None");		
		color = getValue("color.foreground");
		//print(color);
		setColor(255, 0, 0);
		setForegroundColor(255, 0, 0);
		run("Colors...", "foreground=white background=white selection=red");				
		moveTo(XX[1 + MaxCoord * (contornA - 1)], YY[1 + MaxCoord * (contornA - 1)]);

		u = 0;		
		setColor(0, 0, 255);
		for(n = 2; n <= nXX[contornA]; n++) 
		{
			u++;
			if(u == step) setColor(255, 0, 0);
			if(u == (step*2)) setColor(0, 255, 0);
			if(u == (step*3))
			{
				setColor(0, 0, 255);
				u = 0;		
			}
			lineTo(XX[n + MaxCoord * (contornA - 1)], YY[n + MaxCoord * (contornA - 1)]);
		}
		lineTo(XX[1 + MaxCoord * (contornA - 1)], YY[1 + MaxCoord * (contornA - 1)]);
		setForegroundColor(color);
		// waitForUser("");

		// selection to be saved as coordinates
		run("Colors...", "foreground=blue background=yellow selection=red");		
		xcoord=newArray(nXX[contornA]+1);
		ycoord=newArray(nXX[contornA]+1);
		for(n = 1; n <= nXX[contornA]; n++)
		{ 
			xcoord[n-1] = XX[n + MaxCoord * (contornA-1)];
			ycoord[n-1] = YY[n + MaxCoord * (contornA-1)];
		}
		makeSelection("freehand", xcoord, ycoord, nXX[contornA]);
		new_filename = folder + filename_noextension + " -joined_contours.txt";
		saveAs("XY Coordinates", new_filename);
		setTool("hand");		
		
		new_filename = folder + filename_noextension + " -joined_contours.tif";
		// print(new_filename);
		saveAs("Tiff", new_filename);	
		// waitForUser(" ");
	}
}

function delete_all_contours()
{
	for(j = 1; j <= 100; j++)
	{
		new_filename = folder + filename_noextension + " -contour " + toString(j) + ".txt";
		file_esiste = File.exists(new_filename);
		if (file_esiste == 1) h = File.delete(new_filename);
		if(h != 1) exit("error in deleting single contours");
	}				
}

function numero_progressivo()
{
	s = getTime();
	k = d2s(s,0);
	// print(k);
	id = toString(k,0);
	return id;
}

function clear_Roi_manager()
{
	n = roiManager("count");
	// print(n);
	if(n > 0) 
	{
		for(j = n-1; j >=0; j--) 
		{
			roiManager("select", j);
			roiManager("Delete");
		}
	}
}

function create_dialog(stringa)
{
	print(stringa);
	def = "";
	lines = split(stringa, ",");
	n_lines = lengthOf(lines);
	items = newArray(n_lines);
	p=-1;
	for(j = 0; j<n_lines; j++) 
	{
		h = lines[j];
		//print("j =", j,"     h =", h);
		for(k=0; k<n_buttons; k++)
		{	
			if(h == code[k])
			{
				p++;
				items[p] = choice[k];
			}
			if(lengthOf(def) < 1) def = items[p]; // set as default the first option
		}
	}
	scelte = newArray(p);
	D = " ";
	while(lengthOf(D) < 3)
	{		
		Dialog.create(""); Dialog.addRadioButtonGroup("Menù", items, lengthOf(items), 1, def); Dialog.show();
		D = Dialog.getRadioButton();
		print("D = ", D, "   lengthOf(D) =", lengthOf(D));
	}

	for(j=0; j<n_buttons; j++)
	{
		if(D == choice[j]) 
		{
			//print("radial button selected =", j);
			return code[j];
		}
	}
}

function equispac(array1X, array1Y, NumOrigPoints, NumEquiPoints) //===================================================================
{
	RI = newArray(128); NI = newArray(128);
	G = newArray(maxpoints); H = newArray(maxpoints); X = newArray(maxpoints); Y = newArray(maxpoints);
	
	// print("N° of original coordinates: " + toString(NumOrigPoints));
	// print("N° of equispaced coordinates: " + toString(NumEquiPoints));

	// X e Y are base 1
	// array1X e array1Y - come from ROI and are therefore base 0
	for (i = 1; i <= NumOrigPoints; i++)
	{
		// base 1
		X[i]=array1X[i-1];
		Y[i]=array1Y[i-1];
	}

	// to close the contour
	// even if X and Y are base 1, the elements [0] are temporarily set equal to the last element
	X[0]=X[NumOrigPoints];
	Y[0]=X[NumOrigPoints];	
	
	rm = 0; lm = 0; t=0;
	
	G[0] = NumEquiPoints;	
	RI[0] = 0;
	NI[0] = 0;

	t = perimeter(X,Y, NumOrigPoints, 1);
	print("perimeter: " + toString(t));
	
	rri = t/NumEquiPoints*1.3;
	rr = 0;
	flag = 0;
	
	while (flag != 1)
	{
		k = 1;
		G[1] = X[1];
		H[1] = Y[1];
		sx = 0;
		sy = 0;
		i = 0;
		
		// set the virtul cursor in G[1] e H[1]
		for (j = 2; j < NumOrigPoints; j++)
		{
			i = i+1;
	 		sx = sx+X[j]*i;
			sy = sy+Y[j]*i;
			x1 = G[k];
			x2 = sx/((i*(i+1))/2);
			y1 = H[k];
			y2 = sy/((i*(i+1))/2);
			d = sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
	 		if (d >=  rri)
	 		{ 
				k = k+1; 
				G[k] = rri/d*(x2-x1)+x1;
				H[k] = rri/d*(y2-y1)+y1;
				sx = 0;
				sy = 0;
				i = 0;
				j = j-1;
	 		}
		}
		
		d = sqrt((G[k]-X[1])*(G[k]-X[1])+(H[k]-Y[1])*(H[k]-Y[1]));

		while (d >= (rri*1.5))
		{
			x1 = G[k];
			x2 = X[1];
			y1 = H[k];
			y2 = Y[1];
	 		k = k+1;
			G[k] = rri/d*(x2-x1)+x1;
			H[k] = rri/d*(y2-y1)+y1;
			d = sqrt((G[k]-X[1])*(G[k]-X[1])+(H[k]-Y[1])*(H[k]-Y[1]));
		}
		
		if(d<rri*0.5) k = k-1;
		
		if(abs(k-NumEquiPoints)>1)
		{
			rr = rr+1;
			RI[rr] = rri;
			NI[rr] = k;
			if(rr == 1)
			{
				kn = k/NumEquiPoints;
				if (1 > kn) 
					rri = rri*(kn+kn/100);
				else
					rri = rri*(kn-kn/100);
			}
		 	else
			{  // regression 
	 			r1 = 0; r2 = 0; c1 = 0;	c2 = 0; cr = 0; n = 0;
				if (rr >= 5) 
					gg = rr-4;
				else
					gg = 1;
	
				for (j = gg; j <= rr; j++)
				{
		 			r1 = r1+RI[j];
					r2 = r2+RI[j]*RI[j];						// dipendent variable: y - lungh. 
					c1 = c1+NI[j];
					c2 = c2+NI[j]*NI[j];						// indipendente variable:	x - numero
					cr = cr+RI[j]*NI[j];
					n = n+1;
				}
				
				xy = cr-(c1*r1/n);
				xx = c2-(c1*c1/n);				 					
		 		if (xx <=  0.0001)
		 		{
		 			kn = NI[rr]/NumEquiPoints;		
		 			random("seed", 48);					
					if (1 > kn)
			 			rri = rri*(kn+(kn/(3+random())));
					else
			 			rri = rri*(kn-(kn/(3+random())));
		 		}
				else
				{
			 		b = xy/xx;
					a = (r1-b*c1)/n;
					rri = a+b*NumEquiPoints;
				}
			}
		 	
		 	
	 	 	for(j = 1; j<= rr; j++)
			{
				print(NI[j], " segments of length", RI[j]);
				t = perimeter(G,H, k, 1);
				print("perimeter: " + toString(t));
			}
		
		}
		if (abs(k-NumEquiPoints) <= 1) 
		{
			flag = 1;  // exit the cycle
			print("k = ", k, "     -        NumEquiPoints = ", NumEquiPoints); 
		}
	}

	NumEquiPoints = k;
	
	// EQUI... base 0
	// G[] e H[] always base 1
	for (j = 0; j < NumEquiPoints; j++)
	{
		EQUI_ptsX[j] = round(G[j+1]);
		EQUI_ptsY[j] = round(H[j+1]);
	}
		

	print("final number of equally spaced coordinates: " + toString(NumEquiPoints));
	print("start   x =:" + toString(EQUI_ptsX[0]) +  "   y =:" + toString(EQUI_ptsY[0]));
	print("end     x =:" + toString(EQUI_ptsX[NumEquiPoints-1]) +  "   y =:" + toString(EQUI_ptsY[NumEquiPoints-1]));	

	return NumEquiPoints;
}

function perimeter(G, H, N, base)  //===================================================================
{
	perim=0;	

	if (base == 1)
	// from 1 to N
	{	
		for (j=1; j <= (N-1); j++)
		{
			perim=perim+sqrt((G[j]-G[j+1])*(G[j]-G[j+1])+(H[j]-H[j+1])*(H[j]-H[j+1]));
		}
		perim=perim+sqrt((G[N]-G[1])*(G[N]-G[1])+(H[N]-H[1])*(H[N]-H[1]));
	}

	if (base == 0)
	// from 0 to N-1
	{	
		for (j=0; j < (N-1); j++)
		{
			perim=perim+sqrt((G[j]-G[j+1])*(G[j]-G[j+1])+(H[j]-H[j+1])*(H[j]-H[j+1]));
		}
		perim=perim+sqrt((G[N-1]-G[0])*(G[N-1]-G[0])+(H[N-1]-H[0])*(H[N-1]-H[0]));
	}
	return perim;
}