macro "CT helical slice-thickness" {
//-----------------------------------------------------------
// Macro "CT_helical_slice_thickness.txt" written by David Platten
// Started on 13th September 2006
// Most recent amendments 29th April 2010
// Version 2.1
//
// 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.
//
// Calculates the FWHM by finding the bead or disc position then measuring the CT#
// of it. The FWHM is taken as the difference between min and max CT# values.
// Note - this may not be a correct assumption for your scan!
// Note - this will only work well if the contrast between the bead or disc CT no.
// and background is high - use a very small bead / disc roi to help with this.
//

requires("1.33n");

Dialog.create("Helical slice thickness");
Dialog.addMessage("This routine expects a single helical scan run of images contained within a directory.");
Dialog.addMessage("There must be no other files in the directory.");
Dialog.addMessage("Click OK to choose the directory containing the images, or cancel to exit.");
Dialog.addNumber("ROI diameter (mm):", 5);
Dialog.show();
roi_diameter = Dialog.getNumber();

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

// Open the first image in the list
path = dir+list[0];
if (!endsWith(path,"/")) open(path);

// Draw a 5 mm diameter circle on centre of the image
getVoxelSize(pixel_width, pixel_height, pixel_depth, pixel_units);
diameter = roi_diameter / pixel_width;
width = getWidth;
height = getHeight;
makeOval(width/2, height/2, diameter, diameter);

// Ask the user to move the ROI over the centre of the phantom
waitForUser("User action required", "Please move the ROI over the centre of the phantom\nand then click on OK");
getBoundingRect(x1, y1, width, height);

// Prevent the images from being displayed
setBatchMode(true);

// Loop through each of the images in the list recording the mean pixel value of the
// ROI for each image
profile = newArray(list.length * 2); // An array to hold the z-pos and value of the profile
for(i=0; i<list.length; i++) {
   path = dir+list[i];
   showProgress(i, list.length);
   if (!endsWith(path,"/")) open(path);

   if (nImages>=1) {
      // ROI for the bead
      makeOval( x1, y1, width, height);
      //setSelectionName("bead");

      // 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.
      tag = getTag("0020,1041");
      profile[i] = tag;

      // Obtain the CT # of the bead and write the result into the array
      getRawStatistics(area, mean, min, max, std);
      bead = mean;
      profile[i + list.length] = bead;
   }
   // Close the image file
   close();
}

// 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", "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.25, 0.25); 
   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;
}
//-----------------------------------------------------------