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.

Keywords

• mode="nfft" — use "nfft" for Non-uniform FFT (fast) or "dft" for direct DFT (exact)

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::AbstractArray{T,4}, ft::AbstractMatrix, data::AbstractMatrix{<:OIdata}; ...)

Compute the weighted chi-squared of a 4-D image cube x[px, py, wav, epoch] against a Matrix{OIdata} of size (nwav, nepoch).

Returns the raw chi-squared (not divided by ndof). For nwav > 1, cross-channel observables (differential phases, visibility amplitudes, OI_FLUX) are included automatically based on data metadata.

See also: image_to_chi2_fg, image_to_obs.

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.