Basic Fitting Example

This example demonstrates the most common use case: fitting Gaussian blobs with a fixed PSF width using the GaussianXYNB model.

Running the Example

The complete example is available in examples/basic_fitting.jl:

julia --project=. examples/basic_fitting.jl

Step-by-Step Walkthrough

Setup and Imports

using GaussMLE
using Statistics

Creating a Fitter

# Create fitter with default settings
# - GaussianXYNB(0.13f0): fixed PSF width of 130nm
# - Auto backend selection (GPU if available)
# - 20 Newton-Raphson iterations
fitter = GaussMLEConfig()

# Or with explicit configuration
fitter = GaussMLEConfig(
    psf_model = GaussianXYNB(0.13f0),  # sigma = 130nm in microns
    backend = :cpu,                      # Force CPU
    iterations = 20
)

Preparing Data

GaussMLE expects data as a 3D array with dimensions (roi_size, roi_size, n_rois):

# Example: 100 ROIs of 11x11 pixels each
data = rand(Float32, 11, 11, 100)

For real data, each ROI should contain a single fluorescent emitter centered approximately in the ROI.

Fitting the Data

# Perform the fitting
smld, info = fit(data, fitter)

println("Fitted $(length(smld.emitters)) localizations")

Accessing Results

The fit() function returns a tuple (smld, info) where smld is a SMLMData.BasicSMLD containing emitter objects:

# Extract position arrays
x_positions = [e.x for e in smld.emitters]
y_positions = [e.y for e in smld.emitters]

# Extract photometry
photons = [e.photons for e in smld.emitters]
backgrounds = [e.bg for e in smld.emitters]

# Extract uncertainties (CRLB)
precisions_x = [e.σ_x for e in smld.emitters]
precisions_y = [e.σ_y for e in smld.emitters]

# Display statistics
println("Mean position: ($(round(mean(x_positions), digits=2)), $(round(mean(y_positions), digits=2))) microns")
println("Mean photons: $(round(mean(photons), digits=1))")
println("Mean precision: $(round(mean(precisions_x)*1000, digits=1)) nm")

Complete Working Example

using GaussMLE
using Statistics

# ROIBatch typically comes from SMLMBoxer.jl (extracts ROIs from raw movie frames)
# For testing, use generate_roi_batch() or raw arrays:
println("Generating synthetic data...")
data = rand(Float32, 11, 11, 100)

# Create fitter with PSF model (sigma from PSF calibration)
println("Creating fitter with 130nm PSF width...")
fitter = GaussMLEConfig(psf_model = GaussianXYNB(0.13f0))

# Fit
println("Fitting $(size(data, 3)) ROIs...")
smld, info = fit(data, fitter)

# Display results
println("\n=== Results ===")
println("Type: BasicSMLD")
println("Fitted: $(length(smld.emitters)) localizations")

# Extract statistics
x_positions = [e.x for e in smld.emitters]
y_positions = [e.y for e in smld.emitters]
photons = [e.photons for e in smld.emitters]
backgrounds = [e.bg for e in smld.emitters]
precisions_x = [e.σ_x for e in smld.emitters]

println("Mean position: ($(round(mean(x_positions), digits=2)), $(round(mean(y_positions), digits=2))) microns")
println("Mean photons: $(round(mean(photons), digits=1))")
println("Mean background: $(round(mean(backgrounds), digits=1))")
println("Mean precision: $(round(mean(precisions_x)*1000, digits=1)) nm")

Understanding the Results

Emitter2DFitGaussMLE Fields

Each emitter in smld.emitters is an Emitter2DFitGaussMLE containing:

Quality Filtering with @filter

Use SMLMData's @filter macro for quality control:

using GaussMLE

smld, info = fit(data, fitter)

# Filter by precision and photon count
good = @filter(smld, σ_x < 0.020 && photons > 500)

# Filter by multiple criteria
precise = @filter(smld, σ_x < 0.015 && σ_y < 0.015 && bg < 50)

println("Kept $(length(good.emitters)) / $(length(smld.emitters)) localizations")

Using ROIBatch for Real Data

In a typical SMLM pipeline, ROIBatch comes from SMLMBoxer.jl which detects candidates and extracts ROIs from raw movie frames:

Raw Movie → SMLMBoxer.jl → ROIBatch → GaussMLE.fit() → BasicSMLD → Analysis

For testing with proper coordinate handling:

using GaussMLE

# Create camera model (100nm pixels)
camera = IdealCamera(0:1023, 0:1023, 0.1)

# Generate synthetic data with proper camera model
batch = generate_roi_batch(
    camera,
    GaussianXYNB(0.13f0),  # PSF model
    n_rois = 100,
    roi_size = 11
)

# Fit with proper coordinate conversion
fitter = GaussMLEConfig(psf_model = GaussianXYNB(0.13f0))
smld, info = fit(batch, fitter)

# Positions are now in camera coordinates (microns)
x_positions = [e.x for e in smld.emitters]
println("X range: $(extrema(x_positions)) microns")

Performance Optimization

Use Float32

# Float32 is faster and sufficient precision
data = Float32.(your_data)

Batch Processing for Large Datasets

# For GPU: configure batch size based on memory
fitter = GaussMLEConfig(
    backend = :gpu,
    batch_size = 10_000
)

# Fit large dataset
large_data = rand(Float32, 11, 11, 100_000)
@time smld, info = fit(large_data, fitter)

Timing Comparison

using GaussMLE

data = rand(Float32, 11, 11, 10_000)

# CPU timing
fitter_cpu = GaussMLEConfig(backend = :cpu)
t_cpu = @elapsed fit(data, fitter_cpu)
println("CPU: $(round(10_000/t_cpu)) fits/second")

# GPU timing (if available)
fitter_gpu = GaussMLEConfig(backend = :gpu)
t_gpu = @elapsed fit(data, fitter_gpu)
println("GPU: $(round(10_000/t_gpu)) fits/second")

Practical Considerations

Choosing ROI Size

The ROI should be large enough to capture the full PSF:

# Rule of thumb: ROI size >= 4 * PSF_width_in_pixels
# For 130nm PSF with 100nm pixels: PSF = 1.3 pixels
# Minimum ROI: 6 pixels (use 7 or 11 for safety)

Handling Edge Cases

# Filter out failed fits using @filter
good = @filter(smld, photons > 0)
println("Valid fits: $(length(good.emitters)) / $(length(smld.emitters))")

Next Steps

• Try the PSF Width Fitting example for variable PSF scenarios
• Learn about GPU Support for larger datasets
• See the API Reference for all options