// JNCT0.08b
// Initially, an attempt at conversion of Keith Brain's NCT1.6 for NIH image to an ImageJ macro.
// By Rob Amos, August 2006.
// Contact (as of 2009): robert.amos@tessella.com.
// Requires no plugins:
var particle2="Outlines";
var maxPartSize=9999;
var maxPartString="Infinity";
var partsize=0;
var partSResults=0;
var partDResults=1;
var partCResults=1;
var partRResults=0;
var minCircle=0;
var maxCircle=1;
var FlagPretreat=false;
var lowerboundrat;
var upperboundrat;
var cellthreshold;
var GRadius=5;
var rimWidth=5;
var kernel=0;
var analyseOn=1;
var imagepath=0;
var FirstSlices = newArray(1);		// initialised here.
var currentstack;
var autoload=true;
var imagetitle=0;
var refImageTitle=0;
var datapath=0;
var outputFile=0;
var logvar=0;
var extension;
var startstack;
var stackCount;
var preStacks = newArray(1);
var RobLeica=true;		// Variable indicating whether the source images are stacked and named according to the Leica_SP2_Stacker convention.
var scaleDistance=0;    //global scale variables, so they can be accessed by functions loading and treating multiple stacks.
						// Assumes that stacks have the same scale!  
var scaleKnown=0;
var scaleUnit="um";
var myMeasurements=true;
var makeItQuick=true;
var plotTime = 0;
var storepath;
var backgroundval=0;
var excelWrite=1;
var overlay=1;
var include1=0;
var exclude1=1;
var particleLUT1="R1LUT.lut";
var concat1=0;          //after analysis of series, files will be concatenated if this is set to 1;
var concatArray;
var concatcount;
var plotName="DeltaF/Fo Plot";
var referencevar=0;
var messageLevel = 0;  // Flag setting the quantity of messages to display: 0 for normal messages only, 1 for debug.

macro "Load Stack [l]"
{
   StackLoader(2,FirstSlices);		//2 translates to choose file...
}

macro "Load Next Stack [n3]"
{
   StackLoader(1,FirstSlices);  	// 1 translates to next file in series
}

macro "Load Previous Stack [n1]"
{
   StackLoader(-1,FirstSlices);	// -1 translates to previous file in series
}

macro "Reload Current Stack [n2]"
{
   StackLoader(0,FirstSlices);   // 0 translates to reload the current file
}

macro "Set Particle Parameters [1]"
{
  setParticleParameters();
  setOptions();
}

//PartAnal Set of Sets
macro "Analyze Stacks [2]"
{
  time1=getTime();
  FlagPretreat=false;
  PartAnalSetofSets();
  time2=getTime();
  if(concat1==1)
  {
    concatFunc();
  }
  // time is retrieved in milliseconds.
  print("Time taken = " + (time2 - time1) / 1000);
  resetvariables(1);
}

// Allows automatic pretreatment (using Gaussian mean averaging)
macro "Pretreat and Analyze Stacks [3]"
{
  chooseManyStacks();
  currentstack=startstack;
  FlagPretreat=true;
  preStacks=newArray(stackCount+1);
  silentPretreat=true;
  Dialog.create("Pretreatment options")
  Dialog.addCheckbox("Pretreat silently?",silentPretreat);
  Dialog.addCheckbox("Analyse immediately following treatment?",analyseOn);
  Dialog.addNumber("Radius for Guassian Blur (smoothing):\n(Set radius = 0 for no smoothing.)",GRadius);
  Dialog.addNumber("Width of rim to remove (e.g. same as radius):",rimWidth);
  //Dialog.addString("Kernel to use for pretreatment:",kernel);  // #1 If we wished to allow the selection of custom kernels, this code could be uncommented.
  Dialog.show();
  silentPretreat=Dialog.getCheckbox();
  analyseOn=Dialog.getCheckbox();
  GRadius=Dialog.getNumber();
  rimWidth=Dialog.getNumber();
  //kernel=Dialog.getString();		// See #1 above regarding custom kernels.
  if(datapath==0)
  {
    datapath=getDirectory("Select location for output of data files and stacks");
  }
  setParticleParameters();
  setOptions();
  time1=getTime();
  if(silentPretreat==1)	{setBatchMode(1);}
  StackLoader(0,FirstSlices);
  for(i=0;i<=stackCount;i++)
  {
	if(i>0) { StackLoader(1,FirstSlices); }     //load next in list
	if(GRadius>0) {ConvolveStack(GRadius,kernel);}           //This code could be modified to pass in custom kernels. See #1
	MaskandRim(rimWidth);  // We may exclude a set portion of the image on the outer edges 
	setOutput();
	saveAs("TIFF...",datapath+outputFile+"_treated.tif");
	preStacks[i]=outputFile+"_treated.tif";
	print(outputFile+"_treated.tif saved");
	closeAll();
  }
  setBatchMode(0);
  closeAll();
  if(analyseOn==1)
  {
	PartAnalSetofSets();
	if(concat1==1)
	{
		concatFunc();
	}
  }
  time2=getTime();
  print("Time taken = " + (time2 - time1) / 1000);
  resetvariables(1);
}

