Building Custom Phantoms

This guide shows how to create new phantom types from scratch using the geometric primitives provided by GeometricMedicalPhantoms. We'll walk through the key concepts and build a complete example.

Overview: From Primitives to Phantoms

Every phantom in GeometricMedicalPhantoms follows the same pattern:

Define geometry parameters: Create a struct to parameters of your phantom geometry (e.g., centers, radii, etc.) – Optional
Define intensities: Create a struct specifying how to assign intensities to tissues
Define drawing function: Write a draw_*!(ctx, geometry, intensities) function that calls draw_shape!(ctx, shape) for each primitive — shapes are drawn directly onto the context rather than collected in an intermediate array
Render: Construct a DrawContext3D or DrawContext2D and call your drawing function, then return the image

Intensity Assignment: Additive vs. Masking

When primitives overlap, their intensities combine according to one of two modes:

Additive Mode

Intensities add together. If two primitives overlap at a point, the result is the sum of both intensities (or maximum, depending on implementation):

# Pseudo-code
for point in phantom
    intensity = 0
    for primitive in primitives
        if point_in_primitive(point, primitive)
            intensity += primitive.intensity
        end
    end
end

Use when: Creating phantoms where structures have distinct tissue types that don't overlap (e.g., Shepp-Logan brain).

Masking Mode

Later primitives overwrite earlier ones. Order matters:

# Pseudo-code
for point in phantom
    intensity = 0
    for primitive in primitives
        if point_in_primitive(point, primitive)
            intensity = primitive.intensity  # Overwrite
        end
    end
end

Use when: Creating composites where you want inner structures to override outer ones (e.g., heart chambers inside the torso, blood vessels overriding tissue).

Combining Primitives

Complex phantoms are created by combining multiple primitives. The order and intensity assignment mode (Additive vs. Masking) determine the final result.

# Combine ellipsoid and cylinder
# Note: Constructors use positional arguments (cx, cy, cz, rx, ry, rz, intensity)
combined_shapes = [
    Ellipsoid(0.0, 0.0, 0.0, 0.7, 0.6, 0.5, 0.5),
    CylinderZ(0.0, 0.0, 0.0, 0.2, 1.5, 0.9)
]

# When using draw!, the order in the array matters for Masking primitives.
# Additive primitives sum up.

Example: Building a Simple Brain Phantom

Let's create a minimal custom phantom showing a brain with ventricles and a tumor. This example follows the pattern used in the package's built-in phantoms.

Step 1: Define Geometry Structure

First, define the geometry parameters and a function that generates the shapes:

# Custom geometry structure for our brain phantom
struct SimpleBrainGeometry
    brain_radius::Float32      # Radius of brain ellipsoid
    ventricle_radius::Float32  # Central ventricles
    tumor_radius::Float32      # Lesion/tumor size
    tumor_position::NTuple{3, Float32}  # (x, y, z) position
end

# Default constructor
SimpleBrainGeometry() = SimpleBrainGeometry(0.8f0, 0.15f0, 0.1f0, (0.3f0, 0.2f0, 0.0f0))

# Function that draws all shapes directly onto a rendering context
function draw_brain_shapes!(ctx, geom::SimpleBrainGeometry, intensities)
    # Brain ellipsoid (outermost) - using MaskingIntensityValue for masking mode
    GeometricMedicalPhantoms.draw_shape!(ctx, Ellipsoid(
        0.0f0, 0.0f0, 0.0f0,  # center
        geom.brain_radius, geom.brain_radius * 1.1f0, geom.brain_radius * 0.9f0,  # radii
        GeometricMedicalPhantoms.MaskingIntensityValue(Float32(intensities.brain))  # wrapped intensity
    ))

    # Ventricles (overwrites brain interior)
    GeometricMedicalPhantoms.draw_shape!(ctx, Ellipsoid(
        0.0f0, 0.0f0, 0.0f0,
        geom.ventricle_radius, geom.ventricle_radius * 1.5f0, geom.ventricle_radius * 2f0,
        GeometricMedicalPhantoms.MaskingIntensityValue(Float32(intensities.ventricles))
    ))

    # Tumor (overwrites everything at its location)
    tx, ty, tz = geom.tumor_position
    GeometricMedicalPhantoms.draw_shape!(ctx, Ellipsoid(
        tx, ty, tz,
        geom.tumor_radius * 1.2f0, geom.tumor_radius * 1.0f0, geom.tumor_radius * 0.8f0,
        GeometricMedicalPhantoms.MaskingIntensityValue(Float32(intensities.tumor))
    ))
    return nothing
end

Step 2: Define Intensity Structure

This struct holds intensity values. We also define helper functions for both numeric and boolean (mask) intensities:

# Numeric intensities
struct SimpleBrainIntensities
    brain::Float32
    ventricles::Float32
    tumor::Float32
end

# Default constructor
SimpleBrainIntensities() = SimpleBrainIntensities(0.6f0, 0.9f0, 1.0f0)

# Helper function to get intensity by name (for numeric)
function get_intensity(intensities::SimpleBrainIntensities, tissue::Symbol)
    if tissue == :brain
        return intensities.brain
    elseif tissue == :ventricles
        return intensities.ventricles
    elseif tissue == :tumor
        return intensities.tumor
    else
        return 0.0f0
    end
