From 008a47fd66875dd5cd4e13d0b5b3347ff0aa37a2 Mon Sep 17 00:00:00 2001
From: Marius Causemann <mariuscausemann@yahoo.de>
Date: Fri, 19 Sep 2025 18:32:46 +0200
Subject: [PATCH 01/21] replace pyclesperanto with nbmorph

---
 config_files/test.yml                         |  10 +-
 workflow/Snakefile                            |   2 +-
 workflow/envs/nbmorph_env.yml                 |  13 ++
 workflow/scripts/process_image_data.py        |  80 +++----
 workflow/scripts/process_image_data_opencl.py | 217 ++++++++++++++++++
 5 files changed, 270 insertions(+), 52 deletions(-)
 create mode 100644 workflow/envs/nbmorph_env.yml
 create mode 100644 workflow/scripts/process_image_data_opencl.py

diff --git a/config_files/test.yml b/config_files/test.yml
index f4cfb9b..75bcc12 100644
--- a/config_files/test.yml
+++ b/config_files/test.yml
@@ -8,15 +8,15 @@ fixed :
     operation : 
       - "merge labels='[864691136534887842, 864691132131484059, 864691132336155009, 864691132348210676, 864691135791214365, 864691135971392446]'"
       - "removeislands minsize=5000"
-      - "dilate iterations=3 radius=3"
+      - "dilate radius=9"
       - "smooth iterations=3 radius=5"
-      - "erode iterations=1 radius=1"
+      - "erode radius=1"
       - "roigenerate labels='[864691136534887842]'"
-      - "roidilate iterations=4 radius=6"
+      - "roidilate radius=24"
       - "roiapply"
-      - "erode iterations=3 radius=5 labels='[864691136534887842]'"
+      - "erode radius=15 labels='[864691136534887842]'"
       - "removeislands minsize=5000"
-      - "roierode iterations=2 radius=1"
+      - "roierode radius=2"
 raw:
   mip : 1
   size : [5000] 
diff --git a/workflow/Snakefile b/workflow/Snakefile
index 7e3bddf..ba67a0d 100644
--- a/workflow/Snakefile
+++ b/workflow/Snakefile
@@ -91,7 +91,7 @@ rule processImageData:
         #imagestatistic,
         outfile=processed_data_path,
     conda:
-        "envs/pycle_env.yml"
+        "envs/nbmorph_env.yml"
     resources:
         ntasks=10,
         time=30,
diff --git a/workflow/envs/nbmorph_env.yml b/workflow/envs/nbmorph_env.yml
new file mode 100644
index 0000000..f8030ea
--- /dev/null
+++ b/workflow/envs/nbmorph_env.yml
@@ -0,0 +1,13 @@
+name: intel
+channels:
+  - conda-forge
+dependencies:
+  - dask-image
+  - pyvista
+  - ipython
+  - numba
+  - pip
+  - pip:
+      - connected-components-3d
+      - fastremap
+      - nbmorph
\ No newline at end of file
diff --git a/workflow/scripts/process_image_data.py b/workflow/scripts/process_image_data.py
index 621f3db..3b68226 100644
--- a/workflow/scripts/process_image_data.py
+++ b/workflow/scripts/process_image_data.py
@@ -5,19 +5,18 @@
 import time
 import pyvista as pv
 from pathlib import Path
-from utils import get_cell_frequencies, np2pv
+from utils import np2pv
 from dask_image.ndinterp import affine_transform
 import fastremap
-import skimage.morphology as skim
-import pyclesperanto_prototype as cle
-from cle_patch import opening_labels, closing_labels, erode_labels
 import dask.array as da
 import dask
 from functools import partial
-from collections.abc import Iterable
 import cc3d
-cle.set_wait_for_kernel_finish(True)
-dask.config.set({"array.chunk-size": "512 MiB"})
+import numba
+import nbmorph
+
+numba.set_num_threads(1)
+dask.config.set({"array.chunk-size": "1024 MiB"})
 
 def mergecells(img, labels):
     print(f"merging cells: {labels},  ({img.shape})")
@@ -28,30 +27,26 @@ def ncells(img, ncells):
     cell_labels, cell_counts = fastremap.unique(img, return_counts=True)
     cell_labels = cell_labels[np.argsort(cell_counts)]
     cois = list(cell_labels[-ncells :])
-    img = da.where(da.isin(img, cois), img, 0)
+    img = np.where(np.isin(img, cois), img, 0)
     return img
     
-def dilate(img, iterations, radius, labels=None):
+def dilate(img, radius, labels=None):
     print(f"dilating cells,  ({img.shape})")
     if labels is None:
-        for i in range(iterations):
-            img = cle.dilate_labels(img, radius=radius)
-        img = cle.pull(img)
+        img = nbmorph.dilate_labels_spherical(img, radius=radius)
     else:
         vipimg = np.where(np.isin(img, labels), img, 0)
-        vipimg = dilate(vipimg, iterations=iterations, radius=radius)
+        vipimg = dilate(vipimg, radius=radius)
         img = np.where(vipimg, vipimg, img)
     return img
 
-def erode(img, iterations, radius, labels=None):
+def erode(img, radius, labels=None):
     print(f"eroding cells,  ({img.shape})")
     if labels is None:
-        for i in range(iterations):
-            img = erode_labels(img, radius=radius)
-        img = cle.pull(img)
+        img = nbmorph.erode_labels_spherical(img, radius=radius)
     else:
         vipimg = np.where(np.isin(img, labels), img, 0)
-        vipimg = erode(vipimg, iterations=iterations, radius=radius)
+        vipimg = erode(vipimg, radius=radius)
         orig_wo_vips = np.where(np.isin(img, labels), 0, img)
         img = np.where(orig_wo_vips > vipimg, orig_wo_vips, vipimg)
     return img
@@ -59,10 +54,8 @@ def erode(img, iterations, radius, labels=None):
 def smooth(img, iterations, radius, labels=None):
     print(f"smoothing cells,  ({img.shape})")
     if labels is None:
-        for i in range(iterations):
-            img = opening_labels(img, radius=radius)
-            img = closing_labels(img, radius=radius)
-        img = cle.pull(img)
+        img = nbmorph.smooth_labels_spherical(img, radius=radius,
+                                              iterations=iterations, dilate_radius=radius)
     else:
         vipimg = np.where(np.isin(img, labels), img, 0)
         vipimg = smooth(vipimg, iterations=iterations, radius=radius)
@@ -75,13 +68,9 @@ def smooth(img, iterations, radius, labels=None):
 def removeislands(img, minsize):
     return cc3d.dust(img, threshold=minsize, connectivity=6)
 
-    #would be more performant, but relabels all cells..
-    #return cle.exclude_small_labels(img, maximum_size=minsize)
-
 opdict ={"merge": mergecells, "smooth":smooth, "dilate":dilate,
          "erode":erode, "removeislands":removeislands, "ncells":ncells}
 
-
 def _parse_to_dict(values):
     result = {}
     for value in values:
@@ -109,22 +98,19 @@ def parse_operations(ops):
     parser.add_argument('-o','--operation', nargs='+', action='append', help="operations to be performed on the segmented image")
 
     args = parser.parse_args()
-    parallel = args.nworkers
-    print(f"Using {parallel} workers...")
-    print(f"Available  devices {cle.available_device_names()}")
-    cle.select_device("NVIDIA")
-    dev = cle.get_device()
-    print(f"Running processing on {dev}")
-    mem_gb = dev.device.global_mem_size / 1e9
-    print(f"device memory detected: {mem_gb:.1f} GB")
-
+    n_parallel = args.nworkers
+    print(f"Using {n_parallel} workers...")
     start = time.time()
+
+    # read image file
     imggrid = pv.read(args.infile)
     img = imggrid["data"]
     dims = imggrid.dimensions
     resolution = imggrid.spacing
     img = img.reshape(dims - np.array([1, 1, 1]), order="F")
     img = da.from_array(img)
+
+    # get cells labels, and filter by n-largest (if requested)
     cell_labels, cell_counts = fastremap.unique(img, return_counts=True)
     if args.ncells:
         cell_labels = list(cell_labels[np.argsort(cell_counts)])
@@ -135,22 +121,27 @@ def parse_operations(ops):
         cell_labels = list(cell_labels)
         cell_labels.remove(0)
 
+    # remap labels to smaller, sequential ints
     remapping = {int(c):i for i,c in enumerate([0] + cell_labels)}
     remap = lambda ids: [remapping[int(i)] for i in ids if i in remapping.keys()]
     img = img.map_blocks(partial(fastremap.remap, table=remapping), dtype=img.dtype)
     img = img.map_blocks(partial(fastremap.refit, value=len(cell_labels)))
 
+    # interpolate into the specified, isotropic grid with size dx
     dx = args.dx
     scale = np.diag([dx / r for r in resolution] + [1])
     new_dims = [int(d * r / dx) for d,r in zip(dims, resolution)]
     img = affine_transform(img, scale, output_shape=new_dims, order=0,
                            output_chunks=500)
+    img = dask.compute(img, num_workers= args.nworkers)[0]
     print(f"image size: {img.shape}")
     resolution = [dx]*3
     roi = None
-    operations = parse_operations(args.operation)
 
+    # parse user specified operations, and iterate over them:
+    operations = parse_operations(args.operation)
     for op, kwargs in operations:
+        print(op, kwargs)
         if "labels" in kwargs.keys():
             labels = kwargs["labels"]
             if kwargs.get("allexcept", False):
@@ -159,28 +150,25 @@ def parse_operations(ops):
                 kwargs["labels"] = remap(labels)
 
         if op=="roigenerate":
-            roi = da.isin(img, kwargs["labels"])
+            roi = np.isin(img, kwargs["labels"])
             continue
         if op=="roiapply":
