Merge pull request #11 from ledigchr/malpem_v1.3

Updates for MALPEM v1.3 release.

Functional:
* Updated initial N4 bias correction from potentially using OTSU mask to always use complete image domain as foreground
* Perform consistently two N4 bias correction (1st: N4 on full image domain, 2nd: N4 w/ mask) regardless of whether user specified a mask
* Updated installer script to allow user to select among various version

Reproducibility:
* Fix issue #8 : Fixed itk threads to 8 to make N4 bias correction deterministic (depends on number of threads)

Robustness/Compatibility:
* Fix issue #9 : Ensure input image and mask have numerically same header information to avoid N4 failing, don't fail silently in case of another unexpected problem.
* Fix issue #2 : Added 'export PROOT_NO_SECCOMP=1' to avoid Error 'root info: pid XXX: terminated with signal 11' 
* ingore_errors = True when attempting to delete temporary files when 'cleanup' flag is used

Minor:
* removed unused imports, cosmetics

bin/malpem

@@ -213,7 +213,9 @@ def main(argv):
     malpem.mytools.check_ex_dir(tmp_dir)
 
     image_n4 = os.path.join(output_dir, base_file + "_N4.nii.gz")
+    image_n4_initial = os.path.join(tmp_dir, base_file + "_N4_initial.nii.gz")
     image_n4_masked = os.path.join(output_dir, base_file + "_N4_masked.nii.gz")
+    image_full_mask = os.path.join(output_dir, base_file + "_mask_full_image.nii.gz")
 
     if not do_n4:
         image_n4 = os.path.join(output_dir, base_file + "_noN4.nii.gz")
@@ -301,10 +303,10 @@ def main(argv):
 
     # Run N4 Bias Correction (ITK implementation, parameters depend on 1.5T/3T)
     if not os.path.isfile(image_n4):
-        if not os.path.isfile(image_mask):
-            malpem.bias_correction.N4(input_file, image_n4, field_strength, "", output_dir)
-        else:
-            malpem.bias_correction.N4(input_file, image_n4, field_strength, image_mask, output_dir)
+        # perform bias correction with the full image domain as foreground to avoid using unpredictable OTSU mask
+        print "Creating artificial 'full image mask' to perform bias correction on complete image domain"
+        malpem.bias_correction.get_full_image_mask(input_file, image_full_mask)
+        malpem.bias_correction.N4(input_file, image_n4, field_strength, image_full_mask, output_dir)
     else:
         print "Skipping bias correction: file exists / specified by user"
 
@@ -319,10 +321,19 @@ def main(argv):
     if not os.path.isfile(image_n4_masked):
         if not os.path.isfile(image_mask):
             malpem.brain_extraction.pincram(image_n4, image_mask, threads, output_dir)
-            # Note that the MNI transformation above was calculated with the originally N4 corrected image
-            if do_n4_pincram:
-                print "Running a second bias correction with the calculated brain mask"
-                malpem.bias_correction.N4(image_n4, image_n4, field_strength, image_mask, output_dir)
+
+        print "Transforming image/mask to ensure numerically matching headers (workaround)"
+        # Nearest neighbor interpolation since this is merely about numerical adjustment of the header
+        # and should have no effect on any pixel values
+        malpem.registration.transform(image_n4, image_mask, image_mask, "", "nn", output_dir)
+        malpem.registration.transform(image_n4, image_n4, image_n4, "", "nn", output_dir)
+
+        # Note that the MNI transformation above was calculated with the originally N4 corrected image
+        if do_n4_pincram:
+            print "Running a second bias correction with the calculated brain mask"
+            # Move file to make sure file doesn't exist if this N4 fails and MALPEM doesn't continue silently
+            shutil.move(image_n4, image_n4_initial)
+            malpem.bias_correction.N4(image_n4_initial, image_n4, field_strength, image_mask, output_dir)
 
         malpem.brain_extraction.apply_mask(image_n4, image_n4_masked, image_mask)
     else:
