carandraug

// "Macros to count things from nuclei masks"
// Copyright (c) 2011 Carnë Draug <[email protected]>

// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; If not, see <http://www.gnu.org/licenses/>.

macro "Count foci" {

  // Set variables
  extension = ".tif";    // extension of files to analyze (don't forget the dot)

  //////////////////////////////////////////////////////////////////////////////
  // Options menu
  //////////////////////////////////////////////////////////////////////////////
  Dialog.create("Options to count foci");

  Dialog.addChoice("Channel for nuclei", newArray("blue", "red", "green"));
  Dialog.addNumber("Threshold value for nuclei", 30);
  Dialog.addNumber("Minimum particle size", 5000);
  Dialog.addCheckbox("Watershed cell nuclei", true);

  Dialog.addMessage("Aside the channel used to find the nuclei, two other\nare left. Both can be used to count foci if so desired.");
  Dialog.addCheckbox("Count foci on other channel #1", true);
  Dialog.addCheckbox("Count foci on other channel #2", true);
  Dialog.addString("Name for other channel #1", "Channel #1");
  Dialog.addString("Name for other channel #2", "Channel #2");
  Dialog.addChoice("Channel for foci #1", newArray("red", "green", "blue"));
  Dialog.addChoice("Channel for foci #2", newArray("green", "red", "blue"));
  Dialog.addNumber("Threshold value for channel #1", 30);
  Dialog.addNumber("Threshold value for channel #2", 30);
  Dialog.addNumber("Minimum particle size", 2);
  Dialog.addNumber("Minimum particle size", 2);
  Dialog.addCheckbox("Watershed channel #1", true);
  Dialog.addCheckbox("Watershed channel #2", true);

  Dialog.show();

  // get options for nuclei
  nuclei_channel  = get_channel_number (Dialog.getChoice());
  nuclei_thres    = Dialog.getNumber();
  nuclei_min_size = Dialog.getNumber();
  nuclei_water    = Dialog.getCheckbox();

  use_1 = Dialog.getCheckbox();
  use_2 = Dialog.getCheckbox();
  // return error if none of the 2 other channels is selected
  if (!use_1 && !use_2) {
    exit("None of the other two channels was selected to count foci. You must select at least one. Exiting macro now");
  }

  // create array for each option. If there's only 1 channel selected
  // they'll be removed later
  // FUCKING limited macro language forces me to make this gymnastic...
  extra_name      = new_array_from (Dialog.getString(), Dialog.getString());
  extra_channel   = new_array_from (get_channel_number(Dialog.getChoice()), get_channel_number(Dialog.getChoice()));
  extra_thres     = new_array_from (Dialog.getNumber(), Dialog.getNumber());
  extra_min_size  = new_array_from (Dialog.getNumber(), Dialog.getNumber());
  extra_water     = new_array_from (Dialog.getCheckbox(), Dialog.getCheckbox());
  // remove values if one of the channels was not selected
  if (use_1 && !use_2) {
    extra_name      = remove_from_array (extra_name, 1);
    extra_channel   = remove_from_array (extra_channel, 1);
    extra_thres     = remove_from_array (extra_thres, 1);
    extra_min_size  = remove_from_array (extra_min_size, 1);
    extra_water     = remove_from_array (extra_water, 1);
  } else if (!use_1 && use_2) {
    extra_name      = remove_from_array (extra_name, 0);
    extra_channel   = remove_from_array (extra_channel, 0);
    extra_thres     = remove_from_array (extra_thres, 0);
    extra_min_size  = remove_from_array (extra_min_size, 0);
    extra_water     = remove_from_array (extra_water, 0);
  }

  // return error if one of the 2 other channels is the same as nuclei but allow
  // for channel #1 and #2 to be the same to test different threshold values or
  // watershed
  if (nuclei_channel == extra_channel[0]) {
    exit("A selected channel is the same as channel for nuclei");
  // only if both channels are selected does the array extends to index 1 so
  // that must be checked first (short circuit comparison operators, GO!!!)
  // FUCK ImageJ macro language. Short circuit operators don't work
  } else if (use_1 && use_2) {
    if (nuclei_channel == extra_channel[1]) {
      exit("Selected channel #2 is the same as channel for nuclei");
    }
  }

  // Get list of files
  file_dir  = getDirectory("Choose a Directory ");
  file_list = get_cleansed_file_list (file_dir, extension);

  // Go light speed fast in batch mode
//  setBatchMode (true);
  if (lengthOf(file_list) == 1 && file_list[0] == 0) {
    exit("No file was found in the selected directory with extension " + extension);
  }

  // Start analyzing images, one file at a time
  for (file_i = 0; file_i < lengthOf(file_list); file_i++) {
    filename  = file_list[file_i];
    // be nice and tell people where we are
    showStatus("Analyzing image " + (file_i+1) + " of " + lengthOf(file_list));
    showProgress((file_i+1)/lengthOf(file_list));

    print ("Taking care of file "+filename);
    // if file no longer exists, skip and warn user
    if (!File.exists(filename)) {
      print ("    File " + filename + " seems to no longer exist. Skipped.");
    } else {
      prepare_environment();
      image_ID      = prepare_image(filename);
      nuclei_count  = nuclei_processing (nuclei_channel, nuclei_thres, nuclei_min_size, nuclei_water);
      roiManager("save", filename+"nuclei_mask.zip");
      print ("    Found "+nuclei_count+" nuclei");

      selectImage(image_ID);
      foci_1_mask_ID  = foci_processing (extra_channel[0], extra_thres[0], extra_water[0]);
      if (lengthOf(extra_name) == 2) {
        selectImage(image_ID);
        foci_2_mask_ID  = foci_processing (extra_channel[1], extra_thres[1], extra_water[1]);
      }

      for (nuclei_i = 0; nuclei_i < nuclei_count; nuclei_i++) {
        selectImage(foci_1_mask_ID);
        foci_count    = foci_per_nuclei (nuclei_i, nuclei_count, extra_min_size[0]);
        roiManager("save", filename+"_foci_of_nuclei_"+nuclei_i+"_channel"+extra_name[0]+".zip");
        print ("    Found "+foci_count+" foci in nuclei "+nuclei_i+" in channel "+extra_name[0]);
        run("Clear Results");

        for (foci_i = 0; foci_i < foci_count; foci_i++) {
          selectImage(image_ID);
          roi_index = nuclei_count + foci_i;
          measure_foci (roi_index, foci_i, get_channel_name(extra_channel[0]));
        }
        saveAs("results", filename+"_foci_"+extra_name[0]+"_nuclei_"+nuclei_i+".xls");
        run("Clear Results");
        remove_foci_rois (nuclei_count);

        if (lengthOf(extra_name) == 2) {
          selectImage(foci_2_mask_ID);
          foci_count    = foci_per_nuclei (nuclei_i, nuclei_count, extra_min_size[1]);
          roiManager("save", filename+"_foci_of_nuclei_"+nuclei_i+"_channel"+extra_name[1]+".zip");
          print ("    Found "+foci_count+" foci in nuclei "+nuclei_i+" in channel "+extra_name[1]);
          run("Clear Results");

          for (foci_i = 0; foci_i < foci_count; foci_i++) {
            selectImage(image_ID);
            roi_index = nuclei_count + foci_i;
            measure_foci (roi_index, foci_i, get_channel_name(extra_channel[1]));
          }
          saveAs("results", filename+"_foci_"+extra_name[1]+"_nuclei_"+nuclei_i+".xls");
          run("Clear Results");
          remove_foci_rois (nuclei_count);
        }

      }
      close_by_ID(foci_1_mask_ID);
      if (lengthOf(extra_name) == 2) {
        close_by_ID(foci_2_mask_ID);
      }

    }
  }
  showStatus("All done captain!");
}

