Imaging (Setup)

setup_ft(data, nx, pixsize; mode="nfft")

Set up Fourier transform plans for the given data matrix.

data is a Matrix{OIdata} of size (nwav, nepoch) (as returned by readoifits or readoifits_multiepochs). Returns a matching matrix of FT plans, in the same floating-point precision as data (so Float32 data gives Float32 plans).

Keywords

• mode="nfft" — use "nfft" for Non-uniform FFT (fast) or "dft" for direct DFT (exact)
• fftflags=FFT_FLAGS — FFTW planning strategy (default FFTW.MEASURE; see FFT_FLAGS). Pass FFTW.ESTIMATE to skip the one-off measurement at the cost of a slower transform.

ft_info(ft)

Print a summary of a Fourier-transform plan set: mode, image size, channels and uv counts.

julia> ft = setup_ft(data, 128, 0.25);
julia> ft_info(ft)
OITOOLS NFFT Fourier transform plans: 8 wavelength(s) × 1 epoch(s)
  Image size: 128×128 pixels
  UV points: 380–412 (total 3168)

Named rather than a Base.display method: display and AbstractMatrix both belong to Base, so a method on them would change how every matrix in the session displays, including in unrelated packages.

OIft{T,C}

Fourier transform plans for a (nwav, nepoch) grid of channels, as built by setup_ft.

Indexes like the matrix it wraps — ft[w, t] is one channel's NFFTCell or DFTCell — and additionally records what was built:

Deliberately does not hold the OIdata it was built from: that would double the data's lifetime and let the two drift apart once the data is filtered. Callers pass both.

NFFTCell{T}

One channel's NFFT plans: the full uv transform and one per observable index set.

Named access (cell.v2, cell.t3_1) says which transform is meant. Positional access is kept, and AbstractVector is the supertype, so the six plans can still be indexed 1:6 in the order setup_nfft has always produced them.

DFTCell{T}

One channel's dense DFT kernel, nuv × nx².

A DFT covers every uv point in a single matrix, so unlike NFFTCell it has no per-observable sub-transforms: the observables are selected by indexing the result. Asking for cell.v2 therefore raises, rather than returning something that looks usable.

gaussian2d(n, m, sigma; T=Float32)

Unit-flux m×n Gaussian of width sigma pixels. T defaults to Float32 to match the default precision of readoifits (and hence of setup_ft's plans).

image_to_residuals(x, ft, data)

Compute normalised residuals (model - data) / error for each observable type. Returns a NamedTuple with fields v2, t3amp, t3phi, visamp, visphi. Phase residuals (t3phi, visphi) are wrapped to [-180, 180] before dividing by error.

image_to_chi2(x::AbstractMatrix, ft, data::OIdata; ...)

Mono convenience: compute chi-squared for a single 2-D image against a single OIdata.

image_to_chi2_fg(x::AbstractArray{T,4}, g::AbstractArray{T,4},
                 ft::AbstractMatrix, data::AbstractMatrix{<:OIdata}; ...)

Compute chi-squared and its gradient w.r.t. the raw (unnormalised) image pixels.

The gradient includes the flux-normalisation chain-rule correction: for each (w,t) cell, g[:,:,w,t] is corrected so that g = (g_raw .- ⟨y, g_raw⟩) / S where y = x/S and S = sum(x).

Returns the raw chi-squared (not divided by ndof).

image_to_chi2_fg(x::AbstractMatrix, g::AbstractMatrix, ft, data::OIdata; ...)

Mono convenience: compute chi-squared and flux-corrected gradient for a single 2-D image against a single OIdata.

crit_fg(x::AbstractArray{T,4}, g::AbstractArray{T,4},
        ft::AbstractMatrix, data::AbstractMatrix{<:OIdata}; ...)

Unified criterion: (χ² + regularization) / ndof, with gradient.

The flux-normalisation chain-rule correction is applied to the combined (χ² + regularization) gradient per cell, matching the original behaviour.

Returns the criterion value. g is modified in place.

Regularization structure

Three categories of regularizers, each a list of ["name", μ, ...] tuples:

• 2-D / spatial (regularizers): applied independently to each (wavelength, epoch) cell, e.g. [["tv", 1e-3], ["l1l2", 1e-3, 0.01]]
• Transspectral (transspectral_regularizers): cross-wavelength penalties applied per epoch when nwav > 1, e.g. [["transspectral_tv", 0.1]]
• Temporal (temporal_regularizers): cross-epoch penalties applied when nepoch > 1, e.g. [["temporal_tvsq", 0.01]]

Keywords

• weights = [1.0, 1.0, 1.0]: relative weights for (V², T3amp, T3phi)
• regularizers = []: 2-D per-cell regularizers
• transspectral_regularizers = []: cross-wavelength regularizers
• temporal_regularizers = []: cross-epoch regularizers
• epochs_weights = []: per-epoch scaling (default: uniform)
• use_diffphases = false: force differential-phase fitting
• verb = false: verbose output
• vonmises = false: von Mises loss for closure phases

crit_f(x::AbstractArray{T,4}, ft::AbstractMatrix,
       data::AbstractMatrix{<:OIdata}; ...)

Unified criterion: (χ² + regularization) / ndof, forward-only (no gradient).

Skips adjoint/gradient computation for speed. Accepts the same keywords as crit_fg (minus g).

lcurve(x_start, data, ft, reg_name, mu_range;
       reg_args=[], fixed_regularizers=[["centering", 1e3]],
       pixsize=nothing, pixsize_ref=0.4, use_dimless=false,
       weights=[1.0,1.0,1.0], maxiter=500, restarts=3,
       verb=false, vonmises=false)

Sweep a single regularizer's weight over mu_range and return L-curve data.

• reg_args: extra arguments for the regularizer (e.g. [1e-3] for l1l2's α).
• If use_dimless=true and pixsize is given, mu_range values are treated as dimensionless μ̃ and converted via scale_regularizers.

Returns a NamedTuple with fields:

• mu: the input weight values
• chi2r: reduced χ² (divided by ndof) at each point
• reg_val: unweighted regularizer functional for the swept regularizer only
• ndof: number of degrees of freedom
• images: vector of reconstructed images

lcurve_elbow(chi2_vals, reg_vals)

Find the L-curve corner (maximum discrete curvature) on a log-log plot. Returns the index of the elbow point.

lcurve_normalize(chi2_vals, reg_vals; ref_index=nothing)

Normalize L-curve values by dividing each axis by its value at a reference point (Renard+ 2011, Fig. 15). If ref_index is not given, the elbow point is used. Returns (chi2_norm, reg_norm).

scale_regularizers(regularizers, pixsize; pixsize_ref=0.4)

Convert dimensionless regularizer weights μ̃ (calibrated at pixsize_ref mas) to actual weights μ for a given pixsize (mas), using the scaling laws from Renard, Thiébaut & Malbet (2011, A&A 533, A64, Appendix C).

For unit-flux normalized images in 2-D, the optimal μ scales as a power of the pixel size. This function applies μ = μ̃ · (pixsize_ref / pixsize)^α where the exponent α depends on the regularizer type:

The returned list has the same structure as the input, with rescaled weights.

scale_tv_epsilon(pixsize; pixsize_ref=0.4, eps_ref=1e-8)

Scale the TV relaxation parameter ε so that the regularization is pixel-size independent (Renard+ 2011, Appendix C.3). ε should scale linearly with pixsize.