Merge branch 'release-0.1.19' -- which actually is 0.1.20. The branch…

… was misnamed

NEWS

@@ -1,3 +1,10 @@
+New in Version 0.1.20
+=====================
+* rotational_minctracc.py does not keep any intermediate files anymore; everything is deleted
+  during the process, and the final transformation is re-created at the end.
+* rotational_minctracc.py now can also handle ctrl-c (SIGINT) in its cleanup.
+* rotational_minctracc.py can be used for lsq12 instead of lsq6 brute force alignment
+
 New in Version 0.1.19
 =====================
 * added script: make_xfm_for_grid.pl which creates an xfm file for a provided deformation grid

python/rotational_minctracc.py

@@ -162,14 +162,18 @@ def resample_volume(source, target, transform):
                           % (transform, target, source, tmp_resampled)).split())
     return tmp_resampled
 
-def minctracc(source, target, mask, stepsize, wtranslations, simplex):
+def minctracc(source, target, mask, stepsize, wtranslations, simplex, use_lsq12_for_alignment):
     wtrans_decomp = array(wtranslations.split(',')).astype("float")
     tmp_transform = get_tempfile('.xfm')
-    cmd = ("minctracc -identity -lsq6 -xcorr -simplex %s -step %s %s %s %s %s %s -w_translations %s %s %s " 
+    cmd = ("minctracc -identity -xcorr -simplex %s -step %s %s %s %s %s %s -w_translations %s %s %s " 
            % (simplex, stepsize, stepsize, stepsize, source, target, tmp_transform,
            wtrans_decomp[0], wtrans_decomp[1], wtrans_decomp[2]))
     if mask:
         cmd += ("-source_mask %s -model_mask %s " % (mask, mask))
+    if use_lsq12_for_alignment:
+        cmd += "-lsq12"
+    else:
+        cmd += "-lsq6"
     print(cmd)
     subprocess.check_call(cmd.split())
 
@@ -191,11 +195,14 @@ def get_cross_correlation_from_coordinate_pair(source_img, target_img, target_vo
     transform_from_coordinates = create_transform(coordinate_pair[1] - coordinate_pair[0], 0, 0, 0, coordinate_pair[0])
     resampled_source = resample_volume(source_img, target_img, transform_from_coordinates)
     xcorr = compute_xcorr(resampled_source, target_vol, mask)
+    # clean up. We do not need either the MINC file nor the transform anymore
     os.remove(resampled_source)
+    os.remove(transform_from_coordinates)
     return float(xcorr)
 
 def loop_rotations(stepsize, source, target, mask, simplex, start=50, interval=10, 
-                   wtranslations="0.2,0.2,0.2", use_multiple_seeds=True, max_number_seeds=5):
+                   wtranslations="0.2,0.2,0.2", use_multiple_seeds=True, max_number_seeds=5,
+                   use_lsq12_for_alignment=False):
     # load the target and mask volumes
     targetvol = volumeFromFile(target)
     maskvol = volumeFromFile(mask) if mask is not None else None
@@ -259,7 +266,7 @@ def loop_rotations(stepsize, source, target, mask, simplex, start=50, interval=1
             print(list_of_coordinate_pairs)
 
     results = []
-    best_xcorr = 0
+    #best_xcorr = 0
     for coordinates_src_target in list_of_coordinate_pairs:
         coor_src = coordinates_src_target[0]
         coor_trgt = coordinates_src_target[1]
@@ -270,27 +277,44 @@ def loop_rotations(stepsize, source, target, mask, simplex, start=50, interval=1
                     init_transform = create_transform(coor_trgt - coor_src, x, y, z, coor_src)
                     init_resampled = resample_volume(source, target, init_transform)
                     transform = minctracc(init_resampled, target, mask, stepsize=stepsize,
-                                          wtranslations=wtranslations, simplex=simplex)
+                                          wtranslations=wtranslations, simplex=simplex,
+                                          use_lsq12_for_alignment=use_lsq12_for_alignment)
                     resampled = resample_volume(init_resampled, target, transform)
                     conc_transform = concat_transforms(init_transform, transform)
                     xcorr = compute_xcorr(resampled, targetvol, maskvol)
                     if isnan(xcorr):
                         xcorr = 0
-                    results.append({'xcorr': xcorr, 'transform': conc_transform, \
-                                    'resampled': resampled, 'x': x, \
-                                    'y': y, 'z': z})
-                    if xcorr > best_xcorr:
-                        best_xcorr = xcorr
+                    results.append({'xcorr': xcorr, 'transform': conc_transform,
+                                    'resampled': resampled, 'xrot': x,
+                                    'yrot': y, 'zrot': z, 'coor_src' : coor_src, 'coor_trgt' : coor_trgt })
+                    #if xcorr > best_xcorr:
+                    #    best_xcorr = xcorr
                     # had some issues with the resampled file being gone...
                     # we'll just resample the final file only at the end
                     os.remove(resampled)
                     os.remove(init_resampled)
+                    os.remove(conc_transform)
+                    os.remove(init_transform)
+                    os.remove(transform)
                     print("FINISHED: %s %s %s :: %s" % (x,y,z, xcorr))
 
     sort_results(results)
     # resample the best result:
-    final_resampled = resample_volume(source, target, results[-1]["transform"])
+    # TODO this is the same code as the inner loop above -- make a procedure?
+    best = results[-1]
+    best_init_transform = create_transform(best['coor_trgt'] - best['coor_src'],
+                                           best['xrot'], best['yrot'], best['zrot'],
+                                           best['coor_src'])
+    best_init_resampled = resample_volume(source, target, best_init_transform)
+    best_transform = minctracc(best_init_resampled, target, mask, stepsize=stepsize,
+                                          wtranslations=wtranslations, simplex=simplex,
+                                          use_lsq12_for_alignment=use_lsq12_for_alignment)
+    best_resampled = resample_volume(best_init_resampled, target, best_transform)
+    best_conc_transform = concat_transforms(best_init_transform, best_transform)
+    final_resampled = resample_volume(source, target, best_conc_transform)
+    #final_resampled = resample_volume(source, target, results[-1]["transform"])
     results[-1]["resampled"] = final_resampled
+    results[-1]["transform"] = best_conc_transform
     targetvol.closeVolume()
     if mask is not None:
         maskvol.closeVolume()
@@ -311,12 +335,14 @@ def downsample(infile, stepsize):
 # clean up tmp on soft kill signal
 def termtrapper(signum, frame):
     print("got kill signal")
-    os.removedirs("%s/rot_%s" % (os.environ["TMPDIR"], os.getpid()))
+    shutil.rmtree("%s/rot_%s" % (os.environ["TMPDIR"], os.getpid()))
+    print("\n\nWent down gracefully!\n")
     exit(1)
 
 if __name__ == "__main__":
     # handle soft kill signal to clean up tmp
     signal.signal(signal.SIGTERM, termtrapper)
+    signal.signal(signal.SIGINT, termtrapper)
 
     parser = ArgumentParser()
     parser.add_argument("-m", "--mask", dest="mask",
@@ -356,6 +382,13 @@ def termtrapper(signum, frame):
                         "pairs are ordered based on the cross correlation gotten "
                         "from the alignment based on only the translation from the "
                         "seed point. [default = %(default)s]")
+    parser.set_defaults(use_lsq12_for_alignment=False)
+
+    parser.add_argument("--use-lsq12-for-alignment", dest="use_lsq12_for_alignment", action="store_true",
+                        help="Instead of aligning the files using rotations and translations only "
+                             "use 3 scaling parameters as well. [default = %(default)s]")
+    parser.add_argument("--no-use-lsq12-for-alignment", dest="use_lsq12_for_alignment", action="store_false",
+                        help="Opposite of --use-lsq12-for-alignment")
     parser.add_argument("source", help="", type=str, metavar="source.mnc")
     parser.add_argument("target", help="", type=str, metavar="target.mnc")
     parser.add_argument("output_xfm", help="", type=str, metavar="output.xfm")
@@ -391,7 +424,8 @@ def termtrapper(signum, frame):
                              wtranslations=options.wtranslations, 
                              simplex=options.simplex,
                              use_multiple_seeds=options.use_multiple_seeds,
-                             max_number_seeds=options.max_number_seeds)
+                             max_number_seeds=options.max_number_seeds,
+                             use_lsq12_for_alignment=options.use_lsq12_for_alignment)
 
     print(results)
     subprocess.check_call(("cp %s %s" % (results[-1]["transform"], output_xfm)).split())

setup.py

@@ -5,7 +5,7 @@
 import sys
 
 setup(name="python-stuffs",
-      version='0.1.19',
+      version='0.1.20',
       install_requires=['numpy', 'scipy', 'pyminc'],
       scripts=["python/TFCE",
                "python/smooth_vector",