.. toctree::
   :maxdepth: 2

.. _examples-analysis:

################################################################################
Geometry Analysis
################################################################################


**********************
Geometrical properties
**********************

The :mod:`CellComplex <cellcomplex>` library comes with tools to automatically compute a whole range of geometrical properties on the tissue mesh (:func:`compute_topomesh_property <cellcomplex.property_topomesh.analysis.compute_topomesh_property>`) and to propagate properties between elements of different dimensions (:func:`compute_topomesh_vertex_property_from_faces <cellcomplex.property_topomesh.analysis.compute_topomesh_vertex_property_from_faces>` for instance). First, let's create a multicellular complex with random cell shapes using :func:`circle_voronoi_topomesh <cellcomplex.property_topomesh.example_topomesh.circle_voronoi_topomesh>`.

.. code-block:: python
    :linenos:

    from cellcomplex.property_topomesh.example_topomesh import circle_voronoi_topomesh

    topomesh = circle_voronoi_topomesh(3)

.. plot::

    import numpy as np
    np.random.seed(583)

    import matplotlib.pyplot as plt

    from cellcomplex.property_topomesh.example_topomesh import circle_voronoi_topomesh
    from cellcomplex.property_topomesh.utils.matplotlib_tools import mpl_draw_topomesh, mpl_draw_incidence_graph

    topomesh = circle_voronoi_topomesh(3)

    figure = plt.figure(0)
    figure.clf()
    figure.gca().axis('equal')

    mpl_draw_topomesh(topomesh,figure,2,color='g')
    mpl_draw_topomesh(topomesh,figure,1,color='b')
    mpl_draw_topomesh(topomesh,figure,0,color='m')

    figure.gca().axis('off')
    figure.tight_layout()

A straightforward property to compute is for instance the length of all the edges. This computation is performed very simply by calling the :func:`compute_topomesh_property <cellcomplex.property_topomesh.analysis.compute_topomesh_property>` on a ``PropertyTopomesh`` object with the arguments 1 for **dimension** and ``'length'`` for **property name**. Being one of the pre-defined property names, the property ``'length'`` will be added to the ``wisps`` of dimension 1 and the computed length value will be assigned to each ID.

.. note::
    The list of pre-defined (computable) property names can be found at the bottom of this page, or in the documentation of the function :func:`compute_topomesh_property <cellcomplex.property_topomesh.analysis.compute_topomesh_property>`

.. code-block:: python
    :linenos:

    from cellcomplex.property_topomesh.analysis import compute_topomesh_property

    compute_topomesh_property(topomesh,'length',1)

.. plot::

    import numpy as np
    np.random.seed(583)

    import matplotlib.pyplot as plt

    from cellcomplex.property_topomesh.example_topomesh import circle_voronoi_topomesh
    from cellcomplex.property_topomesh.analysis import compute_topomesh_property
    from cellcomplex.property_topomesh.utils.matplotlib_tools import mpl_draw_topomesh, mpl_draw_incidence_graph

    topomesh = circle_voronoi_topomesh(3)
    compute_topomesh_property(topomesh,'length',1)

    figure = plt.figure(0)
    figure.clf()
    figure.gca().axis('equal')

    col = mpl_draw_topomesh(topomesh,figure,1,property_name='length',colormap='viridis',linewidth=3)
    mpl_draw_topomesh(topomesh,figure,0,color='k',size=20)

    figure.colorbar(col)

    figure.gca().axis('off')
    figure.tight_layout()

Similarly, it is possible to compute the area of the faces with just one function call. The method takes advantage of the storage of properties in NumPy based structures to perform computations efficiently on whole arrays of values instead of looping on the mesh elements.

.. code-block:: python
    :linenos:

    compute_topomesh_property(topomesh,'area',2)
    print(topomesh.wisp_property('area',2))

.. code-block:: bash

    {0: 0.9092855414754071,
     1: 0.96649956740289,
     2: 0.6270843768350033,
     3: 1.3496005090046475,
     4: 1.0153916778136718,
     5: 1.425343737120764,
     6: 0.863140023276767,
     7: 0.688484985614696,
     8: 1.0622081816602609,
     9: 0.7054585275907782,
     10: 1.1532886832605143,
     11: 1.5323925705923283,
     12: 1.0469905848920904,
     13: 0.8364178763830168,
     14: 0.9512681918544135,
     15: 1.1161306160901407,
     16: 1.396702833427955,
     17: 0.8403108600291247,
     18: 1.161257331179229}

