Anisotropic stresses#

This notebook demonstrates the calculation of anisotropic stresses of injected oil droplets. The measure properties included the following:

  • \(\sigma^A(\vec{x}_s)\): anisotropic stress on surface

  • \(\sigma^A_T(\vec{x}_s)\): anisotropic stress on fitted ellipsoid (e.g. tissue-level)

import napari
import napari_stress
from napari_stress import measurements, approximation
import matplotlib.pyplot as plt

import numpy as np
napari.manifest -> 'napari' could not be imported: Could not find file 'builtins.yaml' in module 'napari'

config#

We need to set some parameters:

max_degree = 5
gamma = 26.0
# Get some sample data
pointcloud = napari_stress.get_droplet_point_cloud()[0]
viewer = napari.Viewer(ndisplay=3)
napari.manifest -> 'napari' could not be imported: Could not find file 'builtins.yaml' in module 'napari'
layer_raw = viewer.add_layer(napari.layers.Layer.create(data = pointcloud[0], meta=pointcloud[1], layer_type=pointcloud[2]))
napari.utils.nbscreenshot(viewer)

Ellipsoid fit#

First, we fit an ellipsoid to the pointcloud. The ellipsoid represents the tissue-scale stresses.

expander =  approximation.EllipsoidExpander()
expander.fit(pointcloud[0][:, 1:])
fitted_ellipsoid = expander.coefficients_
fitted_ellipse_points = expander.fit_expand(pointcloud[0][:, 1:])

We then measure the curvatures on every point on the surface:

curvature_on_ellipsoid = measurements.curvature_on_ellipsoid(fitted_ellipsoid, fitted_ellipse_points)
H_ellipsoid_major_medial_minor = measurements.mean_curvature_on_ellipse_cardinal_points(fitted_ellipsoid)
H_ellipsoid_major_medial_minor
[0.08211623521117135, 0.06563346684868424, 0.055253239527228264]
layer_raw.visible = False
ellipsoid_layer = viewer.add_vectors(fitted_ellipsoid)
viewer.add_points(fitted_ellipse_points, features = curvature_on_ellipsoid[1]['features'], face_color = 'mean_curvature', size=0.5, face_colormap = 'jet')
napari.utils.nbscreenshot(viewer)

We now do a spherical harmonics expansion of the ellipsoid points. From this, we can then calculate the average mean curvature \(H_0\) and the mean curvature \(H_i\) at every point on the surface.

fitted_ellipse_points_spherical_harmoniocs = napari_stress.fit_spherical_harmonics(fitted_ellipse_points, max_degree=max_degree)
coefficients = fitted_ellipse_points_spherical_harmoniocs[1]['metadata']['spherical_harmonics_coefficients']
quadrature_points, lebedev_fit = napari_stress.lebedev_quadrature(coefficients=coefficients, number_of_quadrature_points=500, use_minimal_point_set=False)
manifold_ellipsoid = napari_stress.create_manifold(quadrature_points, lebedev_fit=lebedev_fit, max_degree=max_degree)

Note: \(H_0\) is calculated using the spherical-harmonics related curvature definition. The point-wise curvature \(H_i\) is derived by projecting the quadrature points on the surface of the ellipsoid and using the respective measurement function:

_, _, H0_ellipsoid = measurements.calculate_mean_curvature_on_manifold(manifold_ellipsoid)

quadrature_points_on_ellipsoid = approximation.expand_points_on_ellipse(fitted_ellipsoid, quadrature_points)
mean_curvature_ellipsoid = measurements.curvature_on_ellipsoid(fitted_ellipsoid, quadrature_points_on_ellipsoid)[1]['features']['mean_curvature']

Spherical harmonics expansion#

We now do a spherical harmonics expansion of the raw pointcloud. This also captures the higher orders of curvature on the surface of the droplet.

fitted_points = napari_stress.fit_spherical_harmonics(pointcloud[0], max_degree=max_degree)
fitted_points_layer = viewer.add_points(fitted_points[0], **fitted_points[1], name='spherical harmonics expansion')
coefficients = fitted_points[1]['metadata']['spherical_harmonics_coefficients']
quadrature_points, lebedev_fit = napari_stress.lebedev_quadrature(coefficients=coefficients, number_of_quadrature_points=500, use_minimal_point_set=False)
manifold_droplet = napari_stress.create_manifold(quadrature_points, lebedev_fit=lebedev_fit, max_degree=max_degree)
mean_curvature_spherical_harmonics, _, H0_surface_integral_spherical_harmonics = measurements.calculate_mean_curvature_on_manifold(manifold_droplet)

Stresses#

We first calculate the tissue stress tensors along the ajor axes of the ellipsoid and the cartesian cardinal directions (x/y/z):

ellptical, cartesian = measurements.tissue_stress_tensor(fitted_ellipsoid, H0_ellipsoid, gamma=gamma)

Secondly, we calculate the anisotropic stress for every point on the surface:

stress, stress_tissue, stress_droplet = measurements.anisotropic_stress(mean_curvature_spherical_harmonics, H0_surface_integral_spherical_harmonics,
                                                                        mean_curvature_ellipsoid, H0_ellipsoid,
                                                                        gamma=gamma)
nbins = 20
fig, axes = plt.subplots(ncols=3, figsize=(15,5))
axes[0].hist(stress, nbins)
axes[0].set_title('Stress')
axes[1].hist(stress_tissue, nbins)
axes[1].set_title('tissue-scale stresses')
axes[2].hist(stress_droplet, nbins)
axes[2].set_title('Cell scale stresses')

for ax in axes:
    ax.set_xlabel('anisotropic stress [mN/m]')
    ax.set_ylabel('Occurrences [#]')
../../_images/286df67084643ff1ad2e2b7aff6208c71772d2624742201b12030ee2935492ef.png

Lastly, we can also calculate the anisotropic tissue stressses \(\sigma^A_T\) from the fitted ellipsoid above:

ansiostropic_tissue_stress = measurements.maximal_tissue_anisotropy(fitted_ellipsoid, gamma=gamma)
ansiostropic_tissue_stress
1.396875775565042