From 67debdccc1f42aecf66ea8030a7aa57fc87164f5 Mon Sep 17 00:00:00 2001 From: Matthijs van Eede Date: Mon, 2 May 2016 10:14:04 -0400 Subject: [PATCH] * binarizing the input files based on the bimodalt value does not always produce the desired result. To ensure we really find the best seeds to start the rotational part of the registration with, we now also look for peaks in the blurred input files. * there have been some instances where rotational_minctracc.py failed because it wasn't able to copy the resampled file belonging to the best transform. Changed the code to remove all files while running the main loop, and generate the desired sampled file only when we've already found the best transform * lowered the number of seed pairs to use from 5 to 3. Now that we are using a potential high number of possible starting points, there's a very high chance that one of the top few already creates a really good starting point --- python/rotational_minctracc.py | 37 +++++++++++++++++++++++++++------- 1 file changed, 30 insertions(+), 7 deletions(-) diff --git a/python/rotational_minctracc.py b/python/rotational_minctracc.py index 1a53a1c..d1df504 100755 --- a/python/rotational_minctracc.py +++ b/python/rotational_minctracc.py @@ -94,6 +94,14 @@ def get_distance_transform_peaks(input_file, peak_distance): subprocess.check_call(("find_peaks -pos_only -min_distance %s %s %s" % (peak_distance, distance_transform, peak_tags)).split()) all_coors = get_coordinates_from_tag_file(peak_tags) return all_coors + +def get_blur_peaks(input_file, blur_kernel, peak_distance): + blurred_input = get_tempfile('_blur.mnc') + subprocess.check_call(("mincblur -no_apo -fwhm %s %s %s" % (blur_kernel, input_file, blurred_input.split('_blur.mnc')[0])).split()) + peak_tags = get_tempfile('.tag') + subprocess.check_call(("find_peaks -pos_only -min_distance %s %s %s" % (peak_distance, blurred_input, peak_tags)).split()) + all_coors = get_coordinates_from_tag_file(peak_tags) + return all_coors def compute_xcorr(sourcefile, targetvol, maskvol): try: @@ -185,16 +193,27 @@ def loop_rotations(stepsize, source, target, mask, simplex, start=50, interval=1 # between peaks is based on the stepsize used for the # registrations if use_multiple_seeds: - list_source_peaks = get_distance_transform_peaks(input_file=source, peak_distance=stepsize*10) - print("\n\nPeaks found in the source image:") + list_source_peaks = get_distance_transform_peaks(input_file=source, peak_distance=stepsize) + print("\n\nPeaks found in the source image (Distance Transform):") for coor_src in list_source_peaks: print(coor_src) + # also add peaks from the blurred version of the input file + blurred_peaks_source = get_blur_peaks(input_file=source, blur_kernel=stepsize, peak_distance=stepsize) + print("\n\nPeaks found in the source image (blurrred image):") + for coor_src in blurred_peaks_source: + print(coor_src) + list_source_peaks.append(coor_src) # also add the center of gravity of the source image list_source_peaks.append(cog_source) - list_target_peaks = get_distance_transform_peaks(input_file=target, peak_distance=stepsize*10) - print("\n\nPeaks found in the target image:") + list_target_peaks = get_distance_transform_peaks(input_file=target, peak_distance=stepsize) + print("\n\nPeaks found in the target image (Distance Transform):") for coor_trgt in list_target_peaks: print(coor_trgt) + blurred_peaks_target = get_blur_peaks(input_file=target, blur_kernel=stepsize, peak_distance=stepsize) + print("\n\nPeaks found in the target image (blurrred image):") + for coor_target in blurred_peaks_target: + print(coor_target) + list_target_peaks.append(coor_target) # same for the target; add the center of gravity: list_target_peaks.append(cog_target) for source_coor in list_source_peaks: @@ -244,12 +263,16 @@ def loop_rotations(stepsize, source, target, mask, simplex, start=50, interval=1 'y': y, 'z': z}) if xcorr > best_xcorr: best_xcorr = xcorr - else: - os.remove(resampled) + # 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) 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"]) + results[-1]["resampled"] = final_resampled targetvol.closeVolume() if mask is not None: maskvol.closeVolume() @@ -308,7 +331,7 @@ def termtrapper(signum, frame): "(of the intensities) of the input files. [default = %(default)s]") parser.add_argument("--no-use-multiple-seeds", dest="use_multiple_seeds", action="store_false", help="Opposite of --use-multiple-seeds") - parser.set_defaults(max_number_seeds=5) + parser.set_defaults(max_number_seeds=3) parser.add_argument("--max-number-seeds", dest="max_number_seeds", type=int, help="Specify the maximum number of seed-pair starting points " "to use for the rotational part of the code. The seed "