Updates for shape injection and surf processing of injected surfaces (#…

…374)

- high-res ref space for hipp, dentate from shape inject (hipp default 0.3, dentate 0.1)
  - required adding a label wildcard to the dsegtissue and upsampled_ref
files
- fixes to get injected dentate/hipp working well (enough) 
- updated template upenn tissue dseg to include DG (only local
currently) 
- updates to greedy registration for shape inject
   - downsampling for greedy template shape reg 
   - enforce det=1 for template shape reg moments (was flipping)
- dilate dentate pd src/sink (not used yet)
- adjust decimation
- use workaround for transforming surfaces (need to look into this bug - perhaps provide an example to HCP devs)
- use two different, new,  sdt threshold approaches
  - percentile-based for AP, and firstminima for PD
  - NOTE: still not perfect and should really implement #372 to make more robust

hippunfold/config/snakebids.yml

@@ -122,12 +122,6 @@ parse_args:
   --inject_template:
     choices:
       - 'upenn'
-      - 'dHCP'
-      - 'MBMv2'
-      - 'MBMv3'
-      - 'CIVM'
-      - 'ABAv3'
-      - 'bigbrain'
     default: 'upenn'
     help: 'Set the template to use for shape injection. (default: %(default)s)'
 
@@ -189,6 +183,16 @@ parse_args:
     help: 'Resample the manual segmentation (for --modality dsegtissue), keeping the bounding box the same, but changing the number of voxels in the image. This can be specified as voxels ("300x300x300"), or a percentage ("20%%") or percentage in each direction ("20x20x40%%"). Note: the segmentation will always be subsequently cropped around the segmentation, with a padding of 5 voxels.  (default: %(default)s)'
     default: null
 
+  --laminar_coords_res_hipp:
+    help: 'Set the resolution that the laminar (IO) coords will be computed with for hippocampus. This is implemented by resampling the segmentation to this resolution prior to generating the coords. (default: %(default)s)'
+    default: '0.3x0.3x0.3mm'
+    dest: "laminar_coords_res.hipp"
+
+  --laminar_coords_res_dentate:
+    help: 'Set the resolution that the laminar (IO) coords will be computed with for the dentate gyrus. This is implemented by resampling the segmentation to this resolution prior to generating the coords. (default: %(default)s)'
+    default: '0.1x0.1x0.1mm'
+    dest: "laminar_coords_res.dentate"
+
 
   --atlas:
     choices:
@@ -426,10 +430,18 @@ template_files:
   upenn:
     T1w: tpl-upenn_desc-hipptissue_dseg.nii.gz
     dseg: tpl-upenn_desc-hipptissue_dseg.nii.gz
+    dseg_dentate: tpl-upenn_desc-hippdentate_dseg.nii.gz # hack: separate files for DG
     coords: tpl-upenn_dir-{dir}_label-{label}_coords.nii.gz
     hemi_wildcards:
       - R
 
+  bbhist:
+    T1w: sub-bbhist_hemi-{hemi}_100um.nii.gz
+    dseg: sub-bbhist_hemi-{hemi}_100um_dseg.nii.gz # aka dsegtissue
+    hemi_wildcards:
+      - L
+      - R
+
   ABAv3:
     T1w: tpl-ABAv3_T1w.nii.gz
     xfm_corobl: tpl-ABAv3_from-native_to-corobl_type-itk_affine.txt
@@ -475,15 +487,15 @@ atlas_files:
     lut: desc-subfields_dseg.tsv
     surf_gii: sub-bbhist_hemi-{hemi}_space-unfold_label-{label}_{surf_type}.surf.gii
     metric_gii: sub-bbhist_hemi-{hemi}_space-corobl_label-{label}_{metric}.shape.gii
-    label_gii: sub-bbhist_hemi-{hemi}_space-corobl_label-{label}_desc-manual_subfields.label.gii
+    label_gii: sub-bbhist_hemi-{hemi}_space-corobl_label-{label}_subfields.label.gii
     label_wildcards:
      - hipp
+     - dentate
     hemi_wildcards:
      - L
      - R
     metric_wildcards:
      - thickness
-     - gyrification
      - curvature
 
   bigbrain:
@@ -570,6 +582,7 @@ resource_urls:
     CIVM: 'files.ca-1.osf.io/v1/resources/v8acf/providers/osfstorage/65395bf62827451220b86e24/?zip='
     ABAv3: 'files.ca-1.osf.io/v1/resources/v8acf/providers/osfstorage/6668855b6b6c8e2cc704ca97/?zip='
     bigbrain: 'files.ca-1.osf.io/v1/resources/v8acf/providers/osfstorage/666b1bc765e1de5972893e59/?zip='
