Initial implementations of the Rhoana's classifier training and GALA's 2D segmentation scripts. Added or modified a number of utility scripts.

Author: denisovg

AgglomarativeSegmentation2D/MS_AS_GalaSegmentation2D.py

@@ -0,0 +1,187 @@
+#! /usr/local/python-2.7.6/bin/python
+#
+# Copyright (C) 2015 by Howard Hughes Medical Institute.
+#
+# Purpose: perform 2D segmentation using GALA
+#          then post-process results
+#
+# ------------------------- imports -------------------------
+
+#from __future__ import absolute_import
+#from __future__ import print_function
+# imports
+import os, sys, re, optparse
+from gala import imio, classify, features, agglo, evaluate as ev
+from six.moves import map
+import numpy as np
+
+import MS_LIB_Dict
+import MS_LIB_Options
+import MS_LIB_IO
+
+from PIL import Image
+
+ms_home = os.environ['MS_HOME']   # where source code is located
+ms_data = os.environ['MS_DATA']   # dir for initial input data and final output
+ms_temp = os.environ['MS_TEMP']   # dir for intermediate data files
+gala_home    = os.environ['GALA_HOME']    # where source code is located
+ilastik_home = os.environ['ILASTIK_HOME'] # where source code is located
+
+# -----------------------------------------------------------------------------
+
+def perform_segmentation(training_raw, groundtruth_annot_dir, \
+                         groundtruth_seg_dir, production_raw, \
+                         options):
+    raw_files = os.listdir(os.path.join(ms_data, training_raw))
+    labels    = os.listdir(os.path.join(ms_data, groundtruth_annot_dir))
+
+    # Assume only one raw and one labels file for now
+    command = "MS_UT_CreateH5LabelsForIlastik.py " + groundtruth_annot_dir + \
+              " groundtruth.h5"
+    print "\nRunning the commmand: " , command
+    os.system(command)
+
+    # Generate Ilastik project file
+    command = "python " + os.path.join(ilastik_home, "ilastik", "bin", "train_headless.py") + \
+              " training.ilp '" + training_raw +  "/*.png' groundtruth.h5/main"
+    print "\nRunning the commmand: " , command
+    os.system(command)
+
+    # Produce the probabilities and oversegmentation files from the inputs
+    if not os.path.isdir('seg_data'):
+        print "\nRunning the commmand: mkdir seg_data"
+        os.mkdir('seg_data')
+    else:
+        print "\nRunning the commmand: rm -f seg_data/*"
+        os.system("rm -rf seg_data/*")
+    if os.path.isfile('STACKED_prediction.h5'):
+        command = "rm -f STACKED_prediction.h5"
+        print "\nRunning the commmand: " , command
+        os.system(command)
+    if not os.path.islink('gala-segmentation-pipeline'):
+        os.symlink('gala-segmentation-pipeline', os.path.join(gala_home, 'gala',\
+                   'bin', 'gala-segmentation-pipeline'))
+    command = "python gala-segmentation-pipeline " + \
+              " . -I '" + training_raw + "/*.png' --ilp-file training.ilp " + \
+              " --enable-gen-supervoxels --enable-gen-agglomeration " + \
+              " --enable-gen-pixel --seed-size 5   --segmentation-thresholds 0.0"
+    print "\nRunning the commmand: " , command
+    os.system(command)
+
+    # Put probabilities in the right format:
+    command = "MS_UT_ProbsNF2GALA.py STACKED_prediction.h5 probabilities_training.h5"
+    print "\nRunning the commmand: ", command
+    os.system(command)
+
+    # Rename the oversegmentation training file
+    for file in os.listdir("seg_data"):
+        if re.search(".h5", file):
+            oversegmentation_file = os.path.join("seg_data", file)
+    command = "mv " + oversegmentation_file + " oversegmentation_training.h5"
+    print "\nRunning the commmand: ", command
+    os.system(command)
+
+    # Create the groundtruuth segmentation labels file
+    command = "MS_UT_CreateH5Groundtruth.py " + groundtruth_seg_dir + " seg_labels_file.h5"
+    print "\nRunning the commmand: " , command
+    os.system(command)
+
+    # Read in training data
+    print "\nReading training data into gala..."
+    gt_train, pr_train, ws_train = (map(imio.read_h5_stack,
+                                    ['seg_labels_file.h5', 'probabilities_training.h5',
+                                     'oversegmentation_training.h5']))
+
+    # create a feature manager
+    fm = features.moments.Manager()
+    fh = features.histogram.Manager()
+    fc = features.base.Composite(children=[fm, fh])
+
+    # create graph and obtain a training dataset
+    g_train = agglo.Rag(ws_train, pr_train, feature_manager=fc)
+    (X, y, w, merges) = g_train.learn_agglomerate(gt_train, fc)[0]
+    y = y[:, 0] # gala has 3 truth labeling schemes, pick the first one
+    print((X.shape, y.shape)) # standard scikit-learn input format
+
+    # train a classifier, scikit-learn syntax
+    rf = classify.DefaultRandomForest().fit(X, y)
+    # a policy is the composition of a feature map and a classifier
+    learned_policy = agglo.classifier_probability(fc, rf)
+
+    if not os.path.isdir('seg_data'):
+        print "\nRunning the commmand: mkdir seg_data"
+        os.mkdir('seg_data')
+    else:
+        command = "rm -f seg_data/*"
+        print "\nRunning the commmand: ", command
+        os.system(command)
+    if os.path.isfile('STACKED_prediction.h5'):
+        command = "rm -f STACKED_prediction.h5"
+        print "\nRunning the commmand: " , command
+        os.system(command)
+
+    # Produce the probabilities and oversegmentation files from the inputs
+    command = "python " + os.path.join(gala_home, "gala", "bin", "gala-segmentation-pipeline") + \
+              " . -I '" + productiobn_raw + "/*.png' --ilp-file training.ilp " + \
+              " --enable-gen-supervoxels --enable-gen-agglomeration " + \
+              " --enable-gen-pixel --seed-size 5   --segmentation-thresholds 0.0"
+    print "\nRunning the commmand: ", command
+    os.system(command)
+
+    # Put probabilities in the right format:
+    if os.path.isdir("seg_data"):
+        print "\nRunning the commmand: rm -rf seg_data"
+        os.system("rm -rf seg_data")
+    command = "MS_UT_ProbsNF2GALA.py STACKED_prediction.h5 probabilities_production.h5"
+    print "\nRunning the commmand: ", command
+    os.system(command)
+
+    # Rename the oversegmentation production file
+    for file in os.listdir("seg_data"):
+        if re.search(".h5", file):
+            oversegmentation_file = os.path.join("seg_data", file)
+    command = "mv " + oversegmentation_file + " oversegmentation_production.h5"
+
+    # get the production data and make a RAG with the trained policy
+    pr_test, ws_test = (map(imio.read_h5_stack,
+                        ['probabilities_production.h5', 'oversegmentation_production.h5']))
+    g_test = agglo.Rag(ws_test, pr_test, learned_policy, feature_manager=fc)
+    g_test.agglomerate(0.5) # best expected segmentation
+    seg_production = g_test.get_segmentation()
+
+    # Output the segmentation file
+    img = Image.fromarray(seg_production)
+    img.save('segmentation_production.h5')
+
+
+# -----------------------------------------------------------------------------
+
+def output_results(ws_test, gt_test, seg_test1):
+    results = np.vstack((
+        ev.split_vi(ws_test, gt_test),
+        ev.split_vi(seg_test1, gt_test),
+        ))
+    print(results)
+
+# -----------------------------------------------------------------------------
+
+if __name__ == "__main__":
+
+    usage = "Usage: \n\
+    %prog <raw_train_data> <gt_labels> <raw_data> [options (-h to list)]"
+
+    parser = optparse.OptionParser(usage=usage, version="%%prog ")
+    parser = MS_LIB_Options.GalaSegmentation2D_command_line_parser(parser)
+    (options, args) = parser.parse_args()
+
+    if len(args) == 4:
+        training_raw    = args[0] 
+        gt_annot_labels = args[1]
+        gt_seg_labels   = args[2]
+        production_raw  = args[3] 
+        perform_segmentation(training_raw, gt_annot_labels, gt_seg_labels, \
+                             production_raw, options)
+    else:
+        parser.print_usage()
+        sys.exit(2)
+ 

