Added method geometry_utils.create_icosphere

Author: adamhusar

CMakeLists.txt

@@ -165,6 +165,8 @@ include_directories(
   "${VTK_DIR}/Filters/Core/"
   "${CMAKE_SOURCE_DIR}/../VTK/Filters/Modeling/"
   "${VTK_DIR}/Filters/Modeling/"
+  "${CMAKE_SOURCE_DIR}/../VTK/Filters/Sources/"
+  "${VTK_DIR}/Filters/Sources/"
   "${CMAKE_SOURCE_DIR}/../VTK/Rendering/Core/"
   "${VTK_DIR}/Rendering/Core/"
   "${CMAKE_SOURCE_DIR}/../VTK/Rendering/Context2D/"

libmcell/api/geometry_utils.cpp

@@ -23,6 +23,9 @@
 #include "generated/gen_geometry_utils.h"
 
 #include "api/geometry_object.h"
+#include <vtkPlatonicSolidSource.h>
+#include <vtkPolyData.h>
+#include <vtkSmartPointer.h>
 
 using namespace std;
 
@@ -72,6 +75,224 @@ std::shared_ptr<GeometryObject> create_box(const std::string& name, const double
   return res;
 }
 
+
+
+// used code from the link below due to absence of any reasonable library
+// https://github.com/caosdoar/spheres/blob/master/src/spheres.cpp
+// published under MIT license
+
+struct Edge
+{
+  uint32_t v0;
+  uint32_t v1;
+
+  Edge(uint32_t v0, uint32_t v1)
+    : v0(v0 < v1 ? v0 : v1)
+    , v1(v0 < v1 ? v1 : v0)
+  {
+  }
+
+  bool operator <(const Edge &rhs) const
+  {
+    return v0 < rhs.v0 || (v0 == rhs.v0 && v1 < rhs.v1);
+  }
+};
+
+
+static Vec3 normalize(const Vec3 &a)
+{
+  const double lrcp = 1.0 / std::sqrt(dot(a, a));
+  return Vec3(a.x * lrcp, a.y * lrcp, a.z * lrcp);
+}
+
+
+static void init_icosphere(
+    vector<Vec3>& vertex_list,
+    vector<vector<int>>& wall_list) {
+
+  vertex_list.clear();
+  wall_list.clear();
+
+  const double t = (1.0 + std::sqrt(5.0)) / 2.0;
+
+  // Vertices
+  vertex_list.push_back(normalize(Vec3(-1.0,  t, 0.0)));
+  vertex_list.push_back(normalize(Vec3( 1.0,  t, 0.0)));
+  vertex_list.push_back(normalize(Vec3(-1.0, -t, 0.0)));
+  vertex_list.push_back(normalize(Vec3( 1.0, -t, 0.0)));
+  vertex_list.push_back(normalize(Vec3(0.0, -1.0,  t)));
+  vertex_list.push_back(normalize(Vec3(0.0,  1.0,  t)));
+  vertex_list.push_back(normalize(Vec3(0.0, -1.0, -t)));
+  vertex_list.push_back(normalize(Vec3(0.0,  1.0, -t)));
+  vertex_list.push_back(normalize(Vec3( t, 0.0, -1.0)));
+  vertex_list.push_back(normalize(Vec3( t, 0.0,  1.0)));
+  vertex_list.push_back(normalize(Vec3(-t, 0.0, -1.0)));
+  vertex_list.push_back(normalize(Vec3(-t, 0.0,  1.0)));
+
+  // Faces
+  wall_list.push_back(vector<int>{0, 11, 5});
+  wall_list.push_back(vector<int>{0, 5, 1});
+  wall_list.push_back(vector<int>{0, 1, 7});
+  wall_list.push_back(vector<int>{0, 7, 10});
+  wall_list.push_back(vector<int>{0, 10, 11});
+  wall_list.push_back(vector<int>{1, 5, 9});
+  wall_list.push_back(vector<int>{5, 11, 4});
+  wall_list.push_back(vector<int>{11, 10, 2});
+  wall_list.push_back(vector<int>{10, 7, 6});
+  wall_list.push_back(vector<int>{7, 1, 8});
+  wall_list.push_back(vector<int>{3, 9, 4});
+  wall_list.push_back(vector<int>{3, 4, 2});
+  wall_list.push_back(vector<int>{3, 2, 6});
+  wall_list.push_back(vector<int>{3, 6, 8});
+  wall_list.push_back(vector<int>{3, 8, 9});
+  wall_list.push_back(vector<int>{4, 9, 5});
+  wall_list.push_back(vector<int>{2, 4, 11});
+  wall_list.push_back(vector<int>{6, 2, 10});
+  wall_list.push_back(vector<int>{8, 6, 7});
+  wall_list.push_back(vector<int>{9, 8, 1});
+}
+
+
+int subdivide_icosphere_edge(
+    int f0, int f1,
+    const Vec3& v0, const Vec3& v1,
+    vector<Vec3>& vertex_list_out,
+    std::map<Edge, uint32_t> &io_divisions
+    )
+{
+  // required to make the mesh watertight
+  const Edge edge(f0, f1);
+  auto it = io_divisions.find(edge);
+  if (it != io_divisions.end()) {
+    return it->second;
+  }
+
+  const Vec3 v = normalize(Vec3(0.5) * (v0 + v1));
+  const uint32_t f = vertex_list_out.size();
+  vertex_list_out.push_back(v);
+  io_divisions[edge] = f;
+  return f;
+}
+
+
+static void subdivide_icosphere_mesh(
+    const vector<Vec3>& vertex_list_in,
+    const vector<vector<int>>& wall_list_in,
+    vector<Vec3>& vertex_list_out,
+    vector<vector<int>>& wall_list_out) {
+
+  vertex_list_out = vertex_list_in;
+
+  std::map<Edge, uint32_t> divisions; // Edge -> new vertex
+
+  for (size_t i = 0; i < wall_list_in.size(); ++i)
+  {
+    const int f0 = wall_list_in[i][0];
+    const int f1 = wall_list_in[i][1];
+    const int f2 = wall_list_in[i][2];
+
+    const Vec3& v0 = vertex_list_in[f0];
+    const Vec3& v1 = vertex_list_in[f1];
+    const Vec3& v2 = vertex_list_in[f2];
+
+    const int f3 = subdivide_icosphere_edge(f0, f1, v0, v1, vertex_list_out, divisions);
+    const int f4 = subdivide_icosphere_edge(f1, f2, v1, v2, vertex_list_out, divisions);
+    const int f5 = subdivide_icosphere_edge(f2, f0, v2, v0, vertex_list_out, divisions);
+
+    wall_list_out.push_back(vector<int>{f0, f3, f5});
+    wall_list_out.push_back(vector<int>{f3, f1, f4});
+    wall_list_out.push_back(vector<int>{f4, f2, f5});
+    wall_list_out.push_back(vector<int>{f3, f4, f5});
+  }
+}
+
+
+std::shared_ptr<GeometryObject> create_icosphere(
+    const std::string& name, const double radius, const int subdivisions) {
+
+  if (subdivisions < 1) {
+    throw ValueError(S("Value of parameter ") + NAME_SUBDIVISIONS + " must be greater or equal to 1.");
+  }
+
+  if (subdivisions > 8) {
+    throw ValueError(S("Value of parameter ") + NAME_SUBDIVISIONS + " must be less or equal to 8.");
+  }
+
+  vector<Vec3> vertex_list;
+  vector<vector<int>> wall_list;
+
+  init_icosphere(vertex_list, wall_list);
+
+  for (int i = 0; i < subdivisions - 1; i++) {
+    vector<Vec3> vertex_list_new;
+    vector<vector<int>> wall_list_new;
+    subdivide_icosphere_mesh(vertex_list, wall_list, vertex_list_new, wall_list_new);
+
+    vertex_list.swap(vertex_list_new);
+    wall_list.swap(wall_list_new);
+  }
+
+  // convert vertices to the required format and scale
+  vector<vector<double>> double_vertices;
+  for (const Vec3& v: vertex_list) {
+    double_vertices.push_back(vector<double>{v.x * radius, v.y * radius, v.z * radius});
+  }
+
+  auto res = make_shared<GeometryObject>(
+      name,
+      double_vertices,
+      wall_list
+  );
+
+  return res;
+}
+
+#if 0
+// keeping this code as an example on how to transform vtkPolyData
+std::shared_ptr<GeometryObject> create_sphere(
+    const std::string& name, const double radius, const int resolution) {
+
+  vtkNew<vtkSphereSource> sphere;
+  sphere->SetCenter(0.0, 0.0, 0.0);
+  sphere->SetRadius(radius);
+  sphere->SetPhiResolution(resolution);
+  sphere->SetThetaResolution(resolution);
+  sphere->Update();
+
+  vtkSmartPointer<vtkPolyData> polydata = sphere->GetOutput();
+
+  vtkPoints* points = polydata->GetPoints();
+
+  vector<vector<double>> vertex_list;
+  int num_vertices = points->GetNumberOfPoints();
+  for (int i = 0; i < num_vertices; i++) {
+    double pt[3];
+    points->GetPoint(i, pt);
+    vertex_list.push_back(vector<double>{pt[0], pt[1], pt[2]});
+  }
+
+  vector<vector<int>> wall_list;
+  int num_walls = polydata->GetNumberOfCells();
+  for (int i = 0; i < num_walls; i++) {
+    vtkCell* cell = polydata->GetCell(i);
+
+    vector<int> wall;
+    for (int k = 0; k < 3; k++) {
+      wall.push_back(cell->GetPointId(k));
+    }
+    wall_list.push_back(wall);
+  }
+
+  auto res = make_shared<GeometryObject>(
+      name,
+      vertex_list,
+      wall_list
+  );
+
+  return res;
+}
+#endif
+
 } // namespace geometry_utils
 
 } // namespace API

