Refactor code into module

.github/workflows/test_conda.yml

@@ -21,7 +21,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        os: [ubuntu-latest, windows-latest, macos-15-intel, macos-15]
+        os: [ubuntu-latest, windows-latest, macos-15]
 
 
     steps:
@@ -57,4 +57,4 @@ jobs:
 
       - name: Run snakemake
         run: |
-          snakemake --cores 2 -p --configfile ./config_files/motta2019.yml
+          snakemake --cores 2 -p --configfile ./config_files/test.yml

pyproject.toml

@@ -0,0 +1,90 @@
+[build-system] # Require setuptool version due to https://github.com/pypa/setuptools/issues/2938
+requires = ["setuptools>=61.0.0", "wheel"]
+
+[project]
+dependencies = [
+    "matplotlib",
+    "pyvista",
+    "scikit-image",
+    "meshio",
+    "dask-image",
+    "seaborn",
+    "pyacvd",
+    "pytetwild",
+    "fastremap",
+    "cloud-volume",
+    "connected-components-3d",
+    "cmocean",
+    "nbmorph",
+]
+name = "EMIMesh"
+version = "1.0.0.dev0"
+description = "Generating high-quality extracellular-membrane-intracellular meshes from imaging data"
+authors = [{ name = "Marius Causemann", email = "mariusca@simula.no" }]
+license = { file = "LICENSE" }
+readme = "readme.md"
+
+
+[project.scripts]
+emi-generate-mesh = "emimesh.clients.generate_mesh:main"
+emi-download-data = "emimesh.clients.download_data:main"
+emi-evaluate-mesh = "emimesh.clients.evaluate_mesh:main"
+emi-extract-surfaces = "emimesh.clients.extract_surfaces:main"
+emi-plot-analysis = "emimesh.clients.analysis_plot:main"
+emi-process-image-data = "emimesh.clients.process_image_data:main"
+emi-generate-screenshot = "emimesh.clients.generate_screenshot:main"
+
+[project.optional-dependencies]
+test = []
+dev = ["pdbpp", "ipython", "mypy", "ruff"]
+webknossos = ["webknossos"]
+docs = [
+    "jupyter-book<2.0.0",
+    "jupytext",
+    "ipykernel<7.0.0",     # Note: Remove once https://github.com/ipython/ipykernel/issues/1450 is in a release
+    "sphinx-codeautolink",
+]
+all = ["EMIMesh[test,dev,docs]"]
+
+[tool.pytest.ini_options]
+addopts = ["--import-mode=importlib"]
+testpaths = ["tests"]
+
+[tool.mypy]
+ignore_missing_imports = true
+# Folders to exclude
+exclude = ["docs/", "build/"]
+# Folder to check with mypy
+files = ["src", "tests"]
+
+[tool.ruff]
+src = ["src", "tests", "docs"]
+line-length = 100
+indent-width = 4
+
+[tool.ruff.lint]
+select = [
+    # Pyflakes
+    "F",
+    # Pycodestyle
+    "E",
+    "W",
+    # isort
+    "I001",
+]
+
+
+[tool.ruff.lint.isort]
+known-first-party = ["EMIMesh"]
+known-third-party = ["numpy", "pytest"]
+section-order = [
+    "future",
+    "standard-library",
+    "mpi",
+    "third-party",
+    "first-party",
+    "local-folder",
+]
+
+[tool.ruff.lint.isort.sections]
+"mpi" = ["mpi4py", "petsc4py"]

src/emimesh/clients/analysis_plot.py

@@ -0,0 +1,22 @@
+import argparse
+import yaml
+from emimesh.generate_analysis_plots import plot_cell_sizes
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        "--infile",
+        help="path to mesh statistics infile",
+        type=str,
+    )
+    parser.add_argument(
+        "--output",
+        help="outfile name",
+        type=str,
+    )
+    args = parser.parse_args()
+
+    with open(args.infile) as infile:
+        mesh_statistic = yaml.load(infile, Loader=yaml.FullLoader)
+
+    plot_cell_sizes(mesh_statistic, args.output)

