--- a
+++ b/gen_unfold_template/Snakefile
@@ -0,0 +1,269 @@
+from snakebids import bids
+import numpy as np
+from snakebids.utils import glob_wildcards
+
+
+#configfile: 'config_dentate.yml'
+
+nbins = config['num_bins']
+targetareas = config['targetareas']
+
+subjects, = glob_wildcards(os.path.join(config['hippunfold_dir'],config['in_cifti']))
+print(subjects)
+
+rule all_tplmesh:
+    input:
+        surf_gii = expand(config['out_final_root']+'/tpl-avg_space-unfold_den-{density}_{surfname}.surf.gii',
+                        density=config['density_to_area'].keys(),
+                        surfname=config['surf_to_io'].keys()),
+        surf_area = expand(config['out_final_root']+'/tpl-avg_space-unfold_den-{density}_surfarea.shape.gii',
+                        density=config['density_to_area'].keys())
+
+
+rule all_area:
+    input:
+        txt = expand(bids(root=config['out_root'],subject='{subject}',hemi='{hemi}',space='corobl',
+                suffix='{surfname}.surfarea.stats.alltargets.txt',nbins='{nbins}'),
+                    subject=subjects[0],
+                    hemi='R',
+                    surfname=config['surf_to_io'].keys(),
+                    nbins=nbins)
+
+rule all_example:
+    input:
+        gii = expand(bids(root='example_transformed',
+                        include_subject_dir=False,
+                        subject='{subject}',
+                        hemi='{hemi}',
+                        space='corobl',
+                        suffix='{surfname}.surf.gii',den='{density}'),
+                    subject='0001',
+                    hemi='R',
+                    density=config['density_to_area'].keys(),
+                    surfname=config['surf_to_io'].keys())
+ 
+
+rule all_subjsurf:
+    input:
+        gii = expand(bids(root=config['out_root'],subject='{subject}',hemi='{hemi}',space='corobl',
+                suffix='{surfname}.surfarea.stats.txt',targetarea='{targetarea}',nbins='{nbins}'),
+                    subject='0001',
+                    hemi='R',
+                    surfname=config['surf_to_io'].keys(),
+                    targetarea=targetareas,
+                    nbins=nbins)
+
+
+def get_average_surfarea_cifti_cmd(wildcards, input,output):
+    cmd = ['wb_command','-cifti-average',f'{output.cifti}','-exclude-outliers 2 2']
+    for cifti in input.cifti:
+        cmd.append(f'-cifti {cifti}')
+    return ' '.join(cmd)
+
+rule create_surfarea_cifti:
+    """this will soon be added to main hippunfold worflow, but using this until then"""
+    input:
+        left_metric = lambda wildcards: os.path.join(config['hippunfold_dir'],config['in_gifti']).format(
+                    hemi='L',
+                    subject=wildcards.subject),
+        right_metric = lambda wildcards: os.path.join(config['hippunfold_dir'],config['in_gifti']).format(
+                    hemi='R',
+                    subject=wildcards.subject),
+ 
+
+    output:
+        cifti = bids(root=config['out_root'],subject='{subject}',suffix='surfarea.dscalar.nii')
+    shell:
+        'wb_command  -cifti-create-dense-scalar {output}'
+        ' -left-metric {input.left_metric}'
+        ' -right-metric {input.right_metric}'
+        
+
+    
+rule average_surfarea_cifti:
+    input: 
+        cifti = expand(bids(root=config['out_root'],subject='{subject}',suffix='surfarea.dscalar.nii'),
+                    subject=subjects)
+    params:
+        cmd = get_average_surfarea_cifti_cmd
+    output:
+        cifti = config['out_root']+'/avg-HCPUR100_surfarea.dscalar.nii'
+    shell:
+        '{params.cmd}'
+
+rule separate_avg_left_right:
+    input:
+        cifti = config['out_root']+'/avg-HCPUR100_surfarea.dscalar.nii'
+    output:
+        gifti_left = config['out_root']+'/avg-HCPUR100_hemi-L_surfarea.shape.gii',
+        gifti_right = config['out_root']+'/avg-HCPUR100_hemi-R_surfarea.shape.gii'
+    shell:
+        'wb_command -cifti-separate {input} COLUMN -metric CORTEX_LEFT {output.gifti_left} -metric CORTEX_RIGHT {output.gifti_right}'
+
+rule avg_left_and_right:
+    """ NOT USING THIS, SINCE LEFT NEEDS TO BE FLIPPED in both dims"""
+    input:
+        left = config['out_root']+'/avg-HCPUR100_hemi-L_surfarea.shape.gii',
+        right = config['out_root']+'/avg-HCPUR100_hemi-R_surfarea.shape.gii'
+    output: config['out_root']+'/avg-HCPUR100_hemi-avg_surfarea.shape.gii' 
+    shell:
+        'wb_command -metric-math "0.5*(left+right)" {output} -var left {input.left}  -var right {input.right}'
+
+
+
+rule gen_adaptive_mesh:
+    input:
+        surfarea_gii = config['out_root']+'/avg-HCPUR100_hemi-R_surfarea.shape.gii',
+    params:
+        targetarea = lambda wildcards: float(wildcards.targetarea),
+        nbins = lambda wildcards: int(wildcards.nbins),
+    output:
+        points_csv = config['out_root']+'/tpl-avg_targetarea-{targetarea}_nbins-{nbins}_points.csv',
+        triangles_csv = config['out_root']+'/tpl-avg_targetarea-{targetarea}_nbins-{nbins}_triangles.csv',
+        grid_png = config['out_root']+'/tpl-avg_targetarea-{targetarea}_nbins-{nbins}_grid.png'
+    script: 'gen_adaptivesurfs.py'
+       
+
+rule gen_uniform_mesh:
+    output:
+        triangles_csv = config['out_root']+'/tpl-avg_targetarea-uniform_nbins-uniform_triangles.csv',
+        points_csv = config['out_root']+'/tpl-avg_targetarea-uniform_nbins-uniform_points.csv',
+        grid_png = config['out_root']+'/tpl-avg_desc-uniform_grid.png'
+    script: 'gen_uniformgrid.py'
+
+
+rule csv_to_gii:
+    """ ideally can replace this with writing directly from nibabel.."""
+    input:
+        points_csv = config['out_root']+'/tpl-avg_targetarea-{targetarea}_nbins-{nbins}_points.csv',
+        triangles_csv = config['out_root']+'/tpl-avg_targetarea-{targetarea}_nbins-{nbins}_triangles.csv',
+    params:
+        io_value = lambda wildcards: config['surf_to_io'][wildcards.surfname]
+    output:
+        surf_gii = config['out_root']+'/tpl-avg_targetarea-{targetarea}_nbins-{nbins}_{surfname}.surf.gii',
+    shell:
+        "matlab -batch \"convert_csv_to_gifti('{input.triangles_csv}','{input.points_csv}','{output.surf_gii}',{params.io_value})\""
+
+#lookup tables for structure:
+hemi_to_structure = {'L': 'CORTEX_LEFT', 'Lflip': 'CORTEX_LEFT', 'R': 'CORTEX_RIGHT'}
+surf_to_secondary_type = {'midthickness': 'MIDTHICKNESS', 'inner': 'PIAL', 'outer': 'GRAY_WHITE'}
+
+
+rule surf_to_final_name:
+    input:
+        surf_gii = lambda wildcards: config['out_root']+'/tpl-avg_targetarea-{targetarea}_nbins-{nbins}_{surfname}.surf.gii'.format(targetarea=config['density_to_area'][wildcards.density], nbins=config['num_bins'],surfname=wildcards.surfname)
+    output:
+        surf_gii = config['out_final_root']+'/tpl-avg_space-unfold_den-{density}_{surfname}.surf.gii',
+    shell:
+        'cp {input} {output}'
+
+
+#calc surface area, hippunfold will need this for the templates to calculate gyrification
+rule template_surfarea:
+    input:
+        surf_gii = config['out_final_root']+'/tpl-avg_space-unfold_den-{density}_midthickness.surf.gii',
+    output:
+        surf_gii = config['out_final_root']+'/tpl-avg_space-unfold_den-{density}_surfarea.shape.gii',
+    shell:
+        'wb_command -surface-vertex-areas {input} {output}'
+
+
+rule examplesurf_to_final_name:
+    input:
+        gii = lambda wildcards: bids(root=config['out_root'],
+                    subject=f'{wildcards.subject}',
+                    hemi=f'{wildcards.hemi}',
+                    space='corobl',
+                    suffix=f'{wildcards.surfname}.surf.gii',
+                    targetarea=config['density_to_area'][wildcards.density],
+                    nbins=config['num_bins']),
+    output:
+        gii = bids(root='example_transformed',include_subject_dir=False,subject='{subject}',hemi='{hemi,R|Lflip}',space='corobl',
+                suffix='{surfname}.surf.gii',den='{density}'),
+    shell:
+        'cp {input} {output}'
+#
+#rule warp_gii_unfoldtemplate2unfold: 
+#    """warp from template space to subj unfolded"""
+#    input: 
+#        warp ='../work/sub-{subject}/sub-{subject}_hemi-{hemi}_space-corobl_desc-cropped_modality-T2w_autotop/Warp_unfoldtemplate2unfold.nii',
+#        gii = config['out_root']+'/tpl-avg_targetarea-{targetarea}_nbins-{nbins}_{surfname}.surf.gii',
+#
+#    params:
+#        structure_type = lambda wildcards: hemi_to_structure[wildcards.hemi],
+#        secondary_type = lambda wildcards: surf_to_secondary_type[wildcards.surfname],
+#        surface_type = 'FLAT'
+#    output:
+#        gii = bids(root=config['out_root'],subject='{subject}',hemi='{hemi,R|Lflip}',space='unfolded',
+#                suffix='{surfname}.surf.gii',targetarea='{targetarea}',nbins='{nbins}'),
+# 
+#    shell:
+#        'wb_command -surface-apply-warpfield {input.gii} {input.warp} {output.gii} && '
+#        'wb_command -set-structure {output.gii} {params.structure_type} -surface-type {params.surface_type}'
+#            ' -surface-secondary-type {params.secondary_type}'
+
+#subj unfolded surf might have a few vertices outside the bounding box.. this constrains all the vertices to the warp bounding box
+#rule constrain_surf_to_bbox:
+#    input:
+#        gii = bids(root=config['out_root'],subject='{subject}',hemi='{hemi,R|Lflip}',space='unfolded',
+#                suffix='{surfname}.surf.gii',targetarea='{targetarea}',nbins='{nbins}'),
+#        ref_nii = '../work/sub-{subject}/sub-{subject}_space-unfold_refvol.nii.gz'
+#    output:
+#        gii = bids(root=config['out_root'],subject='{subject}',hemi='{hemi,R|Lflip}',space='unfolded',
+#                suffix='{surfname}.surf.gii',desc='constrainbbox',targetarea='{targetarea}',nbins='{nbins}'),
+#    group: 'subj'
+#    script: 'constrain_surf_to_bbox.py'
+
+#warp from subj unfolded to corobl
+rule warp_gii_unfold2native: 
+    input: 
+        gii = config['out_root']+'/tpl-avg_targetarea-{targetarea}_nbins-{nbins}_{surfname}.surf.gii',
+        warp = config['in_warp_unfold2corobl']
+    params:
+        structure_type = lambda wildcards: hemi_to_structure[wildcards.hemi],
+        secondary_type = lambda wildcards: surf_to_secondary_type[wildcards.surfname],
+        surface_type = 'ANATOMICAL'
+    output:
+        gii = bids(root=config['out_root'],subject='{subject}',hemi='{hemi,R|Lflip}',space='corobl',
+                suffix='{surfname}.surf.gii',targetarea='{targetarea}',nbins='{nbins}'),
+    shell:
+        'wb_command -surface-apply-warpfield {input.gii} {input.warp} {output.gii} && '
+        'wb_command -set-structure {output.gii} {params.structure_type} -surface-type {params.surface_type}'
+            ' -surface-secondary-type {params.secondary_type}'
+
+
+
+rule calc_surfarea:
+    input:
+        gii = bids(root=config['out_root'],subject='{subject}',hemi='{hemi,R|Lflip}',space='corobl',
+                suffix='{surfname}.surf.gii',targetarea='{targetarea}',nbins='{nbins}'),
+    output:
+        gii = bids(root=config['out_root'],subject='{subject}',hemi='{hemi,R|Lflip}',space='corobl',
+                suffix='{surfname}.surfarea.shape.gii',targetarea='{targetarea}',nbins='{nbins}'),
+    shell: 'wb_command -surface-vertex-areas {input} {output}'
+
+rule write_surfarea_stats:
+    input:
+        gii = bids(root=config['out_root'],subject='{subject}',hemi='{hemi,R|Lflip}',space='corobl',
+                suffix='{surfname}.surfarea.shape.gii',targetarea='{targetarea}',nbins='{nbins}'),
+    output:  
+        txt = bids(root=config['out_root'],subject='{subject}',hemi='{hemi,R|Lflip}',space='corobl',
+                suffix='{surfname}.surfarea.stats.txt',targetarea='{targetarea}',nbins='{nbins}'),
+    shell:
+        'wb_command -file-information {input} | tail -n 3 | head -n 1 > {output}'
+
+rule combine_surfarea_stats:
+    input:
+        txt = expand(bids(root=config['out_root'],subject='{subject}',hemi='{hemi,R|Lflip}',space='corobl',
+                suffix='{surfname}.surfarea.stats.txt',targetarea='{targetarea}',nbins='{nbins}'),
+                    targetarea=targetareas,allow_missing=True)
+    output:
+        txt = bids(root=config['out_root'],subject='{subject}',hemi='{hemi,R|Lflip}',space='corobl',
+                suffix='{surfname}.surfarea.stats.alltargets.txt',nbins='{nbins}'),
+    run:
+        shell('echo "Map   Minimum   Maximum    Mean   Sample Dev   % Positive   % Negative   Inf/NaN   Map Name" >> {output}')
+        for txt in input.txt:
+            shell('echo {txt} >> {output}')
+            shell('cat {txt} >> {output}')
+
+