PL / Raman Mapping

PLMap <: AbstractSpectroscopyData

Photoluminescence intensity map from a CCD raster scan.

A 2D spatial grid where each point has a full CCD spectrum. The intensity field holds the integrated PL signal at each position; the full spectra are stored in spectra for extraction at individual positions.

Fields

• intensity::Matrix{Float64} — Integrated PL intensity (nx, ny)
• spectra::Array{Float64,3} — Raw CCD counts, channel-major (n_pixel, nx, ny) (one pixel's spectrum is the contiguous first axis — see pixel_view)
• x::Vector{Float64} — Spatial x positions (μm)
• y::Vector{Float64} — Spatial y positions (μm)
• pixel::Vector{Float64} — Pixel indices (or wavelength if calibrated)
• metadata::Dict{Symbol,Any} — Symbol-keyed metadata token dict (shared with the rest of the AbstractSpectroscopyData family): provenance (:source_file, :nx, :ny, :step_size, :pixel_range) plus axis tokens (:position_unit for the spatial axes, :yquantity/:yunit for the intensity).

extract_spectrum(m::PLMap, ix::Int, iy::Int) -> NamedTuple

Extract the CCD spectrum at grid index (ix, iy).

Returns (pixel=..., signal=..., x=..., y=...).

extract_spectrum(m::PLMap; x::Real, y::Real) -> NamedTuple

Extract the CCD spectrum at the spatial position nearest to (x, y).

Returns (pixel=..., signal=..., x=..., y=..., ix=..., iy=...).

subtract_background(matrix::TimeResolvedMatrix; t_range=nothing) -> TimeResolvedMatrix

Subtract pre-pump background from a TA matrix by averaging and removing the signal in the baseline region (before pump arrival).

NaN samples inside the baseline window are skipped: each wavelength column's baseline is the mean of its finite values only. This keeps the correction usable after correct_chirp, which NaN-fills samples shifted past the time-axis edges — for negative shifts that lands at the start of the axis, inside typical pre-pump windows. A column with no finite values in the window gets a NaN baseline (the whole column becomes NaN) and a warning is emitted.

subtract_background(m::PLMap; positions=nothing, margin=5) -> PLMap

Subtract a background spectrum from every grid point.

The background is the average CCD spectrum over the reference positions. After subtraction, the intensity map is recomputed from the corrected spectra using the same pixel_range as the original load (if any).

Arguments

• positions: Vector of (x, y) spatial coordinate tuples (μm) for background reference points. These should be off-flake positions with no PL signal.
• margin: Number of grid points from each edge used for auto-detection when positions is not given. Auto mode averages the corners of the bottom half of the map (avoids top-row artifacts). Default: 5. Must be ≥ 1; on maps smaller than the margin it is clamped so every corner point is counted exactly once (no bounds errors, no double-weighted overlap).

Example

# spectra: (n_pixel, nx, ny) CCD counts; x, y: stage positions (μm)
intensity = dropdims(sum(spectra; dims=1); dims=1)
m = PLMap(intensity, spectra, x, y, collect(1.0:size(spectra, 1)))

# Explicit background positions
m_bg = subtract_background(m; positions=[(-40, -40), (40, -40), (-40, -20)])

# Auto mode (bottom corners)
m_bg = subtract_background(m)

integrated_intensity(m::PLMap; pixel_range=nothing) -> Matrix{Float64}

Compute the integrated intensity at each grid point.

Without pixel_range, returns m.intensity (the precomputed full-range sum). With pixel_range, sums m.spectra[p1:p2, :, :] over the given pixel window.

Arguments

• pixel_range: (start, stop) pixel indices to integrate over. Falls back to the pixel_range stored in metadata, or uses m.intensity if unset.

intensity_mask(m::PLMap; pixel_range=nothing, threshold=0.05, exclude=nothing) -> NamedTuple

Compute a boolean mask over the PLMap grid based on integrated intensity.

Grid points where the integrated intensity is below threshold fraction of the intensity range are excluded (false). Spatial regions listed in exclude are also set to false regardless of intensity. This is useful for filtering out off-sample regions and known artifacts (e.g., substrate bands) before batch operations like fitting.

The threshold is computed as: cutoff = min + threshold × (max - min).

Arguments

• pixel_range: (start, stop) pixel indices for computing integrated intensity. If nothing, uses the full precomputed m.intensity.
• threshold: Fraction of the intensity range (0.0 to 1.0). Default 0.05 (5%).
• exclude: Spatial regions to exclude. Each region is a tuple of x and y ranges in spatial coordinates (μm): ((x_min, x_max), (y_min, y_max)). Accepts a single region or a vector of regions. Grid points whose spatial coordinates fall within any excluded region are masked out.

Returns

A named tuple with fields:

• mask::BitMatrix — true for included grid points
• n_included::Int — number of included points
• n_total::Int — total number of grid points
• intensity_min::Float64 — minimum integrated intensity
• intensity_max::Float64 — maximum integrated intensity
• cutoff::Float64 — absolute intensity cutoff value

Examples

m = PLMap(intensity, spectra, x, y, pixel)
m = subtract_background(m)

# Threshold only
result = intensity_mask(m; pixel_range=(950, 1100), threshold=0.1)

# Exclude a substrate band at the top of the map
result = intensity_mask(m; threshold=0.05,
    exclude=((-Inf, Inf), (40.0, Inf)))

# Multiple exclusion regions
result = intensity_mask(m; threshold=0.05,
    exclude=[
        ((-Inf, Inf), (40.0, Inf)),     # top substrate band
        ((-50.0, -40.0), (-50.0, -40.0))  # noisy corner
    ])

peak_centers(m::PLMap; pixel_range=nothing, threshold=0.05) -> Matrix{Float64}

Compute the centroid (intensity-weighted average pixel) at each grid point.

Returns a (nx, ny) matrix of peak center positions in pixel units. Grid points where the PL intensity (m.intensity) is below threshold × the map maximum are set to NaN (renders as transparent in heatmaps with nan_color=:transparent).

Masking uses the PLMap's intensity field — the integrated PL signal already stored in the map. This produces clean masks that match the intensity heatmap: off-flake regions (low PL signal) are transparent, on-flake regions show centroid positions.

Arguments

• pixel_range: (start, stop) pixel range to compute the centroid over. Falls back to the pixel_range stored in metadata, or all pixels if unset.
• threshold: Fraction of the maximum PL intensity below which a grid point is masked as NaN. Default 0.05 (5%). Set to 0 to disable masking.

Example

m = PLMap(intensity, spectra, x, y, pixel)
m = subtract_background(m)
centers = peak_centers(m; pixel_range=(950, 1100))
heatmap(m.x, m.y, centers; colormap=:viridis, nan_color=:transparent)

fit_map(m::PLMap; kwargs...) -> FitMapResult

Fit peaks at every grid point in a PLMap, returning a FitMapResult with per-pixel fit results and summary arrays for heatmap visualization.

Uses reference pixel seeding: the brightest included pixel is fitted first with automatic initial parameter detection, and its converged parameters are used as starting values for all remaining pixels. This dramatically improves convergence for broad or noisy peaks where auto-detection is unreliable, since peak shapes vary smoothly across the spatial map.

Fitting is multi-threaded when Julia is started with multiple threads (julia --threads=auto).

Keyword Arguments

• model::Function = gaussian — Peak model (gaussian, lorentzian, pseudo_voigt)
• n_peaks::Int = 1 — Number of peaks to fit
• region::Union{Tuple{Real,Real}, Nothing} = nothing — Spectral region to fit within
• baseline_order::Int = 1 — Polynomial baseline order
• mask::Union{BitMatrix, Nothing} = nothing — Pre-computed fit mask. If provided, threshold/pixel_range/exclude are ignored.
• threshold::Real = 0.0 — Intensity threshold for intensity_mask (0–1)
• pixel_range::Union{Tuple{Int,Int}, Nothing} = nothing — Pixel range for mask computation
• exclude — Spatial exclusion regions for intensity_mask
• abort::Union{Threads.Atomic{Bool}, Nothing} = nothing — Set to true to stop early
• progress::Union{Function, Nothing} = nothing — Callback (n_done, n_total) -> nothing

Example

m = PLMap(intensity, spectra, x, y, pixel)
m = subtract_background(m)

result = fit_map(m;
    model = lorentzian,
    region = (950, 1050),
    threshold = 0.05)

# Summary arrays for heatmaps
heatmap(m.x, m.y, result.centers)
heatmap(m.x, m.y, result.fwhms)

# Per-pixel result
result[10, 15]  # MultiPeakFitResult or nothing

FitMapResult

Result of batch peak fitting across a spatial map (e.g., PLMap).

Contains per-pixel fit results and pre-extracted summary arrays suitable for heatmap visualization. Pixels that were skipped (masked) or failed to converge have nothing in the results matrix and NaN in summary arrays.

Summary arrays are 3D (nx, ny, n_peaks) for per-peak quantities (centers, fwhms, amplitudes) and 2D (nx, ny) for per-pixel quantities (r_squareds).

Created by fit_map.

Fields

• results::Matrix{Any} — Per-pixel MultiPeakFitResult or nothing
• mask::Union{BitMatrix, Nothing} — Fit mask used (if any)
• centers::Array{Float64, 3} — Peak center per peak (nx, ny, n_peaks)
• fwhms::Array{Float64, 3} — FWHM per peak (nx, ny, n_peaks)
• amplitudes::Array{Float64, 3} — Amplitude per peak (nx, ny, n_peaks)
• r_squareds::Matrix{Float64} — Per-pixel R²
• n_peaks::Int — Number of peaks fitted
• n_converged::Int — Number of successfully fitted pixels
• n_failed::Int — Number of pixels where fitting threw an error
• n_skipped::Int — Number of pixels excluded by the mask
• median_r_squared::Float64 — Median R² across converged pixels

pca_map(m::PLMap; n_components::Int=3, pixel_range=nothing) -> DecompositionResult

Principal Component Analysis of PLMap spectra via truncated SVD.

Each spatial pixel has a full spectrum. PCA finds the orthogonal spectral components that capture the most variance across the map, along with their spatial loading maps.

Arguments

• m::PLMap: Input PL map with channel-major spectra of shape (n_pixel, nx, ny).
• n_components::Int=3: Number of principal components to retain.
• pixel_range: Optional (start, stop) tuple to restrict the spectral range before decomposition.

Returns

A DecompositionResult with:

• loadings: Spatial score maps (nx, ny, n_components) — the weight of each component at each spatial position.
• components: Spectral profiles (n_components, n_spectral) — the principal spectral shapes (eigenvectors of the covariance matrix).
• explained_variance: Fraction of total variance captured by each component.

Example

result = pca_map(m; n_components=5, pixel_range=(900, 1100))
# result.loadings[:, :, 1]  — spatial map of the first principal component
# result.components[1, :]    — spectral shape of the first component
# result.explained_variance  — [0.85, 0.08, 0.03, ...]

nmf_map(m::PLMap; n_components::Int=3, pixel_range=nothing,
        max_iter::Int=200, tol::Float64=1e-4, rng=nothing) -> DecompositionResult

Non-negative Matrix Factorization of PLMap spectra.

Factorizes the non-negative spectra matrix V ≈ W * H where W contains spatial weights and H contains non-negative spectral components. Unlike PCA, NMF enforces non-negativity, producing physically interpretable components that can represent distinct emitters or spectral species.

Uses Lee & Seung multiplicative update rules with Frobenius norm objective.

Arguments

• m::PLMap: Input PL map with channel-major spectra of shape (n_pixel, nx, ny).
• n_components::Int=3: Number of NMF components.
• pixel_range: Optional (start, stop) tuple to restrict the spectral range before decomposition.
• max_iter::Int=200: Maximum number of multiplicative update iterations.
• tol::Float64=1e-4: Convergence tolerance on the relative change in reconstruction error between iterations.
• rng: Optional AbstractRNG for the random W/H initialization. NMF is non-convex, so pass a seeded RNG (e.g. MersenneTwister(42)) for reproducible decompositions; nothing (default) uses the global RNG.

Returns

A DecompositionResult with:

• loadings: Non-negative spatial weight maps (nx, ny, n_components).
• components: Non-negative spectral profiles (n_components, n_spectral).
• explained_variance: Per-component fraction of the total Frobenius norm (‖wₖ hₖ‖² / ‖V‖² for each component k, not a cumulative sum).

Example

result = nmf_map(m; n_components=3, max_iter=500)
# result.loadings[:, :, 2]  — spatial distribution of the second species
# result.components[2, :]    — emission spectrum of the second species

DecompositionResult

Result of PCA or NMF decomposition of a PLMap.

Fields

• loadings::Array{Float64,3} — Spatial maps (nx, ny, n_components)
• components::Matrix{Float64} — Spectral profiles (n_components, n_spectral)
• explained_variance::Vector{Float64} — Per-component fraction of the total: variance (PCA), or Frobenius norm ‖wₖhₖ‖²/‖V‖² (NMF)

normalize_intensity(m::PLMap) -> PLMap

Return a new PLMap with intensity normalized to [0, 1].

intensity(m::PLMap) -> Matrix{Float64}

Return the PL intensity matrix.

pixel_view(m::PLMap, ix::Integer, iy::Integer) -> AbstractVector

View (no copy) of the CCD spectrum at grid index (ix, iy). With the channel-major cube layout this is a unit-stride, allocation-free column.

eachpixel(m::PLMap)

Iterator over per-pixel spectrum views (contiguous columns), in column-major grid order (ix fastest).

spectra_matrix(m::PLMap) -> AbstractMatrix

The cube as a contiguous (n_pixel, nx*ny) matrix view — each column is one pixel's spectrum. Zero-copy.

CosmicRayResult

Result of cosmic ray detection on a 1D signal.

Fields

• indices::Vector{Int} — flagged channel indices
• count::Int — number of flagged channels

CosmicRayMapResult

Result of cosmic ray detection on a PLMap.

Fields

• mask::BitArray{3} — channel-major (n_pixel, nx, ny) like PLMap.spectra, true = cosmic ray
• count::Int — total flagged voxels
• affected_spectra::Int — number of spectra with at least one cosmic ray
• channel_counts::Vector{Int} — cosmic ray count per spectral channel
• msn_index::Matrix{Int} — (nx, ny) Most Similar Neighbor cache from detection: for each pixel, the position of its MSN in the 4-connected neighbor list (0 = no neighbors/unknown). Reused by remove_cosmic_rays to skip recomputing Pearson correlations. The cache reflects the detection pixel_range window.

The convenience constructor CosmicRayMapResult(mask, count, affected_spectra, channel_counts) zero-fills msn_index, so removal recomputes the MSN.

detect_cosmic_rays(signal::AbstractVector; threshold=5.0) -> CosmicRayResult

Detect cosmic ray spikes in a 1D spectrum using modified z-scores on first differences (Whitaker-Hayes method).

Computes first differences Δ[k] = signal[k+1] - signal[k], then flags channels where the modified z-score exceeds threshold. The modified z-score uses the Median Absolute Deviation (MAD) for robust scale estimation.

Arguments

• signal: 1D spectral signal
• threshold: z-score cutoff (default 5.0). Lower values detect more spikes but may flag real features.

Returns

A CosmicRayResult with the flagged channel indices.

Example

result = detect_cosmic_rays(spectrum; threshold=5.0)
println("Found $(result.count) cosmic rays at indices $(result.indices)")

Reference

D. A. Whitaker and K. Hayes, "A simple algorithm for despiking Raman spectra", Chemom. Intell. Lab. Syst. 179, 82 (2018). doi:10.1016/j.chemolab.2018.06.009

detect_cosmic_rays(m::PLMap; threshold=5.0, pixel_range=nothing, max_spike_width=7) -> CosmicRayMapResult

Detect cosmic ray spikes across all spectra in a PLMap using Most Similar Neighbor (MSN) comparison.

For each pixel, selects the 4-connected neighbor with the highest Pearson correlation as the reference spectrum. The residual (pixel minus MSN) is then tested for positive outliers using the MAD-based robust scale estimate. Channels where the residual exceeds threshold × σ above the median residual are flagged.

Using the most-similar neighbor avoids false positives at spatial boundaries: real spectral differences cancel because the MSN is on the same side of the boundary. Only true cosmic ray spikes — absent from all neighbors — remain.

Additional filters:

• Width filter: contiguous flagged runs wider than max_spike_width are unflagged — cosmic rays span 1–5 channels; wider runs are spectral variation.
• Fraction cap: if more than 5% of channels are still flagged after the width filter, all flags for that pixel are cleared.

A global noise floor (median of per-pixel σ estimates) prevents over-sensitivity in spatially homogeneous regions.

Arguments

• m: PLMap with channel-major 3D spectra array (n_pixel, nx, ny)
• threshold: outlier cutoff in MAD-scaled units (default 5.0)
• pixel_range: (start, stop) channel indices to analyze. When set, only this subrange of each spectrum is checked for cosmic rays. Channels outside the range are never flagged. Defaults to the pixel_range in m.metadata, or the full spectrum if unset.
• max_spike_width: maximum width (in channels) of a contiguous flagged run to keep (default 7). Runs wider than this are unflagged as spectral variation.

Returns

A CosmicRayMapResult with a 3D boolean mask and summary statistics.

Example

cr = detect_cosmic_rays(plmap; threshold=5.0)
println("Found $(cr.count) cosmic rays in $(cr.affected_spectra) spectra")

detect_cosmic_rays(m::TimeResolvedMatrix; threshold=5.0, row_fraction_limit=0.25) -> CosmicRayMatrixResult

Detect cosmic ray hits in a 2D time × wavelength image.

Computes the residual of the image over its 3×3 median filter (which smooths isolated spikes while preserving gradual spatial/temporal structure). The residuals are then Poisson-scaled before the noise-floor computation: each residual is divided by sqrt(max(filtered, 1)), so shot noise in bright bands does not inflate the scale estimate. The max(·, 1) floor makes the scaling a no-op on ΔA-scale data (where all values are small), and Poisson-correct for photon-counting data (where filtered ≈ N counts). The noise scale is estimated from the positive scaled residuals only: their median med₊ and MAD MAD₊. A pixel is flagged when

scaled_resid > med₊ + threshold × MAD₊ × 1.4826

Only positive outliers are flagged — cosmic rays deposit charge, so they are always bright.

Arguments

• m: TimeResolvedMatrix with 2D data array (n_time, n_wavelength)
• threshold: outlier cutoff in MAD-scaled units (default 5.0). Lower values detect more spikes but may flag real features.
• row_fraction_limit: rows with more than this fraction of flagged channels are treated as real temporal structure, not cosmic rays (default 0.25). Run detection on unbinned data; on very narrow matrices (few wavelength columns) raise row_fraction_limit, since wide cosmic-ray streaks or hot rows spanning more than that fraction of a row are treated as temporal structure and left unflagged.

Returns

A CosmicRayMatrixResult with the flagged pixel indices.

Limitations

Returns zero detections when positive residuals are absent or have zero spread (e.g. a near-flat background with fewer than ~3 distinct positive residual values) — there is no noise floor to measure against in that regime.

remove_cosmic_rays(signal::AbstractVector, result::CosmicRayResult) -> Vector

Replace flagged cosmic ray channels with linearly interpolated values from the nearest non-flagged neighbors.

Returns a new vector; does not mutate the input.

Example

result = detect_cosmic_rays(spectrum)
cleaned = remove_cosmic_rays(spectrum, result)

remove_cosmic_rays(m::PLMap, result::CosmicRayMapResult; pixel_range=nothing) -> PLMap

Remove cosmic rays from a PLMap using scaled Most Similar Neighbor (MSN) replacement.

Pass the same pixel_range used for detection so the MSN scale factor is computed over the same channel window; it defaults to the pixel_range in m.metadata (or the full spectrum) exactly like detection.

For each affected pixel, uses the 4-connected neighbor with the highest Pearson correlation (the MSN). The MSN index cached in result.msn_index at detection time is reused directly — the correlations are only recomputed for pixels where the cache is 0 (unknown). The cache reflects the detection pixel_range window: calling removal with a different pixel_range than detection keeps the detection-time neighbor choice rather than re-ranking within the new window. A scale factor is computed from non-flagged channels to match the pixel's intensity level. Flagged channels are then replaced with scale × MSN[k]. This preserves the pixel's overall intensity while removing only the spike shape, avoiding the dark-patch artifacts that raw neighbor median replacement produces.

Falls back to spectral interpolation (1D linear from nearest non-flagged channels) for edge pixels with no spatial neighbors.

Returns a new PLMap with recomputed intensity.

Example

cr = detect_cosmic_rays(plmap)
cleaned = remove_cosmic_rays(plmap, cr)

remove_cosmic_rays(m::TimeResolvedMatrix, result::CosmicRayMatrixResult) -> TimeResolvedMatrix

Replace flagged pixels with the median of their non-flagged 3×3 neighbors. If all 8 neighbors of a flagged pixel are themselves flagged, falls back to the global median of the image. Records :cosmic_rays_removed in metadata.

Returns a new TimeResolvedMatrix with cleaned data; does not mutate the input.