-            img = da.where(roi, img, 0) 
+            img = np.where(roi, img, 0) 
             continue
         if op=="ncells":
             img = ncells(img, **kwargs)
             continue
         if op.startswith("roi"):
             roiop = op[3:]
-            roi = roi.map_overlap(partial(opdict[roiop], **kwargs), depth=30, dtype=np.uint8)
+            roi = opdict[roiop](roi, **kwargs)
         else:
-            img = img.map_overlap(partial(opdict[op], **kwargs), depth=30, dtype=img.dtype)
+            img =opdict[op](img, **kwargs)
     
-    chunk_mem_gb = np.prod(img.chunksize) * np.nbytes[img.dtype] / 1e9
-    print(f"chunk size: {chunk_mem_gb} GB")
-    max_workers = int(mem_gb / (chunk_mem_gb*3)) # assume max 3x chunk memory is used
-    img, roi = dask.compute(img, roi, num_workers=min(max_workers, args.nworkers))
+    #img, roi = dask.compute(img, roi, num_workers= args.nworkers)
     print(f"processed! {img.shape}")
     proc_time = time.time()
 
-    img = da.array(img)
+    # remap labels to smaller, sequential ints again, since many labels might have disappeared...
     cell_labels, cell_counts = fastremap.unique(img, return_counts=True)
     cell_labels = list(cell_labels[np.argsort(cell_counts)])
     cell_labels.remove(0)
@@ -190,7 +178,7 @@ def parse_operations(ops):
     remapping2.update({int(c):i + 2 for i,c in enumerate(cell_labels)})
     img = img.map_blocks(partial(fastremap.remap, table=remapping2), dtype=img.dtype)
     img = img.map_blocks(partial(fastremap.refit, value=max(remapping2.values())))
-    img = dask.compute(img, num_workers=min(max_workers, args.nworkers))[0]
+    img = dask.compute(img, num_workers= args.nworkers)[0]
 
     imggrid = np2pv(img, resolution)
     if roi is not None:
diff --git a/workflow/scripts/process_image_data_opencl.py b/workflow/scripts/process_image_data_opencl.py
new file mode 100644
index 0000000..621f3db
--- /dev/null
+++ b/workflow/scripts/process_image_data_opencl.py
@@ -0,0 +1,217 @@
+import numpy as np
+import argparse
+import yaml
+import json
+import time
+import pyvista as pv
+from pathlib import Path
+from utils import get_cell_frequencies, np2pv
+from dask_image.ndinterp import affine_transform
+import fastremap
+import skimage.morphology as skim
+import pyclesperanto_prototype as cle
+from cle_patch import opening_labels, closing_labels, erode_labels
+import dask.array as da
+import dask
+from functools import partial
+from collections.abc import Iterable
+import cc3d
+cle.set_wait_for_kernel_finish(True)
+dask.config.set({"array.chunk-size": "512 MiB"})
+
+def mergecells(img, labels):
+    print(f"merging cells: {labels},  ({img.shape})")
+    img = np.where(np.isin(img, labels), labels[0], img)
+    return img
+
+def ncells(img, ncells):
+    cell_labels, cell_counts = fastremap.unique(img, return_counts=True)
+    cell_labels = cell_labels[np.argsort(cell_counts)]
+    cois = list(cell_labels[-ncells :])
+    img = da.where(da.isin(img, cois), img, 0)
+    return img
+    
+def dilate(img, iterations, radius, labels=None):
+    print(f"dilating cells,  ({img.shape})")
+    if labels is None:
+        for i in range(iterations):
+            img = cle.dilate_labels(img, radius=radius)
+        img = cle.pull(img)
+    else:
+        vipimg = np.where(np.isin(img, labels), img, 0)
+        vipimg = dilate(vipimg, iterations=iterations, radius=radius)
+        img = np.where(vipimg, vipimg, img)
+    return img
+
+def erode(img, iterations, radius, labels=None):
+    print(f"eroding cells,  ({img.shape})")
+    if labels is None:
+        for i in range(iterations):
+            img = erode_labels(img, radius=radius)
+        img = cle.pull(img)
+    else:
+        vipimg = np.where(np.isin(img, labels), img, 0)
+        vipimg = erode(vipimg, iterations=iterations, radius=radius)
+        orig_wo_vips = np.where(np.isin(img, labels), 0, img)
+        img = np.where(orig_wo_vips > vipimg, orig_wo_vips, vipimg)
+    return img
+
+def smooth(img, iterations, radius, labels=None):
+    print(f"smoothing cells,  ({img.shape})")
+    if labels is None:
+        for i in range(iterations):
+            img = opening_labels(img, radius=radius)
+            img = closing_labels(img, radius=radius)
+        img = cle.pull(img)
+    else:
+        vipimg = np.where(np.isin(img, labels), img, 0)
+        vipimg = smooth(vipimg, iterations=iterations, radius=radius)
+        # remove labelled cells from original image
+        orig_wo_vips = np.where(np.isin(img, labels), 0, img)
+        # insert smoothed labeled cells in original (overwrite original)
+        img = np.where(vipimg, vipimg, orig_wo_vips)
+    return img
+
+def removeislands(img, minsize):
+    return cc3d.dust(img, threshold=minsize, connectivity=6)
+
+    #would be more performant, but relabels all cells..
+    #return cle.exclude_small_labels(img, maximum_size=minsize)
+
+opdict ={"merge": mergecells, "smooth":smooth, "dilate":dilate,
+         "erode":erode, "removeislands":removeislands, "ncells":ncells}
+
+
+def _parse_to_dict(values):
+    result = {}
+    for value in values:
+        k, v = value.split('=')
+        result[k.strip()] = yaml.safe_load(v.strip())
+    return result
+
+def parse_operations(ops):
+    parsed = []
+    for op in ops:
+        subargs =  _parse_to_dict(op[1:])
+        parsed.append((op[0], subargs))
+    return parsed
+    
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--infile", help="input data", type=str)
+    parser.add_argument(
+        "--output", help="output filename", type=str, default="processeddata.vtk"
+    )
+    parser.add_argument("--nworkers", help="number of workers", type=int, default=1)
+    parser.add_argument("--dx", help="target resolution", type=int, default=None)
+    parser.add_argument("--ncells", help="number of cells", type=int, default=None)
+    parser.add_argument('-o','--operation', nargs='+', action='append', help="operations to be performed on the segmented image")
+
+    args = parser.parse_args()
+    parallel = args.nworkers
+    print(f"Using {parallel} workers...")
+    print(f"Available  devices {cle.available_device_names()}")
+    cle.select_device("NVIDIA")
+    dev = cle.get_device()
+    print(f"Running processing on {dev}")
+    mem_gb = dev.device.global_mem_size / 1e9
+    print(f"device memory detected: {mem_gb:.1f} GB")
+
+    start = time.time()
+    imggrid = pv.read(args.infile)
+    img = imggrid["data"]
+    dims = imggrid.dimensions
+    resolution = imggrid.spacing
+    img = img.reshape(dims - np.array([1, 1, 1]), order="F")
+    img = da.from_array(img)
+    cell_labels, cell_counts = fastremap.unique(img, return_counts=True)
+    if args.ncells:
+        cell_labels = list(cell_labels[np.argsort(cell_counts)])
+        cell_labels.remove(0)
+        cois = list(cell_labels[-args.ncells :])
+        img = da.where(da.isin(img, cois), img, 0)
+    else:
+        cell_labels = list(cell_labels)
+        cell_labels.remove(0)
+
+    remapping = {int(c):i for i,c in enumerate([0] + cell_labels)}
+    remap = lambda ids: [remapping[int(i)] for i in ids if i in remapping.keys()]
+    img = img.map_blocks(partial(fastremap.remap, table=remapping), dtype=img.dtype)
+    img = img.map_blocks(partial(fastremap.refit, value=len(cell_labels)))
+
+    dx = args.dx
+    scale = np.diag([dx / r for r in resolution] + [1])
+    new_dims = [int(d * r / dx) for d,r in zip(dims, resolution)]
+    img = affine_transform(img, scale, output_shape=new_dims, order=0,
+                           output_chunks=500)
+    print(f"image size: {img.shape}")
+    resolution = [dx]*3
+    roi = None
+    operations = parse_operations(args.operation)
+
+    for op, kwargs in operations:
+        if "labels" in kwargs.keys():
+            labels = kwargs["labels"]
+            if kwargs.get("allexcept", False):
+                kwargs["labels"] = list(set(remapping.values()) - set(remap(labels)))
+            else:
+                kwargs["labels"] = remap(labels)
+
+        if op=="roigenerate":
+            roi = da.isin(img, kwargs["labels"])
+            continue
+        if op=="roiapply":
+            img = da.where(roi, img, 0) 
+            continue
+        if op=="ncells":
+            img = ncells(img, **kwargs)
+            continue
+        if op.startswith("roi"):
+            roiop = op[3:]
+            roi = roi.map_overlap(partial(opdict[roiop], **kwargs), depth=30, dtype=np.uint8)
+        else:
+            img = img.map_overlap(partial(opdict[op], **kwargs), depth=30, dtype=img.dtype)
+    
+    chunk_mem_gb = np.prod(img.chunksize) * np.nbytes[img.dtype] / 1e9
+    print(f"chunk size: {chunk_mem_gb} GB")
+    max_workers = int(mem_gb / (chunk_mem_gb*3)) # assume max 3x chunk memory is used
+    img, roi = dask.compute(img, roi, num_workers=min(max_workers, args.nworkers))
+    print(f"processed! {img.shape}")
+    proc_time = time.time()
+
+    img = da.array(img)
+    cell_labels, cell_counts = fastremap.unique(img, return_counts=True)
+    cell_labels = list(cell_labels[np.argsort(cell_counts)])
+    cell_labels.remove(0)
+
+    img = da.array(img)
+    remapping2 = {0:0}
+    remapping2.update({int(c):i + 2 for i,c in enumerate(cell_labels)})
+    img = img.map_blocks(partial(fastremap.remap, table=remapping2), dtype=img.dtype)
+    img = img.map_blocks(partial(fastremap.refit, value=max(remapping2.values())))
+    img = dask.compute(img, num_workers=min(max_workers, args.nworkers))[0]
+
+    imggrid = np2pv(img, resolution)
+    if roi is not None:
+        imggrid["roimask"] = np.array(roi).flatten(order="F")
+
+    resdir = Path(args.output).parent
+
+    resdir.mkdir(parents=True, exist_ok=True)
+    imggrid.save(args.output)
+    save_time = time.time()
+    print(f"saving time: {save_time - proc_time} s")
+
+
+    combinedmap = {k:remapping2[v] for k,v in remapping.items() if v in remapping2.keys()}
+    mesh_statistics = dict()
+    mesh_statistics["cell_labels"] = cell_labels
+    mesh_statistics["cell_counts"] = cell_counts
+    mesh_statistics["mapping"] = combinedmap
+
+    for k, v in mesh_statistics.items():
+        mesh_statistics[k] = np.array(v).tolist()
+
+    with open(resdir / "imagestatistic.yml", "w") as mesh_stat_file:
+        yaml.dump(mesh_statistics, mesh_stat_file)

