// Register_Coverslip_Nikon macro by Christophe Leterrier
// Calculates new positions based on the original and new positions of two reference points
// Used for repetitive imaging of neurons
// Works with Nikon csv stage positions files: generates a new file from an "original positions" file and a "new positions" file containing only the two new positions of the reference points
// Computed transformation is XY rotation + translation, plus a mean Z offset

// 18/11/2011 Zeiss version
// 20/02/2024 Nikon version

macro Register_Coverslip_Nikon {


	is3D = false;

	//Gets two Nikon position export files : one with the original positions, another with the new position of reference points (2 points)
	OriPath=File.openDialog("Original Positions file?");
	RefPath=File.openDialog("New Reference Positions file?");
	OriDirPath=File.getParent(OriPath);
	OriFileName=File.getName(OriPath);

	NewDirPath=File.getParent(RefPath);
	NewFileName=File.getNameWithoutExtension(OriPath)+"_transfo";

	// Gets the file name, and positions set name
	OriFile=File.openAsString(OriPath);
	OriLine=split(OriFile, "\n");

	print("\n");
	print("*** Register Coverslip Nikon Log ***");
	print("Original Positions Path: "+OriPath);

	// Gets the number of positions to initialize the original positions arrays
	OriLength=OriLine.length;
	OriPosNumber=OriLength;

	// Defines the name, X and Y original positions arrays
	OriPosName=newArray(OriPosNumber);
	OriX=newArray(OriPosNumber);
	OriY=newArray(OriPosNumber);
	if (is3D == true) OriZ=newArray(OriPosNumber);

	// Loops on all positions to feed the original positions array
	for (i=0; i<OriLine.length-1; i++) {
		SplitSemi=split(OriLine[i], ";");
		OriX[i]=parseFloat(SplitSemi[0]);
		OriY[i]=parseFloat(SplitSemi[1]);
		if (is3D == true) OriZ[i]=parseFloat(SplitSemi[2]);

		if (is3D == true)  print("Position #"+(i+1)+" at ("+OriX[i]+" , "+OriY[i]+" , "+OriZ[i]+")");
		else print("Position #"+(i+1)+" at ("+OriX[i]+" , "+OriY[i]+")");
	}

	// Does the same to get the path, set name, and position arrays for the reference positions file (only two firts points)
	RefFile=File.openAsString(RefPath);
	RefLine=split(RefFile, "\n");
	print("\n");
	print("Reference Positions Path: "+RefPath);

	RefPosName=newArray(2);
	RefX=newArray(2);
	RefY=newArray(2);
	if (is3D == true) RefZ=newArray(2);
	for (i=0; i<2; i++) {
		SplitSemi=split(RefLine[i], ";");
		RefX[i]=parseFloat(SplitSemi[0]);
		RefY[i]=parseFloat(SplitSemi[1]);
		if (is3D == true) RefZ[i]=parseFloat(SplitSemi[2]);

		if (is3D == true) print("Reference #"+(i+1)+" at ("+RefX[i]+" , "+RefY[i]+" , "+RefZ[i]+")");
		else print("Reference #"+(i+1)+" at ("+RefX[i]+" , "+RefY[i]+")");

	}

	// Initialize the original and new positions arrays.
	// IMPORTANT: assumes that the original positions of the two reference points are the first two points stored in the original positions file
	
	OriRef=newArray(4);
	OriRef[0]=OriX[0];
	OriRef[1]=OriY[0];
	OriRef[2]=OriX[1];
	OriRef[3]=OriY[1];

	NewRef=newArray(4);
	NewRef[0]=RefX[0];
	NewRef[1]=RefY[0];
	NewRef[2]=RefX[1];
	NewRef[3]=RefY[1];

	// Initializes the new positions arrays: X, Y, Names
//	NewPosName=newArray(OriPosName.length);
	NewX=newArray(OriX.length);
	NewY=newArray(OriY.length);
	if (is3D == true) NewZ=newArray(OriZ.length);

// The transformation (translation and rotation) parameters are calculated by the Affine() function
	AffineTrans=Affine(OriRef, NewRef);

	// Retrieve the parameters from the array returned by Affine()
	TransX=AffineTrans[0];
	TransY=AffineTrans[1];
	CosAngle=AffineTrans[2];
	SinAngle=AffineTrans[3];
	Angle=AffineTrans[4];

	// Processes a simple offset for the Z focus by averaging the Z of the two reference points
	if (is3D == true) {
		OffsetZ1=RefZ[0]-OriZ[0];
		OffsetZ2=RefZ[1]-OriZ[1];
		OffsetZ=(OffsetZ1+OffsetZ2)/2;
	}
	
	print("\n");
	
	if (is3D == true) print("Transformation: TransX="+d2s(AffineTrans[0],1)+", TransY="+d2s(AffineTrans[1],1)+", Angle="+d2s(Angle,4)+"°, Z-Offset="+OffsetZ);
	else print("Transformation: TransX="+d2s(AffineTrans[0],1)+", TransY="+d2s(AffineTrans[1],1)+", Angle="+d2s(Angle,4)+"°");

	// sets the file path and positions set name for the new positions

	print("\n");

	// Loops on all positions to calculate the new coordinates and change their name
	for (i=0; i<OriX.length; i++) {

		// Translation
		TempX=OriX[i]+AffineTrans[0];
		TempY=OriY[i]+AffineTrans[1];

		// Rotation around the origin
		NewX[i]=(TempX)*CosAngle-(TempY)*SinAngle;
		NewY[i]=(TempX)*SinAngle+(TempY)*CosAngle;

		// Offset for Z focus
		if (is3D == true) NewZ[i] = OriZ[i] + OffsetZ;
		
		 // Change the name of the position
//		NewPosName[i]=OriPosName[i]+"_new";
		if (is3D == true) print("New Position #" + (i+1) + " at ("+NewX[i]+" , "+NewY[i] + " , " + NewZ[i] + ")");
		else print("New Position #" + (i+1) + " at (" + NewX[i] + " , " + NewY[i] + ")");
	}



// Generate the new positions file
	// Loop through new positions to store X, and Y Z values
	NewFile = "";
	for (i=0; i<NewX.length; i++) {
		if (is3D == true) NewLine = "" + NewX[i] + ";" + NewY[i] + ";" + NewZ[i];
		else NewLine = "" + NewX[i] + ";" + NewY[i];
		NewFile = NewFile + NewLine + "\n";
	}
	
	// Saves the new file by relpacing the ref file with added new positions
	File.saveString(NewFile, NewDirPath + File.separator + NewFileName + ".csv");

	
	
	
// Creates a plot with the original and new positions	
	Plot.create("Affine Transformation", "x", "y");

	// Gets the boundaries for scaling the graph
	Array.getStatistics(OriX, MinOriX, MaxOriX, mean, stdDev);
	Array.getStatistics(OriY, MinOriY, MaxOriY, mean, stdDev);
	Array.getStatistics(NewX, MinNewX, MaxNewX, mean, stdDev);
	Array.getStatistics(NewY, MinNewY, MaxNewY, mean, stdDev);

	MinX=minOf(MinOriX, MinNewX);
	MinY=minOf(MinOriY, MinNewY);
	MaxX=maxOf(MaxOriX, MaxNewX);
	MaxY=maxOf(MaxOriY, MaxNewY);

	Plot.setLimits(MinX-10, MaxX+10,MinY-10, MaxY+10);

	Plot.setFrameSize(400, 400);

	// This builds the coordinates arrays (original and new) of the two reference points
	OriRefX=newArray(2);
	OriRefX[0]=OriRef[0];
	OriRefX[1]=OriRef[2];

	OriRefY=newArray(2);
	OriRefY[0]=OriRef[1];
	OriRefY[1]=OriRef[3];

	NewRefX=newArray(2);
	NewRefX[0]=NewRef[0];
	NewRefX[1]=NewRef[2];

	NewRefY=newArray(2);
	NewRefY[0]=NewRef[1];
	NewRefY[1]=NewRef[3];

	// Traces a green line for the two original reference points 
	Plot.setColor("green");
	Plot.add("line", OriRefX, OriRefY);
	// Traces a pink line for the two new reference points
	Plot.setColor("pink");
	Plot.add("line", NewRefX, NewRefY);
	// Traces the new positions as red circles
	Plot.setColor("red");
	Plot.add("circles", NewX, NewY);
	// Traces the original positions as blue crosses
	Plot.setColor("blue");
	Plot.add("crosses", OriX, OriY);
	// Displays the Z offset
	if (is3D == true) {
		Plot.setColor("black");
		Plot.addText("Z offset ="+OffsetZ, 0, 0);
	}
	Plot.show;

	// Saves the plot
	PlotPath=NewDirPath+File.separator+NewFileName+"_Plot.tif";
	save(PlotPath);
	// close;
	// Saves the log file
	selectWindow("Log");
	LogPath=NewDirPath+File.separator+NewFileName+"_Log.txt";
	save(LogPath);
	// run("Close");

	showStatus("Register_Coverlsip_Nikon finished");
	print("*** Register Coverslip Nikon End ***");
}