libmcell/definition/submodules.yaml

@@ -6,7 +6,7 @@ geometry_utils:
   type: submodule  # default type is class
   methods:
   - name: create_box
-    doc: Creates a GeometryObject whose center is at (0, 0, 0).
+    doc: Creates a GeometryObject in the shape of a cube whose center is at (0, 0, 0).
     return_type: GeometryObject*
     examples: tests/pymcell4/1400_rel_site_for_each_it/model.py
     params:
@@ -18,6 +18,27 @@ geometry_utils:
       type: float
       doc: Specifies length of each edge of the box in um.
 
+  - name: create_icosphere
+    doc: Creates a GeometryObject in the shape of an icosphere whose center is at (0, 0, 0).
+    return_type: GeometryObject*
+    examples: tests/pymcell4/1110_point_release_w_create_icosphere/model.py
+    params:
+    - name: name
+      type: str
+      doc: Name of the created geometry object.
+
+    - name: radius
+      type: float
+      doc: Specifies radius of the sphere.
+
+    - name: subdivisions
+      type: int
+      min: 1
+      max: 8
+      doc: | 
+         Number of subdivisions from the initial icosphere. 
+         The higher this value will be the smoother the icosphere will be.
+         Allowed range is between 1 and 8.
       
 bngl_utils:
   doc: | 