From a6d51fc9389c4f21ed432da2fd14e9d516f551fd Mon Sep 17 00:00:00 2001
From: Marius Causemann <mariusca@simula.no>
Date: Fri, 3 Oct 2025 15:20:44 +0200
Subject: [PATCH 02/21] add webknossos cache

---
 .gitignore | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/.gitignore b/.gitignore
index 82736c2..a4f6ce5 100644
--- a/.gitignore
+++ b/.gitignore
@@ -8,4 +8,5 @@ fTetWild
 zipped
 zipped_surf
 .snakemake
-.vscode
\ No newline at end of file
+.vscode
+.cache
\ No newline at end of file

From eed9e3f3b7216b5381ce978afc3305294220c349 Mon Sep 17 00:00:00 2001
From: Marius Causemann <mariusca@simula.no>
Date: Fri, 3 Oct 2025 15:21:14 +0200
Subject: [PATCH 03/21] update config

---
 config_files/cappilary1.yml   | 12 ++++++------
 config_files/cortical_mm3.yml |  4 ++--
 config_files/motta2019.yml    | 31 +++++++++++++++++++------------
 config_files/test.yml         |  5 ++---
 4 files changed, 29 insertions(+), 23 deletions(-)

diff --git a/config_files/cappilary1.yml b/config_files/cappilary1.yml
index dd2d34f..601a865 100644
--- a/config_files/cappilary1.yml
+++ b/config_files/cappilary1.yml
@@ -9,17 +9,17 @@ fixed :
     #cells : 864691134581644554 864691135228370458 864691135918336944 864691136475874625 864691136723365629
     operation : 
       - "merge labels='[864691136534887842, 864691135476973891, 864691135751005697, 864691133070457251, 864691136178683094, 864691133070456483, 864691133070450083, 864691133003757151, 864691136475845441, 864691134032036203, 864691131901467286]'"
-      - "dilate iterations=2 radius=2 labels='[864691134581644554, 864691135228370458, 864691135918336944, 864691136475874625, 864691136723365629, 864691134308859492]'"
+      - "dilate radius=4 labels='[864691134581644554, 864691135228370458, 864691135918336944, 864691136475874625, 864691136723365629, 864691134308859492]'"
       - "removeislands minsize=5000"
-      - "dilate iterations=3 radius=3"
+      - "dilate radius=9"
       - "smooth iterations=3 radius=5"
-      - "erode iterations=1 radius=1"
+      - "erode radius=1"
       - "roigenerate labels='[864691136534887842]'"
-      - "roidilate iterations=10 radius=8"
+      - "roidilate radius=80"
       - "roiapply"
-      - "erode iterations=3 radius=5 labels='[864691136534887842]'"
+      - "erode radius=15 labels='[864691136534887842]'"
       - "removeislands minsize=5000"
-      - "roierode iterations=2 radius=2"
+      - "roierode radius=4"
 raw:
   mip : 0
   size : [7500] 
diff --git a/config_files/cortical_mm3.yml b/config_files/cortical_mm3.yml
index 58babc2..7a63f7b 100644
--- a/config_files/cortical_mm3.yml
+++ b/config_files/cortical_mm3.yml
@@ -12,9 +12,9 @@ fixed :
       - "erode iterations=1 radius=1"
       - "removeislands minsize=5000"
 raw:
-  size : [5000, 10000, 20000, 40000,]
+  size : [10000] #[5000, 10000, 20000]
 processing: 
   dx : 20
-  ncells : [5, 10, 50, 100, 200, 500, 1000]
+  ncells : [10] #[5,  50, 200, 500]
 meshing:
   envelopsize : [18]
diff --git a/config_files/motta2019.yml b/config_files/motta2019.yml
index c3c71ed..ad0d9fc 100644
--- a/config_files/motta2019.yml
+++ b/config_files/motta2019.yml
@@ -1,12 +1,19 @@
-dataset : "motta"
-cloudpath : 
-  "motta": "https://webknossos.org/datasets/MPI_Brain_Research/2012-09-28_ex145_07x2_ROI2017_connectome/view/"
-mip : [2]
-position : "2016-3595-2015"
-size : [5000, 10000, 20000] 
-ncells : [5, 10, 50, 100, 200, 1000]
-smoothiter : [0,]
-smoothradius : 0
-expand : [0]
-shrink : 0
-envelopsize : [8]
+name : "motta2019"
+fixed : 
+  raw:
+    cloudpath: "https://webknossos.org/datasets/l4dense_motta_et_al_demo-5f33d5340100001800e8cddc/view"
+    position : "2816-4352-1792"
+  processing:
+    operation : 
+      - "removeislands minsize=8000"
+      - "dilate radius=1"
+      - "smooth iterations=3 radius=6"
+      - "erode radius=3"
+raw:
+  mip : 2
+  size : [4000] 
+processing:
+  dx : 18
+meshing:
+  envelopsize : [12]
+  
diff --git a/config_files/test.yml b/config_files/test.yml
index 75bcc12..1448150 100644
--- a/config_files/test.yml
+++ b/config_files/test.yml
@@ -2,9 +2,8 @@ name : "test"
 fixed : 
   raw:
     cloudpath: "precomputed://gs://iarpa_microns/minnie/minnie65/seg"
-    position : "232042-166129-21467"
+    position : "232000-166129-21467"
   processing:
-    #cells : 864691132723008555 864691134309135460 864691134310007140 864691135609682183 864691135753683149 864691135918336944
     operation : 
       - "merge labels='[864691136534887842, 864691132131484059, 864691132336155009, 864691132348210676, 864691135791214365, 864691135971392446]'"
       - "removeislands minsize=5000"
@@ -19,7 +18,7 @@ fixed :
       - "roierode radius=2"
 raw:
   mip : 1
-  size : [5000] 
+  size : [7000] 
 processing:
   dx : 18
 meshing:

From 2ef76fbd31012ff59e7f4d3f4e0dc508897a42aa Mon Sep 17 00:00:00 2001
From: Marius Causemann <mariusca@simula.no>
Date: Fri, 3 Oct 2025 15:22:18 +0200
Subject: [PATCH 04/21] remove fenics seperate cells

---
 workflow/Snakefile | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/workflow/Snakefile b/workflow/Snakefile
index ba67a0d..163c422 100644
--- a/workflow/Snakefile
+++ b/workflow/Snakefile
@@ -115,7 +115,7 @@ rule extractSurfaces:
         csgtree=csgtree
     resources:
         time=120,
-        ntasks=36
+        ntasks=4
     conda:
         conda_env
     shell:
@@ -143,10 +143,12 @@ rule generateMesh:
         --csgtree {input.csgtree} \
         {params.options} \
         --output {output.outfile} \