.. plot::

    import numpy as np
    np.random.seed(583)

    import matplotlib.pyplot as plt

    from cellcomplex.property_topomesh.example_topomesh import circle_voronoi_topomesh
    from cellcomplex.property_topomesh.analysis import compute_topomesh_property
    from cellcomplex.property_topomesh.utils.matplotlib_tools import mpl_draw_topomesh, mpl_draw_incidence_graph

    topomesh = circle_voronoi_topomesh(3)
    compute_topomesh_property(topomesh,'area',2)

    figure = plt.figure(0)
    figure.clf()
    figure.gca().axis('equal')

    col = mpl_draw_topomesh(topomesh,figure,2,property_name='area',colormap='viridis')
    mpl_draw_topomesh(topomesh,figure,2,color='k',alpha=0,plot_ids=True)
    mpl_draw_topomesh(topomesh,figure,1,color='k',linewidth=0.5)
    mpl_draw_topomesh(topomesh,figure,0,color='k',size=10)

    figure.colorbar(col)

    figure.gca().axis('off')
    figure.tight_layout()

The  :func:`compute_topomesh_property <cellcomplex.property_topomesh.analysis.compute_topomesh_property>` function offers the possibility to compute a great number of numerical properties on each element of the mesh. The following lines show for example how to compute the **valency** of all vertices in the mesh, the **normal vectors** on the contour edges of a 2d shape, or the **inertia tensors** describing the principal directions of the mesh faces.

.. code-block:: python
    :linenos:

    compute_topomesh_property(topomesh,'valency',0)
    compute_topomesh_property(topomesh,'normal',1)
    compute_topomesh_property(topomesh,'inertia_tensor',2)

.. plot::

    import numpy as np
    np.random.seed(583)

    import matplotlib.pyplot as plt

    from cellcomplex.property_topomesh.example_topomesh import circle_voronoi_topomesh
    from cellcomplex.property_topomesh.analysis import compute_topomesh_property, compute_topomesh_vertex_property_from_edges
    from cellcomplex.property_topomesh.utils.matplotlib_tools import mpl_draw_topomesh, mpl_draw_incidence_graph

    topomesh = circle_voronoi_topomesh(3)
    compute_topomesh_property(topomesh,'valency',0)
    compute_topomesh_property(topomesh,'normal',1)
    compute_topomesh_vertex_property_from_edges(topomesh,'normal')
    compute_topomesh_property(topomesh,'inertia_tensor',2)

    figure = plt.figure(0)
    figure.clf()
    figure.gca().axis('equal')

    mpl_draw_topomesh(topomesh,figure,2,property_name='inertia_tensor',color='b',coef=1,linewidth=2)
    mpl_draw_topomesh(topomesh,figure,1,color='k',linewidth=0.5)
    mpl_draw_topomesh(topomesh,figure,1,property_name='normal',color='r',coef=0.5)
    col = mpl_draw_topomesh(topomesh,figure,0,property_name='valency',colormap='jet',color='k',size=40,linewidth=0.5,intensity_range=(0,4))

    cbar = figure.colorbar(col)
    cbar.set_label('Vertex Valency')

    figure.gca().axis('off')
    figure.tight_layout()


**********************
3D Example : Curvature
**********************

Now, let's consider a 3D example. The :mod:`example_topomesh <cellcomplex.property_topomesh.example_topomesh>` module provides several simple 3D shapes that can be represented as a :class:`PropertyTopomesh`. Here, we generate a 3D ellipsoid elongated along the X axis, represented as a triangular mesh, and keep only the upper (Z>0) part of the surface.

.. code-block:: python

    from cellcomplex.property_topomesh.example_topomesh import vtk_ellipsoid_topomesh
    from cellcomplex.property_topomesh.extraction import cut_surface_topomesh
    from cellcomplex.property_topomesh.utils.matplotlib_tools import mpl_draw_topomesh, mpl_draw_incidence_graph

    topomesh = vtk_ellipsoid_topomesh(ellipsoid_radius=2,ellipsoid_scales=[1.5,1,1],subdivisions=2)
    topomesh = cut_surface_topomesh(topomesh,z_offset=0,below=False)

