Cgal (#162)

* updated pybind11

* small setup cleanup

* adds spherical sheet and monolayer generators from CGAL

* bump version number

Author: glyg

.travis.yml

@@ -10,7 +10,6 @@ before_install:
   - export PATH="$HOME/miniconda/bin:$PATH"
   - hash -r
   - conda config --set always_yes yes --set changeps1 no
-  - conda update -q conda
   - conda info -a
   - conda env create -q -f environment.yml
   - source activate tyssue

CMakeLists.txt

@@ -28,10 +28,13 @@ endif()
 include( ${CGAL_USE_FILE})
 
 
-set (SOURCE_DIR "tyssue/collisions/cpp")
-include_directories (${SOURCE_DIR})
+set (COLLISION_SRCE "tyssue/collisions/cpp")
+set (GENERATION_SRCE "tyssue/generation/cpp")
+
+include_directories (${COLLISION_SRCE} ${GENERATION_SRCE})
 
 # Assumes pybind11 is a subdirectory at the project root
 # Make sure you git cloned tyssue recursively
 add_subdirectory(pybind11)
-pybind11_add_module(c_collisions ${SOURCES} "${SOURCE_DIR}/c_collisions.cpp" )
+pybind11_add_module(c_collisions ${SOURCES} "${COLLISION_SRCE}/c_collisions.cpp" )
+pybind11_add_module(mesh_generation ${SOURCES} "${GENERATION_SRCE}/mesh_generation.cpp" )

pybind11

@@ -29,8 +29,8 @@
 ## Thanks to them!
 MAJOR = 0
 MINOR = 6
-MICRO = 0
-ISRELEASED = True
+MICRO = 1
+ISRELEASED = False
 VERSION = "%d.%d.%s" % (MAJOR, MINOR, MICRO)
 
 
@@ -166,8 +166,7 @@ def build_extension(self, ext):
         subprocess.check_call(
             ["cmake", "--build", "."] + build_args, cwd=self.build_temp
         )
-        print(env["CXXFLAGS"])
-        print()  # Add an empty line for cleaner output
+        print(env["CXXFLAGS"], "\n")
 
 
 if __name__ == "__main__":
@@ -199,7 +198,10 @@ def build_extension(self, ext):
         packages=find_packages(),
         package_data={"tyssue": files},
         include_package_data=True,
-        ext_modules=[CMakeExtension("tyssue/collisions/cpp/c_collisions")],
+        ext_modules=[
+            CMakeExtension("tyssue/collisions/cpp/c_collisions"),
+            CMakeExtension("tyssue/generation/cpp/mesh_generation"),
+        ],
         cmdclass=dict(build_ext=CMakeBuild),
         zip_safe=False,
     )

setup.py

@@ -5,7 +5,12 @@
 
 import numpy as np
 from numpy.testing import assert_almost_equal
-from tyssue.generation.shapes import generate_ring, ellipsoid_sheet, spherical_sheet
+from tyssue.generation.shapes import (
+    generate_ring,
+    ellipsoid_sheet,
+    spherical_sheet,
+    spherical_monolayer,
+)
 from tyssue import PlanarGeometry
 
 
@@ -62,4 +67,21 @@ def test_ellipsoid():
 def test_spherical():
     sph = spherical_sheet(1.0, 40)
     rho = np.linalg.norm(sph.vert_df[sph.coords], axis=1)
-    np.testing.assert_array_almost_equal(rho, np.ones(sph.Nv), decimal=1)
+    np.testing.assert_almost_equal(rho.mean(), 1.0, decimal=1)
+
+
+def test_spherical_mono():
+    radius = 1  # µm
+    height = 3  # the height of the cells
+    Nc = 400  # the (approximate) number of cells
+    R_in = radius
+    R_out = radius + height
+
+    organo = spherical_monolayer(R_in, R_out, Nc, apical="in")
+
+    rho = np.linalg.norm(organo.vert_df[organo.coords], axis=1)
+
+    mean_apical = rho[organo.vert_df["segment"] == "apical"].mean()
+    mean_basal = rho[organo.vert_df["segment"] == "basal"].mean()
+    np.testing.assert_almost_equal(mean_apical, radius, decimal=1)
+    np.testing.assert_almost_equal(mean_basal, radius + height, decimal=1)

tests/generation/test_shapes.py

@@ -6,3 +6,4 @@
 from .from_voronoi import *
 from .modifiers import *
 from .shapes import *
+from .cpp import mesh_generation

tyssue/generation/__init__.py