src/emimesh/clients/download_data.py

@@ -0,0 +1,53 @@
+import argparse
+from pathlib import Path
+from emimesh.utils import np2pv
+from pathlib import Path
+import argparse
+from emimesh.download_data import download_webknossos, download_cloudvolume
+
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        "--cloudpath",
+        help="path to cloud data",
+        type=str,
+        default="precomputed://gs://iarpa_microns/minnie/minnie65/seg",
+    )
+    parser.add_argument("--mip", help="resolution (0 is highest)", type=int, default=0)
+    parser.add_argument(
+        "--position",
+        help="point position in x-y-z integer pixel position \
+              (can be copied from neuroglancer)",
+        type=str,
+    )
+    parser.add_argument(
+        "--size",
+        help="cube side length of the volume to be downloaded (in nm)",
+        type=str,
+        default=1000,
+    )
+    parser.add_argument(
+        "--output", help="output filename", type=str, default="data.xdmf"
+    )
+
+    args = parser.parse_args()
+
+    position = args.position.split("-")
+    try:
+        size = [int(args.size)] * 3
+    except ValueError:
+        size = [int(s) for s in args.size.split("-")]
+
+    try:
+        img,res = download_cloudvolume(args.cloudpath, args.mip, position, size)
+    except:
+        img,res = download_webknossos(args.cloudpath, args.mip, position, size)
+
+    print(res)
+    data = np2pv(img, res)
+    Path(args.output).parent.mkdir(exist_ok=True, parents=True)
+    data.save(args.output)
+
+
+if __name__ == "__main__":
+    main()

workflow/scripts/evaluate_mesh.py → src/emimesh/clients/evaluate_mesh.py

@@ -1,63 +1,11 @@
+import argparse
 import pyvista as pv
 import numpy as np
-import argparse
-import yaml
-import matplotlib.pyplot as plt
-import seaborn as sns
-
-lstr = "label"
-ecs_id = 1
-
-
-def compute_local_width(mesh, ecs_id, cell_ids):
-    ecs = mesh.extract_cells(np.isin(mesh.cell_data[lstr], ecs_id))
-    distances = []
-    for cid in cell_ids:
-        cell = mesh.extract_cells(np.isin(mesh.cell_data[lstr], [cid]))
-        dist = ecs.compute_implicit_distance(cell.extract_surface())
-        distances.append(abs(dist.point_data["implicit_distance"]))
-
-    dist = ecs.compute_implicit_distance(mesh.extract_surface())
-    distances.append(abs(dist.point_data["implicit_distance"]))
-
-    distances = np.array(distances).T
-    distances.sort()
-    min_dist = distances[:, 0] + distances[:, 1]
-    ecs["local_width"] = min_dist
-    ecs = ecs.point_data_to_cell_data()
-    ecs = ecs.compute_cell_sizes()
-    return ecs.cell_data["local_width"], abs(ecs.cell_data["Volume"])
 
-def compute_surface_volume(mesh, cell_ids):
-    mesh = mesh.compute_cell_sizes()
-    mesh["Volume"] = np.abs(mesh["Volume"])
-    assert (mesh["Volume"] > 0).all()
-    volumes, surface_areas = [], []
-    for cid in cell_ids:
-        cell = mesh.extract_cells(np.isin(mesh.cell_data[lstr], [cid]))
-        surf = cell.extract_surface()
-        surface_areas.append(float(surf.compute_cell_sizes()["Area"].sum()))
-        volumes.append(float(cell["Volume"].sum()))
-    return volumes, surface_areas
-
-
-def plot_local_width(width, volume, filename):
-    plt.figure(dpi=300)
-    sns.histplot(
-        x=width,
-        weights=volume,
-        bins=40,
-        kde=True,
-        stat="percent",
-        edgecolor="white",
-        kde_kws={"bw_adjust": 1},
-    )
-    plt.xlabel("local width (nm)")
-    plt.ylabel("relative frequency (%)")
-    plt.tight_layout()
-    plt.savefig(filename)
+import yaml
+from emimesh.evaluate_mesh import compute_surface_volume, ecs_id, lstr
 