.. plot::

    import os
    has_x = ('DISPLAY' in os.environ.keys()) and (os.environ['DISPLAY'] is not None)

    import numpy as np

    from cellcomplex.property_topomesh.example_topomesh import vtk_ellipsoid_topomesh
    from cellcomplex.property_topomesh.extraction import cut_surface_topomesh
    from cellcomplex.property_topomesh.utils.matplotlib_tools import mpl_draw_topomesh, mpl_draw_incidence_graph

    from cellcomplex.property_topomesh.visualization.vtk_actor_topomesh import VtkActorTopomesh
    from cellcomplex.property_topomesh.visualization.vtk_tools import vtk_image_actors

    topomesh = vtk_ellipsoid_topomesh(ellipsoid_radius=2,ellipsoid_scales=[1.5,1,1],subdivisions=2)
    topomesh = cut_surface_topomesh(topomesh,z_offset=0,below=False)

    figure = plt.figure(0)
    figure.clf()

    if has_x:
        figure.add_subplot(1,2,1)
    figure.gca().axis('equal')

    mpl_draw_topomesh(topomesh,figure,2,color='g')
    mpl_draw_topomesh(topomesh,figure,1,color='b')
    mpl_draw_topomesh(topomesh,figure,0,color='m')

    if has_x:
        figure.add_subplot(1,2,2)

        actors = []

        face_actor = VtkActorTopomesh(topomesh,2)
        face_actor.update(colormap='YlGn',value_range=(0,0))
        actors += [face_actor.actor]

        edge_actor = VtkActorTopomesh(topomesh,1,line_glyph='tube',glyph_scale=0.01)
        edge_actor.update(colormap='Blues',value_range=(0,0),linewidth=3,opacity=1)
        actors += [edge_actor.actor]

        vertex_actor = VtkActorTopomesh(topomesh,0,glyph_scale=0.05)
        vertex_actor.update(colormap='RdPu',value_range=(0,0),opacity=1)
        actors += [vertex_actor.actor]

        figure.gca().imshow(vtk_image_actors(actors,focal_point=(0,0,1),view_up=(0,1,0)))

        figure.gca().axis('off')

        figure.set_size_inches(10,5)
    figure.tight_layout()

Based on this mesh, we would like to compute the curvature of the surface. To do so, we will use the function :func:`compute_topomesh_property <cellcomplex.property_topomesh.analysis.compute_topomesh_property>` to compute the geometrical properties allowing to obtain this information.

The first property to compute is the **normal** of the faces (elements of degree = 2). Given the implementation as an incidence graph, the order of vertices for each face is not unambiguous, hence the orientation of faces is not contained in the surface. It is necessary to have a way of re-orienting consistently the normals. Here we choose the orientation mode that propagates the orientation infomration in a topologically accurate way on a connected surfacic mesh.

.. code-block:: python
    :linenos:

    compute_topomesh_property(topomesh,'normal',2,normal_method='orientation')


Next, we compute the area of all faces, and use it to compute the normal property on the vertices (elements of degree = 0). This is done by averaging the normals of faces neighboring each vertex using the function :func:`compute_topomesh_vertex_property_from_faces <cellcomplex.property_topomesh.analysis.compute_topomesh_vertex_property_from_faces>`. This function is a generic tool that allows to propagate properties between elements of different dimensions.

.. code-block:: python
    :linenos:

    compute_topomesh_property(topomesh,'area',2)
    compute_topomesh_vertex_property_from_faces(topomesh,'normal',weighting='area',adjacency_sigma=1.2,neighborhood=3)


Once each vertex bears a consistent normal, it is possible to compute the **mean curvature** of the each face, using a discrete curvature estimator. This property can be propagated to the vertices using the same method as for the normals.

.. code-block:: python
    :linenos:

    compute_topomesh_property(topomesh,'mean_curvature',2)
    compute_topomesh_vertex_property_from_faces(topomesh,'mean_curvature',weighting='area',adjacency_sigma=1.2,neighborhood=3)