macro "Set Background value (using value from ROI) [b]"
{
  setBckValue();
}

macro "Calculate DeltaF/Fo stack [6]"
{
  Dialog.create("DeltaF/Fo settings");
  Dialog.addNumber("Number of control frames:",1);
  Dialog.addNumber("Background:",backgroundval);
  Dialog.addNumber("Cell edge threshold:",cellthreshold);
  Dialog.addCheckbox("Stacks named according to Rob's Leica format? (eg. xxx_t000_ch00.tif",RobLeica);
  Dialog.show();
  controlframes=Dialog.getNumber();
  backgroundval=Dialog.getNumber();
  cellthreshold=Dialog.getNumber();
  RobLeica=Dialog.getCheckbox();
  title=getTitle();
  print(title+" will be converted to DeltaF/Fo image.");
  print("Fo is average of first "+controlframes+" frames.");
  print("Background set to "+backgroundval+".");
  print("Cell edge threshold set to "+cellthreshold+".");
  if(datapath==0)
  {
    datapath=getDirectory("Select location for output of data files and stacks");
  }
  RatioFramesToFo(controlframes);
  setOutput();
  saveAs("TIFF...",datapath+outputFile+"_Fo.tif");
  print(outputFile+"_Fo.tif saved");
}

macro "Reset Macro Variables [7]"
{
  resetTempVars();
  resetvariables(1);
}

// Sets the background value for the image, using the mean grey value
// in the currently-selected ROI.
function setBckValue()
{
  checkROI();
  currentFocus = getImageID();      //gets current window, by ID number
  getStatistics(area,mean);
  m1=round(mean);
  run("Select None");
  backgroundval=getNumber("Background value:",m1);
  backgroundtitle = getTitle();
  print("Background "+m1+" set for "+backgroundtitle);
  logvar=1;
  if(isOpen(currentFocus)==true)
  {
    selectImage(currentFocus);      //restore focus to original image
  }
}