-if __name__ == "__main__":
+def main():
     parser = argparse.ArgumentParser()
     parser.add_argument(
         "--infile",
@@ -108,4 +56,8 @@ def plot_local_width(width, volume, filename):
     with open(args.output, "w") as outfile:
         yaml.dump(results, outfile)
 
-    #plot_local_width(width, volume, Path(args.infile).parent / "local_width.png")
+    #plot_local_width(width, volume, Path(args.infile).parent / "local_width.png")
+
+
+if __name__ == "__main__":
+    main()

src/emimesh/clients/extract_surfaces.py

@@ -0,0 +1,72 @@
+import argparse
+from pathlib import Path
+import pyvista as pv
+
+from emimesh.extract_surfaces import extract_surface, clip_closed_box, extract_cell_meshes,create_balanced_csg_tree
+import numpy as np
+import json
+import fastremap
+from emimesh.utils import get_bounding_box
+import argparse
+
+
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        "--infile",
+        help="path to image data",
+        type=str,
+    )
+    parser.add_argument(
+        "--outdir", help="directory for output", type=str, default="output"
+    )
+    parser.add_argument(
+        "--ncpus", help="number oc cores, default 1", type=int, default=1
+    )
+    args = parser.parse_args()
+    outdir = Path(args.outdir)
+    print(f"reading file: {args.infile}")
+    img_grid = pv.read(args.infile)
+    dims = img_grid.dimensions
+    resolution = img_grid.spacing
+    img = img_grid["data"].reshape(dims - np.array([1, 1, 1]), order="F")
+    cell_labels, cell_counts = fastremap.unique(img, return_counts=True)
+    cell_labels = list(cell_labels[np.argsort(cell_counts)])
+    cell_labels.remove(0)
+
+    outerbox = get_bounding_box(img_grid, resolution[0]*5 + img_grid.length*0.002)
+
+    if "roimask" in img_grid.array_names:
+        roimask = img_grid["roimask"].reshape(dims - np.array([1, 1, 1]), 
+                                              order="F")
+        roipadded = np.pad(roimask, 1)
+        grid = pv.ImageData(dimensions=roipadded.shape, spacing=resolution,
+                            origin=(0, 0, 0))
+        roisurf = extract_surface(roipadded, grid, mesh_reduction_factor=10,
+                                  taubin_smooth_iter=5)
+
+        clip_closed_box(roisurf, outerbox)
+    else:
+        roisurf = outerbox
+    roi_file = outdir / "roi.ply"
+
+    surfs = extract_cell_meshes(
+        img,
+        cell_labels[::-1],
+        resolution,
+        mesh_reduction_factor=10,
+        taubin_smooth_iter=5,
+        write_dir=outdir,
+        ncpus=args.ncpus
+    )
+
+    mesh_files = [outdir / f"{cid}.ply" for cid in surfs if cid]
+    roisurf.save(roi_file)
+    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))
+
+if __name__ == "__main__":
+    main()

workflow/scripts/generate_mesh.py → src/emimesh/clients/generate_mesh.py

@@ -2,7 +2,7 @@
 import pyvista as pv
 import argparse
 import numpy as np
-import time
+from emimesh.generate_mesh import mesh_surfaces
 
 def get_values(d):
     for v in d.values():
@@ -11,19 +11,8 @@ def get_values(d):
         else:
             yield v
 
-def mesh_surfaces(csg_tree_path, eps, stop_quality, max_threads):
-    from pytetwild import tetrahedralize_csg
-    start = time.time()
-    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__":
+def main():
     parser = argparse.ArgumentParser()
     parser.add_argument(
         "--csgtree",
@@ -64,3 +53,7 @@ def mesh_surfaces(csg_tree_path, eps, stop_quality, max_threads):
     pv.save_meshio(args.output, volmesh)
     print(volmesh.array_names)
 
+
+
+if __name__ == "__main__":
+    main()