function Affine(OriPoints, NewPoints) {

	// This function is fed with the original and new coordinates of two reference points M1 and M2
	// as two arrays: original (x1,y1,x2,y2) and new (X1,Y1,X2,Y2)
	// It calculates the transformation composed by 
	//  1) a translation
	//  2) a rotation around the origin that transform the original coordinates into the new ones
	// it returns an array containing 4 components :
	// X of translation, Y of translation, cos and sin of the rotation .

	// Calculates M1M2 vectors in original and new coordinates
	OriM1M2=newArray(2);

	OriM1M2[0]=OriPoints[2]-OriPoints[0];
	OriM1M2[1]=OriPoints[3]-OriPoints[1];

	NewM1M2=newArray(2);

	NewM1M2[0]=NewPoints[2]-NewPoints[0];
	NewM1M2[1]=NewPoints[3]-NewPoints[1];

	// Calculates the cosinus of the angle between the original and new vectors : this is the cos of the rotation angle
	ScalT=OriM1M2[0]*NewM1M2[0]+OriM1M2[1]*NewM1M2[1];
	LonT=sqrt(OriM1M2[0]*OriM1M2[0]+OriM1M2[1]*OriM1M2[1])*sqrt(NewM1M2[0]*NewM1M2[0]+NewM1M2[1]*NewM1M2[1]);
	CosT=ScalT/LonT;
	Theta=acos(CosT)*360/(2*PI);
	SinT=sqrt(1-(CosT*CosT));
	Sign=(OriM1M2[0]*NewM1M2[1]-OriM1M2[1]*NewM1M2[0]);	
	if (Sign<0) {
		SinT=-SinT;
		Theta=-Theta;
	}
	
	// Rotates the new position of the first point back by the opposite rotation relative to the origin	
	XTemp=(NewPoints[0])*CosT+(NewPoints[1])*SinT;
	YTemp=-(NewPoints[0])*SinT+(NewPoints[1])*CosT;

	// The coordinates of the back-rotated point relative to the original position of the first point are the coordinates for the translation component.
	XTrans=XTemp-OriPoints[0];
	YTrans=YTemp-OriPoints[1];

	Transfo=newArray(XTrans, YTrans, CosT, SinT, Theta);
	return Transfo;
}