// Allows the user to set parameters for ImageJ's built-in particle analysis algorithm.
function setParticleParameters()
{
  // We initialise the parameters with some defaults if they have not yet been set.
  if (lowerboundrat==0) {lowerboundrat=50;}
  if (upperboundrat==0) {upperboundrat=999;}
  if (partsize==0) {partsize=100;}
  if (backgroundval==0) {backgroundval=1;}
  if (cellthreshold==0) {cellthreshold=10;}
  particleOps=newArray("Nothing","Outlines");

  // Set up a dialog box to gather details from the user.
  Dialog.create("Set Particle Parameters");
  Dialog.addNumber("Threshold for ratio (e.g. 50[/32]):",lowerboundrat);
  Dialog.addNumber("Ratio limit (e.g. 255):",upperboundrat);
  Dialog.addNumber("Min particle size (e.g. 100):",partsize);
  Dialog.addNumber("Max particle size (0 = infinity):",maxPartSize);
  Dialog.addNumber("Threshold for subtracted (e.g. 1)?",backgroundval);
  Dialog.addNumber("Cell Edge Threshold (e.g. 10)?",cellthreshold);
  Dialog.addNumber("Circularity (minimum >= 0.00)",minCircle);
  Dialog.addNumber("Circularity (maximum <= 1.00)",maxCircle);
  Dialog.addCheckbox("Include holes in particle?",include1);
  Dialog.addCheckbox("Exclude particles touching edges?",exclude1);
  Dialog.addMessage("Output Options:");
  Dialog.addChoice("Particle Outline:",particleOps,particle2);
  Dialog.addCheckbox("Overlay particle outline on ratio image?",overlay);
  //Dialog.addCheckbox("Clear Previous Results?",partCResults);
  Dialog.addCheckbox("Summarize results?",partSResults);
  Dialog.addCheckbox("Record XY start coordinates\n(useful for some other plugins)?",partRResults);
  Dialog.addCheckbox("Set standard NCT measurement options?",myMeasurements);
  Dialog.addCheckbox("Batch Mode on (speeds up by not displaying images)",makeItQuick);
  Dialog.show();

  // Read the parameter values from the Dialog.
  lowerboundrat=round(Dialog.getNumber());
  upperboundrat=round(Dialog.getNumber());
  partsize=(Dialog.getNumber());
  maxPartSize=(Dialog.getNumber());
  if(maxPartSize==0)
  {
    maxPartString="Infinity";
  }
  else
  {
    maxPartString=maxPartSize;
  }
  backgroundval=round(Dialog.getNumber());
  cellthreshold=round(Dialog.getNumber());
  minCircle=(Dialog.getNumber());
  maxCircle=(Dialog.getNumber());
  include1=Dialog.getCheckbox();
  exclude1=Dialog.getCheckbox();
  particle2=Dialog.getChoice();
  overlay=Dialog.getCheckbox();
  //partCResults=Dialog.getCheckbox();
  partCResults=1;
  partSResults=Dialog.getCheckbox();
  partRResults=Dialog.getCheckbox();
  myMeasurements=Dialog.getCheckbox();
  makeItQuick=Dialog.getCheckbox();
  if(myMeasurements==true)
  {
    run("Set Measurements...", "area mean standard min centroid circularity slice redirect=None decimal=3");
  }
  print("ratio for threshold = "+lowerboundrat);
  print("ratio limit = "+upperboundrat);
  print("particle size = "+partsize+" - "+maxPartString);
  print("background = "+backgroundval);
  print("cell edge threshold = "+cellthreshold);
  print("circularity = "+minCircle+" - "+maxCircle);
  print("include holes set to "+include1);
  print("exclude on edges set to "+exclude1);
  print("overlay = "+overlay);
}

// Perform particle analysis on a series of stacks.
function PartAnalSetofSets()
{
  if(FlagPretreat==false)
  {
    chooseManyStacks();
    if(datapath==0)
    {
       datapath=getDirectory("Select location for output of data files and stacks");
    }
    setParticleParameters();
    setOptions();
    fileArray=FirstSlices;
    currentstack=startstack;     //currently unnecessary, but might be useful if other functions are introduced.
  }
  else
  {
    fileArray=preStacks;
    storepath=imagepath;  //in case of later functions that would need to find real original image!
    imagepath=datapath;   //ensures that preStacks are found in correct directory.
    currentstack=0;
  }
  if(makeItQuick==true)
  {
    setBatchMode(1);
  }
  StackLoader(0,fileArray);
  for (i=0; i<=stackCount; i++)
  {
    if(i>0){StackLoader(1,fileArray);}
    PartAnalSet(1);
    print(fileArray[(currentstack)]+" analyzed for particles.");
	saveLog();
	//setOutput();
    //saveAs("Measurements", datapath+outputFile+"A.xls");
    //run("Clear Results");
    //saveAs("TIFF...",datapath+outputFile+"A.tif");
    closeAll();
  }
  if(makeItQuick==true)
  {
    setBatchMode(0);
  }
}