RhoanaSegmentation/MS_RS_RhoanaSegmentation2D.py

@@ -32,34 +32,25 @@ def produce_segmentation(input_data, probs_file_list, options):
         "Segment", "segment.py")
     for i in range(0, len(probs_file_list)):
         probs_file = ntpath.basename(probs_file_list[i])
-        output_file = probs_file.split(".")[0] + ".png"
-        output_dir = os.path.join(ms_data, "seg_" + options.output_tag)
+        output_file = probs_file.split(".")[0] + ".h5"
+        if re.search("mitochondria", options.classifier):
+            output_dir = os.path.join(ms_data, "segmentation_2D_bw_mitochondria_" \
+                                      + options.output_tag)
+        else:
+            output_dir = os.path.join(ms_data, "segmentation_2D_bw_membranes_" \
+                                      + options.output_tag)
         if not os.path.isdir(output_dir):
             os.mkdir(output_dir)
         command = "python " + segmentation_script + " " + probs_file_list[i] +\
-                  " " + os.path.join(ms_data, output_file)
+                  " " + os.path.join(output_dir, output_file)
+        if options.verbose:
+            print "\nSegmentation command=", command, "\n"
         os.system(command)
     return 
 
 # -----------------------------------------------------------------------------
 
-def invert_probabilities(h5_infile, h5_outfile):
-    print "h5_outfile=", h5_outfile
-    f   = h5py.File(h5_infile,  'r')
-    f2  = h5py.File(h5_outfile, 'w')
-    keys = f.keys()
-    for k in keys:
-        if not k == "probabilities":
-            f2.create_dataset(k,      data = f[k])
-        else: 
-            probs  = f["probabilities"]
-            inverse_probs = 1. - numpy.array(probs)
-            f2.create_dataset("probabilities", data = inverse_probs)
-    return 
-
-# -----------------------------------------------------------------------------
-
-def compute_membrane_probabilities(input_data, file_list, options):
+def compute_segmented_object_probabilities(input_data, file_list, options):
     classify_image_executable = options.executable
     classifier_file = options.classifier
     output_file_list = []
