Code Examples

This page provides practical code examples for common MedImages.jl operations.

Quick Reference Examples

Loading and Saving

using MedImages

# Load a CT scan
ct = load_image("/path/to/ct.nii.gz", "CT")

# Load a PET scan
pet = load_image("/path/to/pet.nii.gz", "PET")

# Save a MedImage as NIfTI
create_nii_from_medimage(ct, "/output/processed")  # Creates processed.nii.gz

Inspecting Image Properties

# Access voxel data
data = image.voxel_data
println("Dimensions: ", size(data))
println("Data type: ", eltype(data))

# Access spatial metadata
println("Origin (x,y,z): ", image.origin)
println("Spacing (x,y,z): ", image.spacing)
println("Direction: ", image.direction)

# Access image type
println("Modality: ", image.image_type)
println("Subtype: ", image.image_subtype)

Basic Transformations

# Crop: Extract region starting at (10,20,30) with size (100,100,50)
cropped = crop_mi(image, (10, 20, 30), (100, 100, 50), Linear_en)

# Pad: Add 5 voxels on each side with value -1000
padded = pad_mi(image, (5, 5, 5), (5, 5, 5), -1000.0, Linear_en)

# Rotate: 30 degrees around axis 2
rotated = rotate_mi(image, 2, 30.0, Linear_en)

# Translate: Move origin by 10 voxels along axis 1
translated = translate_mi(image, 10, 1, Linear_en)

# Scale: Reduce size by half
scaled = scale_mi(image, 0.5, Linear_en)

Resampling

# Resample to new spacing
isotropic = resample_to_spacing(image, (1.0, 1.0, 1.0), Linear_en)

# Resample moving image to match fixed image geometry
registered = resample_to_image(fixed_image, moving_image, Linear_en)

Orientation Changes

# Change to RAS orientation
ras = change_orientation(image, ORIENTATION_RAS)

# Change to LPS orientation
lps = change_orientation(image, ORIENTATION_LPS)

HDF5 Storage

using HDF5

# Save to HDF5
h5open("data.h5", "w") do f
    uuid = save_med_image(f, "images", image)
    # Store uuid for later retrieval
end

# Load from HDF5
h5open("data.h5", "r") do f
    loaded = load_med_image(f, "images", "stored-uuid-here")
end

---

Complete Workflow Examples

Example 1: CT Preprocessing

using MedImages

function preprocess_ct_for_analysis(input_path::String)
    # Load the image
    ct = load_image(input_path, "CT")

    # Convert to LPS orientation (standard for CT)
    ct = change_orientation(ct, ORIENTATION_LPS)

    # Resample to 1mm isotropic spacing
    ct = resample_to_spacing(ct, (1.0, 1.0, 1.0), Linear_en)

    # Window to soft tissue range (-100 to 300 HU)
    windowed_data = clamp.(ct.voxel_data, -100.0, 300.0)
    ct = update_voxel_data(ct, windowed_data)

    return ct
end

# Usage
preprocessed = preprocess_ct_for_analysis("/data/raw_ct.nii.gz")
create_nii_from_medimage(preprocessed, "/data/preprocessed_ct")

Example 2: PET/CT Fusion Preparation

using MedImages

function prepare_petct_fusion(ct_path::String, pet_path::String)
    # Load both modalities
    ct = load_image(ct_path, "CT")
    pet = load_image(pet_path, "PET")

    # Register PET to CT coordinate space
    pet_registered = resample_to_image(ct, pet, Linear_en, 0.0)

    # Normalize PET SUV values to 0-1 range for fusion
    pet_data = pet_registered.voxel_data
    pet_normalized = (pet_data .- minimum(pet_data)) ./ (maximum(pet_data) - minimum(pet_data))
    pet_registered = update_voxel_data(pet_registered, pet_normalized)

    return ct, pet_registered
end

ct, pet = prepare_petct_fusion("/data/ct.nii.gz", "/data/pet.nii.gz")

Example 3: Batch Processing

using MedImages

function batch_resample_to_isotropic(input_dir::String, output_dir::String, spacing::Float64)
    # Find all NIfTI files
    nifti_files = filter(f -> endswith(f, ".nii.gz") || endswith(f, ".nii"),
                         readdir(input_dir))

    for filename in nifti_files
        input_path = joinpath(input_dir, filename)
        output_path = joinpath(output_dir, replace(filename, ".nii" => "_iso.nii"))

        println("Processing: ", filename)

        # Load image
        image = load_image(input_path, "CT")

        # Resample
        resampled = resample_to_spacing(image, (spacing, spacing, spacing), Linear_en)

        # Save
        create_nii_from_medimage(resampled, output_path)
    end

    println("Done processing ", length(nifti_files), " files")
end

# Usage
batch_resample_to_isotropic("/data/input", "/data/output", 1.0)