// Perform particle analysis on a series of images
function PartAnalSet(arg)
{
  if(makeItQuick==true)
  {
    setBatchMode(1);
  }
  // Create the new stack of ratio images from the current stack.
  RatioAdjacentFrames();            
  Focalstack=getImageID();
  if(overlay!=0 && particle2!="Nothing")
  {
    temptitle=getTitle();
  }
  lowerbound=lowerboundrat / 32;
  upperbound=upperboundrat / 32;
  startResult=0;
  ParticleAnalysisSet(lowerbound,upperbound);
  endResult=nResults;
  particleImage = getImageID();
  particlesFound = endResult-startResult;
  if(particlesFound==0)
  {
	print("No particles detected - check settings?");
	selectWindow(temptitle);
	run("8-bit");
	if(messageLevel > 0)
	{
		setBatchMode(0);
		waitForUser("Focalstack: no particles detected.");
		setBatchMode(1);
	}
  }
  else
  {
	print("Particles found = "+ particlesFound);
	// If we have selected to perform an overlay, we need to find the output window of the particle analysis module
	if(overlay==1&&particle2!="Nothing")
	{
	  expectedTitle = "Drawing of "+temptitle;
	  if(particle2=="Outlines" && isOpen(expectedTitle))
      {
		selectWindow(expectedTitle);
      }
      else if(particle2=="Masks" && isOpen(expectedTitle))
      {
		selectWindow("Mask of "+temptitle);
      }
	  else
	  {
	  // In some cases, the above calculation using nResults is incorrect - in this case, we do not bother with the overlay.
	    particlesFound = 0;
	  }
	  if(particlesFound > 0)
	  {
	  // We rename the particle detection output to give it a consistent name, and store a reference to it as image2
        run("Rename...","Outlines");
		if(messageLevel > 0)
		{
			setBatchMode(0);
			waitForUser("Click OK to continue.");
			setBatchMode(1);
		}
        image2=getImageID();
        selectWindow(temptitle);
        run("8-bit");    //since trying the analysis on the 32-bit image
        run("Min...","stack value=2");
        image1=getImageID();
        time3=getTime();
		// We overlay one image over the other, so that we can see the objects that were actually measured. 
        fastOverlay(image1,image2);
        time4=getTime();
        slowTime=time4-time3;
        print("Time taken for overlay was "+slowTime/1000+" s");
	  }
      selectWindow(temptitle);
      Focalstack=getImageID;
	  if(messageLevel > 0)
	  {	
		setBatchMode(0);
		waitForUser("Focalstack: particles detected.");
		setBatchMode(1);
	  }
	  if(bitDepth() != 8)
	  {
	    run("8-bit");    // In order to paste later, we must ensure that the current image is 8 bit
	  }
    }
  }
  temptitle=imagetitle;
  if(FlagPretreat==true)
  {
    fileArray=preStacks;
  }
  else
  {
    fileArray=FirstSlices;
  }
  // We open either the originals or the pretreated (smoothed) images., and stored them as the "Global stack"
  print("Loading files: " + imagepath + fileArray[(currentstack)]);
  open(imagepath+fileArray[(currentstack)]);
  Globalstack=getImageID();
  print("Setting globalstack as "+ getTitle());
  // We create a new stack, with the ratio images side-by-side with the originals or smoothed originals.
  CombineStacks(Focalstack,Globalstack,1);
  if(bitDepth() != 8)
  {
	 run("8-bit");    // In order to paste later, we must ensure that the current image is 8 bit
  }
  
  // If we have performed the overlay and have found some particles in our images, we set the lookup table.
  if(overlay==1&&particle2!="Nothing"&& particlesFound > 0)
  {
    lutPath=getDirectory("plugins")+"LUT"+File.separator+particleLUT1;
	if(File.exists(lutPath))
	{
      run("LUT... ", "open=["+lutPath+"]");
	}
	else
	{
	  run("Grays");
	}
  }
  // Set the root name for the output files
  setOutput();
  for(i=0;i<endResult-startResult;i++)
  {
    setResult("Original",startResult+i,getResult("Slice",startResult+i)+1);
  }
  if(plotTime!=0)
  {
    for(i=0;i<endResult-startResult;i++)
    {
      setResult("Time",startResult+i,i*plotTime);
    }
  }
  updateResults();
  if(excelWrite==1)
  {
    saveAs("Measurements", datapath+outputFile+"A.xls");
    print("Measurements saved as "+outputFile+"A.xls");
    run("Clear Results");
  }
  saveAs("TIFF...",datapath+outputFile+"A.tif");
  closeAll();
  // We turn batch mode off again if it has been set, so as to display the final image now that we have finished processing.
  if(makeItQuick==true&&arg==0)
  {
    setBatchMode(0);
  }
}

// Performs the particle analysis operation
function ParticleAnalysisSet(lowerbound,upperbound)
{
  checkForStack();
  if(messageLevel > 0)
  {
    print("debug: bit depth of image = " + bitDepth());
  }
  //if(bitDepth() != 8)
  //{
  //  print("debug: converting from " + bitDepth + " to 8 bit before particle analysis.");
  //	run("8-bit");  // convert the image to 8-bit, necessary for particle analysis.
  //}
  if(messageLevel > 0)
  {
    print("debug: applying threshold settings: lower=" + lowerbound + ", upper="+ upperbound);
  }
  setThreshold(lowerbound,upperbound);
  setBatchMode(false);
  if(messageLevel > 0)
  {
    waitForUser("Click OK to continue.");
  }
  setBatchMode(true);
  if(partCResults==1) {clear1="clear ";}
  else                {clear1="";}
  
  if(partDResults==1) {display1="display ";}
  else                {display1="";}
  
  if(partRResults==1) {record1="record ";}
  else                {record1="";}
  
  if(partSResults==1) {summarize1="summarize ";}
  else                {summarize1="";}
  
  if(include1==1)	{include2="include ";}
  else		{include2="";}
  
  if(exclude1==1)	{exclude2="exclude ";}
  else		{exclude2="";}

  if(messageLevel > 0)
  {
    print("debug: Number of open images = " + nImages);
    print("debug: imageID = " + getImageID());
  }
  
  run("Analyze Particles...","size="+partsize+"-"+maxPartString+" circularity="+minCircle+"-"+maxCircle+" show="+particle2+" "+exclude2+display1+clear1+include2+summarize1+record1+"stack");
}