@@ -68,28 +59,19 @@ def compute_membrane_probabilities(input_data, file_list, options):
         print "input_file_path=", input_file_path
         input_file       = ntpath.basename(input_file_path)     
         output_file      = input_file.split(".")[0] + ".h5"
-        output_dir       = os.path.join(ms_data, "probs_"     + input_data)
+        input_tag = input_data
+        if input_data[:4] == "raw_":
+            input_tag = input_data[4:]
+        output_dir    = os.path.join(ms_data, options.seg_type + "_prob_" + input_tag)
         if not os.path.isdir(output_dir):
             os.mkdir(output_dir)
-        if input_data[0:4] == "raw_":
-            output_dir   = os.path.join(ms_data, "probs_"     + input_data[4:])
-            output_dir2  = os.path.join(ms_data, "probs_inv_" + input_data[4:])
         output_file_path = os.path.join(output_dir, output_file)
         output_file_list.append(output_file_path)
         command = classify_image_executable + " " + input_file_path + " " + \
                   classifier_file + " " + output_file_path     
         if options.verbose:
             print "...command=", command
         os.system(command) 
-        # Inverse probabilities
-        if options.inverse_probabilities:
-            output_dir2      = os.path.join(ms_data, "probs_inv_" + input_data)
-            if input_data[0:4] == "raw_":
-               output_dir2  = os.path.join(ms_data, "probs_inv_" + input_data[4:])
-            if not os.path.isdir(output_dir2):
-                os.mkdir(output_dir2)
-            output_file_path2= os.path.join(output_dir2, output_file)
-            invert_probabilities(output_file_path, output_file_path2)
     return output_file_list
 
 # -----------------------------------------------------------------------------
@@ -146,7 +128,7 @@ def process_input_data(input_data, input_type, input_label, options):
 
     if len(input_file_list) > 0:
         if int(options.processing_start) == 1:
-            probs_file_list = compute_membrane_probabilities(input_data, \
+            probs_file_list = compute_segmented_object_probabilities(input_data, \
                                   filtered_input_file_list, options)
 
         if int(options.processing_end) == 2:
@@ -190,6 +172,8 @@ def process_input_data(input_data, input_type, input_label, options):
             sys.exit(2)         
         if options.verbose:
             print "\nProcessing input data ..."
+        if len(options.output_tag) == 0:
+            options.output_tag = input_data
         process_input_data(input_data, input_type, input_label, options)
     else:
         parser.print_usage()