Note

This page was generated from 7_3d_mesh_digitization.ipynb. Interactive online version: Colab badge. Some tutorial content may look better in light mode.

7. Digitization of MESH surface in 3D space

This notebook demonstrates how to use digitize_general on a point network, like MESH.

Load packages

[1]:
import sys
import spateo as st

import numpy as np
import pandas as pd
pd.options.mode.chained_assignment = None

import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.rcParams['pdf.fonttype'] = 42

import pyvista as pv
pv.start_xvfb()
2024-08-17 14:07:29.320758: I tensorflow/core/util/port.cc:111] oneDNN custom operations are on. You may see slightly different numerical results due to floating-point round-off errors from different computation orders. To turn them off, set the environment variable `TF_ENABLE_ONEDNN_OPTS=0`.
2024-08-17 14:07:29.352262: E tensorflow/compiler/xla/stream_executor/cuda/cuda_dnn.cc:9342] Unable to register cuDNN factory: Attempting to register factory for plugin cuDNN when one has already been registered
2024-08-17 14:07:29.352292: E tensorflow/compiler/xla/stream_executor/cuda/cuda_fft.cc:609] Unable to register cuFFT factory: Attempting to register factory for plugin cuFFT when one has already been registered
2024-08-17 14:07:29.352312: E tensorflow/compiler/xla/stream_executor/cuda/cuda_blas.cc:1518] Unable to register cuBLAS factory: Attempting to register factory for plugin cuBLAS when one has already been registered
2024-08-17 14:07:29.358475: I tensorflow/core/platform/cpu_feature_guard.cc:182] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations.
To enable the following instructions: AVX2 AVX512F AVX512_VNNI FMA, in other operations, rebuild TensorFlow with the appropriate compiler flags.
2024-08-17 14:07:30.067520: W tensorflow/compiler/tf2tensorrt/utils/py_utils.cc:38] TF-TRT Warning: Could not find TensorRT
2024-08-17 14:07:31.917660: I tensorflow/core/common_runtime/gpu/gpu_device.cc:1886] Created device /job:localhost/replica:0/task:0/device:GPU:0 with 21386 MB memory:  -> device: 0, name: NVIDIA GeForce RTX 3090, pci bus id: 0000:65:00.0, compute capability: 8.6
/home/ftnws/anaconda3/envs/spateo_github/lib/python3.10/site-packages/POT-0.9.1-py3.10-linux-x86_64.egg/ot/backend.py:2998: UserWarning: To use TensorflowBackend, you need to activate the tensorflow numpy API. You can activate it by running:
from tensorflow.python.ops.numpy_ops import np_config
np_config.enable_numpy_behavior()
  register_backend(TensorflowBackend())

Data source

adata_tube: A manually extracted subset of E11.5 mouse embryo near diencephalic neural tube

[3]:
adata_tube=st.read_h5ad("adata_dien_clip.h5ad")

Construct MESH for an diencephalic subset of the E11.5 mouse embryo

[4]:
pc_model, _ = st.tdr.construct_pc(adata=adata_tube.copy(),spatial_key="spatial",groupby="mapped_celltype", key_added="mapped_celltype",colormap="rainbow")

st.pl.three_d_plot(
    model=pc_model,
    key="mapped_celltype",
    model_style="points",
    model_size=3,
    jupyter="static",opacity=0.2,
)
../../../_images/tutorials_notebooks_5_cluster_digitization_7_3d_mesh_digitization_6_0.png
[5]:
# Subset the point cloud to reduce memory consumption
if pc_model.n_points > 500000:
        np.random.seed(19491001)
        sampling = np.random.choice(
            np.asarray(pc_model.point_data["obs_index"]),
            size=100000, #decrease this if memory limit exceeded
            replace=False,
        )
        downsampling_pc_model = pc_model.extract_points(np.isin(np.asarray(pc_model.point_data["obs_index"]), sampling))
else:
    downsampling_pc_model = pc_model.copy()
[8]:
# Here we use a manually defined ellipsoid to cut off the inner surface of the neural tube.
# Use interactive or any selection methods in your case.

ellipsoid = pv.ParametricEllipsoid(xradius = 190 , yradius = 139 , zradius = 120)
ellipsoid.point_data["mapped_celltype"] = np.asarray(["Ellipsoid"] * ellipsoid.n_points)
ellipsoid.points  = np.asarray(ellipsoid.points) - np.asarray(np.asarray(ellipsoid.center) - np.asarray(pc_model.center) + np.array([-25, 65, -80]))