// Applies a Gaussian averaging filter with the specified radius
// to all images in the stack, to smooth out local variations due to
// noise.  A possible extension would be to allow custom kernels for
// different filters.
function ConvolveStack(GRadius,kernel)
{
  checkForStack();
  if(GRadius!=0)
  {
    run("Gaussian Blur...","radius="+GRadius+" stack");
  }
  else
  {
    if(kernel==0)
    {
      run("Convolve...");
    }
    else
    {
      ;//  If custom kernels were supported, we would add the calls here.
    }
  }
}

// Converts the ratio images into images of DeltaF / initialF, using the specified number of control
// frames to calculate the initialF values. 
function RatioFramesToFo(controlframes)
{
  checkForStack();
  j=nSlices;
  mainstack=getImageID();
  selectImage(mainstack);
  run("Subtract...","stack value="+cellthreshold);
  run("Add...","stack value="+cellthreshold);
  run("Subtract...","stack value="+backgroundval);
  run("Z Project...", "start="+0+" stop="+controlframes+" projection=[Average Intensity]");
  Fo=getImageID();
  imageCalculator("divide 32-bit stack",mainstack,Fo);
  run("Rename...","DeltaF/Fo");
}

// Creates the new stack of ratio images: the pixel-by-pixel ratio of successive images to their predecessors.
function RatioAdjacentFrames()
{
  checkForStack();
  j=nSlices;
  mainstack=getImageID();
  setSlice(1);
  run("Subtract...","value="+cellthreshold);
  run("Add...","value="+cellthreshold);
  run("Subtract...","value="+backgroundval);
  run("Duplicate...","title=image1");      //copies slice, and sets title to '1'
  denominator=getImageID();
  for(i=2;i<=j;i++)
  {
    selectImage(mainstack);
    setSlice(i);
    run("Duplicate...","title=image"+i);
    numerator=getImageID();
    run("Subtract...","value="+cellthreshold);
    run("Add...","value="+cellthreshold);
    run("Subtract...","value="+backgroundval);
    imageCalculator("divide create 32-bit","image"+i,"image"+(i-1));
    selectImage(denominator);
    close();
    selectImage(numerator);
    denominator=getImageID();
  }
  selectImage(denominator);
  close();
  selectImage(mainstack);
  close();
  print("Converting ratio images to stack.");
  run("Convert Images to Stack");
  // imageId = getImageID();  // TODO: check. Apparently required in ImageJ 1.41 to reselect new image?
  if(messageLevel > 0)
  {
    print("debug: images open after stack conversion = " + nImages);
	print("debug: stack title = " + getTitle());
  }
  // selectImage(imageId); 
  rename("ratiostack");
  if(messageLevel > 0)
  {
    setBatchMode(false);
    waitForUser("Ratio Stack.  Click OK to continue.");
    setBatchMode(true);
  }
  //run("Multiply...","stack value=32");
  
  //run("Conversions..."," ");
}