@@ -427,15 +438,15 @@ def main(argv):
     if cleanup:
         task_name_cu = "Cleaning up directories / deleting tmp files"
         start_time_cu = malpem.mytools.start_task(task_name_cu)
-        shutil.rmtree(os.path.join(output_dir, "tmp"))
-        shutil.rmtree(os.path.join(output_dir, "tmp_fusion"))
-        shutil.rmtree(os.path.join(output_dir, "tmp_malpem"))
-        shutil.rmtree(os.path.join(output_dir, "tmp_pincram"))
-        shutil.rmtree(os.path.join(output_dir, "dofs"))
+        shutil.rmtree(os.path.join(output_dir, "tmp"), ignore_errors=True)
+        shutil.rmtree(os.path.join(output_dir, "tmp_fusion"), ignore_errors=True)
+        shutil.rmtree(os.path.join(output_dir, "tmp_malpem"), ignore_errors=True)
+        shutil.rmtree(os.path.join(output_dir, "tmp_pincram"), ignore_errors=True)
+        shutil.rmtree(os.path.join(output_dir, "dofs"), ignore_errors=True)
 
 	malpem.mytools.finished_task(start_time_cu, task_name_cu)
 
     malpem.mytools.finished_task(start_time, task_name)
 
 if __name__ == "__main__":
-    main(sys.argv[1:])
+    main(sys.argv[1:])

bin/malpem-proot

@@ -1,5 +1,7 @@
 #! /bin/bash
 