st.pl.three_d_plot(
    model=st.tdr.collect_models([ellipsoid, mesh_model]),
    key="mapped_celltype",
    model_style=["surface", "surface"],
    colormap=["yellowgreen", "gainsboro"],
    opacity=[0.6, 0.6],
    model_size=3,
    jupyter="static",
)
../../../_images/tutorials_notebooks_5_cluster_digitization_7_3d_mesh_digitization_8_0.png
[14]:
select_mesh_model = mesh_model.select_enclosed_points(ellipsoid, check_surface=False)
inside_mesh_model = select_mesh_model.threshold(0.5, scalars="SelectedPoints").extract_surface()

st.pl.three_d_plot(
    model=st.tdr.collect_models([inside_mesh_model]),
    key="mapped_celltype",
    model_style=["surface", "points"],
    model_size=3,
    jupyter="static",
)
../../../_images/tutorials_notebooks_5_cluster_digitization_7_3d_mesh_digitization_9_0.png

Digitization on MESH

[68]:
# Here we use a coordinate cutoff to demonstrate the selection of digitization boundary.
# Use interactive or any selection methods in your case.

df_mesh_net = pd.DataFrame({
    'x':np.array(inside_mesh_model.points)[:,0].flatten(),
    'y':np.array(inside_mesh_model.points)[:,1].flatten(),
    'z':np.array(inside_mesh_model.points)[:,2].flatten(),
    # Here we use a coordinate cutoff for demonstration. Use interactive selection in your case.
    'point':1+((np.abs(np.array(inside_mesh_model.points)[:,1] - 325)<15).astype(int)-
    2*(np.abs(np.array(inside_mesh_model.points)[:,0] - 1420)<25).astype(int))
})

fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
ax.scatter(
    xs=df_mesh_net.x,
    ys=df_mesh_net.y,
    zs=df_mesh_net.z,
    c=df_mesh_net.point,
    s=2,
    cmap="RdYlBu_r",
)
ax.grid(False)
ax.view_init(elev=80, azim=180)
../../../_images/tutorials_notebooks_5_cluster_digitization_7_3d_mesh_digitization_11_0.png
[27]:
# Here we use a coordinate cutoff for demonstration. Use interactive or any selection methods in your case.
boundary_lower = np.where(np.abs(np.array(inside_mesh_model.points)[:,0] - 1420)<25)
boundary_upper = np.where(np.abs(np.array(inside_mesh_model.points)[:,1] - 325)<15)
[ ]:
# Extract the adjacency matrix from mesh network
adj_mtx = np.zeros((len(inside_mesh_model.points), len(inside_mesh_model.points)))

faces = inside_mesh_model.faces
ind = 0
face_degrees = []

while ind < len(faces):
    adj_mtx[faces[ind+1], faces[ind+faces[ind]]] = 1
    adj_mtx[faces[ind+faces[ind]], faces[ind+1]] = 1

    face_degrees.append(faces[ind])

    for i in range(1, faces[ind]):
        ind = ind + 1
        adj_mtx[faces[ind], faces[ind+1]] = 1
        adj_mtx[faces[ind+1], faces[ind]] = 1

    ind = ind + 2

conv_mtx = np.apply_along_axis(lambda x: x/np.sum(x), axis=0, arr=adj_mtx)
[ ]:
# Perform digitization on the mesh network
grid_field = st.dd.utils.digitize_general(
    pc=inside_mesh_model.points,
    adj_mtx=conv_mtx,
    boundary_lower=boundary_lower,
    boundary_upper=boundary_upper,
)
[75]:
# Visualize the digitized mesh

fig, ax = plt.subplots(nrows=1, ncols=2, subplot_kw={"projection": "3d"})

ax[0].scatter(
    xs=df_mesh_net.x,
    ys=df_mesh_net.y,
    zs=df_mesh_net.z,
    c=grid_field,
    s=2,
    cmap="RdYlBu_r",
)
ax[0].grid(False)
ax[0].view_init(elev=80, azim=180)

ax[1].scatter(
    xs=df_mesh_net.x,
    ys=df_mesh_net.y,
    zs=df_mesh_net.z,
    c=(grid_field - 1) // 20,
    s=2,
    cmap="RdYlBu_r",
)
ax[1].grid(False)
ax[1].view_init(elev=85, azim=180)
../../../_images/tutorials_notebooks_5_cluster_digitization_7_3d_mesh_digitization_15_0.png