Cell Nuclei

In this example, we will use a cellular image from the Allen Cell WTC-11 hiPSC Single-Cell Image Dataset (Viana et al. 2023).

import libcarna
import numpy as np
import scipy.ndimage as ndi

Get the data:

data = libcarna.data.nuclei()
data.shape, data.dtype

((60, 256, 256), dtype('uint16'))

The data is 60 × 256 × 256 pixels (uint16).

Maximum Intensity Projection

In the code below, we use normals=True on the volume node; albeit this doesn’t make any difference for the Maximum Intensity Projection (MIP), it is benefitial for other rendering modes, discussed below.

GEOMETRY_TYPE_VOLUME = 1

# Create and configure frame renderer
mip = libcarna.mip(GEOMETRY_TYPE_VOLUME, cmap='jet', sr=500)
r = libcarna.renderer(600, 450, [mip])

# Create and configure scene
root = libcarna.node()

volume = libcarna.volume(
    GEOMETRY_TYPE_VOLUME,
    data,
    parent=root,
    spacing=(1, 0.5, 0.5),
).rotate('x', -35).rotate('y', 90)

camera = libcarna.camera(
    parent=root,
).frustum(fov=90, z_near=1, z_far=500).translate(z=100)

# Render
libcarna.imshow(r.render(camera), mip.cmap.bar(volume))

66k

49k

33k

16k

0

In the MIP, it can easily be seen that there is one mitotic nucleus in the image.

For an even better visual perception of the 3D data, it is best viewed from different angles, that can be achieved with as a subtle animation:

# Render as animation
libcarna.imshow(
    libcarna.animate(
        libcarna.animate.swing_local(camera, amplitude=22),
        n_frames=50,
    ).render(r, camera),
    mip.cmap.bar(volume),
)

66k

49k

33k

16k

0

This makes the camera swing by 22° to the left and to the right.

Direct Volume Rendering

In a Direct Volume Rendering (DVR), surfaces are rendered by simulation of the absorption of light. This simulation is most realstic, when the spatial orientation of the surfaces can be taken into account, which requires that the normals of the volume have been computed (this is why we used normals=True when we created the volume node).

dvr = libcarna.dvr(GEOMETRY_TYPE_VOLUME, sr=500)

We again use the jet colormap, that, as we have seen in the MIP, employs blueish colors for the nuclear envelope, and reddish colors for the chromatin. In addition, we use a linear ramp function for the colormap, because we want the space between the nuclei to be translucent:

dvr.cmap('jet', ramp=(0.15, 0.25))

libcarna.imshow(
    libcarna.animate(
        libcarna.animate.swing_local(camera, amplitude=22),
        n_frames=50,
    ).render(
        libcarna.renderer(600, 450, [dvr]),
        camera,
    ),
    dvr.cmap.bar(volume),
)

66k

49k

33k

16k

0

The DVR of the nuclei in the image allows for a very natural perception of the 3D scene. On the downside, the mitotic nucleus is harder to identify. This is because the chromatin from the inside of the nucleus (should be reddish due to our colormap) is fully occluded by the nuclear envelope (blueish).

We can overlay a DVR of the nuclear envelope with a MIP for the chromatin:

mip = libcarna.mip(GEOMETRY_TYPE_VOLUME, sr=500)
mip.cmap('jet', ramp=(0.5, 0.7))

libcarna.imshow(
    libcarna.animate(
        libcarna.animate.swing_local(camera, amplitude=22),
        n_frames=50,
    ).render(
        libcarna.renderer(600, 450, [dvr.replicate(), mip]),
        camera,
    ),
    dvr.cmap.bar(volume, label='DVR', tick_labels=False),
    mip.cmap.bar(volume, label='MIP'),
)

DVR

66k

49k

33k

16k

0

MIP

Note that we use dvr.replicate() when adding the previously defined DVR to the renderer. This is because each rendering stage can only be added to one renderer, hence, we replicate it this time. Of course, we could have used the dvr.replicate() method the first time that we added the DVR to a renderer, too, but this is not mandatory. All rendering stages provide such a method.

