Visualizing Volumetric Data and Reconstructions

3D Scalar Fields and Reconstructions

Medical imaging (CT, MRI), computational electromagnetics (FDTD), and radar/sonar produce 3D scalar fields that must be visualized as volumes. Isosurface extraction, volume rendering, and slicing are the key techniques.

Definition:
Isosurface Extraction

An isosurface is the 3D equivalent of a contour line — the surface where a scalar field equals a constant value:

$\{(x,y,z) : f(x,y,z) = c\}$

The Marching Cubes algorithm extracts triangle meshes from voxel data:

from skimage.measure import marching_cubes

verts, faces, normals, values = marching_cubes(volume, level=0.5)

Definition:
Voxel (Volume Element)

A voxel is the 3D analogue of a pixel. A volumetric dataset is a 3D array where each element $V[i,j,k]$ represents a scalar value at grid position $(x_i, y_j, z_k)$ .

# Create a 3D dataset
volume = np.zeros((64, 64, 64))
x, y, z = np.mgrid[-3:3:64j, -3:3:64j, -3:3:64j]
volume = np.exp(-(x**2 + y**2 + z**2) / 2)  # 3D Gaussian

Theorem: Marching Cubes Complexity

For a volumetric grid of $N^3$ voxels, the Marching Cubes algorithm:

Processes each cube (8 neighboring voxels) independently
Has 256 possible cube configurations (stored in a lookup table)
Runs in $O(N^3)$ time and produces $O(N^2)$ triangles

The number of output triangles depends on the isosurface complexity, not the volume size.

Marching Cubes visits every cube once and looks up the triangle pattern in a table. It is embarrassingly parallel and runs fast even on large volumes.

Example: Orthogonal Volume Slicing

Display three orthogonal slices (axial, sagittal, coronal) through a 3D EM field simulation volume.

Solution

Implementation

import numpy as np
import matplotlib.pyplot as plt

# Simulate 3D field
x, y, z = np.mgrid[-3:3:64j, -3:3:64j, -3:3:64j]
E_field = np.sinc(np.sqrt(x**2 + y**2 + z**2)) * np.exp(-(x**2+y**2+z**2)/4)

fig, axes = plt.subplots(1, 3, figsize=(14, 4))

# Axial slice (z=0)
axes[0].imshow(E_field[:, :, 32], cmap='RdBu_r', origin='lower')
axes[0].set_title('Axial (z=0)')

# Sagittal slice (x=0)
axes[1].imshow(E_field[32, :, :].T, cmap='RdBu_r', origin='lower')
axes[1].set_title('Sagittal (x=0)')

# Coronal slice (y=0)
axes[2].imshow(E_field[:, 32, :].T, cmap='RdBu_r', origin='lower')
axes[2].set_title('Coronal (y=0)')

for ax in axes:
    ax.set_aspect('equal')
fig.suptitle('EM Field — Orthogonal Slices', fontweight='bold')

3D Volume Slicer

Adjust the slice position through a 3D volume to explore the internal structure of a scalar field.

Parameters

Common Mistake: 3D Arrays Consume Enormous Memory

Mistake:

Creating a (1000, 1000, 1000) float64 array — that is 8 GB of RAM.

Correction:

Use coarser grids for visualization ( $64^3$ to $256^3$ is usually sufficient). Use float32 instead of float64 to halve memory. For production volumes, use chunked I/O (HDF5, Zarr).

Isosurface

A 3D surface where a scalar field equals a specified constant value, analogous to contour lines in 2D.

Voxel

A volume element — the smallest unit in a 3D grid, analogous to a pixel in 2D.

Visualizing Antenna Arrays and Beam Patterns Mesh Processing with Open3D and Trimesh