-        --max_threads {resources.ntasks} && \
-        python3 workflow/scripts/seperate_touching_cells.py --infile {output.outfile} --output {output.outfile}
+        --max_threads {resources.ntasks} \
         """
 
+#   python3 workflow/scripts/seperate_touching_cells.py --infile {output.outfile} --output {output.outfile}
+
+
 rule generateSurfaceTags:
     input:
         meshfile=mesh_path

From db6040eab85a30fd337a6021bc141378490c797e Mon Sep 17 00:00:00 2001
From: Marius Causemann <mariusca@simula.no>
Date: Fri, 3 Oct 2025 15:22:45 +0200
Subject: [PATCH 05/21] exchange wildmeshing with pytetwild wrapper

---
 workflow/envs/environment.yml | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/workflow/envs/environment.yml b/workflow/envs/environment.yml
index cbfa34e..cd55ccf 100644
--- a/workflow/envs/environment.yml
+++ b/workflow/envs/environment.yml
@@ -2,7 +2,7 @@ name: emimesh
 channels:
   - conda-forge
 dependencies:
-  - python
+  - python=3.10
   - numpy
   - vtk=*=*mesa*
   - matplotlib
@@ -20,10 +20,9 @@ dependencies:
   - pip
   - pip:
     - pyacvd>=0.3
-    - wildmeshing
+    - pytetwild
     - fastremap
     - cloud-volume
     - webknossos
-    - cattrs
     - connected-components-3d
     - cmocean

From de338440b4057aac8f948a37da1f9da5dca493eb Mon Sep 17 00:00:00 2001
From: Marius Causemann <mariusca@simula.no>
Date: Fri, 3 Oct 2025 15:23:06 +0200
Subject: [PATCH 06/21] clean up download script

---
 workflow/scripts/download_data.py | 14 ++++++--------
 1 file changed, 6 insertions(+), 8 deletions(-)

diff --git a/workflow/scripts/download_data.py b/workflow/scripts/download_data.py
index ee0653c..a188d93 100644
--- a/workflow/scripts/download_data.py
+++ b/workflow/scripts/download_data.py
@@ -1,14 +1,11 @@
-from cloudvolume import CloudVolume
-import webknossos as wk
 import numpy as np
 import pyvista as pv
 import argparse
-import pyvista as pv
 from utils import np2pv
-import dask.array as da
 from pathlib import Path
 
 def download_webknossos(cloud_path, mip, pos, physical_size):
+    import webknossos as wk
     target_mag=wk.Mag([2,2,1])
     voxel_size = np.array([11.24,11.24, 28])
     mag_voxel_size = (voxel_size * target_mag.to_np())
@@ -16,17 +13,18 @@ def download_webknossos(cloud_path, mip, pos, physical_size):
     size = [int(ps / vs) for ps, vs in zip(physical_size, voxel_size)]
     bbox = wk.BoundingBox(pos, size=size)
     bbox.align_with_mag(target_mag)
-    ds = wk.Dataset.download(cloud_path, mags=[mag], path=f"data/{mip}_{physical_size}", exist_ok=True,
+    ds = wk.Dataset.download(cloud_path, mags=[mag], path=f".cache/webknossos/{mip}_{physical_size}",
                              bbox=bbox, layers="segmentation")
     layer = ds.get_layer("segmentation")
     layer.downsample_mag(from_mag=mag, target_mag=target_mag, allow_overwrite=True)
     mag_view = layer.get_mag(target_mag)
     img = mag_view.read().squeeze()
-    assert img.max() > 0
+    assert img.sum() > 0, "dataset empty!"
     return img, mag_voxel_size
 
 
 def download_cloudvolume(cloud_path, mip, pos, physical_size):
+    from cloudvolume import CloudVolume
     vol = CloudVolume(
         cloud_path, use_https=True, parallel=8, progress=True, mip=mip, cache=True, bounded=True
     )
@@ -47,7 +45,7 @@ def download_cloudvolume(cloud_path, mip, pos, physical_size):
         "--cloudpath",
         help="path to cloud data",
         type=str,
-        #default="precomputed://gs://iarpa_microns/minnie/minnie65/seg",
+        default="precomputed://gs://iarpa_microns/minnie/minnie65/seg",
     )
     parser.add_argument("--mip", help="resolution (0 is highest)", type=int, default=0)
     parser.add_argument(
@@ -81,5 +79,5 @@ def download_cloudvolume(cloud_path, mip, pos, physical_size):
         
     print(res)
     data = np2pv(img, res)
-    Path(args.output).parent.mkdir(exist_ok=True)
+    Path(args.output).parent.mkdir(exist_ok=True, parents=True)
     data.save(args.output)

From bd16378d387c6c590c6a1830fd9bc4f5a97c4840 Mon Sep 17 00:00:00 2001
From: Marius Causemann <mariusca@simula.no>
Date: Fri, 3 Oct 2025 15:23:21 +0200
Subject: [PATCH 07/21] make yaml human readable

---
 workflow/scripts/evaluate_mesh.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/workflow/scripts/evaluate_mesh.py b/workflow/scripts/evaluate_mesh.py
index d9610ea..c9b57f1 100644
--- a/workflow/scripts/evaluate_mesh.py
+++ b/workflow/scripts/evaluate_mesh.py
@@ -2,7 +2,6 @@
 import numpy as np
 import argparse
 import yaml
-from pathlib import Path
 import matplotlib.pyplot as plt
 import seaborn as sns
 
@@ -37,8 +36,8 @@ def compute_surface_volume(mesh, cell_ids):
     for cid in cell_ids:
         cell = mesh.extract_cells(np.isin(mesh.cell_data[lstr], [cid]))
         surf = cell.extract_surface()
-        surface_areas.append(surf.compute_cell_sizes()["Area"].sum())
-        volumes.append(cell["Volume"].sum())
+        surface_areas.append(float(surf.compute_cell_sizes()["Area"].sum()))
+        volumes.append(float(cell["Volume"].sum()))
     return volumes, surface_areas
 
 
@@ -87,6 +86,7 @@ def plot_local_width(width, volume, filename):
     cell_boundary_mesh = boundary_mesh.extract_cells(boundary_mesh[lstr]>ecs_id)
     ecs_boundary_mesh = boundary_mesh.extract_cells(boundary_mesh[lstr]==ecs_id)
     n_ecs_boundary_points = boundary_mesh.number_of_points - cell_boundary_mesh.number_of_points
+
     results = dict(npoints=mesh.number_of_points,
                   ncompcells=mesh.number_of_cells,
                   ecs_volume=ecs_volume, ecs_surface=ecs_surface,

From ea6adf7c39951d6bfbfcab156c2ac77da824f4ad Mon Sep 17 00:00:00 2001
From: Marius Causemann <mariusca@simula.no>
Date: Fri, 3 Oct 2025 15:23:47 +0200
Subject: [PATCH 08/21] make sure surfaces are nan-free

---
 workflow/scripts/extract_surfaces.py | 75 ++++++++++++++++++++++------
 1 file changed, 61 insertions(+), 14 deletions(-)

diff --git a/workflow/scripts/extract_surfaces.py b/workflow/scripts/extract_surfaces.py
index 37d8beb..0215e1a 100644
--- a/workflow/scripts/extract_surfaces.py
+++ b/workflow/scripts/extract_surfaces.py
@@ -39,23 +39,70 @@ def clip_closed_box(surf, box):
     centers = box.cell_centers().points
     for midp, n in zip(centers, box.cell_normals):
         clip_closed_surface(surf, normal=-n, origin=midp, inplace=True)
-        
+
+def clean_mesh_nan_points(grid: pv.PolyData):
+    """
+    Replaces NaN point coordinates in a PyVista grid with the mean of their
+    connected neighbors using the efficient `point_neighbors` method.
+
+    Args:
+        grid: A PyVista PolyData object that may contain NaN values in its points.
+
+    """
+    grid = grid
+    points = grid.points
+
+    # Find the indices of points where any coordinate is NaN
+    nan_point_indices = np.where(np.isnan(points).any(axis=1))[0]
+
+    if nan_point_indices.size == 0:
+        return
+
+    print(f"Found {len(nan_point_indices)} points with NaN coordinates. Cleaning them...")
+
+    # Iterate through each point with NaN coordinates
+    for point_idx in nan_point_indices:
+        # Directly get the indices of the neighboring points
+        neighbor_indices = grid.point_neighbors(point_idx)
+
+        if not neighbor_indices:
+            # As a fallback, if the point has no neighbors, replace with origin.
+            points[point_idx] = [0, 0, 0]
+            continue
+
+        # Get the coordinates of the neighbors
+        neighbor_coords = points[neighbor_indices]
+
+        # Filter out any neighbors that are also NaN
+        valid_neighbors = neighbor_coords[~np.isnan(neighbor_coords).any(axis=1)]
+
+        if valid_neighbors.shape[0] > 0:
+            # Calculate the mean of the valid neighbors and replace the NaN point
+            points[point_idx] = np.mean(valid_neighbors, axis=0)
+        else:
+            # Fallback if all neighbors are also NaN
+            points[point_idx] = [0, 0, 0]
+
 def n_point_target(n, mesh_reduction_factor):
     return int(50 + n / mesh_reduction_factor)
 
 
 def extract_surface(mask, grid, mesh_reduction_factor, taubin_smooth_iter, filename=None):
     mesh = grid.contour([0.5], mask.flatten(order="F"), method="marching_cubes")
-    surf = mesh.extract_geometry()
-    surf.clear_data()
-    n_points = surf.number_of_points
-    clus = pyacvd.Clustering(surf)
+    origsurf = mesh.extract_geometry()
+    origsurf.clear_data()
+    n_points = origsurf.number_of_points
+    if n_points < 10: return False
+    clus = pyacvd.Clustering(origsurf)
     clus.cluster(n_point_target(n_points, mesh_reduction_factor))
     surf = clus.create_mesh()
     surf.smooth_taubin(n_iter=taubin_smooth_iter, inplace=True)
+    clean_mesh_nan_points(surf)
+    assert np.isnan(surf.points).any() == False
     if surf.number_of_points > 10 and filename is not None:
         print(f"saving : {filename}")
         pv.save_meshio(filename, surf)
+    
     return surf
 
 def extract_surf_id(obj_id, mesh_reduction_factor, taubin_smooth_iter,
@@ -71,7 +118,7 @@ def extract_cell_meshes(
     cell_labels,
     resolution,
     mesh_reduction_factor=10,
-    taubin_smooth_iter=0,
+    taubin_smooth_iter=5,
     write_dir=None,
     ncpus=1
 ):
@@ -92,11 +139,11 @@ def extract_cell_meshes(
         pool.join()
     return surfaces
 
-def create_balanced_csg_json_tree(surface_files):
+def create_balanced_csg_tree(surface_files):
     n = len(surface_files)
     if  n >= 2:
-        return {"operation": "union", "left": create_balanced_csg_json_tree(surface_files[:int(n/2)])
-                                    , "right": create_balanced_csg_json_tree(surface_files[int(n/2):])}
+        return {"operation": "union", "left": create_balanced_csg_tree(surface_files[:int(n/2)])
+                                    , "right": create_balanced_csg_tree(surface_files[int(n/2):])}
     return surface_files[0]
 
 if __name__ == "__main__":
@@ -141,7 +188,7 @@ def create_balanced_csg_json_tree(surface_files):
 
     surfs = extract_cell_meshes(
         img,
-        cell_labels,
+        cell_labels[::-1],
         resolution,
         mesh_reduction_factor=10,
         taubin_smooth_iter=5,
@@ -149,10 +196,10 @@ def create_balanced_csg_json_tree(surface_files):
         ncpus=args.ncpus
     )
 
-    mesh_files = [outdir / f"{cid}.ply" for cid in surfs]
+    mesh_files = [outdir / f"{cid}.ply" for cid in surfs if cid]
     roisurf.save(roi_file)
-    csg_tree = create_balanced_csg_json_tree([str(f) for f in mesh_files])
-    csg_tree = create_balanced_csg_json_tree([str(roi_file), csg_tree])
-
+    csg_tree = create_balanced_csg_tree([str(f) for f in mesh_files])
+    csg_tree = create_balanced_csg_tree([str(roi_file), csg_tree])
+    csg_tree = {"operation":"intersection","right":csg_tree, "left":str(roi_file)}
     with open(outdir / "csgtree.json", "w") as outfile:
         outfile.write(json.dumps(csg_tree))

From 9e5724a7ed6fcaea02e3d76ffa60040bd4535756 Mon Sep 17 00:00:00 2001
From: Marius Causemann <mariusca@simula.no>
Date: Fri, 3 Oct 2025 15:24:02 +0200
Subject: [PATCH 09/21] use pytetwild

---
 workflow/scripts/generate_mesh.py | 42 +++++++++++--------------------
 1 file changed, 15 insertions(+), 27 deletions(-)

diff --git a/workflow/scripts/generate_mesh.py b/workflow/scripts/generate_mesh.py
index a0a37e0..5e690db 100644
--- a/workflow/scripts/generate_mesh.py
+++ b/workflow/scripts/generate_mesh.py
@@ -1,5 +1,4 @@
 import json
-import meshio
 import pyvista as pv
 import argparse
 import numpy as np
@@ -12,29 +11,17 @@ def get_values(d):
         else:
             yield v
 
-def mesh_surfaces(csg_tree, eps, stop_quality, max_threads):
-    import wildmeshing as wm
-    tetra = wm.Tetrahedralizer(
-        epsilon=eps,
-        edge_length_r=eps * 50,
-        coarsen=True,
-        stop_quality=stop_quality,
-        max_threads=max_threads,
-        skip_simplify=False,
-    )
-    tetra.load_csg_tree(json.dumps(csg_tree))
+def mesh_surfaces(csg_tree_path, eps, stop_quality, max_threads):
+    from pytetwild import tetrahedralize_csg
     start = time.time()
-    tetra.tetrahedralize()
-    point_array, cell_array, marker = tetra.get_tet_mesh()
-
-    volmesh = pv.from_meshio(
-        meshio.Mesh(
-            point_array, [("tetra", cell_array)], cell_data={"label": [marker.ravel()]}
-        )
-    )
-    volmesh.field_data['runtime'] = time.time() - start
-    volmesh.field_data['threads'] = max_threads
-    return volmesh
+    mesh = tetrahedralize_csg(csg_tree_path, epsilon=eps, edge_length_r=eps*50, 
+                              coarsen=True, stop_energy=stop_quality,
+                              num_threads=max_threads).clean()
+    print("meshing finished!")
+    mesh["label"] = mesh["marker"]
+    mesh.field_data['runtime'] = time.time() - start
+    mesh.field_data['threads'] = max_threads
+    return mesh
 
 if __name__ == "__main__":
     parser = argparse.ArgumentParser()
@@ -67,12 +54,13 @@ def mesh_surfaces(csg_tree, eps, stop_quality, max_threads):
     roifile = [s for s in surfs if "roi.ply" in s][0]
     roi = pv.read(roifile)
     diag = np.sqrt(3) * roi.volume ** (1 / 3)
-    es = args.envelopsize
-    csgtree = {"operation":"intersection","right":csgtree, "left":roifile}
+    abs_eps = args.envelopsize    
     volmesh = mesh_surfaces(
-        csgtree,
-        eps=es / diag,
+        args.csgtree,
+        eps=abs_eps / diag,
         stop_quality=args.stopquality,
         max_threads=args.max_threads,
     )
     pv.save_meshio(args.output, volmesh)
+    print(volmesh.array_names)
+

From 933d21d5ceca38605c97ed77880600bae2f0c819 Mon Sep 17 00:00:00 2001
From: Marius Causemann <mariusca@simula.no>
Date: Fri, 3 Oct 2025 15:24:29 +0200
Subject: [PATCH 10/21] make remove 0 conditional

---
 workflow/scripts/process_image_data.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/workflow/scripts/process_image_data.py b/workflow/scripts/process_image_data.py
index 3b68226..a404d8f 100644
--- a/workflow/scripts/process_image_data.py
+++ b/workflow/scripts/process_image_data.py
@@ -119,7 +119,7 @@ def parse_operations(ops):
         img = da.where(da.isin(img, cois), img, 0)
     else:
         cell_labels = list(cell_labels)
-        cell_labels.remove(0)
+        if 0 in cell_labels: cell_labels.remove(0)
 
     # remap labels to smaller, sequential ints
     remapping = {int(c):i for i,c in enumerate([0] + cell_labels)}

From 87d5c689f12eb7333ade0004be1bbc96485f8325 Mon Sep 17 00:00:00 2001
From: Marius Causemann <mariusca@simula.no>
Date: Tue, 2 Dec 2025 10:52:17 +0100
Subject: [PATCH 11/21] update config files

---
 config_files/astrocyte1.yml   | 14 ++++++--------
 config_files/cappilary1.yml   |  7 ++-----
 config_files/cappilary2.yml   | 28 ++++++++++++++++------------
 config_files/cappilary3.yml   | 26 +++++++++++++++-----------
 config_files/cappilary4.yml   | 26 +++++++++++++++-----------
 config_files/cappilary5.yml   | 13 ++++++-------
 config_files/cappilary6.yml   | 24 ++++++++++++++----------
 config_files/cappilary7.yml   | 14 +++++++-------
 config_files/cortical_mm3.yml |  4 ++--
 config_files/test.yml         |  4 ++--
 10 files changed, 85 insertions(+), 75 deletions(-)

diff --git a/config_files/astrocyte1.yml b/config_files/astrocyte1.yml
index 42ca52f..11b9715 100644
--- a/config_files/astrocyte1.yml
+++ b/config_files/astrocyte1.yml
@@ -3,17 +3,15 @@ fixed :
   raw:
     cloudpath: "precomputed://gs://iarpa_microns/minnie/minnie65/seg"
     position : "189255-171091-22461" # http://tinyurl.com/ycytrm6t
+    mip : 1
   processing:
-      cells : 864691136194301772
+      operation : 
+        - "ncells ncells=1 keep_cell_labels='[864691136194301772]'"
+        - "smooth iterations=1 radius=2"
+        - "erode radius=1"
 raw:
-  mip : 1
-  size : [2500, 5000, 10000] 
+  size : [2500] #, 5000, 10000] 
 processing:
-  ncells : [ 1 ]
-  smoothiter : [1]
-  smoothradius : [2]
-  expand : [1]
-  shrink : 0
   dx : 20
 meshing:
   envelopsize : [8]
\ No newline at end of file
diff --git a/config_files/cappilary1.yml b/config_files/cappilary1.yml
index 601a865..c0bd4fa 100644
--- a/config_files/cappilary1.yml
+++ b/config_files/cappilary1.yml
@@ -4,15 +4,12 @@ fixed :
     cloudpath: "precomputed://gs://iarpa_microns/minnie/minnie65/seg"
     position : "232915-165336-20525"
   processing:
-    #roi : 864691136534887842
-    #merge :  864691136534887842 864691135476973891 864691135751005697 864691133070457251 864691136178683094 864691133070456483 864691133070450083 864691133003757151 864691136475845441 864691134032036203
-    #cells : 864691134581644554 864691135228370458 864691135918336944 864691136475874625 864691136723365629
     operation : 
       - "merge labels='[864691136534887842, 864691135476973891, 864691135751005697, 864691133070457251, 864691136178683094, 864691133070456483, 864691133070450083, 864691133003757151, 864691136475845441, 864691134032036203, 864691131901467286]'"
       - "dilate radius=4 labels='[864691134581644554, 864691135228370458, 864691135918336944, 864691136475874625, 864691136723365629, 864691134308859492]'"
       - "removeislands minsize=5000"
-      - "dilate radius=9"
-      - "smooth iterations=3 radius=5"
+      - "dilate radius=5"
+      - "smooth iterations=2 radius=4"
       - "erode radius=1"
       - "roigenerate labels='[864691136534887842]'"
       - "roidilate radius=80"
diff --git a/config_files/cappilary2.yml b/config_files/cappilary2.yml
index dfef846..b596c3b 100644
--- a/config_files/cappilary2.yml
+++ b/config_files/cappilary2.yml
@@ -4,19 +4,23 @@ fixed :
     cloudpath: "precomputed://gs://iarpa_microns/minnie/minnie65/seg"
     position : "337400-138900-22519"
   processing:
-    roi : 864691136534887842
-    merge : "864691136534887842 864691133995183345 864691134924323841 864691135197974302 864691135501607490 864691135686365079 864691135954940163 864691132805951425 864691132361462832 864691132805954753 864691132805951681 864691134924324353 864691132805952961 864691132805951169"
-    cells : 864691135211452360 864691135212672968 864691135213024712  864691135592508823 864691136215403699 #864691133186201551
+    operation:
+      - "merge labels='[864691136534887842, 864691133995183345, 864691134924323841, 864691135197974302, 864691135501607490, 864691135686365079, 864691135954940163, 864691132805951425, 864691132361462832, 864691132805954753, 864691132805951681, 864691134924324353, 864691132805952961, 864691132805951169]'"
+      - "dilate radius=4 labels='[864691135211452360, 864691135212672968, 864691135213024712,  864691135592508823, 864691136215403699]'"
+      - "removeislands minsize=5000"
+      - "dilate radius=5"
+      - "smooth iterations=2 radius=4"
+      - "erode radius=1"
+      - "roigenerate labels='[864691136534887842]'"
+      - "roidilate radius=64"
+      - "roiapply"
+      - "erode radius=15 labels='[864691136534887842]'"
+      - "removeislands minsize=5000"
+      - "roierode radius=4"
 raw:
-  mip : 1
+  mip : 0
   size : [7500] 
 processing:
-  ncells : [ 1000 ]
-  smoothiter : [3]
-  smoothradius : [3]
-  expand : [3]
-  shrink : 1
-  dx : 20
-  roidilate: "6-6"
+  dx : 8
 meshing:
-  envelopsize : [8]
+  envelopsize : [5]
diff --git a/config_files/cappilary3.yml b/config_files/cappilary3.yml
index ea6771e..8b38c82 100644
--- a/config_files/cappilary3.yml
+++ b/config_files/cappilary3.yml
@@ -4,18 +4,22 @@ fixed :
     cloudpath: "precomputed://gs://iarpa_microns/minnie/minnie65/seg"
     position : "307339-124630-24240"
   processing:
-    roi : 864691136534887842
-    merge : "864691136534887842 864691135941292148 864691132616222812" # 864691135322955164"
+      operation : 
+      - "merge labels='[864691136534887842, 864691135941292148, 864691132616222812]'"
+      - "removeislands minsize=5000"
+      - "dilate radius=5"
+      - "smooth iterations=2 radius=4"
+      - "erode radius=1"
+      - "roigenerate labels='[864691136534887842]'"
+      - "roidilate radius=64"
+      - "roiapply"
+      - "erode radius=15 labels='[864691136534887842]'"
+      - "removeislands minsize=5000"
+      - "roierode radius=4"
 raw:
-  mip : 1
+  mip : 0
   size : [7500] 
 processing:
-  ncells : [ 1000 ]
-  smoothiter : [3]
-  smoothradius : [3]
-  expand : [3]
-  shrink : 1
-  dx : 20
-  roidilate: "6-6"
+  dx : 8
 meshing:
-  envelopsize : [8]
\ No newline at end of file
+  envelopsize : [5]
\ No newline at end of file
diff --git a/config_files/cappilary4.yml b/config_files/cappilary4.yml
index 77da4c7..178b9c6 100644
--- a/config_files/cappilary4.yml
+++ b/config_files/cappilary4.yml
@@ -4,18 +4,22 @@ fixed :
     cloudpath: "precomputed://gs://iarpa_microns/minnie/minnie65/seg"
     position : "306009-143434-23494"
   processing:
-    roi : 864691136534887842
-    merge : "864691136534887842 864691133311938372 864691133711003153 864691135468246461" #864691132824819233 864691133334749831
+    operation : 
+      - "merge labels='[864691136534887842, 864691133311938372, 864691133711003153, 864691135468246461]'"
+      - "removeislands minsize=5000"
+      - "dilate radius=5"
+      - "smooth iterations=2 radius=4"
+      - "erode radius=1"
+      - "roigenerate labels='[864691136534887842]'"
+      - "roidilate radius=64"
+      - "roiapply"
+      - "erode radius=15 labels='[864691136534887842]'"
+      - "removeislands minsize=5000"
+      - "roierode radius=4"
 raw:
-  mip : 1
+  mip : 0
   size : [7500] 
 processing:
-  ncells : [ 1000 ]
-  smoothiter : [3]
-  smoothradius : [3]
-  expand : [3]
-  shrink : 1
-  dx : 20
-  roidilate: "6-6"
+  dx : 8
 meshing:
-  envelopsize : [8]
+  envelopsize : [5]
diff --git a/config_files/cappilary5.yml b/config_files/cappilary5.yml
index a5c0631..e33236c 100644
--- a/config_files/cappilary5.yml
+++ b/config_files/cappilary5.yml
@@ -4,19 +4,18 @@ fixed :
     cloudpath: "precomputed://gs://iarpa_microns/minnie/minnie65/seg"
     position : "232042-166129-21467"
   processing:
-    #cells : 864691132723008555 864691134309135460 864691134310007140 864691135609682183 864691135753683149 864691135918336944
     operation : 
       - "merge labels='[864691136534887842, 864691132131484059, 864691132336155009, 864691132348210676, 864691135791214365, 864691135971392446]'"
       - "removeislands minsize=5000"
-      - "dilate iterations=3 radius=3"
-      - "smooth iterations=3 radius=5"
-      - "erode iterations=1 radius=1"
+      - "dilate radius=5"
+      - "smooth iterations=2 radius=4"
+      - "erode radius=1"
       - "roigenerate labels='[864691136534887842]'"
-      - "roidilate iterations=8 radius=8"
+      - "roidilate radius=64"
       - "roiapply"
-      - "erode iterations=3 radius=5 labels='[864691136534887842]'"
+      - "erode radius=15 labels='[864691136534887842]'"
       - "removeislands minsize=5000"
-      - "roierode iterations=2 radius=2"
+      - "roierode radius=4"
 raw:
   mip : 0
   size : [7500] 
diff --git a/config_files/cappilary6.yml b/config_files/cappilary6.yml
index fd911f1..ff3e187 100644
--- a/config_files/cappilary6.yml
+++ b/config_files/cappilary6.yml
@@ -4,19 +4,23 @@ fixed :
     cloudpath: "precomputed://gs://iarpa_microns/minnie/minnie65/seg"
     position : "230671-173653-22233"
   processing:
-    roi : 864691136534887842
-    merge : 864691135776576352 864691135807087133 864691136534887842
-    cells : 864691135000029429 864691135000202997 864691135609682183 864691135855609134
+    operation : 
+      - "merge labels='[864691136534887842, 864691135776576352, 864691135807087133]'"
+      - "dilate radius=4 labels='[864691135000029429, 864691135000202997, 864691135609682183, 864691135855609134]'"
+      - "removeislands minsize=5000"
+      - "dilate radius=5"
+      - "smooth iterations=2 radius=4"
+      - "erode radius=1"
+      - "roigenerate labels='[864691136534887842]'"
+      - "roidilate radius=80"
+      - "roiapply"
+      - "erode radius=15 labels='[864691136534887842]'"
+      - "removeislands minsize=5000"
+      - "roierode radius=4"
 raw:
   mip : 0
   size : [7500] 
 processing:
-  ncells : [ 1000 ]
-  smoothiter : [3]
-  smoothradius : [3]
-  expand : [3]
-  shrink : 1
   dx : 8
-  roidilate: "8-8"
 meshing:
-  envelopsize : [8]
\ No newline at end of file
+  envelopsize : [5]
\ No newline at end of file
diff --git a/config_files/cappilary7.yml b/config_files/cappilary7.yml
index 4dcbe4a..644294a 100644
--- a/config_files/cappilary7.yml
+++ b/config_files/cappilary7.yml
@@ -6,17 +6,17 @@ fixed :
   processing:
     operation : 
       - "merge labels='[864691136534887842, 864691135662897092, 864691132540602862]'"
-      - "dilate iterations=2 radius=2 labels='[864691135097319093, 864691135097324469, 864691135181838082, 864691135472006194, 864691135777250237, 864691136194177256]'"
+      - "dilate radius=4 labels='[864691135097319093, 864691135097324469, 864691135181838082, 864691135472006194, 864691135777250237, 864691136194177256]'"
       - "removeislands minsize=5000"
-      - "dilate iterations=3 radius=3"
-      - "smooth iterations=3 radius=5"
-      - "erode iterations=2 radius=1"
+      - "dilate radius=5"
+      - "smooth iterations=2 radius=4"
+      - "erode radius=2"
       - "roigenerate labels='[864691136534887842]'"
-      - "roidilate iterations=8 radius=8"
+      - "roidilate radius=64"
       - "roiapply"
-      - "erode iterations=3 radius=5 labels='[864691136534887842]'"
+      - "erode radius=15 labels='[864691136534887842]'"
       - "removeislands minsize=5000"
-      - "roierode iterations=2 radius=2"
+      - "roierode radius=4"
 raw:
   mip : 0
   size : "7500-7500-5000"
diff --git a/config_files/cortical_mm3.yml b/config_files/cortical_mm3.yml
index 7a63f7b..2feabbd 100644
--- a/config_files/cortical_mm3.yml
+++ b/config_files/cortical_mm3.yml
@@ -7,9 +7,9 @@ fixed :
   processing:
     operation : 
       - "removeislands minsize=5000"
-      - "dilate iterations=1 radius=1"
+      - "dilate radius=1"
       - "smooth iterations=1 radius=1"
-      - "erode iterations=1 radius=1"
+      - "erode radius=1"
       - "removeislands minsize=5000"
 raw:
   size : [10000] #[5000, 10000, 20000]
diff --git a/config_files/test.yml b/config_files/test.yml
index 1448150..0ea8f35 100644
--- a/config_files/test.yml
+++ b/config_files/test.yml
@@ -7,8 +7,8 @@ fixed :
     operation : 
       - "merge labels='[864691136534887842, 864691132131484059, 864691132336155009, 864691132348210676, 864691135791214365, 864691135971392446]'"
       - "removeislands minsize=5000"
-      - "dilate radius=9"
-      - "smooth iterations=3 radius=5"
+      - "dilate radius=5"
+      - "smooth iterations=2 radius=4"
       - "erode radius=1"
       - "roigenerate labels='[864691136534887842]'"
       - "roidilate radius=24"

From cd4781dce9c17f8a1de1ac391e1f58c17f016913 Mon Sep 17 00:00:00 2001
From: Marius Causemann <mariusca@simula.no>
Date: Tue, 2 Dec 2025 10:52:38 +0100
Subject: [PATCH 12/21] remove surface tags

---
 workflow/Snakefile | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/workflow/Snakefile b/workflow/Snakefile
index 163c422..64b0559 100644
--- a/workflow/Snakefile
+++ b/workflow/Snakefile
@@ -67,7 +67,7 @@ rule all:
     input:
         expand(mesh_img, rawdata=rd,processing=pd, meshing=md, name=name),  
         expand(meshstatistic, rawdata=rd,processing=pd, meshing=md, name=name),
-        expand(mesh_tag_path, rawdata=rd,processing=pd, meshing=md, name=name),
+        #expand(mesh_tag_path, rawdata=rd,processing=pd, meshing=md, name=name),
 
 rule downloadImageData:
     output:
@@ -91,11 +91,10 @@ rule processImageData:
         #imagestatistic,
         outfile=processed_data_path,
     conda:
-        "envs/nbmorph_env.yml"
+        "envs/environment.yml"
     resources:
-        ntasks=10,
+        ntasks=4,
         time=30,
-        gpus=1,
     params:
         options=lambda wildcards: filename2options(wildcards.processing),
         fopt=lambda wildcards: name2options(wildcards.name, "processing")

From fcf764e6bf80172bbfcbceef4dae5a080bc19a52 Mon Sep 17 00:00:00 2001
From: Marius Causemann <mariusca@simula.no>
Date: Tue, 2 Dec 2025 10:52:49 +0100
Subject: [PATCH 13/21] cleanup env

---
 workflow/envs/environment.yml | 11 ++---------
 workflow/envs/nbmorph_env.yml | 13 -------------
 2 files changed, 2 insertions(+), 22 deletions(-)
 delete mode 100644 workflow/envs/nbmorph_env.yml

diff --git a/workflow/envs/environment.yml b/workflow/envs/environment.yml
index cd55ccf..4a61be7 100644
--- a/workflow/envs/environment.yml
+++ b/workflow/envs/environment.yml
@@ -2,27 +2,20 @@ name: emimesh
 channels:
   - conda-forge
 dependencies:
-  - python=3.10
-  - numpy
-  - vtk=*=*mesa*
   - matplotlib
   - pyvista
-  - pyacvd
   - scikit-image
-  - fenics
   - meshio
-  - dask
-  - k3d
-  - flask
   - dask-image
   - seaborn
+  - pyacvd
   - zip
   - pip
   - pip:
-    - pyacvd>=0.3
     - pytetwild
     - fastremap
     - cloud-volume
     - webknossos
     - connected-components-3d
     - cmocean
+    - nbmorph
diff --git a/workflow/envs/nbmorph_env.yml b/workflow/envs/nbmorph_env.yml
deleted file mode 100644
index f8030ea..0000000
--- a/workflow/envs/nbmorph_env.yml
+++ /dev/null
@@ -1,13 +0,0 @@
-name: intel
-channels:
-  - conda-forge
-dependencies:
-  - dask-image
-  - pyvista
-  - ipython
-  - numba
-  - pip
-  - pip:
-      - connected-components-3d
-      - fastremap
-      - nbmorph
\ No newline at end of file

From 74c89e2b08b8669fb7a6de118bf0f4d0de5b4bc7 Mon Sep 17 00:00:00 2001
From: Marius Causemann <mariusca@simula.no>
Date: Tue, 2 Dec 2025 10:53:02 +0100
Subject: [PATCH 14/21] remive k3d

---
 workflow/scripts/generate_screenshot.py | 22 +++++++++++-----------
 1 file changed, 11 insertions(+), 11 deletions(-)

diff --git a/workflow/scripts/generate_screenshot.py b/workflow/scripts/generate_screenshot.py
index 81bad7b..5f4479d 100644
--- a/workflow/scripts/generate_screenshot.py
+++ b/workflow/scripts/generate_screenshot.py
@@ -1,6 +1,5 @@
 import argparse
 import pyvista as pv
-import k3d
 from plot_utils import get_screenshot
 import numpy as np
 import matplotlib
@@ -48,19 +47,20 @@
     
     newcmap = cmocean.tools.crop_by_percent(cmocean.cm.curl_r, 5, which='max', N=None)
     get_screenshot(cells, args.output, scalar="label", cmap="curl")
+    
+    #import k3d
+    #ecs = mesh.extract_cells((mesh["label"]==1))
 
-    ecs = mesh.extract_cells((mesh["label"]==1))
-
-    color_map =k3d.paraview_color_maps.Rainbow_Uniform
-    color_map = k3d.paraview_color_maps.Linear_Green_Gr4L
-    cells = cells.cell_data_to_point_data()
-    cells_k3d = k3d.vtk_poly_data(cells.extract_surface(), side="double",
+    #color_map =k3d.paraview_color_maps.Rainbow_Uniform
+    #color_map = k3d.paraview_color_maps.Linear_Green_Gr4L
+    #cells = cells.cell_data_to_point_data()
+    #cells_k3d = k3d.vtk_poly_data(cells.extract_surface(), side="double",
                               #color_attribute=("label", 0, float(cells["label"].max())), 
                               #color_map=color_map, name="cells"
-                              color=hexcolor("limegreen"))
-    ecs_k3d = k3d.vtk_poly_data(ecs.extract_surface(), side="double", color=hexcolor("white"),
-                                opacity=0.8,
-                                name="ecs", wireframe=True)
+    #                          color=hexcolor("limegreen"))
+    #ecs_k3d = k3d.vtk_poly_data(ecs.extract_surface(), side="double", color=hexcolor("white"),
+    #                            opacity=0.8,
+    #                            name="ecs", wireframe=True)
 
     #generate_screenshots([cells_k3d, ecs_k3d], args.output, fov=32)
 

From f62e8973d33b52ab231172519d15316b936ce4cf Mon Sep 17 00:00:00 2001
From: Marius Causemann <mariusca@simula.no>
Date: Tue, 2 Dec 2025 10:53:33 +0100
Subject: [PATCH 15/21] add keep_cell_labels to ncells

---
 workflow/scripts/process_image_data.py | 28 +++++++++++++-------------
 1 file changed, 14 insertions(+), 14 deletions(-)

diff --git a/workflow/scripts/process_image_data.py b/workflow/scripts/process_image_data.py
index a404d8f..ca684fe 100644
--- a/workflow/scripts/process_image_data.py
+++ b/workflow/scripts/process_image_data.py
@@ -1,7 +1,6 @@
 import numpy as np
 import argparse
 import yaml
-import json
 import time
 import pyvista as pv
 from pathlib import Path
@@ -15,7 +14,6 @@
 import numba
 import nbmorph
 
-numba.set_num_threads(1)
 dask.config.set({"array.chunk-size": "1024 MiB"})
 
 def mergecells(img, labels):
@@ -23,10 +21,14 @@ def mergecells(img, labels):
     img = np.where(np.isin(img, labels), labels[0], img)
     return img
 
-def ncells(img, ncells):
+def ncells(img, ncells, keep_cell_labels=None):
     cell_labels, cell_counts = fastremap.unique(img, return_counts=True)
     cell_labels = cell_labels[np.argsort(cell_counts)]
-    cois = list(cell_labels[-ncells :])
+    if keep_cell_labels is None: cois =[]
+    cois = set(keep_cell_labels)
+    for cid in cell_labels:
+        if len(cois) >= ncells: break
+        cois.add(cid)
     img = np.where(np.isin(img, cois), img, 0)
     return img
     
@@ -99,6 +101,7 @@ def parse_operations(ops):
 
     args = parser.parse_args()
     n_parallel = args.nworkers
+    numba.set_num_threads(n_parallel)
     print(f"Using {n_parallel} workers...")
     start = time.time()
 
@@ -142,12 +145,13 @@ def parse_operations(ops):
     operations = parse_operations(args.operation)
     for op, kwargs in operations:
         print(op, kwargs)
-        if "labels" in kwargs.keys():
-            labels = kwargs["labels"]
-            if kwargs.get("allexcept", False):
-                kwargs["labels"] = list(set(remapping.values()) - set(remap(labels)))
-            else:
-                kwargs["labels"] = remap(labels)
+        for k in kwargs.keys():
+            if "label" in k:
+                labels = kwargs[k]
+                if kwargs.get("allexcept", False):
+                    kwargs[k] = list(set(remapping.values()) - set(remap(labels)))
+                else:
+                    kwargs[k] = remap(labels)
 
         if op=="roigenerate":
             roi = np.isin(img, kwargs["labels"])
@@ -155,16 +159,12 @@ def parse_operations(ops):
         if op=="roiapply":
             img = np.where(roi, img, 0) 
             continue
-        if op=="ncells":
-            img = ncells(img, **kwargs)
-            continue
         if op.startswith("roi"):
             roiop = op[3:]
             roi = opdict[roiop](roi, **kwargs)
         else:
             img =opdict[op](img, **kwargs)
     
-    #img, roi = dask.compute(img, roi, num_workers= args.nworkers)
     print(f"processed! {img.shape}")
     proc_time = time.time()
 

From b2f80d85ec78d60ec7673aa38788d54ff62854b9 Mon Sep 17 00:00:00 2001
From: Marius Causemann <mariusca@simula.no>
Date: Tue, 2 Dec 2025 11:06:23 +0100
Subject: [PATCH 16/21] fix typo

---
 workflow/scripts/download_data.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/workflow/scripts/download_data.py b/workflow/scripts/download_data.py
index a188d93..832eb7a 100644
--- a/workflow/scripts/download_data.py
+++ b/workflow/scripts/download_data.py
@@ -28,7 +28,7 @@ def download_cloudvolume(cloud_path, mip, pos, physical_size):
     vol = CloudVolume(
         cloud_path, use_https=True, parallel=8, progress=True, mip=mip, cache=True, bounded=True
     )
-    print(f"data resoltion: {vol.resolution}")
+    print(f"data resolution: {vol.resolution}")
     size = [ps / res for ps, res in zip(physical_size, vol.resolution)]
     size = np.array(size).astype("uint64")
 

From 96795876110957ea788b6801c0acf6c80c10a052 Mon Sep 17 00:00:00 2001
From: Marius Causemann <mariusca@simula.no>
Date: Tue, 2 Dec 2025 11:06:44 +0100
Subject: [PATCH 17/21] update CI

---
 .github/workflows/test_conda.yml | 20 ++++++++------------
 1 file changed, 8 insertions(+), 12 deletions(-)

diff --git a/.github/workflows/test_conda.yml b/.github/workflows/test_conda.yml
index 37e5d0e..7a9bd3e 100644
--- a/.github/workflows/test_conda.yml
+++ b/.github/workflows/test_conda.yml
@@ -16,20 +16,16 @@ defaults:
 
 jobs:
   test_scripts:
-    runs-on: ubuntu-latest
-    container: ubuntu:latest
+    name: test on ${{ matrix.os }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        os: [ubuntu-latest, windows-latest, macos-15-intel, macos-15]
+
 
     steps:
       - uses: actions/checkout@v4
-
-      #   - uses: actions/setup-python@v5
-      #     with:
-      #       python-version: "3.11"
-      - name: Install deps
-        run: |
-          DEBIAN_FRONTEND=noninteractive apt-get update
-          DEBIAN_FRONTEND=noninteractive apt-get install -y gcc unzip make cmake build-essential git libgmp-dev xvfb intel-opencl-icd pocl-opencl-icd
-      #       python -m pip install snakemake
       - name: Setup conda-forge
         uses: conda-incubator/setup-miniconda@v3
         with:
@@ -43,4 +39,4 @@ jobs:
           conda install snakemake snakemake-storage-plugin-http snakemake-executor-plugin-cluster-generic
       - name: Run snakemake
         run: |
-          export PYOPENCL_COMPILER_OUTPUT=1 && snakemake --use-conda --cores 2 -p --configfile ./config_files/test.yml
+          snakemake --use-conda --cores 2 -p --configfile ./config_files/test.yml

From 6022b63854220a65c1cb01d0e8b11f6f7b361db8 Mon Sep 17 00:00:00 2001
From: Marius Causemann <mariusca@simula.no>
Date: Tue, 2 Dec 2025 11:13:48 +0100
Subject: [PATCH 18/21] remove zip, update ntasks

---
 workflow/Snakefile            | 6 ++----
 workflow/envs/environment.yml | 1 -
 2 files changed, 2 insertions(+), 5 deletions(-)

diff --git a/workflow/Snakefile b/workflow/Snakefile
index 64b0559..a4e3e15 100644
--- a/workflow/Snakefile
+++ b/workflow/Snakefile
@@ -30,8 +30,6 @@ mesh_path = mesh_dir + "mesh.xdmf"
 mesh_tag_path = mesh_dir + "facets.xdmf"
 meshstatistic = mesh_dir + "meshstatistic.yml"
 mesh_img = mesh_dir + "mesh.png"
-zipped = "zipped/{name}/{meshing}/volmesh_{rawdata}-{processing}.zip"
-zipped_surf = "zipped_surf/{name}/{meshing}/surfaces_{rawdata}-{processing}.zip"
 surfaces_dir = directory(surf_dir)
 
 
@@ -93,7 +91,7 @@ rule processImageData:
     conda:
         "envs/environment.yml"
     resources:
-        ntasks=4,
+        ntasks=lambda wildcards: size_to_ntasks(int(wildcards.rawdata.split("size+")[-1].split("_")[0]) / int(wildcards.processing.split("dx+")[-1].split("_")[0])),
         time=30,
     params:
         options=lambda wildcards: filename2options(wildcards.processing),
@@ -114,7 +112,7 @@ rule extractSurfaces:
         csgtree=csgtree
     resources:
         time=120,
-        ntasks=4
+        ntasks=lambda wildcards: size_to_ntasks(int(wildcards.rawdata.split("size+")[-1].split("_")[0]) / int(wildcards.processing.split("dx+")[-1].split("_")[0])),
     conda:
         conda_env
     shell:
diff --git a/workflow/envs/environment.yml b/workflow/envs/environment.yml
index 4a61be7..b540a9f 100644
--- a/workflow/envs/environment.yml
+++ b/workflow/envs/environment.yml
@@ -9,7 +9,6 @@ dependencies:
   - dask-image
   - seaborn
   - pyacvd
-  - zip
   - pip
   - pip:
     - pytetwild

From a48765998e19694d854a3d0901c61e89692f5807 Mon Sep 17 00:00:00 2001
From: Marius Causemann <mariusca@simula.no>
Date: Tue, 2 Dec 2025 11:52:25 +0100
Subject: [PATCH 19/21] install osmesa

---
 .github/workflows/test_conda.yml | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/.github/workflows/test_conda.yml b/.github/workflows/test_conda.yml
index 7a9bd3e..687bd33 100644
--- a/.github/workflows/test_conda.yml
+++ b/.github/workflows/test_conda.yml
@@ -26,6 +26,13 @@ jobs:
 
     steps:
       - uses: actions/checkout@v4
+
+      - name: Install Headless Plotting Libs
+        if: runner.os == 'Linux'
+        run: |
+          sudo apt-get update
+          sudo apt-get install -y libosmesa6 libgl1-mesa-glx
+
       - name: Setup conda-forge
         uses: conda-incubator/setup-miniconda@v3
         with:

From 5ec7480acff27364a455a904352fe53b41a8c6e5 Mon Sep 17 00:00:00 2001
From: Marius Causemann <mariusca@simula.no>
Date: Tue, 2 Dec 2025 11:52:48 +0100
Subject: [PATCH 20/21] python3 -> python

---
 workflow/Snakefile | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/workflow/Snakefile b/workflow/Snakefile
index a4e3e15..8a27aae 100644
--- a/workflow/Snakefile
+++ b/workflow/Snakefile
@@ -77,7 +77,7 @@ rule downloadImageData:
         fopt=lambda wildcards: name2options(wildcards.name, "raw")
     shell:
         """
