macro "CT helical slice-thickness" {
//-----------------------------------------------------------
// Macro "CT_helical_slice_thickness.txt" written by David Platten
// Started on 13th September 2006
// Version 2.1 29th April 2010
// Version 2.2 23rd April 2013
// Version 2.3 24th April 2013
//
// As of version 2.3 the routine works in the same way as the ImPACT
// z-sensitivity IDL routine.
//
// Expects a directory of DICOM images from a helical scan run through a
// thin metal disc or a very small metal bead. The reconstruction interval
// should be around 1/10 of the slice thickness. For example, if the reconstructed
// slice thickness is 3 mm then you need to set the reconstruction interval to
// 0.3 mm.
//
// The phantom can be a thin gold disc embedded in a Perspex rod, similar to the
// one that ImPACT uses. Alternatively the method may work by scanning through a
// small metal bead, provided that the bead is much smaller than the reconstructed
// slice thickness. However, the routine may fail with this type of test object as
// it tries to locate the centre of the test object using a centre-of-mass
// calculation; this may result in the incorrect positioning of the ROI if the
// test object occupies the whole image.
//
// For each image in the series the centre of mass of pixels that have a value of
// >= max / 2 is calculated. A circular ROI is then centred at this position and
// the mean pixel value calculated. The z-position of the image is recorded. Once
// all images have been examined, the background is subtracted from every point
// (actually the minimum mean ROI value), and then the values are normaised to the
// maximum value. The full-width at half-maximum of the profile is then calculated.
//
// There can be ".txt", ".csv" or "dcminf" files in the directory with your images,
// but nothing else. Any other file in the directory will cause the routine to fail.
//

requires("1.33n");

// Switch display updating off
setBatchMode(true);

Dialog.create("Helical slice thickness");
Dialog.addMessage("This routine is intended to analyse a helical CT scan of a Perspex rod that contains a thin gold disc.\nThe location of the rod is found, and then a 6 mm region of interest placed on the centre of the rod.\nThe mean pixel value of this ROI is measured, and stored against it's z-axis location. This is repeated\nfor every image in the scan. The full-width at half-maximum of these mean values is then calculated,\nand represents the image slice thickness of the scan. The method relies on good z-axis sampling\nof the gold disc, so it is suggested that the image reconstruction increment is set to be 1/10 of the \nnominal slice thickness.");
Dialog.addMessage("The routine requires the images from a single helical scan to be placed in a directory. There can be\nno other images in the directory, but files ending in .txt, .csv or dcminf are allowed.");
Dialog.addMessage("Click OK to choose the directory containing the images, or cancel to exit.");
Dialog.addNumber("       ROI diameter (mm):", 6);
Dialog.show();
ROIdiameter = Dialog.getNumber();

dir = getDirectory("Choose directory with helical slice thickness images");
list = getFileList(dir);

profile = newArray(list.length * 2); // An array to hold the z-pos and value of the profile

imageNumber = 0; // This is used to count the number of images that are loaded

minProfileValue =  99999999; // This is used to store the minimum value of the ROI mean
maxProfileValue = -99999999; // This is used to store the maxumum value of the ROI mean

for(i=0; i<list.length; i++) {
   path = dir+list[i];
   showProgress(i, list.length);

   // If the current file ends in ".txt", ".csv", "dcminf" or "/" then it's ignored.
   if (!endsWith(path,"/") && !endsWith(path,".txt") && !endsWith(path,".csv") && !endsWith(path,"dcminf")) open(path);

   if (nImages>=1) {
      run("Select All");
      getStatistics(area, mean, min, max, std);

      // Determine the centre of mass of all pixels with a value greater than or equal to max / 2.
      width  = getWidth();
      height = getHeight();
      totalMass = 0;
      totalYweight = 0;
      totalXweight = 0;
      threshold = max / 2;
      for (x=0; x<width; x++) {
         for (y=0; y<height; y++) {
            currentValue = calibrate(getPixel(x, y));
	    if (currentValue >= threshold) {
	       totalMass = totalMass + currentValue;
	       totalYweight = totalYweight + (currentValue * y);
	       totalXweight = totalXweight + (currentValue * x);
	    }
	 }
      }
      xCom = totalXweight / totalMass;
      yCom = totalYweight / totalMass;

      // Draw a circular ROI, centred on the centre of mass and obtain the statistics for it
      getVoxelSize(pixel_width, pixel_height, pixel_depth, pixel_units);
      ROIpixelDiameter = ROIdiameter / pixel_width;
      makeOval(xCom-ROIpixelDiameter/2, yCom-ROIpixelDiameter/2, ROIpixelDiameter, ROIpixelDiameter);
      getRawStatistics(area, mean, min, max, std);

      // Set the current profile value to be the mean of the ROI
      profile[imageNumber + list.length] = mean;

      // Check to see if the mean is lower than minProfileValue
      if (mean < minProfileValue) minProfileValue = mean;

      // Check to see if the mean is higher than maxProfileValue
      if (mean > maxProfileValue) maxProfileValue = mean;

      // Find the slice location from the DICOM tag and write this into the results.
      // "0020, 1041" appears to be the tag ID for the "Slice location" tag.
      profile[imageNumber] = parseFloat(getTag("0020,1041"));

      // Increment the image number
      imageNumber++;

      // Close the image file
      close();
   }
}

// See if imageNumber is less than list.length, and if it is then crop the
// profile to remove the empty values. This will occur if there are any .txt,
// .csv or dcminf files in the directory with the images that the for loop
// above skips. Also subtract the background value from each.
if (imageNumber < list.length) {
   croppedProfile = newArray(imageNumber * 2);
   for(i=0; i<imageNumber; i++) {
      croppedProfile[i] = profile[i];

      // Also subtract the background from profile (use the minimum mean ROi value for this) and then
      // normalise to the maximum value, maxProfileValue - minProfileValue.
      croppedProfile[i + imageNumber] = (profile[i + list.length] - minProfileValue) / (maxProfileValue - minProfileValue);
   }
   profile = croppedProfile;
}

// Calculate the FWHM from the profile array
fwhm(profile);

// Select the log file, which contains the fwhm result, into the same directory
// as the images.
selectWindow("Log");
results_name = dir + "helical_slice_thickness_log" + list[0] + ".txt";
saveAs("Text", results_name);

// Switch display updating back on
setBatchMode(false);

// End of the main program
//-----------------------------------------------------------
}