.. plot::

    import os
    has_x = ('DISPLAY' in os.environ.keys()) and (os.environ['DISPLAY'] is not None)

    import numpy as np
    import matplotlib.pyplot as plt

    from cellcomplex.property_topomesh.example_topomesh import vtk_ellipsoid_topomesh
    from cellcomplex.property_topomesh.extraction import cut_surface_topomesh
    from cellcomplex.property_topomesh.analysis import compute_topomesh_property, compute_topomesh_vertex_property_from_faces
    from cellcomplex.property_topomesh.utils.matplotlib_tools import mpl_draw_topomesh, mpl_draw_incidence_graph

    from cellcomplex.property_topomesh.visualization.vtk_actor_topomesh import VtkActorTopomesh
    from cellcomplex.property_topomesh.visualization.vtk_tools import vtk_image_actors

    topomesh = vtk_ellipsoid_topomesh(ellipsoid_radius=2,ellipsoid_scales=[1.5,1,1],subdivisions=2)
    topomesh = cut_surface_topomesh(topomesh,z_offset=0,below=False)

    compute_topomesh_property(topomesh,'normal',2,normal_method='orientation')
    compute_topomesh_property(topomesh,'area',2)
    compute_topomesh_vertex_property_from_faces(topomesh,'normal',weighting='area',adjacency_sigma=1.2,neighborhood=3)
    compute_topomesh_property(topomesh,'mean_curvature',2)

    figure = plt.figure(0)
    figure.clf()

    if has_x:
        figure.add_subplot(1,2,1)
    figure.gca().axis('equal')

    col = mpl_draw_topomesh(topomesh,figure,2,property_name='mean_curvature',colormap='Reds',intensity_range=(0.25,1))
    mpl_draw_topomesh(topomesh,figure,0,property_name='normal',color='g',coef=0.5,linewidth=1)
    mpl_draw_topomesh(topomesh,figure,1,color='k',linewidth=0.5)
    mpl_draw_topomesh(topomesh,figure,0,color='k',size=10)

    cbar = figure.colorbar(col)
    cbar.set_label('Face Mean Curvature')

    if has_x:
        figure.add_subplot(1,2,2)

        actors = []

        curvature_actor = VtkActorTopomesh(topomesh,2,property_name='mean_curvature')
        curvature_actor.update(colormap='Reds',value_range=(0.25,1))
        actors += [curvature_actor.actor]

        edge_actor = VtkActorTopomesh(topomesh,1,line_glyph='tube',glyph_scale=0.01)
        edge_actor.update(colormap='Greys',value_range=(0,0),linewidth=3,opacity=1)
        actors += [edge_actor.actor]

        vertex_actor = VtkActorTopomesh(topomesh,0,glyph_scale=0.05)
        vertex_actor.update(colormap='Greys',value_range=(0,0),opacity=1)
        actors += [vertex_actor.actor]

        normal_actor = VtkActorTopomesh(topomesh,0,property_name='normal',vector_glyph='arrow',glyph_scale=0.5)
        normal_actor.update(colormap='YlGn')
        actors += [normal_actor.actor]

        figure.gca().imshow(vtk_image_actors(actors,focal_point=(0,0,1),view_up=(0,1,0)))

        figure.gca().axis('off')

        figure.set_size_inches(10,5)
    figure.tight_layout()

The numerical property  defined on vertices can then be visualized on the triangles of the mesh, creating a smooth interpolation of curvature values.

While estimating the **mean curvature** of the faces, the algorithm also has filled other properties such as the 3D **curvature tensor** containing all the information on the principal curvatures, in particular mean, gaussian, min and max curvatures. This can be useful for instance to visualize principal directions of curvature on the faces. As a tensor, it is stored as a 2d NumPy array for each ``wisp`` of dimension 2.

.. code-block:: python
    :linenos:

    print("Face 0 : ",topomesh.wisp_property('principal_curvature_tensor',2)[0])

.. code-block:: bash

    Face 0 :  [[ 2.41008840e-01 -6.01676202e-17 -2.43084355e-17]
               [-1.23938667e-17  3.13548659e-01  2.67055929e-01]
               [ 2.88344294e-18  2.67088173e-01  2.27484565e-01]]