-        python3 workflow/scripts/download_data.py \
+        python workflow/scripts/download_data.py \
         {params.options} {params.fopt} \
         --output {output.rawdata}
         """
@@ -98,7 +98,7 @@ rule processImageData:
         fopt=lambda wildcards: name2options(wildcards.name, "processing")
     shell:
         """
-        python3 workflow/scripts/process_image_data.py --infile {input.rawdata} \
+        python workflow/scripts/process_image_data.py --infile {input.rawdata} \
         {params.options} {params.fopt} \
         --output {output.outfile} \
         --nworkers {resources.ntasks}
@@ -117,7 +117,7 @@ rule extractSurfaces:
         conda_env
     shell:
         """
-        python3 workflow/scripts/extract_surfaces.py --infile {input.processeddata} --ncpus {resources.ntasks} \
+        python workflow/scripts/extract_surfaces.py --infile {input.processeddata} --ncpus {resources.ntasks} \
         --outdir {output.outdir}
         """
 
@@ -136,14 +136,14 @@ rule generateMesh:
         options=lambda wildcards: filename2options(wildcards.meshing)
     shell:
         """
-        python3 workflow/scripts/generate_mesh.py \
+        python workflow/scripts/generate_mesh.py \
         --csgtree {input.csgtree} \
         {params.options} \
         --output {output.outfile} \
         --max_threads {resources.ntasks} \
         """
 
