WIP - Laynii native

hippunfold/config/snakebids.yml

@@ -231,10 +231,10 @@ parse_args:
   --laminar_coords_method:
     choices:
       - 'laplace'
-      - 'equivolume'
-    default:
-      - 'equivolume'
-    help: 'Method to use for laminar coordinates. Equivolume uses equivolumetric layering from Waehnert et al 2014 (Nighres implementation). (default: %(default)s)' 
+      - 'equivol'
+      - 'equidist'
+    default: 'equivol'
+    help: 'Method to use for laminar coordinates. Equivol and equidist are from the laynii LN2_LAYERS tool. (default: %(default)s)' 
 
   --autotop_labels:
     choices:

hippunfold/workflow/rules/autotop.smk

@@ -126,98 +126,152 @@ rule laplace_coords_dentate:
     shell:
         "cp {input} {output}"
 
-
-rule prep_equivolume_coords:
+rule morphclose_dg:
     input:
-        get_labels_for_laplace,
-    params:
-        src_labels=lambda wildcards: config["laplace_labels"][wildcards.dir]["src"],
+        dseg_tissue=get_labels_for_laplace,
     output:
-        outerbin=bids(
+        dseg_tissue=bids(
             root=work,
-            datatype="coords",
-            dir="{dir}",
-            desc="all",
-            suffix="mask.nii.gz",
+            datatype="anat",
+            **inputs.subj_wildcards,
+            suffix="dseg.nii.gz",
+            desc="closeDG",
             space="corobl",
-            hemi="{hemi}",
-            **inputs.subj_wildcards
+            hemi="{hemi}"
         ),
-        innerbin=bids(
+    group: 'subj'
+    container:
+        config["singularity"]["autotop"]
+    shell: 
+        "c3d {input} -as DSEG -retain-labels 8 -binarize -dilate 1 3x3x3vox -erode 1 3x3x3vox -scale 100 -push DSEG -max -replace 100 8 -o {output}"
+
+rule prep_dseg_for_laynii_hipp:
+    input:
+        dseg_tissue=bids(
             root=work,
-            datatype="coords",
-            dir="{dir}",
-            desc="SRLM",
-            suffix="mask.nii.gz",
+            datatype="anat",
+            **inputs.subj_wildcards,
+            suffix="dseg.nii.gz",
+            desc="closeDG",
             space="corobl",
-            hemi="{hemi}",
-            **inputs.subj_wildcards
+            hemi="{hemi}"
         ),
-    log:
-        bids(
-            root="logs",
+    params:
+        gm_labels=lambda wildcards: " ".join(
+            [str(lbl) for lbl in config["laplace_labels"][wildcards.dir]["gm_noDG"]]
+        ),
+        src_labels=lambda wildcards: " ".join(
+            [str(lbl) for lbl in config["laplace_labels"][wildcards.dir]["src"]]
+        ),
+        sink_labels=lambda wildcards: " ".join(
+            [str(lbl) for lbl in config["laplace_labels"][wildcards.dir]["sink"]]
+        ),
+    output:
+        dseg_rim=bids(
+            root=work,
+            datatype="anat",
             **inputs.subj_wildcards,
-            dir="{dir}",
-            hemi="{hemi}",
-            suffix="binarize.txt"
+            suffix="dseg.nii.gz",
+            dir="{dir,IO}",
+            desc="laynii",
+            label="{autotop,hipp}",
+            space="corobl",
+            hemi="{hemi}"
         ),
+    container:
+        config["singularity"]["autotop"]
     group:
         "subj"
+    shell:
+        "c3d -background -1 {input} -as DSEG -retain-labels {params.gm_labels} -binarize -scale 3 -popas GM -push DSEG -retain-labels {params.src_labels} -binarize -scale 2 -popas WM -push DSEG -retain-labels {params.sink_labels} -binarize -scale 1 -popas PIAL -push GM -push WM -add -push PIAL -add -o {output}"
+
+
+rule prep_dseg_for_laynii_dentate:
+    input:
+        dseg_tissue=bids(
+            root=work,
+            datatype="anat",
+            **inputs.subj_wildcards,
+            suffix="dseg.nii.gz",
+            desc="closeDG",
+            space="corobl",
+            hemi="{hemi}"
+        ),
+    params:
+        gm_labels=lambda wildcards: " ".join(
+            [str(lbl) for lbl in [8]]
+        ),
+        src_labels=lambda wildcards: " ".join(
+            [str(lbl) for lbl in [2,4,7,0] ]
+        ),
+        sink_labels=lambda wildcards: " ".join(
+            [str(lbl) for lbl in [1]]
+        ),
+    output:
+        dseg_rim=bids(
+            root=work,
+            datatype="anat",
+            **inputs.subj_wildcards,
+            suffix="dseg.nii.gz",
+            dir="{dir,IO}",
+            desc="laynii",
+            label="{autotop,dentate}",
+            space="corobl",
+            hemi="{hemi}"
+        ),
     container:
         config["singularity"]["autotop"]
-    script:
-        "../scripts/prep_equivolume_coords.py"
+    group:
+        "subj"
+    shell:
+        "c3d -background -1 {input} -as DSEG -retain-labels {params.gm_labels} -binarize -scale 3 -popas GM -push DSEG -retain-labels {params.src_labels} -binarize -scale 2 -popas WM -push DSEG -retain-labels {params.sink_labels} -binarize -scale 1 -popas PIAL -push GM -push WM -add -push PIAL -add -o {output}"
 
 
-rule equivolume_coords:
+
+rule laynii_layers:
     input:
-        outerbin=bids(
+        dseg_rim=bids(
             root=work,
-            datatype="coords",
+            datatype="anat",
+            **inputs.subj_wildcards,
+            suffix="dseg.nii.gz",
             dir="{dir}",
-            desc="all",
-            suffix="mask.nii.gz",
+            desc="laynii",
+            label="{autotop}",
             space="corobl",
-            hemi="{hemi}",
-            **inputs.subj_wildcards
+            hemi="{hemi}"
         ),
-        innerbin=bids(
+    output:
+        equivol=bids(
             root=work,
             datatype="coords",
-            dir="{dir}",
-            desc="SRLM",
-            suffix="mask.nii.gz",
+            dir="{dir,IO}",
+            label="{autotop}",
+            suffix="coords.nii.gz",
+            desc="equivol",
             space="corobl",
             hemi="{hemi}",
             **inputs.subj_wildcards
         ),
-    params:
-        script=os.path.join(workflow.basedir, "scripts/equivolume_coords.py"),
-    output:
-        coords=bids(
+        equidist=bids(
             root=work,
             datatype="coords",
-            dir="{dir}",
-            label="hipp",
+            dir="{dir,IO}",
+            label="{autotop}",
             suffix="coords.nii.gz",
-            desc="equivol",
+            desc="equidist",
             space="corobl",
             hemi="{hemi}",
             **inputs.subj_wildcards
         ),
+    shadow:
+        "minimal"
+    container:
+        config["singularity"]["autotop"] 
     group:
         "subj"
-    resources:
-        time=30,
-    log:
-        bids(
-            root="logs",
-            **inputs.subj_wildcards,
-            dir="{dir}",
-            hemi="{hemi}",
-            suffix="equivolume.txt"
-        ),
-    container:
-        config["singularity"]["autotop"]
     shell:
-        "python {params.script} {resources.tmpdir} {input.innerbin} {input.outerbin} {output.coords} &> {log}"
+        "cp {input} dseg.nii.gz && "
+        "LN2_LAYERS  -rim dseg.nii.gz -equivol && "
+        "cp dseg_metric_equidist.nii.gz {output.equidist} && "
+        "cp dseg_metric_equivol.nii.gz {output.equivol}"

hippunfold/workflow/rules/common.smk

@@ -123,11 +123,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(
@@ -144,7 +139,7 @@ def get_final_coords():
                 **inputs.subj_wildcards,
             ),
             desc="laplace",
-            dir=["AP", "PD", "IO"],
+            dir=["AP", "PD"],
             autotop=config["autotop_labels"],
             hemi=config["hemi"],
             space=crop_ref_spaces,
@@ -164,7 +159,7 @@ def get_final_coords():
                 label="hipp",
                 **inputs.subj_wildcards,
             ),
-            desc=[desc_io],
+            desc=config['laminar_coords_method'],
             dir=["IO"],
             hemi=config["hemi"],
             space=crop_ref_spaces,

hippunfold/workflow/rules/native_surf.smk

@@ -30,7 +30,7 @@ nan_labels = {
 }
 
 desc_io = {
-    "hipp": "equivol" if "equivolume" in config["laminar_coords_method"] else "laplace",
+    "hipp": config["laminar_coords_method"],
     "dentate": "laplace",
 }
 

hippunfold/workflow/rules/warps.smk

@@ -286,31 +286,18 @@ rule create_unfold_coord_map:
 
 
 def get_laminar_coords(wildcards):
-    if "laplace" in config["laminar_coords_method"]:
-        coords_io = bids(
+
+    return bids(
             root=work,
             datatype="coords",
             dir="IO",
             label="hipp",
             suffix="coords.nii.gz",
-            desc="laplace",
-            space="corobl",
-            hemi="{hemi}",
-            **inputs.subj_wildcards
-        )
-    elif "equivolume" in config["laminar_coords_method"]:
-        coords_io = bids(
-            root=work,
-            datatype="coords",
-            dir="IO",
-            label="hipp",
-            suffix="coords.nii.gz",
-            desc="equivol",
+            desc=config["laminar_coords_method"],
             space="corobl",
             hemi="{hemi}",
             **inputs.subj_wildcards
         )
-    return coords_io
 
 
 rule create_warps_hipp: