feat: swc to graph or surface and tree morphological precision, recal…

…l, f1, iou via `recut truth.swc --test test.swc`. Graph to surface via `python3 scripts/visualize-graph.py neuron.graph

fix: close #30

flake.nix

@@ -4,7 +4,7 @@
     # pick latest commit from stable branch and test it so no suprises
     nixpkgs.url = "github:NixOS/nixpkgs/d53978239b265066804a45b7607b010b9cb4c50c";
     openvdb.url = "github:UCLA-VAST/openvdb?ref=feat/reachable-resurfacing-nix";
-    gel.url = "github:UCLA-VAST/GEL";
+    gel.url = "github:UCLA-VAST/GEL?ref=experimental2";
 
     # pin nix package manager versions to exact match to recut
     openvdb.inputs.nixpkgs.follows = "nixpkgs";

scripts/recut_interface.py

@@ -0,0 +1,23 @@
+import subprocess
+
+def valid(voxel_size):
+    return type(voxel_size) == list and len(voxel_size) == 3
+
+def call_recut(**kwargs):
+    if 'inferenced_path' in kwargs:
+        run_dir = kwargs['inferenced_path']
+    else:
+        # whitelist certain arguments to pass directly to recut
+        include = ['min_radius', 'max_radius', 'open_steps', 'close_steps', 'fg_percent',
+                   'preserve_topology', 'seed_action']
+        args = "".join([f"--{k} {v} ".replace('_', '-') for k, v in kwargs.items() if k in include])
+        cmd = f"/home/kdmarrett/recut/result/bin/recut {kwargs['image']} --seeds {kwargs['image']}/somas {args} --voxel-size {kwargs['voxel_size_x']} {kwargs['voxel_size_y']} " \
+              f"{kwargs['voxel_size_z']}"
+        print('  ')
+        print(cmd)
+
+        output = subprocess.check_output(cmd.split()).strip().decode().split('\n')
+        run_dir = [v.split()[-1] for v in output if "written to:" in v][0]
+        kwargs['inferenced_path'] = f"{run_dir}/seeds"
+    return run_dir;
+

scripts/recut_qa.py

@@ -1,16 +1,14 @@
-import subprocess
 import argparse
 import re
 # import pandas as pd
 from test_precision_recall import precision_recall
 from plots import get_hash
 from TMD_classify import filter_dir_by_model
-
+from recut_interface import call_recut
 
 def parse_range(string):
     if string.contains('-'):
         m = re.match(r'(\d+)(?:-(\d+))?$', string)
-        # ^ (or use .split('-'). anyway you like.)
         if not m:
             raise argparse.ArgumentTypeError(
                 "'" + string + "' is not a range of number. Expected forms like '0-5' or '2'.")
@@ -21,25 +19,9 @@ def parse_range(string):
         return list(int(string))
 
 
-def call_recut(**kwargs):
-    if kwargs['inferenced_path'] is None:
-        # whitelist certain arguments to pass directly to recut
-        include = ['min_radius', 'max_radius', 'open_steps', 'close_steps', 'fg_percent']
-        args = "".join([f"--{k} {v} ".replace('_', '-') for k, v in kwargs.items() if k in include if v])
-        cmd = f"recut {kwargs['image']} {args} --voxel-size {kwargs['voxel_size_x']} {kwargs['voxel_size_y']} " \
-              f"{kwargs['voxel_size_z']}"
-        print(cmd)
-
-        output = subprocess.check_output(cmd.split()).strip().decode().split('\n')
-        run_dir = [v.split()[-1] for v in output if "written to:" in v][0]
-        kwargs['inferenced_path'] = f"{run_dir}/seeds"
-    else:
-        run_dir = kwargs['inferenced_path']
-
+def persistence_precision_recall(run_dir):
     #precision_recall(**kwargs)
 
-    git_hash = get_hash()
-
     true_neuron_count, junk_neuron_count = filter_dir_by_model(run_dir, kwargs['model'])
     neuron_count = true_neuron_count + junk_neuron_count
     fraction_true_neurons = true_neuron_count / neuron_count
@@ -79,7 +61,8 @@ def main():
                         help=f"path to the TMD classifier model; defaults to the MSN model: {default_model}")
     args = parser.parse_args()
 
-    call_recut(**vars(args))
+    run_dir = call_recut(**vars(args))
+    persistence_precision_recall(run_dir)
 
 if __name__ == "__main__":
     main()

scripts/run-diadem.py

@@ -0,0 +1,86 @@
+import argparse
+import shutil
+import os, glob
+from recut_interface import call_recut, valid
+
+# 40x lens with 1.5 zoom
+fgpct = .5
+lab_args = {'OlfactoryProjectFibers' : [1, 1, 1]  }
+
+# extract voxel size per lab
+def extract_voxel_size(lab, paths, show=False):
+    lab_voxel_size = 'invalid'
+    if lab in lab_args: 
+        lab_voxel_size = lab_args[lab]
+    else:
+        txt = list(filter(lambda path: os.path.isfile(path) and ("txt" in path), paths))
+        if len(txt) > 0:
+            temp = file_to_voxel_size(txt[0])
+            if len(temp) == 3:
+                lab_voxel_size = temp
+
+    if show:
+        if type(lab_voxel_size) == list:
+            print('  ', end="") # offset
+            print(lab_voxel_size)
+        elif type(lab_voxel_size) == dict:
+            print(lab_voxel_size)
+    return lab_voxel_size
+
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument('dir')
+    # recut arguments
+    parser.add_argument('--fg-percent', type=float)
+    args = parser.parse_args()
+    show = False
+    delete_run_dirs = True
+    if delete_run_dirs:
+        for d in glob.glob("run-*"):
+            shutil.rmtree(d)
+
+    for lab in os.listdir(args.dir):
+        lab_path = os.path.join(args.dir, lab)
+        if os.path.isfile(lab_path):
+            continue;
+        print(lab)
+        image_path = os.path.join(lab_path, 'ImageStacks')
+        subsubdirs = [image for image in os.listdir(image_path)]
+        paths = list(map(lambda image: os.path.join(image_path, image), subsubdirs))
+
+        # TODO match the voxel size to each image of the lab
+        lab_voxel_size = extract_voxel_size(lab, paths)
+
+        # TODO extract soma from each proofread swc
+        # TODO write seed
+        # TODO convert file format
+
+        # extract individual images
+        image_paths = list(filter(lambda dirp: "run" not in dirp, filter(lambda image: os.path.isdir(image), paths)))
+        image_paths.sort()
+        for image in image_paths:
+            #TODO extract voxel size from lab's dict
+            image_voxel_size = 'invalid'
+            if type(lab_voxel_size) == dict:
+                _, tail = os.path.split(image)
+                image_voxel_size = lab_voxel_size[tail]
+            elif type(lab_voxel_size) == list:
+                image_voxel_size = lab_voxel_size
+            else:
+                for file in os.listdir(image):
+                    if "txt" in file:
+                        fpath = os.path.join(image, file)
+                        image_voxel_size = file_to_voxel_size(fpath)
+            if valid(image_voxel_size):
+                if show:
+                    print('    ', end="")
+                    print(image_voxel_size)
+                # build args
+                args = { 'image': image, 'voxel_size_x' : image_voxel_size[0],
+                        'voxel_size_y' : image_voxel_size[1], 'voxel_size_z' : image_voxel_size[2],
+                        'seed_action' : 'find', 'preserve_topology' : '',
+                        'fg_percent' : fgpct}
+                call_recut(**args)
+
+if __name__ == "__main__":
+    main()

scripts/visualize-graph.py

@@ -0,0 +1,25 @@
+#!/opt/local/bin/python
+from pygel3d import hmesh, graph, gl_display as gl
+from os import getcwd
+import argparse
+
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument('graph', help="*.graph file")
+    args = parser.parse_args()
+
+    viewer = gl.Viewer()
+
+    print('Building FEQ...')
+    s = graph.load(args.graph)
+    # viewer.display(s, reset_view=True)
+    m_skel = hmesh.skeleton_to_feq(s)#, [5.0]*len(s.nodes()))
+    hmesh.cc_split(m_skel)
+    hmesh.cc_smooth(m_skel)
+    manifold = hmesh.Manifold(m_skel)
+
+    print("Displaying. HIT ESC IN GRAPHICS WINDOW TO PROCEED...")
+    viewer.display(m_skel, reset_view=True)
+
+if __name__ == "__main__":
+    main()