Pointwise Annotations

Visual validation of detections, for example, requires visualization of those detections within the spatial context of the original image data. Since LibCarna permits combining different renderers with great flexibility, there is nothing in the way of rendering some opaque markers on top of the DVR.

First, lets detect the chromatin spots in the 3D image:

data_denoised = ndi.gaussian_filter(data, 1)
data_max = ndi.maximum_filter(data_denoised, size=5)
detections = np.where(
    np.logical_and(
        data_denoised == data_max,
        data_denoised >= 30_000,
    )
)

We now visualize the detected chromatin spots by marking each detection with a red ball. By using parent=volume, we attach those markers to the volume node. To correctly position the markers, we take advantage of the transform_from_voxels_into method of the volume, that maps the voxel coordinates of the original data to the coordinates of the volume node:

GEOMETRY_TYPE_OPAQUE = 2

ball = libcarna.meshes.create_ball(5)
red = libcarna.material('solid', color=libcarna.color.RED)

for xyz in zip(*detections):
    libcarna.geometry(
        GEOMETRY_TYPE_OPAQUE,
        parent=volume,
        mesh=ball,
        material=red,
    ).translate(
        *volume.transform_from_voxels_into(volume).point(xyz)
    )

libcarna.imshow(
    libcarna.animate(
        libcarna.animate.swing_local(camera, amplitude=22),
        n_frames=50,
    ).render(
        libcarna.renderer(600, 450, [
            dvr.replicate(),
            libcarna.opaque_renderer(GEOMETRY_TYPE_OPAQUE),
        ]),
        camera,
    ),
)

Mask Visualization

Segmentation is the identification of image regions that correspond to individual objects. To first obtain a segmentation of the image, we apply an intensity threshold, followed by a connected component analysis to identify the individual nuclei:

data_seg = ndi.label(data_denoised > 10_000)[0]

The data_seg array is a labeled mask, where 0 corresponds to the image background, and each spatially connected component of nuclei has a unique label (intensity value).

We can now use the Mask Renderer to visualize the segmentation results by rendering their outlines on top of the DVR:

GEOMETRY_TYPE_MASK = 3

libcarna.volume(
    GEOMETRY_TYPE_MASK,
    data_seg,
    parent=volume,
    spacing=volume.spacing,
)

libcarna.imshow(
    libcarna.animate(
        libcarna.animate.swing_local(camera, amplitude=22),
        n_frames=50,
    ).render(
        libcarna.renderer(600, 450, [
            dvr.replicate(),
            libcarna.mask_renderer(GEOMETRY_TYPE_MASK, sr=800),
        ]),
        camera,
    ),
)

The 3D visualization clearly shows that most nuclei are properly segmented, but there are several occasions of falsely detected (e.g., bottom left) as well as falsely merged nuclei (e.g., bottom right).

Track Visualization

Use line strips to visualize tracks. First, we simulate a track of Brownian motion:

np.random.seed(2)
track = np.cumsum(np.random.randn(10, 3), axis=0)

Then, visualize the track as a red line strip:

GEOMETRY_TYPE_TRACK = 4

libcarna.geometry(
    GEOMETRY_TYPE_TRACK,
    parent=volume,
    mesh=libcarna.meshes.create_line_strip(track),
    material=libcarna.material(
        'unshaded', color=libcarna.color.RED, lw=4,
    ),
).translate(y=2).scale(5)

dvr.sample_rate = 800
camera.translate(y=-10, z=-60)

libcarna.imshow(
    libcarna.animate(
        libcarna.animate.rotate_local(camera),
        n_frames=100,
    ).render(
        libcarna.renderer(600, 450, [
            dvr.replicate(),
            libcarna.opaque_renderer(GEOMETRY_TYPE_TRACK),
        ]),
        camera,
    ),
)