-#   python3 workflow/scripts/seperate_touching_cells.py --infile {output.outfile} --output {output.outfile}
+#   python workflow/scripts/seperate_touching_cells.py --infile {output.outfile} --output {output.outfile}
 
 
 rule generateSurfaceTags:
@@ -157,7 +157,7 @@ rule generateSurfaceTags:
         time="120:00"
     shell:
         """
-        python3 workflow/scripts/generate_surface_tags.py \
+        python workflow/scripts/generate_surface_tags.py \
         --infile {input.meshfile} --output {output}
         """
 
@@ -170,7 +170,7 @@ rule evaluateMesh:
         conda_env
     shell:
         """
-        python3 workflow/scripts/evaluate_mesh.py \
+        python workflow/scripts/evaluate_mesh.py \
         --infile {input} --output {output}
         """
 
@@ -183,7 +183,7 @@ rule takeScreenshot:
         conda_env
     shell:
         """
-        python3 workflow/scripts/generate_screenshot.py \
+        python workflow/scripts/generate_screenshot.py \
         --infile {input} --output {output}
         """
 
@@ -196,7 +196,7 @@ rule generateAnalysisPlot:
         conda_env
     shell:
         """
-        python3 workflow/scripts/generate_analysis_plots.py \
+        python workflow/scripts/generate_analysis_plots.py \
         --infile {input.infile} --output {output.plotfile}
         """
 

From 7a8d6b006e5530ef2369bf7a0dfff48a71ab7779 Mon Sep 17 00:00:00 2001
From: Marius Causemann <mariusca@simula.no>
Date: Tue, 2 Dec 2025 11:54:19 +0100
Subject: [PATCH 21/21] fix libgl

---
 .github/workflows/test_conda.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/test_conda.yml b/.github/workflows/test_conda.yml
index 687bd33..433dafb 100644
--- a/.github/workflows/test_conda.yml
+++ b/.github/workflows/test_conda.yml
@@ -31,7 +31,7 @@ jobs:
         if: runner.os == 'Linux'
         run: |
           sudo apt-get update
-          sudo apt-get install -y libosmesa6 libgl1-mesa-glx
+          sudo apt-get install -y libosmesa6 libgl1
 
       - name: Setup conda-forge
         uses: conda-incubator/setup-miniconda@v3