// Combines the reference and ratio image stacks into a single stack.
function CombineStacks(stack1,stack2,userShift)  //was MergeFocalGlobal
{
  if(!isOpen(stack1))
  {
	print("Stack 1 is missing.");
    waitForUser("Unable to combine stacks: an invalid stack has been specified.")
	return;
  }
  if(!isOpen(stack2))
  {
	print("Stack 2 is missing.");
    waitForUser("Unable to combine stacks: an invalid stack has been specified.")
	return;
  }
  
  selectImage(stack1);
  if(messageLevel > 0)
  {
	setBatchMode(0);
	waitForUser("Stack1. Click to Continue");
	setBatchMode(1);
  }
  w1=getWidth();
  h1=getHeight();
  d1=nSlices;
  stop1=0;
  selectImage(stack2);
  if(messageLevel > 0)
  {
	setBatchMode(0);
	waitForUser("Stack2.  Click to Continue");
	setBatchMode(1);
  }
  w2=getWidth();
  h2=getHeight();
  d2=nSlices;
  stop2=0;
  mod=d1-d2;
  if(mod<=0)
  {
    d3=d1;
  }
  else
  {
    d3=d2;
  }
  w3=w1+w2;
  if(h1>=h2)
  {
    h3=h1;
  }
  else 
  {
    h3=h2;
  }
  newImage("Merged","8-bit Black", w3, h3, d3);
  newstack=getImageID();
  for(sliceNumber=1; sliceNumber<=d3; sliceNumber++)
  {
    if(userShift>0)
    {
      userShift--;
    }
    else
    {
      selectImage(stack1);
      if(nSlices==1&&stop1<2)
      {
        stop1=1;
      }
      if(stop1<2)
      {
        setSlice(sliceNumber);
        run("Select All");
        run("Copy");
        selectImage(newstack);
        setSlice(sliceNumber);
        makeRectangle(0,0,w1,h1);
        run("Paste");
        if(stop1==1)
        {
          stop1=2;
        }
      }
    }
    selectImage(stack2);
    if(nSlices==1&&stop2<2)
    {
      stop2=1;
    }
      if(stop2<2)
      {
        setSlice(sliceNumber);
        run("Select All");
        run("Copy");
        selectImage(newstack);
        setSlice(sliceNumber);
        makeRectangle(w1,0,w2,h2);
        run("Paste");
        if(stop2==1)
        {
          stop2=2;
        }
      }
   }
   if(isOpen(stack1)==true)
   {
     selectImage(stack1);
     close();
   }
   if(isOpen(stack2)==true)
   {
     selectImage(stack2);
     close();
   }
   // We leave the newly created combined stack as the active image
   selectImage(newstack);
}

// Closes all open image and plot windows
function closeAll()
{
  if(referencevar==1)
  {
    if(messageLevel > 0)
    {
      print("debug: saving reference image.");
	}
    saveRefWindow();
    refImageTitle=getTitle();
    if(isOpen(refImageTitle)==true)
    {
      selectWindow(refImageTitle);
      run("Close");
    }
  }
  for(i=1;i<=nImages;i++)
  {
    if(isOpen(i)==true)
    {
      selectImage(i);
      run("Close");
    }
  }
  if(isOpen(plotName))
  {
    selectWindow("DeltaF/Fo Plot");
    run("Close");
  }
  if(isOpen("profile"))
  {
    selectWindow("profile");
    run("Close");
  }
  if(messageLevel > 0)
  {
    print("debug: windows closed.");
  }
}

// Selects the image, less a specified border at the edge (since we wish to exclude objects that are not completely contained within the image)
function MaskandRim(rimWidth)
{
  setBatchMode(1);
  checkForStack();
  mainimage=getImageID();
  edge=rimWidth-1;	//to take into account the drawn border
  setLineWidth(1);
  setColor(255);
  width=getWidth();
  height=getHeight();
  lengthToPretreat=nSlices;
  for(i=1;i<lengthToPretreat+1;i++)
  {
    selectImage(mainimage);
    setSlice(i);
    makeRectangle(edge,edge,width-2*edge,height-2*edge);
    run("Duplicate...","title="+i);
    run("Select All");
    run("Draw");
  }
  selectImage(mainimage);
  close();
  run("Convert Images to Stack");
}


//Checks for stack.  If one image in stack only, returns an error
function checkForStack()  
{
   if(nSlices<2)
   {
      exit("This procedure requires a stack.");
   }
}

// Saves current windows, and loads up a new file from the supplied list of files
function StackLoader(setting,fileArray)
{
    if(imagepath==0) {imagepath=getDirectory("Location (Dir) of image files?");}
	if(fileArray[0]==0) {fileArray=filestore(imagepath,extension);}
	closeAll();
	if(referencevar==1)
    {
       saveRefWindow();
    }
    if(logvar==1)
    {
       saveLog();
    }
    resetTempVars();                      //Resets imagetitle, refImage but NOT backgroundval
    if(setting<2)				//if using sequential functions
	{
       currentstack=currentstack+setting;
	   if(currentstack<0 || currentstack==fileArray.length)
	   {
		 currentstack=currentstack-setting;
		 exit("No more files in this direction");
	   }
	   imagetitle=fileArray[(currentstack)];
	   open(imagepath+imagetitle);
	   showStatus("file "+currentstack+" of "+fileArray.length);
	}
	if(setting==2)
    {
	   currentstack=chooseStack();
	   imagetitle=fileArray[(currentstack)];
	   open(imagepath+imagetitle);
	   showStatus("file "+currentstack+" of "+fileArray.length);
    }
	print(imagetitle+" loaded.");
	if(scaleDistance!=0&&scaleKnown!=0)
    {
       run("Set Scale...","distance="+scaleDistance+" known="+scaleKnown+" pixel=1 unit="+scaleUnit);
    }
    if(plotTime>0)
    {
       run("Properties...", "interval="+plotTime+" frames="+nSlices + " slices=1");
       print("Aquisition rate = "+1/plotTime);
    }
}