Example 4: ROI Extraction

using MedImages

function extract_roi_around_point(image::MedImage,
                                  center_voxel::Tuple{Int,Int,Int},
                                  roi_size::Tuple{Int,Int,Int})
    # Calculate crop start position (ensure non-negative)
    half_size = roi_size ./ 2 .|> floor .|> Int

    start_pos = max.(center_voxel .- half_size, 0)

    # Adjust size if ROI extends beyond image boundaries
    image_size = size(image.voxel_data)
    actual_size = min.(roi_size, image_size .- start_pos)

    # Extract ROI
    roi = crop_mi(image, start_pos, actual_size, Linear_en)

    return roi
end

# Usage
image = load_image("/data/scan.nii.gz", "CT")
roi = extract_roi_around_point(image, (100, 200, 50), (64, 64, 64))

Example 5: Multi-Resolution Pyramid

using MedImages

function create_image_pyramid(image::MedImage, levels::Int)
    pyramid = MedImage[]
    push!(pyramid, image)

    current = image
    for level in 2:levels
        # Downscale by factor of 2
        downscaled = scale_mi(current, 0.5, Linear_en)
        push!(pyramid, downscaled)
        current = downscaled

        println("Level $level: ", size(downscaled.voxel_data))
    end

    return pyramid
end

# Usage
image = load_image("/data/highres.nii.gz", "CT")
pyramid = create_image_pyramid(image, 4)

# Access different resolution levels
full_res = pyramid[1]
half_res = pyramid[2]
quarter_res = pyramid[3]

---

Working with Metadata

Accessing Clinical Data

# Set clinical data when creating/modifying images
image.clinical_data["age"] = 65
image.clinical_data["sex"] = "M"
image.clinical_data["diagnosis"] = "Normal"

# Access clinical data
println("Patient age: ", get(image.clinical_data, "age", "Unknown"))

Custom Metadata

# Store custom processing metadata
image.metadata["preprocessing_version"] = "1.0"
image.metadata["windowing_applied"] = true
image.metadata["window_center"] = 40
image.metadata["window_width"] = 400

# Save with metadata preserved
h5open("data.h5", "w") do f
    save_med_image(f, "processed", image)
end

---

Interpolation Methods Comparison

using MedImages

function compare_interpolators(image::MedImage, scale_factor::Float64)
    # Nearest neighbor - fastest, preserves values
    nn = scale_mi(image, scale_factor, Nearest_neighbour_en)

    # Linear - balanced quality and speed
    linear = scale_mi(image, scale_factor, Linear_en)

    # B-spline - smoothest results
    bspline = scale_mi(image, scale_factor, B_spline_en)

    return nn, linear, bspline
end

# Usage
image = load_image("/data/image.nii.gz", "CT")
nn, linear, bspline = compare_interpolators(image, 2.0)

# Save for visual comparison
create_nii_from_medimage(nn, "/output/scaled_nn")
create_nii_from_medimage(linear, "/output/scaled_linear")
create_nii_from_medimage(bspline, "/output/scaled_bspline")

When to use each interpolator:

• Nearest_neighbour_en: Segmentation masks, label maps

• Linear_en: General purpose, CT/MRI intensity images

• B_spline_en: When smooth gradients are critical

---

Error Handling Patterns

using MedImages

function safe_load_image(path::String, modality::String)
    if !isfile(path)
        error("File not found: $path")
    end

    try
        image = load_image(path, modality)
        return image
    catch e
        @warn "Failed to load image: $e"
        return nothing
    end
end

function safe_resample(image::MedImage, target_spacing::Tuple{Float64,Float64,Float64})
    # Validate spacing values
    if any(s <= 0 for s in target_spacing)
        error("Spacing values must be positive")
    end

    # Check for extremely small or large resampling factors
    original_spacing = image.spacing
    factors = original_spacing ./ target_spacing

    if any(f > 10 for f in factors)
        @warn "Large upsampling factor detected. This may produce artifacts."
    end

    return resample_to_spacing(image, target_spacing, Linear_en)
end

---

Performance Tips

Memory-Efficient Processing

# Use views instead of copies when possible
cropped_view = @view image.voxel_data[10:100, 20:150, 30:80]

# Process large datasets in chunks
function process_in_slices(image::MedImage, process_func::Function)
    data = image.voxel_data
    result = similar(data)

    for z in axes(data, 1)
        result[z, :, :] = process_func(data[z, :, :])
    end

    return update_voxel_data(image, result)
end

Type Stability

# Ensure consistent types for better performance
function normalize_intensity(image::MedImage)
    data = Float64.(image.voxel_data)
    min_val, max_val = extrema(data)

    if min_val != max_val
        normalized = (data .- min_val) ./ (max_val - min_val)
    else
        normalized = zeros(size(data))
    end

    return update_voxel_data(image, normalized)
end