+    bbhist: 'placeholder.url'
 #to get hash, see https://github.com/CenterForOpenScience/osf.io/issues/8256#issuecomment-379833911
 
 
@@ -710,6 +723,14 @@ laplace_labels:
         - 9
         - 10
         - 0
+
+
+tight_crop_labels: #defines what labels to use for tight cropping (for upsampled ref)
+  hipp:
+    - 1
+  dentate:
+    - 8
+
 output_spaces:
   - native
 
@@ -737,3 +758,8 @@ workdir: null
 use_template_seg: False
 template_seg_smoothing_factor: 2
 resample_dsegtissue: null
+laminar_coords_res:
+  hipp: 0.3x0.3x0.3mm
+  dentate: 0.1x0.1x0.1mm
+resample_dseg_for_templatereg: 0.3x0.3x0.3mm
+

hippunfold/workflow/rules/common.smk

@@ -83,6 +83,7 @@ def get_final_subfields():
             space="{space}",
             hemi="{hemi}",
             atlas="{atlas}",
+            label="hipp",
             **inputs.subj_wildcards,
         ),
         hemi=config["hemi"],
@@ -331,3 +332,26 @@ def get_create_template_output():
             )
         )
     return files
+
+
+def get_input_for_shape_inject(wildcards):
+    if config["modality"] == "dsegtissue":
+        seg = bids(
+            root=work,
+            datatype="anat",
+            **inputs.subj_wildcards,
+            suffix="dseg.nii.gz",
+            space="corobl",
+            hemi="{hemi}",
+        ).format(**wildcards)
+    else:
+        seg = bids(
+            root=work,
+            datatype="anat",
+            **inputs.subj_wildcards,
+            suffix="dseg.nii.gz",
+            desc="nnunet",
+            space="corobl",
+            hemi="{hemi}",
+        ).format(**wildcards)
+    return seg

hippunfold/workflow/rules/coords.smk

@@ -13,6 +13,7 @@ def get_labels_for_laplace(wildcards):
             desc="postproc",
             space="corobl",
             hemi="{hemi}",
+            label="{label}",
         ).format(**wildcards)
     return seg
 
@@ -185,6 +186,36 @@ rule get_nan_mask:
         "c3d {input} -background -1 -retain-labels {params} -binarize {output}"
 
 
+rule create_upsampled_coords_ref:
+    input:
+        seg=get_input_for_shape_inject,
+    params:
+        tight_crop_labels=lambda wildcards: config["tight_crop_labels"][wildcards.label],
+        resample_res=lambda wildcards: config["laminar_coords_res"][wildcards.label],
+        trim_padding="3mm",
+    output:
+        upsampled_ref=bids(
+            root=work,
+            datatype="anat",
+            **inputs.subj_wildcards,
+            suffix="ref.nii.gz",
+            desc="resampled",
+            space="corobl",
+            label="{label}",
+            hemi="{hemi}",
+        ),
+    container:
+        config["singularity"]["autotop"]
+    conda:
+        "../envs/c3d.yaml"
+    group:
+        "subj"
+    container:
+        config["singularity"]["autotop"]
+    shell:
+        "c3d {input} -retain-labels {params.tight_crop_labels} -trim {params.trim_padding} -resample-mm {params.resample_res} -o {output}"
+
+
 rule prep_dseg_for_laynii:
     input:
         dseg_tissue=get_labels_for_laplace,

hippunfold/workflow/rules/download.smk

@@ -87,6 +87,43 @@ rule import_template_dseg:
         "{params.copy_or_flip_cmd} {output.template_seg}"
 
 
+rule import_template_dseg_dentate:
+    input:
+        template_dir=Path(download_dir) / "template" / config["inject_template"],
+    params:
+        template_seg=lambda wildcards: Path(download_dir)
+        / "template"
+        / config["inject_template"]
+        / config["template_files"][config["inject_template"]]["dseg_dentate"].format(
+            **wildcards
+        ),
+        copy_or_flip_cmd=lambda wildcards: copy_or_flip(
+            wildcards,
+            Path(download_dir)
+            / "template"
+            / config["inject_template"]
+            / config["template_files"][config["inject_template"]][
+                "dseg_dentate"
+            ].format(**wildcards),
+        ),
+    output:
+        template_seg=bids(
+            root=work,
+            datatype="anat",
+            space="template",
+            **inputs.subj_wildcards,
+            desc="dentatetissue",
+            hemi="{hemi}",
+            suffix="dseg.nii.gz",
+        ),
+    group:
+        "subj"
+    container:
+        config["singularity"]["autotop"]
+    shell:
+        "{params.copy_or_flip_cmd} {output.template_seg}"
+
+
 rule import_template_coords:
     input:
         template_dir=Path(download_dir) / "template" / config["inject_template"],

