Merge branch 'release-0.1.15'

NEWS

@@ -1,3 +1,16 @@
+New in Version 0.1.15
+=====================
+* fairly big overhaul of rotational_minctracc.py: it turns out that the centre of gravity (cog)
+  of the two input files to rotational_minctracc.py do not always reflect the same part of the 
+  object. In those cases, the alignment most likely will fail because a mask is used, and 
+  representative areas of the input files might not even be compared. For this reason 
+  rotation_minctracc.py now also runs a "translational" part. Multiple centers of gravity
+  are determined which serve as seeds for the initial translational alignment. All possible
+  pairs of seeds are tried and by default the best 3 pairs are used for the rotational part
+  of the code.
+* some more Python3 compatibility related fixes
+* added options for minc_label_ops
+
 New in Version 0.1.14
 =====================
 * completely remove the multiprocessing machinery from pmincaverage since this caused extra processes

python/minc_label_ops

@@ -19,6 +19,10 @@ if __name__ == "__main__":
     parser.add_option("--merge", dest="merge",
                       help="Second label file to whose labels are to be merged with the input_labels",
                       type="string")
+    parser.add_option("--replace-all-different", dest="replace_all_different",
+                      help="Similar to the --merge option, but in this case, when the given "
+                      "file is zero, and the other input is not, it will be set to zero. "
+                      "The idea is that the fixes include voxels that were actually background.")
     parser.add_option("--remap", dest="remap",
                       help="Label numbers to replace with new label numbers. For example, to turn label 3 into label 11, it would be --remap 3:11. More than one can be specified at once, separated by commas, i.e. --remap 3:11,21:55. It is also possible to specify a range of values to be remapped as follows: to map 3,4,5,6,7, and 8 onto the value 11, specify --remap 3-8:11. This can also be combined: --remap 3-8:11,66:98.",
                       type="string")
@@ -31,18 +35,28 @@ if __name__ == "__main__":
     parser.add_option("--mask", dest="mask",
                       help="Erase all labels where mask >0.5",
                       type="string")
+    parser.add_option("--diff", help="Output all voxels (binary) that differ between the two input files (difference labels are taken from the first input file)", action="store_true", default=False)
     parser.add_option("--binarize", help="threshold >0.5 and binarize output", action="store_true", default=False)
     parser.add_option("-v", "--verbose", help='verbose output', action='store_true', default=False)
 
     (options, args) = parser.parse_args()
 
-    if len(args) != 2:
+    if not options.diff and len(args) != 2:
         parser.error("Incorrect number of arguments")
 
 
-    inlabels = volumeFromFile(args[0], dtype='ushort')
-    outlabels = volumeLikeFile(args[0], args[1], dtype='ushort', 
-                               volumeType="ushort", labels=True)
+    if options.diff and len(args) != 3:
+        parser.error("Incorrect number of arguments. Please specify two input label files and an output label file to hold the difference between the two input files")
+
+    if not options.diff:
+        inlabels = volumeFromFile(args[0], dtype='ushort')
+        outlabels = volumeLikeFile(args[0], args[1], dtype='ushort', 
+                                   volumeType="ushort", labels=True)
+    else:
+        inlabels_1 = volumeFromFile(args[0], dtype='ushort')
+        inlabels_2 = volumeFromFile(args[1], dtype='ushort')
+        outlabels  = volumeLikeFile(args[0], args[2], dtype='ushort', 
+                                    volumeType="ushort", labels=True)
 
     if options.verbose:
         print("Labels found in input: %s" % unique(inlabels.data[::]))
@@ -98,6 +112,11 @@ if __name__ == "__main__":
         outlabels.data[::] = inlabels.data[::] + addlabels.data[::]
         addlabels.closeVolume()
 
+    if options.replace_all_different:
+        addlabels = volumeFromFile(options.replace_all_different, dtype='ushort', labels=True)
+        outlabels.data = where(addlabels.data != inlabels.data, addlabels.data, inlabels.data)
+        addlabels.closeVolume()
+
     if options.mask:
         maskforlabels = volumeFromFile(options.mask, dtype='ushort')
         # initially copy over input labels
@@ -106,11 +125,19 @@ if __name__ == "__main__":
         outlabels.data[maskforlabels.data > 0.5] = 0
         maskforlabels.closeVolume()
 
+    if options.diff:
+        outlabels.data = where(inlabels_1.data != inlabels_2.data, 1, 0)
+
+
     if options.binarize:
         outlabels.data[inlabels.data > 0.5] = 1
 
     # write to file
     outlabels.writeFile()
     outlabels.closeVolume()
-    inlabels.closeVolume()
+    if not options.diff:
+        inlabels.closeVolume()
+    else:
+        inlabels_1.closeVolume()
+        inlabels_2.closeVolume()
 

python/pmincaverage

@@ -5,7 +5,10 @@ from numpy import *
 from argparse import ArgumentParser
 import os.path
 import sys