////////////////////////////////////////////////////////////////////////////////
// Functions from this point on
////////////////////////////////////////////////////////////////////////////////

// takes the name of a channel (blue, red or green) and converts to the channel
// number on composite images
function get_channel_number (name) {
  if (name == "red"){
    number = 1;
  } else if (name == "green") {
    number = 2;
  } else if (name == "blue") {
    number = 3;
  } else {
    exit ("Some error occurred when trying to guess the number of the channel");
  }
  return number;
}

function get_channel_name (number) {
  if (number == 1){
    name = "red";
  } else if (number == 2) {
    name = "green";
  } else if (number == 3) {
    name = "blue";
  } else {
    exit ("Some error occurred when trying to guess the name of the channel");
  }
  return name;
}

// because the newArray function can't accept the return value of fucntions such
// as Dialog.getString as values for the array, this new function should be able
// to do it
function new_array_from (value1, value2) {
  array     = newArray(2);
  array[0]  = value1;
  array[1]  = value2;
  return array;
}

// this function prepares the environment for the analysis. Closes other windows
// inclusive the ROI manager and results windows. It also removes all other ROI
function prepare_environment() {
  if (isOpen("ROI Manager")) {
     selectWindow("ROI Manager");
     run("Close");
  }
  roiManager("reset");
}