//-----------------------------------------------------------
// Function to sort a 2D array of coordinates into ascending order of x,
// where the array is array[x, y]
function sort2d(unsorted) {
   b = newArray(unsorted.length/2);
   for (i=0; i<b.length; i++) {
      b[i] = (unsorted[i]) * 1.0; // The * 1.0 forces ImageJ to recognise the variable as a number rather than text. The () are also important
   }
   b = sort(b);

   sorted = newArray(unsorted.length);
   for (i=0; i<b.length; i++) {
      for (j=0; j<b.length; j++) {
         if (b[i] == unsorted[j]) {
            sorted[i] = unsorted[j];
            sorted[i+sorted.length/2] = unsorted[j+unsorted.length/2];
         }
      }
   }
   return sorted;
}
//-----------------------------------------------------------


//-----------------------------------------------------------
// The sort and quickSort functions below are from the sort demo at
// http://rsb.info.nih.gov/ij/macros/
function sort(a) {
   a = quickSort(a, 0, lengthOf(a)-1);
   return a;
}

function quickSort(a, from, to) {
   i = from; j = to;
   center = a[(from+to)/2];
   do {
       while (i<to && center>a[i]) i++;
       while (j>from && center<a[j]) j--;
       if (i<j) {temp=a[i]; a[i]=a[j]; a[j]=temp;}
       if (i<=j) {i++; j--;}
   } while(i<=j);
   if (from<j) quickSort(a, from, j);
   if (i<to) quickSort(a, i, to);
   return a;
}
//-----------------------------------------------------------