// Gets a list of images and stacks in a specified directory
function filestore(dir,fileExt)
{
   requires("1.32f");
   for (j=0; j<FirstSlices.length; j++)	//clears any previously stored array
   {
     FirstSlices[j]=0;
   }
   currentstack=0;
   tempstore = newArray (1000);      // We assume no more than 1000 stacks in one folder.
   j=0;
   list=getFileList(dir);
     for(i=0; i<list.length; i++)
          {
          if(endsWith(list[i],"tif")||endsWith(list[i],"tiff")||endsWith(list[i],"TIF")||endsWith(list[i],"TIFF"))
            {
               tempstore[j]=list[i];
               j++;
            }
          }
    templength=findend(tempstore);           //gets number of valid entries in tempstore array
    FirstSlices=newArray(templength);       //sets global var FirstSlices to contain this number of terms
    for(i=0; i<templength; i++)
    {
	  FirstSlices[i]=tempstore[i];
	}      //copies valid terms from tempstore
    return FirstSlices;
}

// Saves the reference window, a frame or average showing the original image and particle locations on it.
function saveRefWindow()
{
  if(refImageTitle==0)
  {
    if(referencevar==1)
    {
      referencevar=0;
    }
    exit("No reference image registered.");
  }
  else
  {
	currentFocus = getImageID();
	if(isOpen(refImageTitle))
    {
      selectWindow(refImageTitle);
      run("Tiff...",datapath+outputFile+"Ref.tif");
      print("Reference image saved as "+outputFile+".tif");
      refImageTitle=getTitle();
      logvar=1;
      referencevar = 2;
      if(isOpen(currentFocus)==true)
      {
	    selectImage(currentFocus);      //restore focus to original
      }
    }
	else
    {
    referencevar = 0;
    }
  }
}

// Saves the current contents of the log window, and clears it.
function saveLog()
{
  if(isOpen("Log")==true)
  {
    selectWindow("Log");
    saveAs("Text",datapath+outputFile+"Log.txt");
    logvar = 2;
    run("Close");
    showStatus("Log saved as "+outputFile+"Log.txt");
  }
  else
  {
    showMessage("No log window found.");
  }
}

// Creates a dialogue, allowing a single sequence and returns number of choice in FirstSlices array. 
function chooseStack()
{						
   Dialog.create("Choose sequence or file to load");
   Dialog.addChoice("Select:",FirstSlices,currentstack);		
   Dialog.show();						
   choicetemp=Dialog.getChoice();
   for (i=0; i<FirstSlices.length; i++)
   {
	 if(choicetemp==FirstSlices[i]){modifier=i;}
   }
   currentstack=0;		//resets currentstack before returning the modifier
   return modifier;
}

// Creates a dialogue, allowing sequences to be chosen from a drop-down list
function chooseManyStacks()
{
  if(imagepath==0) {imagepath=getDirectory("Location (Dir) of image files?");}
  if(FlagPretreat==false)
   {
     if(FirstSlices[0]==0) {dummyvar=filestore(imagepath,extension);}
   }
   if(FlagPretreat==true)
   {
     if(FirstSlices[0]==0) {dummyvar=filestore(imagepath,"_treated.tif");}
   }
   Dialog.create("Choose stacks to process");
   Dialog.addChoice("Select first stack:",FirstSlices,startstack);
   Dialog.addChoice("Select last stack:",FirstSlices,startstack+stackCount);
   Dialog.show();						
   choicetemp=Dialog.getChoice();
   choicetemp2=Dialog.getChoice();
   for (i=0; i<FirstSlices.length; i++)
   {
      if(choicetemp==FirstSlices[i]){firstvalue=i;}
      if(choicetemp2==FirstSlices[i]){secondvalue=i;}
   }
   startstack=firstvalue;		//sets currentstack before returning stackCount
   stackCount=secondvalue-firstvalue;
   concatArray=newArray(stackCount);
   if(stackCount<0)
   {
     exit("Macro aborted: First stack must precede last stack on the list!");
   }
}

