Wrl processing

preprocessor/cellstar_preprocessor/flows/common.py

@@ -480,12 +480,9 @@ def prepare_ometiff_for_writing(
     missing_dims = []
 
     if len(img_array.shape) != 5:
-        local_d = {"T": 0, "Z": 1, "C": 2, "Y": 3, "X": 4}
         missing_dims = _get_missing_dims(metadata["Sizes BF"])
         for missing_dim in missing_dims:
-            img_array = da.expand_dims(img_array, axis=local_d[missing_dim])
-
-        d = local_d
+            img_array = da.expand_dims(img_array, axis=d[missing_dim])
 
     CORRECT_ORDER = "TCXYZ"
     reorder_tuple = _create_reorder_tuple(d, CORRECT_ORDER)

preprocessor/cellstar_preprocessor/flows/segmentation/helper_methods.py

@@ -3,7 +3,8 @@
 import math
 import zlib
 from pathlib import Path
-
+from vtk import vtkPolyData
+from cellstar_preprocessor.tools.open_wrl_file.mesh import create_mesh_from_vertices_and_triangles, decimate_mesh, get_vertices_and_triangles_from_vtk_mesh
 import h5py
 import numcodecs
 import numpy as np
@@ -251,6 +252,12 @@ def compute_vertex_density(mesh_list_group: zarr.hierarchy.group, mode="area"):
         return total_vertex_count / total_volume
 
 
+def _convert_mesh_to_vtk_mesh(mesh_from_zarr: zarr.Group):
+    vertices = mesh_from_zarr.vertices[...]
+    triangles = mesh_from_zarr.triangles[...]
+    vtk_mesh = create_mesh_from_vertices_and_triangles(triangles, vertices)
+    return vtk_mesh
+
 def _convert_mesh_to_vedo_obj(mesh_from_zarr):
     vertices = mesh_from_zarr.vertices[...]
     triangles = mesh_from_zarr.triangles[...]
@@ -259,8 +266,18 @@ def _convert_mesh_to_vedo_obj(mesh_from_zarr):
 
 
 def _decimate_vedo_obj(vedo_obj, ratio):
+    # TODO: breaks here while trying to decimate mesh
+    # would it break with empiar 10070?
+    # no
     return vedo_obj.decimate(fraction=ratio)
 
+def _decimate_vtk_mesh(vtk_mesh: vtkPolyData, fraction: float):
+    return decimate_mesh(vtk_mesh, fraction)
+
+def _get_mesh_data_from_vtk_mesh(vtk_mesh: vtkPolyData):
+   triangles, verts  = get_vertices_and_triangles_from_vtk_mesh(vtk_mesh)
+   return triangles, verts
+
 
 def _get_mesh_data_from_vedo_obj(vedo_obj):
     d = {"arrays": {}, "attrs": {}}
@@ -315,7 +332,7 @@ def store_mesh_data_in_zarr(
 
 
 def simplify_meshes(
-    mesh_list_group: zarr.hierarchy.Group, ratio: float, segment_id: int
+    mesh_list_group: zarr.Group, ratio: float
 ):
     """Returns dict with mesh data for each mesh in mesh list"""
     # for each mesh
@@ -324,8 +341,12 @@ def simplify_meshes(
     # get vertices and triangles back
     d = {}
     for mesh_id, mesh in mesh_list_group.groups():
-        vedo_obj = _convert_mesh_to_vedo_obj(mesh)
-        decimated_vedo_obj = _decimate_vedo_obj(vedo_obj, ratio)
-        mesh_data = _get_mesh_data_from_vedo_obj(decimated_vedo_obj)
+        vtk_mesh = _convert_mesh_to_vtk_mesh(mesh)
+        decimated_vtk_mesh = _decimate_vtk_mesh(vtk_mesh, ratio)
+        mesh_data = _get_mesh_data_from_vtk_mesh(decimated_vtk_mesh)
+        print(mesh_data)
+        # vedo_obj = _convert_mesh_to_vedo_obj(mesh)
+        # decimated_vedo_obj = _decimate_vedo_obj(vedo_obj, ratio)
+        # mesh_data = _get_mesh_data_from_vedo_obj(decimated_vedo_obj)
         d[mesh_id] = mesh_data
     return d

preprocessor/cellstar_preprocessor/flows/segmentation/segmentation_downsampling.py

@@ -110,12 +110,12 @@ def sff_segmentation_downsampling(internal_segmentation: InternalSegmentation):
                             <= density_threshold
                         ):
                             break
-                        if fraction == 1:
+                        if fraction == 1.0:
                             continue  # original data, don't need to compute anything
+                        # breaks here
                         mesh_data_dict = simplify_meshes(
                             original_detail_lvl_mesh_list_group,
-                            ratio=fraction,
-                            segment_id=segment_name_id,
+                            ratio=fraction
                         )
                         # TODO: potentially simplify meshes may output mesh with 0 vertices, normals, triangles
                         # it should not be stored?
@@ -147,7 +147,9 @@ def sff_segmentation_downsampling(internal_segmentation: InternalSegmentation):
         if internal_segmentation.downsampling_parameters.remove_original_resolution:
             # del internal_segmentation.simplification_curve[1]
             internal_segmentation.simplification_curve.pop(1, None)
-
+    else:
+        # Breaks
+        raise Exception(f"Primary descriptor {internal_segmentation.primary_descriptor} is not supported")
     print("Segmentation downsampled")
 
 

preprocessor/cellstar_preprocessor/flows/segmentation/sff_preprocessing.py

@@ -34,7 +34,9 @@ def sff_preprocessing(internal_segmentation: InternalSegmentation):
     # PLAN:
     # 1. Convert hff to intermediate zarr structure
     # 2. Process it with one of 2 methods (3d volume segmentation, mesh segmentation)
-    if zarr_structure.primary_descriptor[0] == b"three_d_volume":
+    # TODO: check primary d
+    primary_descriptor = zarr_structure.primary_descriptor[0]
+    if primary_descriptor == b"three_d_volume":
         segm_data_gr: zarr.Group = zarr_structure.create_group(
             LATTICE_SEGMENTATION_DATA_GROUPNAME
         )
@@ -47,7 +49,7 @@ def sff_preprocessing(internal_segmentation: InternalSegmentation):
         _process_three_d_volume_segmentation_data(
             segm_data_gr, zarr_structure, internal_segmentation=internal_segmentation
         )
-    elif zarr_structure.primary_descriptor[0] == b"mesh_list":
+    elif primary_descriptor == b"mesh_list":
         segm_data_gr: zarr.Group = zarr_structure.create_group(
             MESH_SEGMENTATION_DATA_GROUPNAME
         )
@@ -65,7 +67,8 @@ def sff_preprocessing(internal_segmentation: InternalSegmentation):
         _process_mesh_segmentation_data(
             timeframe_gr, zarr_structure, internal_segmentation=internal_segmentation
         )
-
+    else:
+        raise Exception(f"Primary descriptor {primary_descriptor} is not supported")
     print("Segmentation processed")
 
 
@@ -105,6 +108,7 @@ def _process_mesh_segmentation_data(
 ):
     params_for_storing = internal_segmentation.params_for_storing
 
+    # single segment for wrl
     for segment_name, segment in zarr_structure.segment_list.groups():
         segment_id = int(segment.id[...])
         single_segment_group = timeframe_gr.create_group(segment_id)
@@ -131,3 +135,5 @@ def _process_mesh_segmentation_data(
                 )
                 single_mesh_group.attrs["area"] = vedo_mesh_obj.area()
                 # single_mesh_group.attrs['volume'] = vedo_mesh_obj.volume()
+        else:
+            raise Exception(f"No mesh list is provided for segment with name: {segment_name}")

preprocessor/cellstar_preprocessor/tools/open_wrl_file/mesh.py

@@ -0,0 +1,111 @@
+from pathlib import Path
+import numpy as np
+from vtk import vtkSphereSource, vtkPolyData, vtkDecimatePro, vtkIdList, vtkPoints, vtkFloatArray, vtkCellArray
+from vtkmodules.util.numpy_support import vtk_to_numpy
+
+# create mesh
+def mkVtkIdList(it):
+  vil = vtkIdList()
+  for i in it:
+    vil.InsertNextId(int(i))
+  return vil
+
+def CreateMesh(modelNode, arrayVertices, arrayTriangles):
+  # modelNode : a vtkMRMLModelNode in the Slicer scene which will hold the mesh
+  # arrayVertices : list of triples [[x1,y1,z2], [x2,y2,z2], ... ,[xn,yn,zn]] of vertex coordinates
+  # arrayVertexNormals : list of triples [[nx1,ny1,nz2], [nx2,ny2,nz2], ... ] of vertex normals
+  # arrayTriangles : list of triples of 0-based indices defining triangles
+  # labelsScalars : list of strings such as ["bipolar", "unipolar"] to label the individual scalars data sets
+  # arrayScalars : an array of n rows for n vertices and m colums for m inidividual scalar sets
+
+  # create the building blocks of polydata including data attributes.
+  mesh    = vtkPolyData()
+  points  = vtkPoints()
+#   normals = vtkFloatArray()
+  polys   = vtkCellArray()
+
+  # load the array data into the respective VTK data structures
+  for i in range(len(arrayVertices)):
+    points.InsertPoint(i, arrayVertices[i])
+
+  for i in range(len(arrayTriangles)):
+    polys.InsertNextCell(mkVtkIdList(arrayTriangles[i]))
+
+#   for i in range(len(arrayVertexNormals)):
+#     normals.InsertTuple3(i, arrayVertexNormals[i][0], arrayVertexNormals[i][1], arrayVertexNormals[i][2])
+
+  # put together the mesh object
+  mesh.SetPoints(points)
+  mesh.SetPolys(polys)
+#   if(len(arrayVertexNormals) == len(arrayVertices)):
+#     mesh.GetPointData().SetNormals(normals)
+
+  # Add scalars
+
+  modelNode.SetAndObservePolyData(mesh)
+
+
+def decimation(inputPoly: vtkPolyData):
+    sphereS = vtkSphereSource()
+    sphereS.Update()
+
+    # inputPoly = vtkPolyData()
+
+    # inputPoly.ShallowCopy(sphereS.GetOutput())
+
+    print("Before decimation\n"
+          "-----------------\n"
+          "There are " + str(inputPoly.GetNumberOfPoints()) + "points.\n"
+          "There are " + str(inputPoly.GetNumberOfPolys()) + "polygons.\n")
+
+    decimate = vtkDecimatePro()
+    decimate.SetInputData(inputPoly)
+    decimate.SetTargetReduction(.10)
+    decimate.Update()
+
+    decimatedPoly = vtkPolyData()
+    decimatedPoly.ShallowCopy(decimate.GetOutput())
+
+    print("After decimation \n"
+          "-----------------\n"
+          "There are " + str(decimatedPoly.GetNumberOfPoints()) + "points.\n"
+          "There are " + str(decimatedPoly.GetNumberOfPolys()) + "polygons.\n")
+
+
+def create_mesh_from_vertices_and_triangles(triangles: np.ndarray, vertices: np.ndarray):
+    mesh    = vtkPolyData()
+    points  = vtkPoints()
+#   normals = vtkFloatArray()
+    polys   = vtkCellArray()
+
+    for i in range(len(vertices)):
+        points.InsertPoint(i, vertices[i])
+
+    for i in range(len(triangles)):
+        polys.InsertNextCell(mkVtkIdList(triangles[i]))
+
+
+    mesh.SetPoints(points)
+    mesh.SetPolys(polys)
+    return mesh
+
+def get_vertices_and_triangles_from_vtk_mesh(mesh: vtkPolyData):
+    triangles: vtkCellArray = mesh.GetPolys()
+    vertices: vtkPoints = mesh.GetPoints()
+    # TODO: convert to np arrs
+    # https://discourse.vtk.org/t/convert-vtk-array-to-numpy-array/3152
+    # numpy_interface
+    # print(len(triangles), len(vertices))
+    return (vtk_to_numpy(triangles), vtk_to_numpy(vertices))
+
+def decimate_mesh(mesh: vtkPolyData, fraction: float):
+    assert fraction <= 1.0
+  # TODO: quadratic?
+    decimate = vtkDecimatePro()
+    decimate.SetInputData(mesh)
+    decimate.SetTargetReduction(fraction)
+    decimate.Update()
+
+    decimatedPoly = vtkPolyData()
+    decimatedPoly.ShallowCopy(decimate.GetOutput())
+    return decimatedPoly

preprocessor/cellstar_preprocessor/tools/open_wrl_file/open_wrl_file.py

@@ -0,0 +1,20 @@
+
+from pathlib import Path
+from cellstar_preprocessor.flows.segmentation.helper_methods import open_hdf5_as_segmentation_object
+from cellstar_preprocessor.tools.convert_app_specific_segm_to_sff.convert_app_specific_segm_to_sff import convert_app_specific_segm_to_sff
+import pymeshlab
+
+
+def convert_wrl_to_sff(input_file: Path, working_folder: Path):
+    ms = pymeshlab.MeshSet()
+# You can load, save meshes and apply MeshLab filters:
+    stl_path = working_folder / 'temp_stl_dir' / 'temp_stl.stl'
+    ms.load_new_mesh(str(input_file.resolve()))
+    ms.save_current_mesh(str(stl_path.resolve()))
+    sff_path = convert_app_specific_segm_to_sff(stl_path)
+    # sff_obj = open_hdf5_as_segmentation_object(sff)
+    return sff_path
+
+
+if __name__ == "__main__":
+    convert_wrl_to_sff(Path('preprocessor/cellstar_preprocessor/tools/open_wrl_file/actin.wrl'), Path('preprocessor/cellstar_preprocessor/tools/open_wrl_file/'))