{ "cells": [ { "cell_type": "markdown", "id": "3cdf478e", "metadata": {}, "source": [ "# Cell Nuclei\n", "\n", "In this example, we will use a cellular image from the Allen Cell WTC-11 hiPSC Single-Cell Image Dataset ([Viana et al. 2023](https://doi.org/10.1038/s41586-022-05563-7))." ] }, { "cell_type": "code", "execution_count": 1, "id": "59e8a2c3", "metadata": { "execution": { "iopub.execute_input": "2025-05-13T11:30:25.120269Z", "iopub.status.busy": "2025-05-13T11:30:25.120166Z", "iopub.status.idle": "2025-05-13T11:30:25.533939Z", "shell.execute_reply": "2025-05-13T11:30:25.533481Z" }, "vscode": { "languageId": "plaintext" } }, "outputs": [], "source": [ "import libcarna\n", "import numpy as np\n", "import scipy.ndimage as ndi" ] }, { "cell_type": "markdown", "id": "a3157707", "metadata": {}, "source": [ "Get the data:" ] }, { "cell_type": "code", "execution_count": 2, "id": "97c2c4ed", "metadata": { "execution": { "iopub.execute_input": "2025-05-13T11:30:25.535447Z", "iopub.status.busy": "2025-05-13T11:30:25.535270Z", "iopub.status.idle": "2025-05-13T11:30:25.596255Z", "shell.execute_reply": "2025-05-13T11:30:25.595909Z" }, "vscode": { "languageId": "plaintext" } }, "outputs": [ { "data": { "text/plain": [ "((60, 256, 256), dtype('uint16'))" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data = libcarna.data.nuclei()\n", "data.shape, data.dtype" ] }, { "cell_type": "markdown", "id": "fff4fac6", "metadata": {}, "source": [ "The data is 60 × 256 × 256 pixels (uint16).\n", "\n", "## Maximum Intensity Projection\n", "\n", "In the code below, we use `normals=True` on the `volume` node; albeit this doesn't make any difference for the *Maximum\n", "Intensity Projection* (MIP), it is benefitial for other rendering modes, [discussed below](#Direct-Volume-Rendering)." ] }, { "cell_type": "code", "execution_count": 3, "id": "ed16156a", "metadata": { "execution": { "iopub.execute_input": "2025-05-13T11:30:25.597390Z", "iopub.status.busy": "2025-05-13T11:30:25.597224Z", "iopub.status.idle": "2025-05-13T11:30:26.469406Z", "shell.execute_reply": "2025-05-13T11:30:26.469030Z" }, "vscode": { "languageId": "plaintext" } }, "outputs": [ { "data": { "text/html": [ "\n", "
\n", "
\n", " \n", "
\n", " \n", "
\n", "
\n", "
\n", "\n", "
\n", "
\n", "
\n", " \n", "
\n", " 66k\n", " \n", "
\n", "\n", "
\n", " 49k\n", " \n", "
\n", "\n", "
\n", " 33k\n", " \n", "
\n", "\n", "
\n", " 16k\n", " \n", "
\n", "\n", "
\n", " 0\n", " \n", "
\n", "
\n", "
\n", " \n", "
\n", "
\n", " " ], "text/plain": [ "" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "GEOMETRY_TYPE_VOLUME = 1\n", "\n", "# Create and configure frame renderer\n", "mip = libcarna.mip(GEOMETRY_TYPE_VOLUME, cmap='jet', sr=500)\n", "r = libcarna.renderer(600, 450, [mip])\n", "\n", "# Create and configure scene\n", "root = libcarna.node()\n", "\n", "volume = libcarna.volume(\n", " GEOMETRY_TYPE_VOLUME,\n", " data,\n", " parent=root,\n", " spacing=(1, 0.5, 0.5),\n", ").rotate('x', -35).rotate('y', 90)\n", "\n", "camera = libcarna.camera(\n", " parent=root,\n", ").frustum(fov=90, z_near=1, z_far=500).translate(z=100)\n", "\n", "# Render\n", "libcarna.imshow(r.render(camera), mip.cmap.bar(volume))" ] }, { "cell_type": "markdown", "id": "465034ee", "metadata": {}, "source": [ "In the MIP, it can easily be seen that there is one mitotic nucleus in the image.\n", "\n", "For an even better visual perception of the 3D data, it is best viewed from different angles, that can be achieved with\n", "as a subtle animation:" ] }, { "cell_type": "code", "execution_count": 4, "id": "d35452de", "metadata": { "execution": { "iopub.execute_input": "2025-05-13T11:30:26.472115Z", "iopub.status.busy": "2025-05-13T11:30:26.471993Z", "iopub.status.idle": "2025-05-13T11:30:26.889987Z", "shell.execute_reply": "2025-05-13T11:30:26.889554Z" }, "vscode": { "languageId": "plaintext" } }, "outputs": [ { "data": { "text/html": [ "\n", "
\n", "
\n", " \n", "
\n", " \n", "
\n", "
\n", "
\n", "\n", "
\n", "
\n", "
\n", " \n", "
\n", " 66k\n", " \n", "
\n", "\n", "
\n", " 49k\n", " \n", "
\n", "\n", "
\n", " 33k\n", " \n", "
\n", "\n", "
\n", " 16k\n", " \n", "
\n", "\n", "
\n", " 0\n", " \n", "
\n", "
\n", "
\n", " \n", "
\n", "
\n", " " ], "text/plain": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Render as animation\n", "libcarna.imshow(\n", " libcarna.animate(\n", " libcarna.animate.swing_local(camera, amplitude=22),\n", " n_frames=50,\n", " ).render(r, camera),\n", " mip.cmap.bar(volume),\n", ")" ] }, { "cell_type": "markdown", "id": "17ab03b3", "metadata": {}, "source": [ "This makes the camera swing by 22° to the left and to the right." ] }, { "cell_type": "markdown", "id": "ab9dc2e1", "metadata": {}, "source": [ "## Direct Volume Rendering\n", "\n", "In a *Direct Volume Rendering* (DVR), surfaces are rendered by simulation of the absorption of light. This simulation\n", "is most realstic, when the spatial orientation of the surfaces can be taken into account, which requires that the\n", "normals of the volume have been computed (this is why we used `normals=True` when we created the `volume` node)." ] }, { "cell_type": "code", "execution_count": 5, "id": "c33d38bc", "metadata": { "execution": { "iopub.execute_input": "2025-05-13T11:30:26.894851Z", "iopub.status.busy": "2025-05-13T11:30:26.894647Z", "iopub.status.idle": "2025-05-13T11:30:26.899321Z", "shell.execute_reply": "2025-05-13T11:30:26.898913Z" }, "vscode": { "languageId": "plaintext" } }, "outputs": [], "source": [ "dvr = libcarna.dvr(GEOMETRY_TYPE_VOLUME, sr=500)" ] }, { "cell_type": "markdown", "id": "0db629af", "metadata": {}, "source": [ "We again use the `jet` colormap, that, as we have seen in the MIP, employs blueish colors for the nuclear envelope, and\n", "reddish colors for the chromatin. In addition, we use a linear `ramp` function for the colormap, because we want the\n", "space between the nuclei to be translucent:" ] }, { "cell_type": "code", "execution_count": 6, "id": "ea6ed2dd", "metadata": { "execution": { "iopub.execute_input": "2025-05-13T11:30:26.901128Z", "iopub.status.busy": "2025-05-13T11:30:26.900781Z", "iopub.status.idle": "2025-05-13T11:30:27.319418Z", "shell.execute_reply": "2025-05-13T11:30:27.318973Z" }, "vscode": { "languageId": "plaintext" } }, "outputs": [ { "data": { "text/html": [ "\n", "
\n", "
\n", " \n", "
\n", " \n", "
\n", "
\n", "
\n", "\n", "
\n", "
\n", "
\n", " \n", "
\n", " 66k\n", " \n", "
\n", "\n", "
\n", " 49k\n", " \n", "
\n", "\n", "
\n", " 33k\n", " \n", "
\n", "\n", "
\n", " 16k\n", " \n", "
\n", "\n", "
\n", " 0\n", " \n", "
\n", "
\n", "
\n", " \n", "
\n", "
\n", " " ], "text/plain": [ "" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dvr.cmap('jet', ramp=(0.15, 0.25))\n", "\n", "libcarna.imshow(\n", " libcarna.animate(\n", " libcarna.animate.swing_local(camera, amplitude=22),\n", " n_frames=50,\n", " ).render(\n", " libcarna.renderer(600, 450, [dvr]),\n", " camera,\n", " ),\n", " dvr.cmap.bar(volume),\n", ")" ] }, { "cell_type": "markdown", "id": "74821fb1", "metadata": {}, "source": [ "The DVR of the nuclei in the image allows for a very natural perception of the 3D scene. On the downside, the mitotic\n", "nucleus is harder to identify. This is because the chromatin from the inside of the nucleus (should be reddish due to\n", "our colormap) is fully occluded by the nuclear envelope (blueish).\n", "\n", "We can overlay a DVR of the nuclear envelope with a MIP for the chromatin:" ] }, { "cell_type": "code", "execution_count": 7, "id": "b00a9d85", "metadata": { "execution": { "iopub.execute_input": "2025-05-13T11:30:27.322333Z", "iopub.status.busy": "2025-05-13T11:30:27.322209Z", "iopub.status.idle": "2025-05-13T11:30:27.855893Z", "shell.execute_reply": "2025-05-13T11:30:27.855442Z" }, "vscode": { "languageId": "plaintext" } }, "outputs": [ { "data": { "text/html": [ "\n", "
\n", "
\n", " \n", "
\n", " \n", "
\n", "
\n", "
\n", "\n", "
\n", "\n", "
\n", "
\n", "
\n", " \n", "
\n", " 66k\n", " \n", "
\n", "\n", "
\n", " 49k\n", " \n", "
\n", "\n", "
\n", " 33k\n", " \n", "
\n", "\n", "
\n", " 16k\n", " \n", "
\n", "\n", "
\n", " 0\n", " \n", "
\n", "
\n", "
\n", " MIP\n", "
\n", "
\n", " " ], "text/plain": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mip = libcarna.mip(GEOMETRY_TYPE_VOLUME, sr=500)\n", "mip.cmap('jet', ramp=(0.5, 0.7))\n", "\n", "libcarna.imshow(\n", " libcarna.animate(\n", " libcarna.animate.swing_local(camera, amplitude=22),\n", " n_frames=50,\n", " ).render(\n", " libcarna.renderer(600, 450, [dvr.replicate(), mip]),\n", " camera,\n", " ),\n", " dvr.cmap.bar(volume, label='DVR', tick_labels=False),\n", " mip.cmap.bar(volume, label='MIP'),\n", ")" ] }, { "cell_type": "markdown", "id": "8d5f0cf4", "metadata": {}, "source": [ "Note that we use `dvr.replicate()` when adding the previously defined DVR to the renderer. This is because each\n", "rendering stage can only be added to one renderer, hence, we replicate it this time. Of course, we could have used the\n", "`dvr.replicate()` method the first time that we added the DVR to a renderer, too, but this is not mandatory. All\n", "rendering stages provide such a method.\n", "\n", "## Pointwise Annotations\n", "\n", "Visual validation of detections, for example, requires visualization of those detections within the spatial context of\n", "the original image data. Since LibCarna permits combining different renderers with great flexibility, there is nothing\n", "in the way of rendering some opaque markers on top of the DVR.\n", "\n", "First, lets detect the chromatin spots in the 3D image:" ] }, { "cell_type": "code", "execution_count": 8, "id": "9931ba22", "metadata": { "execution": { "iopub.execute_input": "2025-05-13T11:30:27.859332Z", "iopub.status.busy": "2025-05-13T11:30:27.858992Z", "iopub.status.idle": "2025-05-13T11:30:28.082316Z", "shell.execute_reply": "2025-05-13T11:30:28.081718Z" }, "vscode": { "languageId": "plaintext" } }, "outputs": [], "source": [ "data_denoised = ndi.gaussian_filter(data, 1)\n", "data_max = ndi.maximum_filter(data_denoised, size=5)\n", "detections = np.where(\n", " np.logical_and(\n", " data_denoised == data_max,\n", " data_denoised >= 30_000,\n", " )\n", ")" ] }, { "cell_type": "markdown", "id": "f61c72df", "metadata": {}, "source": [ "We now visualize the detected chromatin spots by marking each detection with a red ball. By using `parent=volume`, we\n", "attach those markers to the `volume` node. To correctly position the markers, we take advantage of the\n", "`transform_from_voxels_into` method of the `volume`, that maps the voxel coordinates of the original data to the\n", "coordinates of the `volume` node:" ] }, { "cell_type": "code", "execution_count": 9, "id": "537cfdd2", "metadata": { "execution": { "iopub.execute_input": "2025-05-13T11:30:28.084053Z", "iopub.status.busy": "2025-05-13T11:30:28.083694Z", "iopub.status.idle": "2025-05-13T11:30:28.472777Z", "shell.execute_reply": "2025-05-13T11:30:28.472296Z" }, "vscode": { "languageId": "plaintext" } }, "outputs": [ { "data": { "text/html": [ "\n", "
\n", "
\n", " \n", "
\n", " \n", "
\n", " " ], "text/plain": [ "" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "GEOMETRY_TYPE_OPAQUE = 2\n", "\n", "ball = libcarna.meshes.create_ball(5)\n", "red = libcarna.material('solid', color=libcarna.color.RED)\n", "\n", "for xyz in zip(*detections):\n", " libcarna.geometry(\n", " GEOMETRY_TYPE_OPAQUE,\n", " parent=volume,\n", " mesh=ball,\n", " material=red,\n", " ).translate(\n", " *volume.transform_from_voxels_into(volume).point(xyz)\n", " )\n", "\n", "libcarna.imshow(\n", " libcarna.animate(\n", " libcarna.animate.swing_local(camera, amplitude=22),\n", " n_frames=50,\n", " ).render(\n", " libcarna.renderer(600, 450, [\n", " dvr.replicate(),\n", " libcarna.opaque_renderer(GEOMETRY_TYPE_OPAQUE),\n", " ]),\n", " camera,\n", " ),\n", ")" ] }, { "cell_type": "markdown", "id": "bd7785f1", "metadata": {}, "source": [ "## Mask Visualization\n", "\n", "Segmentation is the identification of image regions that correspond to individual objects. To first obtain a\n", "segmentation of the image, we apply an intensity threshold, followed by a connected component analysis to identify the\n", "individual nuclei:" ] }, { "cell_type": "code", "execution_count": 10, "id": "417a2589", "metadata": { "execution": { "iopub.execute_input": "2025-05-13T11:30:28.475480Z", "iopub.status.busy": "2025-05-13T11:30:28.475278Z", "iopub.status.idle": "2025-05-13T11:30:28.500487Z", "shell.execute_reply": "2025-05-13T11:30:28.499892Z" }, "vscode": { "languageId": "plaintext" } }, "outputs": [], "source": [ "data_seg = ndi.label(data_denoised > 10_000)[0]" ] }, { "cell_type": "markdown", "id": "df9688f5", "metadata": {}, "source": [ "The `data_seg` array is a *labeled mask*, where 0 corresponds to the image background, and each spatially connected\n", "component of nuclei has a unique *label* (intensity value).\n", "\n", "We can now use the *Mask Renderer* to visualize the segmentation results by rendering their outlines on top of the DVR:" ] }, { "cell_type": "code", "execution_count": 11, "id": "b8af663f", "metadata": { "execution": { "iopub.execute_input": "2025-05-13T11:30:28.501942Z", "iopub.status.busy": "2025-05-13T11:30:28.501831Z", "iopub.status.idle": "2025-05-13T11:30:29.434996Z", "shell.execute_reply": "2025-05-13T11:30:29.434505Z" }, "vscode": { "languageId": "plaintext" } }, "outputs": [ { "data": { "text/html": [ "\n", "
\n", "
\n", " \n", "
\n", " \n", "
\n", " " ], "text/plain": [ "" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "GEOMETRY_TYPE_MASK = 3\n", "\n", "libcarna.volume(\n", " GEOMETRY_TYPE_MASK,\n", " data_seg,\n", " parent=volume,\n", " spacing=volume.spacing,\n", ")\n", "\n", "libcarna.imshow(\n", " libcarna.animate(\n", " libcarna.animate.swing_local(camera, amplitude=22),\n", " n_frames=50,\n", " ).render(\n", " libcarna.renderer(600, 450, [\n", " dvr.replicate(),\n", " libcarna.mask_renderer(GEOMETRY_TYPE_MASK, sr=800),\n", " ]),\n", " camera,\n", " ),\n", ")" ] }, { "cell_type": "markdown", "id": "963a2341", "metadata": {}, "source": [ "The 3D visualization clearly shows that most nuclei are properly segmented, but there are several occasions of falsely\n", "detected (e.g., bottom left) as well as falsely merged nuclei (e.g., bottom right).\n", "\n", "## Track Visualization\n", "\n", "Use *line strips* to visualize tracks. First, we simulate a track of Brownian motion:" ] }, { "cell_type": "code", "execution_count": 12, "id": "1e2cbf14", "metadata": { "execution": { "iopub.execute_input": "2025-05-13T11:30:29.442424Z", "iopub.status.busy": "2025-05-13T11:30:29.442310Z", "iopub.status.idle": "2025-05-13T11:30:29.444747Z", "shell.execute_reply": "2025-05-13T11:30:29.444365Z" }, "vscode": { "languageId": "plaintext" } }, "outputs": [], "source": [ "np.random.seed(2)\n", "track = np.cumsum(np.random.randn(10, 3), axis=0)" ] }, { "cell_type": "markdown", "id": "b6ac2c7c", "metadata": {}, "source": [ "Then, visualize the track as a red line strip:" ] }, { "cell_type": "code", "execution_count": 13, "id": "dab0513a", "metadata": { "execution": { "iopub.execute_input": "2025-05-13T11:30:29.446255Z", "iopub.status.busy": "2025-05-13T11:30:29.446056Z", "iopub.status.idle": "2025-05-13T11:30:30.235444Z", "shell.execute_reply": "2025-05-13T11:30:30.234995Z" }, "vscode": { "languageId": "plaintext" } }, "outputs": [ { "data": { "text/html": [ "\n", "
\n", "
\n", " \n", "
\n", " \n", "
\n", " " ], "text/plain": [ "" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "GEOMETRY_TYPE_TRACK = 4\n", "\n", "libcarna.geometry(\n", " GEOMETRY_TYPE_TRACK,\n", " parent=volume,\n", " mesh=libcarna.meshes.create_line_strip(track),\n", " material=libcarna.material(\n", " 'unshaded', color=libcarna.color.RED, lw=4,\n", " ),\n", ").translate(y=2).scale(5)\n", "\n", "dvr.sample_rate = 800\n", "camera.translate(y=-10, z=-60)\n", "\n", "libcarna.imshow(\n", " libcarna.animate(\n", " libcarna.animate.rotate_local(camera),\n", " n_frames=100,\n", " ).render(\n", " libcarna.renderer(600, 450, [\n", " dvr.replicate(),\n", " libcarna.opaque_renderer(GEOMETRY_TYPE_TRACK),\n", " ]),\n", " camera,\n", " ),\n", ")" ] } ], "metadata": { "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.10" } }, "nbformat": 4, "nbformat_minor": 5 }