-import Queue
+try:
+    import queue
+except:
+    import Queue as queue
 
 def getslice(volhandle, slice, q,nslices):
     """collects hyperslabs for a particular volume handle. Designed to
@@ -94,7 +97,7 @@ if __name__ == "__main__":
 
         # without multiprocessing, this queueing stuff probably doesn't make sense -- simplify
         print("SLICE: %i" % i)
-        q = Queue.Queue()
+        q = queue.Queue()
         for j in range(nfiles):
             t = getslice(j,i,q,nslices)
 

python/replace_label_with_nearest_valid_label

@@ -3,13 +3,15 @@
 
 from pyminc.volumes.factory import *
 import numpy as np
+import scipy.stats as scipy_stats
 from optparse import OptionParser
 import os
 import os.path
-import cython_functions
 import datetime as datetime
 import sys
 
+SDTYPE = np.uint16
+
 ################################################################################
 ################################################################################
 ################################################################################
@@ -27,6 +29,75 @@ def is_minc_2_file(filename):
 ################################################################################
 ################################################################################
 ################################################################################
+
+#
+# This function takes an input data structure which
+# contains label values and a list of unwanted/undesired
+# labels. The output data structure contains the input
+# labels where the unwanted/undesired labels have been
+# replaced by the closest valid label as follows:
+#
+# For an unwanted label:
+#
+# 1) all 26 neighbors of a voxel are considered to have
+#    the same distance. If 1 or more of these 26 
+#    neighbors has a valid label (not 0, nor one of the
+#    values in the unwanted list), it will be used
+#    to re-label the unwanted label (using the 
+#    most occuring label)
+#
+# Input:  inlabels,       the current set of labels
+#         outlables,      a copy of the current set of labels, to be edited
+#         unwantedlabels, an array containing the labels numbers that need to be replaced
+#
+# Output: the function will return the number of labels that could not be replaced,
+#         and outlabels will be updated as described above
+# 
+def substitute_labels_using_26_nearest_neighbors(labels,
+                                                 outlabels,
+                                                 unwantedlabels):
+
+  # make sure that the input files have the appropriate types:
+  assert inlabels.dtype == SDTYPE and outlabels.dtype == SDTYPE and unwantedlabels.dtype == SDTYPE
+
+  # get dimension information from the minc volume
+  nv0 = inlabels.shape[0]
+  nv1 = inlabels.shape[1]
+  nv2 = inlabels.shape[2]
+
+  # the number of unwanted labels
+  nlabels = unwantedlabels.shape[0]
+
+  # number of unwanted labels that could not be replaced using its 26 neighbors
+  nunresolved = 0
+
+  # from each dimension exclude the first and last element, because these have
+  # a number of undefined neighbors: the ones that fall outside of the file
+  for v0 in range(1, nv0 - 1):
+    for v1 in range(1, nv1 - 1):
+      for v2 in range(1, nv2 - 1):
+        # first check whether this label value should be replaced
+        if( np.rint(inlabels[v0,v1,v2]) in unwantedlabels):
+          # create array of possible substitues
+          possible_labels = []
+          # we can safely look at all 27 voxels in the 3*3 block, because
+          # the voxel in question will not be added to the list of 
+          # possible labels
+          for i0 in (v0-1, v0, v0+1):
+            for i1 in (v1-1, v1, v1+1):
+              for i2 in (v2-1, v2, v2+1):
+                if( not(np.rint(inlabels[i0,i1,i2]) == 0) and not(np.rint(inlabels[i0,i1,i2]) in unwantedlabels)):
+                  possible_labels.append(np.rint(inlabels[i0,i1,i2]))
+          # potentially we did not find any valid substitues
+          if len(possible_labels) == 0:
+            nunresolved += 1
+          else: # we have candidates, simply take the most occuring label number
+            outlabels[v0,v1,v2] = np.rint(scipy_stats.mode(possible_labels)[0])
+
+  return nunresolved
+
+
+
 if __name__ == "__main__":
 
   usage = "usage: %prog [options] input_labels.mnc updated_labels.mnc\n"
@@ -94,13 +165,13 @@ only desired label numbers.\n"""
   array_undesired_labels = np.fromstring(options.undesired_labels, dtype=np.uint16, sep=',')
 
   previous_unresolved = -1
-  current_unresolved = cython_functions.substitute_labels_using_26_nearest_neighbors(inlabels.data, outlabels.data, array_undesired_labels)
+  current_unresolved = substitute_labels_using_26_nearest_neighbors(inlabels.data, outlabels.data, array_undesired_labels)
   if(options.verbose):
     print("Number of unresolved labels: %f \n" % (current_unresolved))
   while( current_unresolved > 0 and not(previous_unresolved == current_unresolved)):
     previous_unresolved = current_unresolved
     copy_of_data = outlabels.data[::]
-    current_unresolved = cython_functions.substitute_labels_using_26_nearest_neighbors(copy_of_data, outlabels.data, array_undesired_labels)
+    current_unresolved = substitute_labels_using_26_nearest_neighbors(copy_of_data, outlabels.data, array_undesired_labels)
     if(options.verbose):
       print("Number of unresolved labels: %f \n" % (current_unresolved))
 
@@ -181,4 +252,4 @@ only desired label numbers.\n"""
       #####ndx = distance.argsort()                                                # indirect sort
       #####print("Current point <%s> " % (point))
       #####print("closest point <%s> " % (label_coordinates[ndx[:1]]))
-    #####label_iterator_main.iternext()
+    #####label_iterator_main.iternext()