// Sets the root title to use for output files, based on the current stack name
function setOutput()
{
   if(indexOf(imagetitle,"_")!=-1&&RobLeica==true)
   {
     outputFile = substring(imagetitle,0,indexOf(imagetitle,"_",indexOf(imagetitle,"_")+1));
   }
   else
   {
     if(indexOf(imagetitle,".tif")!=-1)
     {
       outputFile=substring(imagetitle,0,indexOf(imagetitle,".tif"));
     }
     else
     {
       outputFile=imagetitle;
     }
   }
}

function findend(dir)			//finds last meaningful (non-zero) value in array
{
  i=0;
  j=0;
  do
  {
    if(dir[i]==0){j=i;}
    i++;
  } while(j==0);
  return j;
}

function resetTempVars()
{
  refImageTitle = 0;
  extractedvar = 0;
  referencevar = 0;
  logvar = 0;
  specResults = 0;
}

function resetvariables(arg)
{
  currentstack=0;
  stackCount=0;
  imagetitle=0;
  imagepath=0;
  datapath=0;
  FirstSlices=newArray(1);
  concatcount=0;
}

function fastOverlay(image1,image2)
{
  setBatchMode(1);
  selectImage(image1);
  im1=getTitle();
  selectImage(image2);
  im2=getTitle();
  if(particle2=="Outlines")
  {
    if(messageLevel > 0)
    {
      setBatchMode(0);
	  waitForUser("CLick to continue");
	  setBatchMode(1);
    }
	// Outlines are typically a white background (255) with dark particles (0), so we run a min calculation so that the particles will come through as an overlay.
    imageCalculator("min stack",image1,image2);
  }
  if(particle2=="Masks")
  {
    if(messageLevel > 0)
	{
      setBatchMode(0);
	  waitForUser("CLick to continue");
	  setBatchMode(1);
	}
    // Masks are typically a bright particle shapes on a dark background, so we run a max calculation so that the particles will come through.
    imageCalculator("max stack",image1,image2);
    if(isActive(image1)!=1)
    {
      selectImage(image1);
    }
    for(i=0;i<nResults;i++)
    {
	// We label each particle site with a number, so that the image can be correlated easily with the result outputs.
      xparticle=getResult("X",i);
      yparticle=getResult("Y",i);
      zparticle=getResult("Slice",i);
      setSlice(zparticle);
      setColor(1);
      drawString(i+1,xparticle,yparticle);
    }
  }
}

// Preserved for posterity, though currently not used.
function slowOverlay(image1,image2)
{
  setBatchMode(1);
  selectImage(image2);
  width = getWidth();
  height = getHeight();
  for(a=1;a<=nSlices;a++)
  {
    for(b=0;b<height;b++)
    {
      for(i=0;i<width;i++)
      {
        setSlice(a);
        p=getPixel(i,b);
        if(p<2)
        {
          if(p==0){v=255;}
          if(p==1){v=254;}
          selectImage(image1);
          setSlice(a);
          setPixel(i,b,v);
          selectImage(image2);
        }
      }
    }
  }
}

// Provides a splash dialog allowing many parameters to be set.
function setOptions()
{
  Dialog.create("Set File Options:");
  Dialog.addCheckbox("Save results to excel file and clear after each stack measurement",excelWrite);
  Dialog.addCheckbox("Are stacks named in Rob's Leica Stacker format?",RobLeica);
  Dialog.addNumber("Time interval (s) between frames for xyt series (0 for off)",plotTime);
  Dialog.addMessage("NB:  Currently, if a scale is set, particle sizes will be calculated accordingly:\n e.g. min particle size 50 will become 50 um2.");
  Dialog.addNumber("Scale settings: distance",scaleDistance);
  Dialog.addNumber("known:",scaleKnown);
  Dialog.addString("unit:",scaleUnit);
  Dialog.show();
  excelWrite = Dialog.getCheckbox();
  RobLeica = Dialog.getCheckbox();
  plotTime=Dialog.getNumber();
  scaleDistance=Dialog.getNumber();
  scaleKnown=Dialog.getNumber();
  scaleUnit=Dialog.getString();
  if (datapath==0)
  {
    datapath=getDirectory("Select location for output of data files and stacks");
  }
}