// Prepares the image for everything and returns the image_ID. Opens it, makes it
// composite and may make any adjustement necessary
function prepare_image (filepath) {
  open(filepath);
  run("Make Composite");
  Stack.setDisplayMode("color");
  // composite creates a new image, the ID is different
  image_ID = getImageID();
  return image_ID;
}

// Makes all the morphological processing of the image to get boundaries of the
// cell nuclei.
function nuclei_processing (channel, thres, min_size, water) {
  // morphological processing of the mask must end in a dilate or the borders will
  // be black and "exclude on edges" won't work
  Stack.setChannel(channel);
  run("Gaussian Blur...", "sigma=5 slice");
  setThreshold(thres, 255);
  run("Create Mask");
  ID = getImageID();
  run("Erode");
  run("Dilate");
  run("Fill Holes");
  run("Erode");
  run("Erode");
  run("Dilate");
  run("Dilate");
  run("Fill Holes");
  if (water) {
    run("Watershed");
  }
  // Run analyze particles. Consider particles larger than min-size until infinity
  // and circularity between 0 and 1 (all). It doesn't show anything, clears
  // the results window, excludes on edges, includes holes, and add selections
  // to ROI manager
  roiManager("reset")
  run("Analyze Particles...", "size="+min_size+"-Infinity circularity=0.00-1.00 show=Nothing exclude clear include add");
  count = roiManager("count");
  close_by_ID(ID);
  return count;
}

function foci_processing (channel, thres, water) {
  Stack.setChannel(channel);
  setThreshold(thres, 255);
  run("Create Mask");
  ID = getImageID();
  if (water) {
    run("Watershed");
  }
  return ID;
}

function foci_per_nuclei (nuclei_i, nuclei_count, min_size) {
  roiManager("Select", nuclei_i);
  // do NOT clear the roi manager now
  run("Analyze Particles...", "size="+min_size+"-Infinity circularity=0.00-1.00 show=Nothing exclude add");
  count = roiManager("count");
  count = count - nuclei_count;
  return count;
}

function measure_foci (roi_index, foci_i, channel_name) {
  roiManager('select', roi_index);
  getStatistics(area, mean, min, max);
  sum = get_intensity_sum();

  setResult("Area", foci_i, area);
  setResult("Mean", foci_i, mean);
  setResult("Min", foci_i, min);
  setResult("Max", foci_i, max);
  setResult("Sum", foci_i, sum);
}

function remove_foci_rois (nuclei_count) {
  count = roiManager("count");
  for (i=(count-1); i>=nuclei_count; i--) {
    roiManager("Select", i);
    roiManager("Delete");
  }
}

// Calculates the sum of the intensities of each pixel in the selection
function get_intensity_sum() {
  nBins = 256;
  total = 0;
  getHistogram(values, counts, nBins);
  for (i=values[0]; i<values[nBins-1]; i++) {
    total = total + (values[i]*counts[i]);
  }
  return total;
}