end

# Boolean mask structure
struct SimpleBrainMask
    brain::Bool
    ventricles::Bool
    tumor::Bool
end

# Default: all false
SimpleBrainMask(; brain=false, ventricles=false, tumor=false) =
    SimpleBrainMask(brain, ventricles, tumor)

# Helper function for mask (returns 1.0 if enabled, 0.0 otherwise)
function get_intensity(mask::SimpleBrainMask, tissue::Symbol)
    if tissue == :brain && mask.brain
        return 1.0f0
    elseif tissue == :ventricles && mask.ventricles
        return 1.0f0
    elseif tissue == :tumor && mask.tumor
        return 1.0f0
    else
        return 0.0f0
    end
end

Step 3: Create Rendering Function

Now create the main function that ties everything together:

function create_simple_brain_phantom(
    nx::Int, ny::Int, nz::Int;
    fov::NTuple{3, Float32} = (2.0f0, 2.0f0, 2.0f0),
    geometry::SimpleBrainGeometry = SimpleBrainGeometry(),
    intensities = SimpleBrainIntensities(),
    eltype::Type = Float32
)
    # Initialize phantom array
    phantom = zeros(eltype, nx, ny, nz)

    # Create coordinate axes in physical space
    ax_x = collect(range(-fov[1]/2, fov[1]/2, length=nx))
    ax_y = collect(range(-fov[2]/2, fov[2]/2, length=ny))
    ax_z = collect(range(-fov[3]/2, fov[3]/2, length=nz))

    # Create drawing context and draw all shapes onto the phantom
    ctx = GeometricMedicalPhantoms.DrawContext3D(phantom, ax_x, ax_y, ax_z)
    draw_brain_shapes!(ctx, geometry, intensities)

    return phantom
end

# Create the phantom
brain_phantom = create_simple_brain_phantom(128, 128, 128)
println("Brain phantom created with size: $(size(brain_phantom))")

Brain phantom created with size: (128, 128, 128)

Step 4: Visualize Your Custom Result

# Extract center slices for all three planes
axial_slice = brain_phantom[:, :, div(128, 2)]
coronal_slice = brain_phantom[:, div(128, 2), :]
sagittal_slice = brain_phantom[div(128, 2), :, :]

# Create a single row plot with three images
fig = plot(layout=(1, 3), size=(1200, 400))
heatmap!(fig[1], axial_slice', title="Axial", aspect_ratio=:equal, color=:grays, legend=false, axis=false)
heatmap!(fig[2], coronal_slice', title="Coronal", aspect_ratio=:equal, color=:grays, legend=false, axis=false)
heatmap!(fig[3], sagittal_slice', title="Sagittal", aspect_ratio=:equal, color=:grays, legend=false, axis=false)

GKS: cannot open display - headless operation mode active

You can also create masks to extract specific structures:

# Create a mask showing only the tumor
tumor_mask = create_simple_brain_phantom(
    128, 128, 128;
    intensities=SimpleBrainMask(tumor=true)
)

jim(tumor_mask[:, :, div(128, 2)]; title="Tumor Mask")

Key Design Principles

When building custom phantoms, follow these practices:

1. Separate Concerns

Geometry parameters: Define a struct to hold all shape parameters (radii, centers, etc.)
Intensity structure: Holds tissue properties (signal values or boolean masks)
Drawing function: draw_*!(ctx, geometry, intensities) draws shapes directly onto the context — no intermediate shape arrays
Rendering function: Creates axes, allocates the phantom, builds the DrawContext, calls the drawing function, and returns the result

2. Coordinate System

Always work in physical space (cm), not voxel indices:

x_coords = range(-fov[1]/2, fov[1]/2, length=nx)
# Avoids confusion and makes phantoms independent of resolution

3. Rendering Order

In masking mode, render in order from:

Large outer structures
Medium interior structures
Small details/pathology last

This ensures later structures overwrite earlier ones correctly.

4. Documentation

Include docstrings explaining:

Purpose of each parameter
Typical values
Biological meaning

Example:

"""
    create_custom_phantom(nx, ny, nz; fov=(20,20,20), geometry=..., intensities=...)

Create a custom phantom.

# Arguments
- `nx, ny, nz`: Grid dimensions
- `fov`: Field of view in cm (default: (20, 20, 20))
- `geometry`: Geometry structure with shape parameters
- `intensities`: Intensity structure with tissue values

# Returns
- `phantom::Array{Float32, 3}`: Generated phantom array
"""

5. Parameter Validation

Validate inputs to catch errors early:

function create_phantom(nx, ny, nz; ...)
    nx > 0 || throw(ArgumentError("nx must be positive"))
    ny > 0 || throw(ArgumentError("ny must be positive"))
    nz > 0 || throw(ArgumentError("nz must be positive"))
    # ... rest of function
end

Integration into the Package

Once you've defined a custom phantom, you can integrate it with the package by:

Creating a module file in src/
Defining your geometry and intensity structures
Implementing a public create_*_phantom function
Adding tests to test/
Documenting in the guide (this page)

The package structure supports arbitrary custom phantoms while maintaining consistency with built-in phantoms.