//-----------------------------------------------------------
// This function calculates the full width at half maximum of
// a set of (x, y) coordinates.
function fwhm(coordinates) {
   setBatchMode(false);

   // Create a new array containing the data and sort it into ascending x-order
   sorted_coordinates = newArray(coordinates.length);
   sorted_coordinates = sort2d(coordinates);

   // find the y max in the data
   max = sorted_coordinates[sorted_coordinates.length/2];
   y_start = sorted_coordinates.length/2;
   y_end = sorted_coordinates.length-1;
   for (i=y_start; i<y_end; i++) {
      current_y = sorted_coordinates[i];
      if (current_y > max) {
         max = current_y;
      }
   }
   print ("Maximum is " + max);

   // find the y min in the data
   min = sorted_coordinates[sorted_coordinates.length/2];
   y_start = sorted_coordinates.length/2;
   y_end = sorted_coordinates.length-1;
   for (i=y_start; i<y_end; i++) {
      current_y = sorted_coordinates[i];
      if (current_y < min) {
         min = current_y;
      }
   }
   print ("Minimum is " + min);

   // Work out y value that corresponds to half way
   fwhm_level = min + ((max - min) / 2);
   print ("FWHM level is " + fwhm_level);

   // find the points to the left and right of the peak where half-max is exceeded
   // first the left-hand side - start with the last point and work backwards
   found = 0;
   i = y_start;
   while(found==0 && i<y_end) {
      current_y = sorted_coordinates[i];
      if (current_y > fwhm_level) {
         left_half_max_index = i-sorted_coordinates.length/2;
         found = 1;
      }
      i++;
   }

   // now the right-hand side
   // start with the first point and work forwards
   found = 0;
   i = y_end;
   while(found==0 && i>y_start) {
      current_y = sorted_coordinates[i];
      if (current_y > fwhm_level) {
         right_half_max_index = i-sorted_coordinates.length/2;
         found = 1;
      }
      i--;
   }

   // Calc the exact x-coord of the left-hand side of fwhm using y = mx + c
   // The * 1.0 is to force ImageJ to treat the object as a number rather than text
   x1 = (sorted_coordinates[left_half_max_index-1]) * 1.0;
   y1 = (sorted_coordinates[left_half_max_index-1 + sorted_coordinates.length/2]) * 1.0;
   x2 = (sorted_coordinates[left_half_max_index]) * 1.0;
   y2 = (sorted_coordinates[left_half_max_index + sorted_coordinates.length/2]) * 1.0;
   m = (y1-y2) / (x1-x2);
   c = y1 - (m*x1);
   x_left_fwhm = (fwhm_level - c) / m;
   print ("Left-hand side of fwhm is at " + x_left_fwhm);

   // and now for the right-hand side
   x1 = (sorted_coordinates[right_half_max_index]) * 1.0;
   y1 = (sorted_coordinates[right_half_max_index + sorted_coordinates.length/2]) * 1.0;
   x2 = (sorted_coordinates[right_half_max_index-1]) * 1.0;
   y2 = (sorted_coordinates[right_half_max_index-1 + sorted_coordinates.length/2]) * 1.0;
   m = (y1-y2) / (x1-x2);
   c = y1 - (m*x1);
   x_right_fwhm = (fwhm_level - c) / m;
   print ("Right-hand side of fwhm is at " + x_right_fwhm);

   // can now calc fwhm
   fwhm_value = x_right_fwhm - x_left_fwhm ;
   print ("FWHM is " + fwhm_value);

   // now plot the profile and include the fwhm on it
   num_points = coordinates.length / 2;
   x_coords = newArray(num_points);
   y_coords = newArray(num_points);
   for (i=0; i<num_points; i++) {
      x_coords[i] = sorted_coordinates[i];
      y_coords[i] = sorted_coordinates[i + num_points];
      print(x_coords[i], y_coords[i]);
   }

   // plot the data
   Plot.create("Profile", "z-position (mm)", "Normalised (CT# - bgd)", x_coords, y_coords);
   fwhm_xpoints = newArray(2);
   fwhm_ypoints = newArray(2);
   fwhm_xpoints[0] = x_left_fwhm;
   fwhm_xpoints[1] = x_right_fwhm;
   fwhm_ypoints[0] = fwhm_level;
   fwhm_ypoints[1] = fwhm_level;
   Plot.setColor("red");
   Plot.add("line", fwhm_xpoints, fwhm_ypoints);
   Plot.addText("FWHM is " + fwhm_value + " mm", 0.15, 0.15); 
   Plot.show();
}
//-----------------------------------------------------------


//-----------------------------------------------------------
// This function returns the value of the specified 
// tag  (e.g., "0010,0010") as a string. Returns "" 
// if the tag is not found.
function getTag(tag) {
   info = getImageInfo();
   index1 = indexOf(info, tag);
   if (index1==-1) return "";
   index1 = indexOf(info, ":", index1);
   if (index1==-1) return "";
   index2 = indexOf(info, "\n", index1);
   value = substring(info, index1+1, index2);
   return value;
}
//-----------------------------------------------------------


//-----------------------------------------------------------
function centreOfMass(xStart, yStart, width, height) {

	totalMass = 0;
	totalYweight = 0;
	totalXweight = 0;

	for (x=xStart; x<(xStart+width); x++) {
		showProgress(x, xStart+width);

		for (y=yStart; y<(yStart+height); y++) {

			pixelValue = getPixel(x, y);
			totalMass = totalMass + pixelValue;
			totalYweight = totalYweight + (pixelValue * y);
			totalXweight = totalXweight + (pixelValue * x);
		}
	}

	xCom = totalXweight / totalMass;
	yCom = totalYweight / totalMass;

//	print("Total mass is:", totalMass);
//	print("Total y weight is:", totalYweight);
//	print("Y centre of mass is:", totalYweight / totalMass);
//	print("X centre of mass is:", totalXweight / totalMass);

//	setResult("X com", nResults, totalXweight / totalMass);
//	setResult("Y com", nResults-1, totalYweight / totalMass);
//	updateResults();

	returnValue = newArray(2);
	returnValue[0] = xCom;
	returnValue[1] = yCom;
	return returnValue;
}
//-----------------------------------------------------------