merge upstream

hippunfold/config/snakebids.yml

@@ -257,11 +257,6 @@ parse_args:
     nargs: '+'
     help: 'Run hipp (CA + subiculum) alone or include dentate (default: %(default)s)'
 
-  --keep_work:
-    help: 'Keep work folder intact instead of archiving it for each subject (default: %(default)s)'
-    default: False
-    action: 'store_true'
-
   --force_nnunet_model:
     help: 'Force nnunet model to use (expert option). (default: %(default)s)'
     default: False
@@ -287,6 +282,14 @@ parse_args:
     action: "store_false"
     default: True
 
+  --workdir:
+    help: |
+      Folder for storing working files. If not specified, a temporary folder
+      will be made in your system's temp directory (e.g. /tmp). You can also
+      use environment variables when setting the workdir, e.g. --workdir '$SLURM_TMPDIR'.
+    default: null
+    type: str
+
 
   --version:
     help: 'Print the version of HippUnfold'
@@ -608,11 +611,6 @@ unfold_vol_ref:
       - '2.5'
     orient: RPI
 
-unfold_crop_epsilon_fractions:
-    - 0
-    - 0
-    - 1 #in unfolded space voxels
-
 
 # space for uniform unfolded grid:
 #  currently only used for interpolating hipp subfields on surface
@@ -695,5 +693,6 @@ force_nnunet_model: False
 t1_reg_template: False
 generate_myelin_map: False
 root: results
+workdir: null
 use_template_seg: False
 template_seg_smoothing_factor: 2

hippunfold/workflow/Snakefile

@@ -2,6 +2,7 @@
 
 import snakebids
 from snakebids import bids, set_bids_spec
+from tempfile import mkdtemp
 
 
 configfile: "config/snakebids.yml"
@@ -90,8 +91,10 @@ wildcard_constraints:
     template="[a-zA-Z0-9]+",
 
 
-root = os.path.normpath(os.path.join(config["root"], "hippunfold"))
-work = os.path.normpath(os.path.join(config["root"], "work"))
+if config["workdir"] == None:
+    config["workdir"] = mkdtemp()
+work = os.path.expandvars(config["workdir"])
+root = os.path.expandvars(config["root"])
 
 
 include: "rules/common.smk"

hippunfold/workflow/rules/common.smk

@@ -271,7 +271,7 @@ def get_final_qc():
     return qc
 
 
-def get_final_subj_output():
+def get_final_output():
     subj_output = []
     subj_output.extend(get_final_spec())
     subj_output.extend(get_final_subfields())
@@ -280,16 +280,6 @@ def get_final_subj_output():
     subj_output.extend(get_final_qc())
     return subj_output
 
-
-def get_final_output():
-    if config["keep_work"]:
-        subj_output = get_final_subj_output()
-    else:
-        subj_output = get_final_work_tar()
-
-    return subj_output
-
-
 if "corobl" in ref_spaces:
 
     rule copy_coords_to_results:
@@ -323,19 +313,6 @@ if "corobl" in ref_spaces:
             "cp {input} {output}"
 
 
-def get_final_work_tar():
-    bids = bids_factory(specs.v0_0_0(subject_dir=False, session_dir=False))
-    return bids(root=work, suffix="work.tar.gz", **inputs.subj_wildcards)
-
-
-def get_work_dir(wildcards):
-    folder_with_file = inputs[config["modality"]].expand(
-        bids(root=work, **inputs.subj_wildcards), **wildcards
-    )
-    folder_without_file = os.path.dirname(folder_with_file[0])
-    return folder_without_file
-
-
 def get_download_dir():
     if "HIPPUNFOLD_CACHE_DIR" in os.environ.keys():
         download_dir = os.environ["HIPPUNFOLD_CACHE_DIR"]
@@ -344,21 +321,3 @@ def get_download_dir():
         dirs = AppDirs("hippunfold", "khanlab")
         download_dir = dirs.user_cache_dir
     return download_dir
-
-
-rule archive_work_after_final:
-    input:
-        get_final_subj_output(),
-    params:
-        work_dir=get_work_dir,
-    output:
-        get_final_work_tar(),
-    group:
-        "subj"
-    shell:
-        #exit code 0 or 1 is acceptable (2 is fatal)
-        "tar -czf {output} {params.work_dir}; "
-        "if [ $? -le 1 ]; then "
-        "  rm -rf {params.work_dir}; "
-        "else exit 1; "
-        "fi"

hippunfold/workflow/rules/native_surf.smk

@@ -338,6 +338,8 @@ rule compute_halfthick_mask:
         ),
     group:
         "subj"
+    container:
+        config["singularity"]["autotop"]
     shell:
         "c3d {input.coords} -threshold {params.threshold_tofrom} 1 0 {input.mask} -multiply -o {output}"
 
@@ -381,6 +383,8 @@ rule register_midthickness:
         ),
     group:
         "subj"
+    container:
+        config["singularity"]["autotop"]
     shell:
         "greedy -d 3 -i {input.fixed} {input.moving} -n 30 -o {output.warp}"
 
@@ -435,6 +439,8 @@ rule apply_halfsurf_warp_to_img:
         ),
     group:
         "subj"
+    container:
+        config["singularity"]["autotop"]
     shell:
         "greedy -d 3  -rf {input.fixed} -rm {input.moving} {output.warped}  -r {input.warp} "
 
@@ -469,6 +475,8 @@ rule convert_inout_warp_from_itk_to_world:
         ),
     group:
         "subj"
+    container:
+        config["singularity"]["autotop"]
     shell:
         "wb_command -convert-warpfield -from-itk {input} -to-world {output}"
 

hippunfold/workflow/scripts/create_warps.py

@@ -149,25 +149,15 @@ def summary(name, array):
 # we have points defined by coord_flat_{ap,pd,io}, and corresponding value as native_coords_phys[:,i]
 # and we want to interpolate on a grid in the unfolded space
 
-# unnormalize and add a bit of noise so points don't ever perfectly overlap
-coord_flat_ap_unnorm = coord_flat_ap * unfold_dims[0]
-coord_flat_pd_unnorm = coord_flat_pd * unfold_dims[1]
-coord_flat_io_unnorm = coord_flat_io * unfold_dims[2]
-
-# get unfolded grid
-
-points = np.stack(
-    [
-        coord_flat_ap * unfold_dims[0]
-        + (np.random.rand(coord_flat_ap.shape[0]) - 0.5) * 1e-6,
-        coord_flat_pd * unfold_dims[1]
-        + (np.random.rand(coord_flat_ap.shape[0]) - 0.5) * 1e-6,
-        coord_flat_io * unfold_dims[2]
-        + (np.random.rand(coord_flat_ap.shape[0]) - 0.5) * 1e-6,
-    ],
-    axis=1,
+# add some noise to avoid perfectly overlapping datapoints!
+points = (
+    coord_flat_ap * unfold_dims[0]
+    + (np.random.rand(coord_flat_ap.shape[0]) - 0.5) * 1e-6,
+    coord_flat_pd * unfold_dims[1]
+    + (np.random.rand(coord_flat_ap.shape[0]) - 0.5) * 1e-6,
+    coord_flat_io * unfold_dims[2]
+    + (np.random.rand(coord_flat_ap.shape[0]) - 0.5) * 1e-6,
 )
-summary("points", points)
 
 
 # perform the interpolation