+export PROOT_NO_SECCOMP=1
+
 function ensureDir {
         if [ ! -d $1 ]; then
                 mkdir $1;
@@ -193,6 +195,7 @@ $PROOT_NO_SECCOMP \
 'TEXTDOMAINDIR=/usr/share/locale/'  \
 'COLORTERM=gnome-terminal'  \
 'XAUTHORITY=/tmp/.mdmCHXAXX'  \
+'ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=8' \
 "${PROOT}" \
 -b "/dev"  \
 -b "/proc"  \

installer/malpem-install

@@ -36,7 +36,34 @@ import errno
 import platform
 
 # MALPEM version and source URL
-version = "1.2"
+
+DEFAULT_VERSION="1.3"
+VERSIONS = {
+    "1.2": "Primary version since 2015",
+    "1.3": "(default) Updated version since 05/2021 w/ "
+           "robustness, reproducibility improvements over v1.2"
+}
+
+def get_version():
+    input_text = "\nChoose Version (for details: https://github.com/ledigchr/MALPEM)\n"
+    for k, v in VERSIONS.items():
+        input_text += "[" + k + "]: " + v + "\n"
+
+    version = raw_input(input_text)
+
+    if version == "":
+        version = DEFAULT_VERSION
+        print "Using default version " + DEFAULT_VERSION
+    else:
+        print "Using version " + version
+
+    msg = "version not in " + ','.join(VERSIONS.keys())
+    assert version in VERSIONS.keys(), msg
+
+    return version
+
+
+version = get_version()
 base_url = "http://www.christianledig.com/Material/MALPEM/"
 
 

lib/malpem/bias_correction.py

@@ -4,8 +4,6 @@
 #
 #         see license file in project root directory
 
-import time
-from subprocess import call
 import os.path
 import malpem.mytools
 
@@ -26,28 +24,42 @@ def N4(input_file, output_file, field_strength, input_mask, output_dir):
     if input_mask == "":
         if field_strength == '1.5T':
             parameters_N4 = " -d 3 -i " + input_file + " -o " + output_file + " -s 2 -c [50x40x30x20,0.0000001] " \
-                                                                              "-b [200,3,0.0,0.5] -t [0.15,0.01,200] 2"
+                                                                              "-b [200,3,0.0,0.5] -t [0.15,0.01,200]"
         elif field_strength == '3T':
             parameters_N4 = " -d 3 -i " + input_file + " -o " + output_file + " -s 2 -c [50x40x30x20,0.0000001] " \
-                                                                              "-b [75,3,0.0,0.5] -t [0.15,0.01,200] 2"
+                                                                              "-b [75,3,0.0,0.5] -t [0.15,0.01,200]"
         else:
             print "Warning N4: Unknown field strength defaulting to 1.5T parameters"
             parameters_N4 = " -d 3 -i " + input_file + " -o " + output_file + " -s 2 -c [50x40x30x20,0.0000001] " \
-                                                                              "-b [200,3,0.0,0.5] -t [0.15,0.01,200] 2"
+                                                                              "-b [200,3,0.0,0.5] -t [0.15,0.01,200]"
     else:
         if field_strength == '1.5T':
             parameters_N4 = " -d 3 -i " + input_file + " -x " + input_mask + " -o " + output_file + " -s 2 " \
-                            "-c [50x40x30x20,0.0000001] -b [200,3,0.0,0.5] -t [0.15,0.01,200] 2"
+                            "-c [50x40x30x20,0.0000001] -b [200,3,0.0,0.5] -t [0.15,0.01,200]"
         elif field_strength == '3T':
             parameters_N4 = " -d 3 -i " + input_file + " -x " + input_mask + " -o " + output_file + " -s 2 " \
-                            "-c [50x40x30x20,0.0000001] -b [75,3,0.0,0.5] -t [0.15,0.01,200] 2"
+                            "-c [50x40x30x20,0.0000001] -b [75,3,0.0,0.5] -t [0.15,0.01,200]"
         else:
             print "Warning N4: Unknown field strength defaulting to 3T parameters"
             parameters_N4 = " -d 3 -i " + input_file + " -x " + input_mask + " -o " + output_file + " -s 2 " \
-                            "-c [50x40x30x20,0.0000001] -b [75,3,0.0,0.5] -t [0.15,0.01,200] 2"
+                            "-c [50x40x30x20,0.0000001] -b [75,3,0.0,0.5] -t [0.15,0.01,200]"
 
     malpem.mytools.execute_cmd(binary_N4, parameters_N4, logfile)
     malpem.mytools.ensure_file(output_file, "")
 
     malpem.mytools.finished_task(start_time, task_name)
     return True
+
+
+def get_full_image_mask(input_file, output_file):
+# DEFINITIONS
+    binary_fsl = os.path.join(malpem.mytools.__malpem_path__, "lib", "niftyseg", "seg_maths")
+# END
+    task_name = "getting mask for full image"
+    start_time = malpem.mytools.start_task(task_name)
+
+    parameters_binarise = input_file + " -bin -add 1 -bin " + output_file
+    malpem.mytools.execute_cmd(binary_fsl, parameters_binarise, "")
+    malpem.mytools.ensure_file(output_file, "")
+
+    malpem.mytools.finished_task(start_time, task_name)

lib/malpem/brain_extraction.py

@@ -74,8 +74,8 @@ def binarise(input_file, output_file):
     task_name = "binarising image for brain mask"
     start_time = malpem.mytools.start_task(task_name)
 
-    parameters_biarise = input_file + " -bin " + output_file
-    malpem.mytools.execute_cmd(binary_fsl, parameters_biarise, "")
+    parameters_binarise = input_file + " -bin " + output_file
+    malpem.mytools.execute_cmd(binary_fsl, parameters_binarise, "")
     malpem.mytools.ensure_file(output_file, "")
 
     malpem.mytools.finished_task(start_time, task_name)

lib/malpem/mytools.py

@@ -43,7 +43,7 @@ def ensure_file(filename, message):
 
 def execute_cmd(cmd, parameters, logfile):
     if logfile == "":
-        final_cmd = final_cmd = cmd + " " + parameters
+        final_cmd = cmd + " " + parameters
     else:
         final_cmd = cmd + " " + parameters + " >> " + logfile + " 2>&1"
         f = open(logfile, 'w')