////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
/////////                                                      ///////// 
/////////                  OPTIMISED SCAN PLUGIN               ///////// 
/////////                                                      /////////
/////////      Created by Hector Dejea, hector.dejea@psi.ch    /////////
/////////                     February, 2018                   /////////
/////////                                                      /////////
////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////

/* Optimised Scan

Opt_Scan.ijm is an ImageJ plugin to find the optimal scanning grid for 
irregular objects. Binary and non binary images are accepted. The option
to have a constant grid on the whole object or independent grids for 
each of the different Y layers is given. The plugin generates a .txt file
that can be used as input for the beamline macro.

Last modified: 14/10/2019 */ 

//Check if there is a file containing the OptScan settings
settings_path = getDirectory("home") + "/.optscan_settings.txt";

if (File.exists(settings_path)){

	//Set default parameters specified in the file
	parameters_file = File.openAsString(settings_path);
	lines=split(parameters_file,"\n");

	LowResimagePixelSize = lines[3];
	ZZcenterLowRes = lines[5];
	XXcenterLowRes = lines[7];
	YposLowRes = lines[9];
	LowResFOVY = lines[11];
	cameraFOVX = lines[13];
	cameraFOVY = lines[15];
	HighResimagePixelSize = lines[17];
	CameraMoveX = lines[19];
	CameraMoveY = lines[21];
	overlapXZ = lines[23];
	overlapY = lines[25];

} else {

	//Otherwise set these default values
	LowResimagePixelSize = 6.5;
	ZZcenterLowRes = 0;
	XXcenterLowRes = 0;
	YposLowRes = 0;
	LowResFOVY = 590;
	cameraFOVX = 2560;
	cameraFOVY = 2160;
	HighResimagePixelSize = 0.65;
	CameraMoveX = 0;
	CameraMoveY = 0;
	overlapXZ = cameraFOVX*0.3;
	overlapY = 50;

}

