Merge pull request #33 from AllenCellModeling/julie_dev_variance_data

Julie dev variance data

Author: jcass11

.gitignore

@@ -2,6 +2,7 @@
 /local_staging
 /data
 /dask-worker-space
+/.dask_logs
 
 # notebooks bcz eww (force add them if you must?)
 *.ipynb

mti_nma/bin/all.py

@@ -11,6 +11,7 @@
 from datetime import datetime
 from pathlib import Path
 
+import dask
 from dask_jobqueue import SLURMCluster
 from distributed import LocalCluster
 from prefect import Flow
@@ -33,22 +34,22 @@ def __init__(self):
         """
         step_list = [
             [
-                steps.Singlecell(step_name=f"single_{name}"),
+                steps.LoadData(),
                 steps.Shparam(step_name=f"shparam_{name}"),
                 steps.Avgshape(step_name=f"avgshape_{name}"),
                 steps.Nma(step_name=f"nma_{name}")
             ]
-            for name in ["nuc", "cell"]
+            for name in ["nuc"]
         ]
         self.step_list = [step for sublist in step_list for step in sublist]
         self.step_list.append(steps.CompareNucCell())
 
     def run(
         self,
-        distributed: bool = True,
+        distributed: bool = False,
         clean: bool = False,
         debug: bool = False,
-        structs: list = ["Nuc", "Cell"],
+        structs: list = ["Nuc"],
         flow_viz: bool = False,
         **kwargs,
     ):
@@ -83,7 +84,7 @@ def run(
 
         # Initalize steps
         if "Nuc" in structs:
-            single_nuc = steps.Singlecell(step_name="single_nuc")
+            loaddata_nuc = steps.LoadData()
             shparam_nuc = steps.Shparam(step_name="shparam_nuc")
             avgshape_nuc = steps.Avgshape(step_name="avgshape_nuc")
             nma_nuc = steps.Nma(step_name="nma_nuc")
@@ -165,15 +166,15 @@ def run(
                 if "Nuc" in structs:
                     struct = "Nuc"
 
-                    sc_nuc_df = single_nuc(
+                    ld_nuc_df = loaddata_nuc(
                         distributed_executor_address=distributed_executor_address,
                         clean=clean,
                         debug=debug,
                         struct=struct,
                         **kwargs
                     )
                     sh_nuc_df = shparam_nuc(
-                        sc_df=sc_nuc_df,
+                        sc_df=ld_nuc_df,
                         distributed_executor_address=distributed_executor_address,
                         clean=clean,
                         debug=debug,

mti_nma/bin/make_9grid_mode_video.py

@@ -0,0 +1,29 @@
+import imageio
+import os
+
+# Steps to generate GIF of first nine modes
+# Edit filepaths in mti_nma/bin/view_modes.pvsm to match your own filepaths
+# In Paraview: File> Load State
+# Select view_modes.pvsm located in mti_nma/bin/
+# At Load State Data File Options "Use file names from state"
+# You should now see the 9-grid of modes
+# To save pngs of frames: File> Save Animation
+# Save in whatever directory you'd like, and use this as "pngdir" below
+# In the base mti_nma directory for this repo, run "python make_9grid_mode_video.png"
+# The gif should now be saved in your specified save_file_path below
+
+
+def create_gif_from_png_dir(
+        basedir="../../local_staging/nma/nma_Nuc/nma_data/9grid_frames/",
+        save_file_path="../../local_staging/nma/nma_Nuc/nma_data/9grid_videos/"):
+
+    # Create empty array
+    for dirname in os.listdir(basedir):
+        images = []
+        pngdir = basedir + dirname
+        for filename in os.listdir(pngdir):
+            images.append(imageio.imread(pngdir + "/" + filename))
+        imageio.mimsave(save_file_path + f"9grid_{dirname}_video.gif" , images)
+
+
+create_gif_from_png_dir()

mti_nma/bin/view_modes.pvsm

@@ -5,11 +5,13 @@
 from .avgshape import Avgshape
 from .nma import Nma
 from .compare_nuc_cell import CompareNucCell
+from .load_data import LoadData
 
 
 __all__ = [
     "Singlecell",
     "Shparam",
     "Avgshape",
     "Nma",
-    "CompareNucCell"]
+    "CompareNucCell",
+    "LoadData"]

mti_nma/steps/__init__.py

@@ -4,10 +4,13 @@
 
 import numpy as np
 import pandas as pd
+import os
 from aicsshparam import shtools
 
 from datastep import Step, log_run_params
 from .avgshape_utils import run_shcoeffs_analysis, save_mesh_as_stl
+from .avgshape_utils import get_smooth_and_coarse_mesh_from_voxelization
+from .avgshape_utils import save_mesh_as_obj, save_voxelization, save_displacement_map
 
 ###############################################################################
 
@@ -50,13 +53,10 @@ def run(self, sh_df=None, struct="Nuc", **kwargs):
         # If no dataframe is passed in, load manifest from previous step
         if sh_df is None:
             sh_df = pd.read_csv(
-                self.step_local_staging_dir.parent / "shparam_"
-                f"{struct}" / "manifest.csv"
+                self.step_local_staging_dir.parent / "shparam" /
+                f"shparam_{struct}" / "manifest.csv", index_col="CellId"
             )
 
-        # Fix filepaths and use cell id as dataframe index
-        sh_df = sh_df.set_index("CellId", drop=True)
-
         # Load sh coefficients of all samples in manifest
         coeffs_df = pd.DataFrame([])
         for CellId in sh_df.index:
@@ -66,9 +66,20 @@ def run(self, sh_df=None, struct="Nuc", **kwargs):
             )
 
         # Create directory to hold results from this step
-        avg_data_dir = self.step_local_staging_dir / "avgshape_data"
+        struct_dir = self.step_local_staging_dir / f"avgshape_{struct}"
+        struct_dir.mkdir(parents=True, exist_ok=True)
+
+        avg_data_dir = struct_dir / f"avgshape_data"
         avg_data_dir.mkdir(parents=True, exist_ok=True)
 
+        # Move init and run parameters to structure dir to avoid overwriting
+        for filetype in ["init", "run"]:
+            filename = f"{filetype}_parameters.json"
+            os.rename(
+                self.step_local_staging_dir / filename,
+                struct_dir / filename
+            )
+
         # Perform some per-cell analysis
         run_shcoeffs_analysis(df=coeffs_df, savedir=avg_data_dir, struct=struct)
 
@@ -82,18 +93,42 @@ def run(self, sh_df=None, struct="Nuc", **kwargs):
         coeffs_avg = coeffs_avg.reshape(-2, lmax, lmax)
 
         # Here we use the new meshing implementation for a more evenly distributed mesh
+        '''
         mesh_avg, _ = shtools.get_even_reconstruction_from_coeffs(
             coeffs=coeffs_avg,
             npoints=1024
         )
+        '''
+
+        mesh_avg, grid_avg = shtools.get_reconstruction_from_coeffs(
+            coeffs=coeffs_avg,
+        )
 
         shtools.save_polydata(
             mesh=mesh_avg,
             filename=str(avg_data_dir / f"avgshape_{struct}.ply")
         )
 
+        # Save mesh as obj
+        save_mesh_as_obj(mesh_avg, avg_data_dir / f"avgshape_{struct}.obj")
+
+        # Save displacement map
+        save_displacement_map(grid_avg, avg_data_dir / f"avgshape_dmap_{struct}.tif")
+
+        # Save mesh as image
+        domain = save_voxelization(mesh_avg, avg_data_dir / f"avgshape_{struct}.tif")
+
         # Save mesh as stl file for blender import
-        save_mesh_as_stl(mesh_avg, str(avg_data_dir / f"avgshape_{struct}.stl"))
+        save_mesh_as_stl(mesh_avg, avg_data_dir / f"avgshape_{struct}.stl")
+
+        # Remesh voxelization
+        remesh_avg = get_smooth_and_coarse_mesh_from_voxelization(domain, sigma=3, npoints=2000)
+
+        # Save remesh as PLY
+        shtools.save_polydata(
+            mesh=remesh_avg,
+            filename=str(avg_data_dir / f"avgshape_remesh_{struct}.ply")
+        )
 
         # Save avg coeffs to csv file
         coeffs_df_avg.to_csv(
@@ -104,12 +139,16 @@ def run(self, sh_df=None, struct="Nuc", **kwargs):
         self.manifest = pd.DataFrame({
             "Label": "Average_mesh",
             "AvgShapeFilePath": avg_data_dir / f"avgshape_{struct}.ply",
+            "AvgShapeRemeshFilePath": avg_data_dir / f"avgshape_remesh_{struct}.ply",
             "AvgShapeFilePathStl": avg_data_dir / f"avgshape_{struct}.stl",
+            "AvgShapeFilePathObj": avg_data_dir / f"avgshape_{struct}.obj",
+            "AvgShapeFilePathTif": avg_data_dir / f"avgshape_{struct}.tif",
+            "AvgShapeDMapFilePathTif": avg_data_dir / f"avgshape_dmap_{struct}.tif",
             "Structure": struct,
         }, index=[0])
 
         # Save manifest as csv
         self.manifest.to_csv(
-            self.step_local_staging_dir / Path(f"manifest.csv"), index=False
+            struct_dir / Path(f"manifest.csv"), index=False
         )
         return self.manifest

mti_nma/steps/avgshape/avgshape.py

@@ -1,9 +1,16 @@
 import logging
 from pathlib import Path
 
+import vtk
 import matplotlib.pyplot as plt
 import numpy as np
 from stl import mesh
+from aicscytoparam import cytoparam
+from aicsimageio import writers
+from skimage import filters as skfilters
+from vtk.util import numpy_support as vtknp
+import pyvista as pv
+import pyacvd
 
 ###############################################################################
 
@@ -131,3 +138,89 @@ def save_mesh_as_stl(polydata, fname):
 
     # Write the mesh to file"
     nuc_mesh.save(fname)
+
+
+def save_mesh_as_obj(polydata, fname):
+
+    writer = vtk.vtkOBJWriter()
+    writer.SetInputData(polydata)
+    writer.SetFileName(str(fname))
+    writer.Write()
+
+
+def save_voxelization(polydata, fname):
+
+    domain, _ = cytoparam.voxelize_meshes([polydata])
+    with writers.ome_tiff_writer.OmeTiffWriter(fname, overwrite_file=True) as writer:
+        writer.save(
+            255 * domain,
+            dimension_order='ZYX',
+            image_name=fname.stem
+        )
+    return domain
+
+
+def save_displacement_map(grid, fname):
+
+    grid = grid.reshape(1, *grid.shape).astype(np.float32)
+
+    with writers.ome_tiff_writer.OmeTiffWriter(fname, overwrite_file=True) as writer:
+        writer.save(
+            grid,
+            dimension_order='ZYX',
+            image_name=fname.stem
+        )
+
+
+def get_smooth_and_coarse_mesh_from_voxelization(img, sigma, npoints):
+    """
+    Converts an image into a triangle mesh with even distributed points.
+    First we use a Gaussian kernel with size (sigma**3) to smooth the
+    input image. Next we apply marching cubes (vtkContourFilter) to obtain
+    a first mesh, which is used as input to a Voronoi-based clustering
+    that is responsible for remeshing. Details can be found here:
+    https://github.com/pyvista/pyacvd
+
+    Parameters
+    ----------
+    img: np.array
+        Input image corresponding to the voxelized version of the original
+        average mesh.
+    sigma: float
+        Gaussian kernel size.
+    npoints: int
+        Number of points used to create the Voronoi clustering. The larger
+        this value the more points the final mesh will have.
+    Returns
+    -------
+    remesh_vtk: vtkPolyData
+        Triangle with even distirbuted points.
+    """
+
+
+    rad = 5
+    img = np.pad(img, ((rad, rad), (rad, rad), (rad, rad)))
+    d, h, w = img.shape
+    img = skfilters.gaussian(img > 0, sigma=sigma, preserve_range=True)
+    imagedata = vtk.vtkImageData()
+    imagedata.SetDimensions([w, h, d])
+    imagedata.SetExtent(0, w - 1, 0, h - 1, 0, d - 1)
+    imagedata.AllocateScalars(vtk.VTK_UNSIGNED_CHAR, 1)
+    values = (255 * img).ravel().astype(np.uint8)
+    values = vtknp.numpy_to_vtk(values, deep=True, array_type=vtk.VTK_UNSIGNED_CHAR)
+    imagedata.GetPointData().SetScalars(values)
+    cf = vtk.vtkContourFilter()
+    cf.SetInputData(imagedata)
+    cf.SetValue(0, 255.0 / np.exp(1.0))
+    cf.Update()
+    mesh = cf.GetOutput()
+
+    pv_temp = pv.PolyData(mesh)
+    cluster = pyacvd.Clustering(pv_temp)
+    cluster.cluster(npoints)
+    remesh = cluster.create_mesh()
+    remesh_vtk = vtk.vtkPolyData()
+    remesh_vtk.SetPoints(remesh.GetPoints())
+    remesh_vtk.SetVerts(remesh.GetVerts())
+    remesh_vtk.SetPolys(remesh.GetPolys())
+    return remesh_vtk

mti_nma/steps/avgshape/avgshape_utils.py

@@ -0,0 +1,5 @@
+# -*- coding: utf-8 -*-
+
+from .load_data import LoadData  # noqa: F401
+
+__all__ = ["LoadData"]