Merge pull request #353 from khanlab/speeed

Laplace-beltrami surface-based approach for generating coords

hippunfold/workflow/Snakefile

@@ -124,7 +124,7 @@ if config["modality"] == "hippb500":
     include: "rules/preproc_hippb500.smk"
 
 
-include: "rules/autotop.smk"
+include: "rules/coords.smk"
 
 
 include: "rules/warps.smk"

hippunfold/workflow/rules/common.smk

@@ -92,58 +92,6 @@ def get_final_subfields():
     )
 
 
-def get_final_coords():
-    if "laplace" in config["laminar_coords_method"]:
-        desc_io = "laplace"
-    elif "equivolume" in config["laminar_coords_method"]:
-        desc_io = "equivol"
-
-    coords = []
-    # compute all laplace coords by default (incl IO)
-    coords.extend(
-        inputs[config["modality"]].expand(
-            bids(
-                root=root,
-                datatype="coords",
-                dir="{dir}",
-                suffix="coords.nii.gz",
-                desc="{desc}",
-                space="{space}",
-                hemi="{hemi}",
-                label="{autotop}",
-                **inputs.subj_wildcards,
-            ),
-            desc="laplace",
-            dir=["AP", "PD", "IO"],
-            autotop=config["autotop_labels"],
-            hemi=config["hemi"],
-            space=crop_ref_spaces,
-            allow_missing=True,
-        )
-    )
-    coords.extend(
-        inputs[config["modality"]].expand(
-            bids(
-                root=root,
-                datatype="coords",
-                dir="{dir}",
-                suffix="coords.nii.gz",
-                desc="{desc}",
-                space="{space}",
-                hemi="{hemi}",
-                label="hipp",
-                **inputs.subj_wildcards,
-            ),
-            desc=[desc_io],
-            dir=["IO"],
-            hemi=config["hemi"],
-            space=crop_ref_spaces,
-            allow_missing=True,
-        )
-    )
-    return coords
-
-
 def get_final_anat():
     anat = []
     if "T1w" in ref_spaces or "T2w" in ref_spaces:
@@ -264,15 +212,14 @@ def get_final_output():
     subj_output = []
     subj_output.extend(get_final_spec())
     subj_output.extend(get_final_subfields())
-    subj_output.extend(get_final_coords())
     subj_output.extend(get_final_anat())
     subj_output.extend(get_final_qc())
     return subj_output
 
 
 if "corobl" in ref_spaces:
 
-    ruleorder: equivolume_coords > laplace_coords_hipp > laplace_coords_dentate > copy_coords_to_results
+    ruleorder: equivolume_coords > laplace_coords_dentate > copy_coords_to_results
 
     rule copy_coords_to_results:
         input:

hippunfold/workflow/rules/autotop.smk → hippunfold/workflow/rules/coords.smk

@@ -1,14 +1,6 @@
 import os
 
 
-def get_cmd_laplace_coords():
-    if config["skip_inject_template_labels"]:
-        cmd = "../scripts/laplace_coords.py"
-    else:
-        cmd = "../scripts/laplace_coords_withinit.py"
-    return cmd
-
-
 def get_labels_for_laplace(wildcards):
     if config["skip_inject_template_labels"]:
         seg = get_input_for_shape_inject(wildcards)
@@ -45,45 +37,6 @@ def get_inputs_laplace(wildcards):
     return files
 
 
-rule laplace_coords_hipp:
-    input:
-        unpack(get_inputs_laplace),
-    params:
-        gm_labels=lambda wildcards: config["laplace_labels"][wildcards.dir]["gm"],
-        src_labels=lambda wildcards: config["laplace_labels"][wildcards.dir]["src"],
-        sink_labels=lambda wildcards: config["laplace_labels"][wildcards.dir]["sink"],
-        convergence_threshold=1e-5,
-        max_iters=10000,
-    output:
-        coords=bids(
-            root=work,
-            datatype="coords",
-            dir="{dir}",
-            label="hipp",
-            suffix="coords.nii.gz",
-            desc="laplace",
-            space="corobl",
-            hemi="{hemi}",
-            **inputs.subj_wildcards
-        ),
-    group:
-        "subj"
-    resources:
-        time=30,
-    log:
-        bids(
-            root="logs",
-            **inputs.subj_wildcards,
-            dir="{dir}",
-            hemi="{hemi}",
-            suffix="laplace-hipp.txt"
-        ),
-    container:
-        config["singularity"]["autotop"]
-    script:
-        get_cmd_laplace_coords()
-
-
 rule laplace_coords_dentate:
     input:
         coords=bids(

hippunfold/workflow/rules/myelin_map.smk

@@ -37,38 +37,6 @@ rule divide_t1_by_t2:
         "c3d {input.t1} {input.t2} -divide -replace inf 1000 -inf -1000 NaN 0 -o {output}"
 
 
-rule create_ribbon:
-    input:
-        ap_coord=bids(
-            root=work,
-            datatype="coords",
-            dir="AP",
-            label="{autotop}",
-            suffix="coords.nii.gz",
-            space="{space}",
-            desc="laplace",
-            hemi="{hemi}",
-            **inputs.subj_wildcards
-        ),
-    output:
-        ribbon=bids(
-            root=work,
-            datatype="anat",
-            label="{autotop}",
-            suffix="mask.nii.gz",
-            space="{space}",
-            desc="ribbon",
-            hemi="{hemi}",
-            **inputs.subj_wildcards
-        ),
-    group:
-        "subj"
-    container:
-        config["singularity"]["autotop"]
-    shell:
-        "c3d {input} -binarize -o {output}"
-
-
 # sample on hipp & dg midthickness surfaces
 rule sample_myelin_map_surf:
     """ samples myelin map on surf using corobl space """
@@ -109,16 +77,6 @@ rule sample_myelin_map_surf:
             label="{autotop}",
             **inputs.subj_wildcards
         ),
-        ribbon=bids(
-            root=work,
-            datatype="anat",
-            label="{autotop}",
-            suffix="mask.nii.gz",
-            space="corobl",
-            desc="ribbon",
-            hemi="{hemi}",
-            **inputs.subj_wildcards
-        ),
     output:
         metric=bids(
             root=root,
@@ -134,4 +92,4 @@ rule sample_myelin_map_surf:
     container:
         config["singularity"]["autotop"]
     shell:
-        "wb_command -volume-to-surface-mapping {input.vol} {input.mid} {output.metric} -ribbon-constrained {input.outer} {input.inner} -volume-roi {input.ribbon}"
+        "wb_command -volume-to-surface-mapping {input.vol} {input.mid} {output.metric} -ribbon-constrained {input.outer} {input.inner}"