repeatFlag = 1;
while (repeatFlag == 1){
	print("Welcome to the Scan Optimisation plugin!\n");

	//Ask the user for LowRes and HighRes acquisition parameters
	Dialog.create("Magnification and Camera Settings");
	Dialog.addNumber("LowRes imagePixelSize:", LowResimagePixelSize);
	Dialog.addNumber("LowRes XX center:", XXcenterLowRes);
	Dialog.addNumber("LowRes ZZ center:", ZZcenterLowRes);
	Dialog.addNumber("LowRes Y position (B1):", YposLowRes); //Y of highest volume
	Dialog.addNumber("LowRes FOV Y:", LowResFOVY);
	Dialog.addNumber("HighRes FOV X:", cameraFOVX);
	Dialog.addNumber("HighRes FOV Y:", cameraFOVY);
	Dialog.addNumber("HighRes imagePixelSize:", HighResimagePixelSize);
	Dialog.addNumber("Camera movement X:", CameraMoveX);
	Dialog.addNumber("Camera movement Y:", CameraMoveY);
	Dialog.addNumber("Overlap in X and Z (pixels):", overlapXZ);
	Dialog.addNumber("Overlap in Y (pixels):", overlapY);
	Dialog.addChoice("Projection Method:", newArray("MaxIP" ,"MinIP"), "MaxIP");

	Dialog.show();

	LowResimagePixelSize = Dialog.getNumber();
	XXcenterLowRes = Dialog.getNumber();
	ZZcenterLowRes = Dialog.getNumber();
	YposLowRes = Dialog.getNumber();
	LowResFOVY = Dialog.getNumber();
	cameraFOVX = Dialog.getNumber();
	cameraFOVY = Dialog.getNumber();
	HighResimagePixelSize = Dialog.getNumber();
	CameraMoveX = Dialog.getNumber();
	CameraMoveY = Dialog.getNumber();
	overlapXZ = Dialog.getNumber();
	overlapY = Dialog.getNumber();
	projMethod = Dialog.getChoice();

	// Real magnification
	real_magnification = LowResimagePixelSize/HighResimagePixelSize;

	// Set the starting slice
	waitForUser("Need for action", "Scroll to the desired starting slice and click OK when you are done!\n");
	StartSlice = getSliceNumber();

	//Set the ending slice
	waitForUser("Need for action", "Scroll to the desired ending slice and click OK when you are done!\n");
	EndSlice = getSliceNumber(); 

	//Compute number of acquisition layers in Y and layer size
	blockSizeY = cameraFOVY-overlapY;
	print("Volume Size in Y is " + (cameraFOVY/real_magnification) + " pixels");
	print("Block Size in Y is " + (blockSizeY/real_magnification) + " pixels");
	nLayersY = computeNumberBlocks("Y",blockSizeY,cameraFOVY,StartSlice,EndSlice,real_magnification);

	//Compute blocksize XZ
	blockSizeXZ = cameraFOVX-overlapXZ;
	print("Volume Size in XZ is " + (cameraFOVX/real_magnification) + " pixels");
	print("Block Size in XZ is " + (blockSizeXZ/real_magnification) + " pixels");

	//Check if image is binary
	binaryFlag = checkBinaryImg();

	getDimensions(imagewidth, imageheight, stackchannels, stackslices, stackframes);

	auxInit = StartSlice;
	xVals = newArray(nLayersY);
	yVals = newArray(nLayersY);
	xEndVals = newArray(nLayersY);
	yEndVals = newArray(nLayersY);

	imgList = getList("image.titles"); //Get name of main stack
	imgName = imgList[0];

	if (binaryFlag) {
		//For binary images
		//Generate maximum/minimum intensity projection images
		for (i = 0;  i < nLayersY; i++) {
			auxEnd = auxInit+(cameraFOVY-overlapY)/real_magnification;
			if (auxEnd > EndSlice) {
				auxEnd = EndSlice;
			}
			selectWindow(imgName);
			if (projMethod == "MaxIP") {
				run("Z Project...", "start=" + auxInit + " stop=" + auxEnd + " projection=[Max Intensity]");
			} else {
				run("Z Project...", "start=" + auxInit + " stop=" + auxEnd + " projection=[Min Intensity]");	
			}
			auxInit = auxEnd + 1;

			if (i > 0) { //Show previous rectangle
				makeRectangle(xVals[i-1], yVals[i-1], xEndVals[i-1]-xVals[i-1], yEndVals[i-1]-yVals[i-1]);
			}

			waitForUser("Need for action", "Draw a rectangle around the object and click OK when you are done!\n");
			getSelectionBounds(x, y, width, height);

			//Save definition of all rectangles
			xVals[i] = x;
			yVals[i] = y;
			xEndVals[i] = x + width -1;
			yEndVals[i] = y + height -1;	
		}

		//GENERATE GRID
		centresXXZZY = generateYIndependentGrid(blockSizeXZ,blockSizeY,cameraFOVX,cameraFOVY,nLayersY,xVals,
			yVals,xEndVals,yEndVals,real_magnification,XXcenterLowRes,ZZcenterLowRes,LowResimagePixelSize,
			HighResimagePixelSize,LowResFOVY,CameraMoveX,CameraMoveY,imagewidth,imageheight,StartSlice);

		//COMPLETE FOR BINARY IMAGES -> ELIMINATE EMPTY VOLUMES



	} else {
		//For non binary images
		//Generate maximum/minimum intensity projection images and surround object with bounding box
		for (i = 0;  i < nLayersY; i++) {
			auxEnd = auxInit+(cameraFOVY-overlapY)/real_magnification;
			if (auxEnd > EndSlice) {
				auxEnd = EndSlice;
			}
			selectWindow(imgName);
			if (projMethod == "MaxIP") {
				run("Z Project...", "start=" + auxInit + " stop=" + auxEnd + " projection=[Max Intensity]");
			} else {
				run("Z Project...", "start=" + auxInit + " stop=" + auxEnd + " projection=[Min Intensity]");	
			}
			auxInit = auxEnd + 1;

			if (i > 0) { //Show previous rectangle
				makeRectangle(xVals[i-1], yVals[i-1], xEndVals[i-1]-xVals[i-1], yEndVals[i-1]-yVals[i-1]);
			}

			waitForUser("Need for action", "Draw a rectangle around the object and click OK when you are done!\n");
			getSelectionBounds(x, y, width, height);

			//Save definition of all rectangles
			xVals[i] = x;
			yVals[i] = y;
			xEndVals[i] = x + width -1;
			yEndVals[i] = y + height -1;
		}

		//GENERATE GRID
		centresXXZZY = generateYIndependentGrid(blockSizeXZ,blockSizeY,cameraFOVX,cameraFOVY,nLayersY,xVals,
			yVals,xEndVals,yEndVals,real_magnification,XXcenterLowRes,ZZcenterLowRes,LowResimagePixelSize,
			HighResimagePixelSize,LowResFOVY,CameraMoveX,CameraMoveY,imagewidth,imageheight,StartSlice);

	}

	//Print .txt file with the generated coordinates
	writeCoordsFile(centresXXZZY);
	//Last dialog and restart flag
	repeatFlag = finalDialog(centresXXZZY, nLayersY);

	//Close all windows except stack
	if (repeatFlag) {
		closeGeneratedImgs();
	}

	//Save txt file with the parameters used, so that they are default for next use.
	settings_file = File.open(settings_path);
	settings_array = "File containing the last settings used in Opt_Scan.ijm\n" +
	"\nLowResimagePixelSize\n" + LowResimagePixelSize +
	"\nZZcenterLowRes\n" + ZZcenterLowRes +
	"\nXXcenterLowRes\n" + XXcenterLowRes +
	"\nYposLowRes\n" + YposLowRes +
	"\nLowResFOVY\n" + LowResFOVY +
	"\ncameraFOVX\n" + cameraFOVX +
	"\ncameraFOVY\n" + cameraFOVY +
	"\nHighResimagePixelSize\n" + HighResimagePixelSize +
	"\nCameraMoveX\n" + CameraMoveX +
	"\nCameraMoveY\n" + CameraMoveY +
	"\noverlapXZ\n" + overlapXZ +
	"\noverlapY\n" + overlapY;
	print(settings_file, settings_array);
	File.close(settings_file);
  
} //END


  
  ////////////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////////////
  //////////////////////////////// FUNCTIONS /////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////////////

  function checkBinaryImg() {
	//Function that checks if the current image is binary
	//Use bitdepth() to adapt to 16bit
	getHistogram(values, counts, 256); 
	gvs = 0;
	binary_flag = true;
	for (i=0;i<=255;i++){
		if (counts[i]>0) {
			gvs++;
		}
		if (gvs>2) {
			binary_flag = false;
			break;
		}
	} 
	return binary_flag;
  }

  ////////////////////////////////////////
  
  function computeNumberBlocks(axis,blockSize,cameraFOV,StartPoint,EndPoint,real_magnification) {
  	//Function that computes the number of blocks in the given axis direction.	
  	delta = (EndPoint-StartPoint)*real_magnification;

	//Include the fact that the first Y layer is complete
	if (delta <= cameraFOV) {
		nBlocks = 1; //Number of heights in the new scan
	}else{
		ratio = (delta - cameraFOV)/blockSize;
		nBlocks = floor(ratio) + 2; //Number of heights in the new scan
	}
	print("Number of " + axis + " blocks: " + nBlocks);

	return nBlocks;
  }

  ////////////////////////////////////////

  function computeYCoords(nLayersY,YposLowRes,LowResFOVY,StartSlice,LowResimagePixelSize,CameraMoveY,cameraFOVY,
  	HighResimagePixelSize,blockSizeY) {
    //Function that computes the Y positions in stage coordinates

    //Y position of center of first volume in stage coords
	centresYreal  = newArray(nLayersY);
	centresYreal[0]  = YposLowRes - (LowResFOVY/2 - StartSlice) * LowResimagePixelSize + CameraMoveY  + 
		(cameraFOVY*HighResimagePixelSize)/2;
  	//Y position of each of the layers in stage coords
  	for (l = 1; l < nLayersY; l++) {
  		centresYreal[l] = centresYreal[l-1] + (blockSizeY * HighResimagePixelSize);
  	}
  	return centresYreal;
  }

  ////////////////////////////////////////

  function lowResPx2StageCoords (numBlocksXX,numBlocksZZ,totalBlocks,nLayersY,yVals,xVals,cameraFOVX,
  	real_magnification,imagewidth,imageheight,XXcenterLowRes,ZZcenterLowRes,LowResimagePixelSize,
  	HighResimagePixelSize,CameraMoveX,blockSizeXZ,centresYreal) {
  	//Function that transforms low resolution pixel coordinates to stage coordinates.

  	//Grid in stage coords giving (XX,ZZ,Y) for each volume, from top to bottom, left to right and up to down.
	centresXXZZY = newArray(totalBlocks*3);
  	//Counter for filling up the last array
  	cnt = 0;

  	for (ll = 0; ll<nLayersY; ll++) {

	  	//Position of center of first volume in LowRes pixels
	  	centre1XXpix = yVals[ll] + cameraFOVX/2/real_magnification;
	  	centre1ZZpix = xVals[ll] + cameraFOVX/2/real_magnification;

	  	im_center_x = imagewidth/2;
	  	im_center_y = imageheight/2;

	  	centresXXreal = newArray(numBlocksXX[ll]);
	  	centresZZreal = newArray(numBlocksZZ[ll]);

	  	//Position of center of first volume in stage coords
	  	centresXXreal[0] = XXcenterLowRes - (im_center_y - centre1XXpix) * LowResimagePixelSize + CameraMoveX;
	  	centresZZreal[0] = ZZcenterLowRes + (im_center_x - centre1ZZpix) * LowResimagePixelSize;

	  	//XX position of each of the other volumes in stage coords
	  	for (j = 1; j < numBlocksXX[ll]; j++) {
	  		centresXXreal[j] = centresXXreal[j-1] + (blockSizeXZ * HighResimagePixelSize);
	  	}
	  	//ZZ position of each of the other volumes in stage coords
	  	for (k = 1; k < numBlocksZZ[ll]; k++) {
	  		centresZZreal[k] = centresZZreal[k-1] - (blockSizeXZ * HighResimagePixelSize);
	  	}

	  	//Fill grid in stage coords giving (XX,ZZ,Y) for each volume, from top to bottom, left to right 
	  	//and up to down.
	  	for (jj = 0; jj<numBlocksXX[ll]; jj++) {
	  		for (kk = 0; kk<numBlocksZZ[ll]; kk++) {
	  			centresXXZZY[cnt] = centresXXreal[jj];
	  			centresXXZZY[cnt+1] = centresZZreal[kk];
	  			centresXXZZY[cnt+2] = centresYreal[ll];
	  			cnt = cnt+3;
	  		}
	  	}
  	}
  	return centresXXZZY;
  }

  ////////////////////////////////////////

  function generateYIndependentGrid(blockSizeXZ,blockSizeY,cameraFOVX,cameraFOVY,nLayersY,xVals,yVals,xEndVals,
  	yEndVals,real_magnification,XXcenterLowRes,ZZcenterLowRes,LowResimagePixelSize,HighResimagePixelSize,
  	LowResFOVY,CameraMoveX,CameraMoveY,imagewidth,imageheight,StartSlice) {
  	//Function that generates independent grids for each Y
  	
  	//Compute Y positions in stage coords
    centresYreal = computeYCoords(nLayersY,YposLowRes,LowResFOVY,StartSlice,LowResimagePixelSize,CameraMoveY,
    	cameraFOVY,HighResimagePixelSize,blockSizeY);

  	//Compute number of blocks per layer and total number
  	numBlocksXX = newArray(nLayersY);
  	numBlocksZZ = newArray(nLayersY);
  	totalBlocksLayer = 0;
  	totalBlocks = 0;
  	for (i = 0; i < nLayersY; i++ ) {
  		print("\nY layer" + (i+1) + ":");
  		numBlocksXX[i] = computeNumberBlocks("XX",blockSizeXZ,cameraFOVX,yVals[i],yEndVals[i],real_magnification);
	  	numBlocksZZ[i] = computeNumberBlocks("ZZ",blockSizeXZ,cameraFOVX,xVals[i],xEndVals[i],real_magnification);

	  	//Total number of XX and ZZ blocks
	  	totalBlocksLayer = numBlocksXX[i]*numBlocksZZ[i];
	  	totalBlocks = totalBlocks + totalBlocksLayer;
  	}

	//Grid in stage coords giving (XX,ZZ,Y) for each volume, from top to bottom, left to right and up to down.
	centresXXZZY = lowResPx2StageCoords(numBlocksXX,numBlocksZZ,totalBlocks,nLayersY,yVals,xVals,cameraFOVX,
		real_magnification,imagewidth,imageheight,XXcenterLowRes,ZZcenterLowRes,LowResimagePixelSize,
		HighResimagePixelSize,CameraMoveX,blockSizeXZ,centresYreal);
  	return centresXXZZY;
  }

  ////////////////////////////////////////

  function writeCoordsFile(centresXXZZY) {
  	  //Function that writes the computed coordinates to a .txt file and saves it.
	  file = File.open("");
	  for (i = 0; i <  lengthOf(centresXXZZY); i=i+3) {
	  	print(file, centresXXZZY[i] + "," + centresXXZZY[i+1] + "," + centresXXZZY[i+2]);
	  }
	  File.close(file);
  }

  ////////////////////////////////////////

  function finalDialog(centresXXZZY, nLayersY) {
  	  //Function that shows the last dialog with a summary and the possibility of restart.
	  info = "The motor positions of " + lengthOf(centresXXZZY)/3 + " scans have been computed (" + 
	  	nLayersY + " layer(s)). \n"; 

	  Dialog.create("MY WORK HERE IS DONE")
	  Dialog.addMessage("Calculations finished and saved. Use the txt file as input for your acquisition.  \n" +
	  	" \n\n" + info + " \n\n");
	  Dialog.addCheckbox("Good luck! Would you like to restart?",0);
	  Dialog.show();

	  repeatFlag = Dialog.getCheckbox();
	  return repeatFlag;
  }

  ////////////////////////////////////////

  function closeGeneratedImgs() {
  	//Function that closes all MIP images generated.
  	imgList = getList("image.titles");
	for (i=1; i<lengthOf(imgList); i++) {
		selectWindow(imgList[i]); 
		run("Close");
	}
  }

  ////////////////////////////////////////

  ////////////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////////////