// Register_Coverslip_Zeiss 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 Axiovision 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 First commented version

macro Register_Coverlip_Zeiss {

	//Gets two Zeiss 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.getName(RefPath);

	// Gets the file name, and positions set name
	OriFile=File.openAsString(OriPath);
	OriLine=split(OriFile, "\n");
	OriLine3=split(OriLine[2], "\"");
	OriName=OriLine3[0];

	print("\n");
	print("*** Register Coverlip Zeiss Log ***");
	print("Original Positions Path: "+OriPath);
	print("Original Positions Name: "+OriName);

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

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

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

		print("Position #"+(i-6)+": "+OriPosName[i-7]+" at ("+OriX[i-7]+" , "+OriY[i-7]+" , "+OriZ[i-7]+")");
	}

	// 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");
	RefLine3=split(RefLine[2], "\"");
	RefName=RefLine3[0];
	print("\n");
	print("Reference Positions Path: "+RefPath);
	print("Reference Positions Name: "+RefName);

	RefPosName=newArray(2);
	RefX=newArray(2);
	RefY=newArray(2);
	RefZ=newArray(2);
	for (i=7; i<9; i++) {
		SplitSemi=split(RefLine[i], ";");
		SplitApo=split(SplitSemi[0], "\"");
		RefPosName[i-7]=SplitApo[0];
		RefX[i-7]=parseFloat(SplitSemi[1]);
		RefY[i-7]=parseFloat(SplitSemi[2]);
		RefZ[i-7]=parseFloat(SplitSemi[3]);

		print("Reference #"+(i-6)+": "+RefPosName[i-7]+" at ("+RefX[i-7]+" , "+RefY[i-7]+" , "+RefZ[i-7]+")");

	}

	// Initialize de 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);
	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
	OffsetZ1=RefZ[0]-OriZ[0];
	OffsetZ2=RefZ[1]-OriZ[1];
	OffsetZ=(OffsetZ1+OffsetZ2)/2;
	
	print("\n");
	print("Transformation: TransX="+d2s(AffineTrans[0],1)+", TransY="+d2s(AffineTrans[1],1)+", Angle="+d2s(Angle,4)+"°, Z-Offset="+OffsetZ);

	// 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
		NewZ[i]=OriZ[i]+OffsetZ;
		
		 // Change the name of the position
		NewPosName[i]=OriPosName[i]+"_new";
		print("New Position #"+(i+1)+": "+NewPosName[i]+" at ("+NewX[i]+" , "+NewY[i]+" , "+NewZ[i]+")");
	}



// Modifies the original positions file to generate the new positions file
	// Loops through positions to change the name, X and Y values
	for (i=7; i<OriLine.length-1; i++) {
		OriLine[i]=replace(OriLine[i],"\";"+OriX[i-7]+";"+OriY[i-7]+";"+OriZ[i-7],"_new\";"+d2s(NewX[i-7],1)+";"+d2s(NewY[i-7],1)+";"+d2s(NewZ[i-7],3));
		OriLine[i]=replace(OriLine[i],"\"Blue\";\"blue\"","\"Red\";\"red\"");
	}

	// Fills the new file	
	NewFile=RefLine[0];
	for (i=1; i<9; i++) {
		NewFile=NewFile+"\n"+RefLine[i];
	}
	for (i=7; i<OriLine.length; i++) {
		NewFile=NewFile+"\n"+OriLine[i];
	}
	NewFile=NewFile+"\n";
	// Saves the new file by relpacing the ref file with added new positions
	File.saveString(NewFile, RefPath);

	
	

// 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-1000, MaxX+1000,MinY-1000, MaxY+1000);

	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
	Plot.setColor("black");
	Plot.addText("Z offset ="+OffsetZ, 0, 0);
	Plot.show;

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

	showStatus("Register_Coverlsip_Zeiss finished");
	print("*** Register Coverlip Zeiss 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;
}