libmcell/documentation/generated/submodules.rst

@@ -49,11 +49,32 @@ Methods:
    * | return type: GeometryObject
 
 
-  | Creates a GeometryObject whose center is at (0, 0, 0).
+  | Creates a GeometryObject in the shape of a cube whose center is at (0, 0, 0).
 
   | Example: `1400_rel_site_for_each_it/model.py <https://github.com/mcellteam/mcell_tests/tree/mcell4_dev/tests/pymcell4/1400_rel_site_for_each_it/model.py>`_ 
 
 
+* | **create_icosphere**
+
+   * | name: str
+     | Name of the created geometry object.
+
+   * | radius: float
+     | Specifies radius of the sphere.
+
+   * | subdivisions: int
+     | Number of subdivisions from the initial icosphere. 
+     | The higher this value will be the smoother the icosphere will be.
+     | Allowed range is between 1 and 8.
+
+   * | return type: GeometryObject
+
+
+  | Creates a GeometryObject in the shape of an icosphere whose center is at (0, 0, 0).
+
+  | Example: `1110_point_release_w_create_icosphere/model.py <https://github.com/mcellteam/mcell_tests/tree/mcell4_dev/tests/pymcell4/1110_point_release_w_create_icosphere/model.py>`_ 
+
+
 
 run_utils
 =========