.. plot::

    import os
    has_x = ('DISPLAY' in os.environ.keys()) and (os.environ['DISPLAY'] is not None)

    import numpy as np
    import matplotlib.pyplot as plt

    from cellcomplex.property_topomesh.example_topomesh import vtk_ellipsoid_topomesh
    from cellcomplex.property_topomesh.extraction import cut_surface_topomesh
    from cellcomplex.property_topomesh.analysis import compute_topomesh_property, compute_topomesh_vertex_property_from_faces

    from cellcomplex.property_topomesh.utils.matplotlib_tools import mpl_draw_topomesh, mpl_draw_incidence_graph
    from cellcomplex.property_topomesh.visualization.vtk_actor_topomesh import VtkActorTopomesh
    from cellcomplex.property_topomesh.visualization.vtk_tools import vtk_image_actors

    topomesh = vtk_ellipsoid_topomesh(ellipsoid_radius=2,ellipsoid_scales=[1.5,1,1],subdivisions=2)
    topomesh = cut_surface_topomesh(topomesh,z_offset=0,below=False)

    compute_topomesh_property(topomesh,'normal',2,normal_method='orientation')
    compute_topomesh_property(topomesh,'area',2)
    compute_topomesh_vertex_property_from_faces(topomesh,'normal',weighting='area',adjacency_sigma=1.2,neighborhood=3)
    compute_topomesh_property(topomesh,'mean_curvature',2)
    compute_topomesh_vertex_property_from_faces(topomesh,'mean_curvature',weighting='area',adjacency_sigma=1.2,neighborhood=3)

    figure = plt.figure(0)
    figure.clf()

    if has_x:
        figure.add_subplot(1,2,1)
    figure.gca().axis('equal')

    col = mpl_draw_topomesh(topomesh,figure,2,property_name='mean_curvature',property_degree=0,colormap='Reds',intensity_range=(0.25,1))
    mpl_draw_topomesh(topomesh,figure,1,color='k',linewidth=0.5)
    mpl_draw_topomesh(topomesh,figure,2,property_name='principal_curvature_tensor',color='darkred',coef=0.25,linewidth=1)
    mpl_draw_topomesh(topomesh,figure,0,property_name='mean_curvature',color='k',colormap='Reds',intensity_range=(0.25,1),size=20,linewidth=0.5)

    cbar = figure.colorbar(col)
    cbar.set_label('Vertex Mean Curvature')


    if has_x:
        figure.add_subplot(1,2,2)

        actors = []

        curvature_actor = VtkActorTopomesh(topomesh,2,property_name='mean_curvature',property_degree=0)
        curvature_actor.update(colormap='Reds',value_range=(0.25,1))
        actors += [curvature_actor.actor]

        edge_actor = VtkActorTopomesh(topomesh,1,line_glyph='tube',glyph_scale=0.01)
        edge_actor.update(colormap='Greys',value_range=(0,0),linewidth=3,opacity=1)
        actors += [edge_actor.actor]

        vertex_actor = VtkActorTopomesh(topomesh,0,property_name='mean_curvature',glyph_scale=0.05)
        vertex_actor.update(colormap='Reds',value_range=(0.25,1))
        actors += [vertex_actor.actor]

        tensor_actor = VtkActorTopomesh(topomesh,2,property_name='principal_curvature_tensor',tensor_glyph='crosshair',glyph_scale=0.25)
        tensor_actor.update(colormap='Reds_r',linewidth=5,opacity=1)
        actors += [tensor_actor.actor]

        figure.gca().imshow(vtk_image_actors(actors,focal_point=(0,0,1),view_up=(0,1,0)))
        figure.gca().axis('off')

        figure.set_size_inches(10,5)
    figure.tight_layout()


*********************************
Computable geometrical properties
*********************************

A list of **geometrical** properties that can be passed as `property_name` argument to the function :func:`compute_topomesh_property <cellcomplex.property_topomesh.analysis.compute_topomesh_property>`, along with the `wisp` dimension for which they can be computed:


+---------------------+---+---+---+---+----------+
| property_name       | 0 | 1 | 2 | 3 |   type   |
+=====================+===+===+===+===+==========+
|``barycenter``       | o | x | x | x |  vector  |
+---------------------+---+---+---+---+----------+
|``length``           |   | x |   |   |  scalar  |
+---------------------+---+---+---+---+----------+
|``area``             |   |   | x |   |  scalar  |
+---------------------+---+---+---+---+----------+
|``volume``           |   |   |   | x |  scalar  |
+---------------------+---+---+---+---+----------+
|``normal``           | x | x | x |   |  vector  |
+---------------------+---+---+---+---+----------+
|``eccentricity``     |   |   | x |   |  scalar  |
+---------------------+---+---+---+---+----------+
|``skewness``         |   |   | x |   |  scalar  |
+---------------------+---+---+---+---+----------+
|``roughness``        | x |   |   |   |  scalar  |
+---------------------+---+---+---+---+----------+
|``valency``          | x |   |   |   |  scalar  |
+---------------------+---+---+---+---+----------+
|``mean_curvature``   |   |   | x |   |  scalar  |
+---------------------+---+---+---+---+----------+
|``inertia_axes``     |   |   | x |   |  tensor  |
+---------------------+---+---+---+---+----------+
|``inertia_moments``  |   |   | x |   |  vector  |
+---------------------+---+---+---+---+----------+