@@ -0,0 +1,68 @@
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/complex.h>
+#include <pybind11/numpy.h>
+#include <boost/lexical_cast.hpp>
+
+#include <fstream>
+#include <vector>
+#include <string>
+#include <sstream>
+#include <algorithm>
+#include <iostream>
+
+
+#include <CGAL/Surface_mesh_default_triangulation_3.h>
+#include <CGAL/Complex_2_in_triangulation_3.h>
+#include <CGAL/make_surface_mesh.h>
+#include <CGAL/Implicit_surface_3.h>
+
+// default triangulation for Surface_mesher
+typedef CGAL::Surface_mesh_default_triangulation_3 Tr;
+// c2t3
+typedef CGAL::Complex_2_in_triangulation_3<Tr> C2t3;
+typedef Tr::Geom_traits GT;
+typedef GT::Sphere_3 Sphere_3;
+typedef GT::Point_3 Point_3;
+typedef GT::FT FT;
+typedef FT (*Function)(Point_3);
+typedef CGAL::Implicit_surface_3<GT, Function> Surface_3;
+//typedef C2t3::Vertex_handle vertex_descriptor;
+//using vertex_descriptor       = boost::graph_traits<C2t3>::vertex_descriptor;
+
+
+FT sphere_function (Point_3 p) {
+  const FT x2=p.x()*p.x(), y2=p.y()*p.y(), z2=p.z()*p.z();
+  return x2+y2+z2-1;
+}
+
+namespace py = pybind11;
+
+std::vector<std::tuple<double, double, double>> make_spherical(double num_points) {
+  Tr tr;            // 3D-Delaunay triangulation
+  C2t3 c2t3 (tr);   // 2D-complex in 3D-Delaunay triangulation
+  // defining the surface
+  Surface_3 surface(sphere_function,             // pointer to function
+                    Sphere_3(CGAL::ORIGIN, 2.)); // bounding sphere
+  // Note that "2." above is the *squared* radius of the bounding sphere!
+  // defining meshing criteria
+  double r_max = std::sqrt(1 / (num_points * 0.12));
+  CGAL::Surface_mesh_default_criteria_3<Tr> criteria(30, r_max, r_max);
+  // meshing surface
+  CGAL::make_surface_mesh(c2t3, surface, criteria, CGAL::Non_manifold_tag());
+  std::vector<std::tuple<double, double, double>> points;
+  Tr::Finite_vertices_iterator it;
+  for (it = tr.finite_vertices_begin(); it != tr.finite_vertices_end(); it++)
+    {
+      Point_3 p = tr.point(it);
+      std::tuple<double, double, double> tmp = std::make_tuple(p.x(), p.y(), p.z());
+      points.push_back(tmp);
+    }
+  return points;
+}
+
+
+PYBIND11_MODULE(mesh_generation, m)
+{
+    m.def("make_spherical", &make_spherical);
+}

tyssue/generation/cpp/__init__.py

@@ -9,12 +9,19 @@
 from ..core.sheet import Sheet, get_outer_sheet
 from ..core.objects import get_prev_edges
 from ..core.objects import Epithelium
+from ..core.monolayer import Monolayer
+
 from ..topology import type1_transition
 from .from_voronoi import from_3d_voronoi
-from ..geometry.bulk_geometry import BulkGeometry
-from ..geometry.sheet_geometry import EllipsoidGeometry, SheetGeometry
-
-from ..utils import single_cell
+from ..geometry.bulk_geometry import BulkGeometry, ClosedMonolayerGeometry
+from ..geometry.sheet_geometry import (
+    EllipsoidGeometry,
+    SheetGeometry,
+    ClosedSheetGeometry,
+)
+from .cpp import mesh_generation
+from .modifiers import extrude
+from ..utils import single_cell, swap_apico_basal
 
 
 class AnnularSheet(Sheet):
@@ -272,13 +279,32 @@ def spherical_sheet(radius, Nf, **kwargs):
     the given number of cells
     """
 
-    n_zs = int(np.ceil(np.roots([2, 1.0, -Nf])[-1]))  # determined experimentaly ;p
-    eptm = ellipsoid_sheet(radius, radius, radius, n_zs, **kwargs)
-    eptm.settings.pop("abc")
-    eptm.settings["radius"] = radius
+    centers = np.array(mesh_generation.make_spherical(Nf))
+    eptm = sheet_from_cell_centers(centers, **kwargs)
+
+    rhos = (eptm.vert_df[eptm.coords] ** 2).sum(axis=1).mean()
+    ClosedSheetGeometry.scale(eptm, radius / rhos, eptm.coords)
+
+    ClosedSheetGeometry.update_all(eptm)
     return eptm
 
 
+def spherical_monolayer(R_in, R_out, Nc, apical="out"):
+    """Returns a spherical monolayer with the given inner and
+    outer radii, and approximately the gieven number of cells.
+
+    The `apical` argument can be 'in' out 'out' to specify wether
+    the apical face of the cells faces inward or outward, reespectively.
+    """
+    sheet = spherical_sheet(R_in, Nc)
+    delta_R = R_out - R_in
+    mono = Monolayer("mono", extrude(sheet.datasets, method="normals", scale=-delta_R))
+    if apical == "out":
+        swap_apico_basal(mono)
+    ClosedMonolayerGeometry.update_all(mono)
+    return mono
+
+
 def sheet_from_cell_centers(points, noise=0):
     """Returns a Sheet object from the Voronoï tessalation
     of the cell centers.

tyssue/generation/cpp/mesh_generation.cpp

@@ -314,3 +314,14 @@ def get_next(eptm):
     )
     next_ = fs_indexed.loc[ft_index, "edge"].values
     return next_
+
+
+## small utlity to swap apical and basal segments
+def swap_apico_basal(organo):
+    """Swap apical and basal segments of an organoid
+    """
+    for elem in ["vert", "face", "edge"]:
+        swaped = organo.datasets[elem]["segment"].copy()
+        swaped.loc[organo.segment_index("apical", elem)] = "basal"
+        swaped.loc[organo.segment_index("basal", elem)] = "apical"
+        organo.datasets[elem]["segment"] = swaped