libmcell/generated/gen_geometry_utils.cpp

@@ -31,7 +31,8 @@ namespace API {
 
 void define_pybinding_geometry_utils(py::module& m) {
   m.def_submodule("geometry_utils")
-      .def("create_box", &geometry_utils::create_box, py::arg("name"), py::arg("edge_length"), "Creates a GeometryObject whose center is at (0, 0, 0).\n- name: Name of the created geometry object.\n\n- edge_length: Specifies length of each edge of the box in um.\n\n")
+      .def("create_box", &geometry_utils::create_box, py::arg("name"), py::arg("edge_length"), "Creates a GeometryObject in the shape of a cube whose center is at (0, 0, 0).\n- name: Name of the created geometry object.\n\n- edge_length: Specifies length of each edge of the box in um.\n\n")
+      .def("create_icosphere", &geometry_utils::create_icosphere, py::arg("name"), py::arg("radius"), py::arg("subdivisions"), "Creates a GeometryObject in the shape of an icosphere whose center is at (0, 0, 0).\n- name: Name of the created geometry object.\n\n- radius: Specifies radius of the sphere.\n\n- subdivisions: Number of subdivisions from the initial icosphere. \nThe higher this value will be the smoother the icosphere will be.\nAllowed range is between 1 and 8.\n\n\n")
     ;
 }
 

libmcell/generated/gen_geometry_utils.h

@@ -34,6 +34,7 @@ class PythonExportContext;
 namespace geometry_utils {
 
 std::shared_ptr<GeometryObject> create_box(const std::string& name, const double edge_length);
+std::shared_ptr<GeometryObject> create_icosphere(const std::string& name, const double radius, const int subdivisions);
 
 } // namespace geometry_utils
 

libmcell/generated/gen_names.h

@@ -114,6 +114,7 @@ const char* const NAME_COUNT = "count";
 const char* const NAME_COUNT_EXPRESSION = "count_expression";
 const char* const NAME_COUNTS = "counts";
 const char* const NAME_CREATE_BOX = "create_box";
+const char* const NAME_CREATE_ICOSPHERE = "create_icosphere";
 const char* const NAME_CUSTOM_DIR = "custom_dir";
 const char* const NAME_CUSTOM_SPACE_STEP = "custom_space_step";
 const char* const NAME_CUSTOM_TIME_STEP = "custom_time_step";
@@ -224,6 +225,7 @@ const char* const NAME_PRODUCT_IDS = "product_ids";
 const char* const NAME_PRODUCTS = "products";
 const char* const NAME_PROPERTIES = "properties";
 const char* const NAME_R = "r";
+const char* const NAME_RADIUS = "radius";
 const char* const NAME_RANDCNT = "randcnt";
 const char* const NAME_RANDSLR = "randslr";
 const char* const NAME_REACTANT_IDS = "reactant_ids";
@@ -270,6 +272,7 @@ const char* const NAME_SPECIES_LIST = "species_list";
 const char* const NAME_SPECIES_PATTERN = "species_pattern";
 const char* const NAME_STATE = "state";
 const char* const NAME_STATES = "states";
+const char* const NAME_SUBDIVISIONS = "subdivisions";
 const char* const NAME_SUBPARTITION_DIMENSION = "subpartition_dimension";
 const char* const NAME_SUBSYSTEM = "subsystem";
 const char* const NAME_SURFACE_CLASS = "surface_class";

libmcell/generated/mcell.pyi

@@ -1563,6 +1563,14 @@ class geometry_utils():
         ) -> 'GeometryObject':
         pass
 
+    def create_icosphere(
+            self,
+            name : str,
+            radius : float,
+            subdivisions : int
+        ) -> 'GeometryObject':
+        pass
+
 class run_utils():
     def __init__(
             self,