hippunfold/workflow/rules/native_surf.smk

@@ -5,6 +5,7 @@
 
 surf_thresholds = {"inner": 0, "outer": 1, "midthickness": 0.5}
 
+
 unfoldreg_method = "greedy"  # choices: ["greedy","SyN"]
 
 unfoldreg_padding = "64x64x0vox"
@@ -66,7 +67,7 @@ rule gen_native_mesh:
     params:
         threshold=lambda wildcards: surf_thresholds[wildcards.surfname],
         decimate_opts={
-            "reduction": 0.9,
+            "reduction": 0.7,
             "feature_angle": 25,
             "preserve_topology": True,
         },
@@ -292,9 +293,12 @@ rule laplace_beltrami:
             **inputs.subj_wildcards,
         ),
     params:
-        min_dist_threshold=0.3,
-        max_dist_threshold=1,
-        min_terminal_vertices=5,
+        min_dist_percentile=1,
+        max_dist_percentile=10,
+        min_terminal_vertices=lambda wildcards: 5 if wildcards.dir == "AP" else 100,  #TODO, instead of # of vertices, we should compute the total length of the segment
+        threshold_method=lambda wildcards: (
+            "percentile" if wildcards.dir == "AP" else "firstminima"
+        ),
     output:
         coords=bids(
             root=work,
@@ -692,12 +696,15 @@ rule warp_midthickness_to_inout:
         config["singularity"]["autotop"]
     conda:
         "../envs/workbench.yaml"
+    shadow:
+        "minimal"
     group:
         "subj"
     shell:
-        "wb_command -surface-apply-warpfield {input.surf_gii} {input.warp} {output.surf_gii} && "
-        "wb_command -set-structure {output.surf_gii} {params.structure_type} -surface-type {params.surface_type}"
-        " -surface-secondary-type {params.secondary_type}"
+        "wb_command -volume-to-surface-mapping {input.warp} {input.surf_gii} warp.shape.gii -trilinear && "
+        "wb_command -surface-coordinates-to-metric {input.surf_gii} coords.shape.gii && "
+        "wb_command -metric-math 'COORDS + WARP' warpedcoords.shape.gii -var COORDS coords.shape.gii -var WARP warp.shape.gii && "
+        "wb_command -surface-set-coordinates  {input.surf_gii} warpedcoords.shape.gii {output.surf_gii}"
 
 
 # --- affine transforming anatomical surfaces from corobl to other (T1w, T2w) spaces

hippunfold/workflow/rules/qc.smk

@@ -51,6 +51,7 @@ rule get_subfield_vols_subj:
                 space="{crop_ref_spaces}",
                 desc="subfields",
                 atlas="{atlas}",
+                label="hipp",
                 suffix="dseg.nii.gz",
             ),
             hemi=config["hemi"],
@@ -173,6 +174,7 @@ rule qc_subfield:
             space="{space}",
             hemi="{hemi}",
             atlas="{atlas}",
+            label="hipp",
             **inputs.subj_wildcards,
         ),
     output:

hippunfold/workflow/rules/resample_final_to_crop_native.smk

@@ -11,6 +11,7 @@ rule create_native_crop_ref:
                 space="{native_modality}",
                 hemi="{hemi}",
                 atlas="{atlas}",
+                label="hipp",
                 **inputs.subj_wildcards,
             ),
             atlas=config["atlas"],
@@ -98,6 +99,7 @@ rule resample_postproc_native_crop:
             desc="postproc",
             space="corobl",
             hemi="{hemi}",
+            label=config["autotop_labels"][-1],
         ),
         xfm=bids(
             root=work,
@@ -148,6 +150,7 @@ rule resample_subfields_native_crop:
             space="corobl",
             hemi="{hemi}",
             atlas="{atlas}",
+            label="{label}",
             **inputs.subj_wildcards,
         ),
         xfm=bids(
@@ -177,6 +180,7 @@ rule resample_subfields_native_crop:
             space="crop{native_modality}",
             hemi="{hemi}",
             atlas="{atlas}",
+            label="{label,hipp}",
             **inputs.subj_wildcards,
         ),
     container: