virgil.imaging

The starlet (isotropic undecimated wavelet) transform of an image.

The à-trous algorithm with a B3-spline kernel (Starck, Murtagh & Fadili 2010): the image is smoothed scales times, the kernel's taps twice as far apart each time, and each detail plane is the difference between successive smoothings. Detail plane j holds structure about 2**j pixels across. Beyond the edges the image is zero.

Parameters:

Returns:

Build an image from components, one at a time: gradient CLEAN.

Högbom's CLEAN repeatedly finds the peak of the residual dirty image and adds a fraction (the loop gain) of a point source there. For data that are linear in the image, the residual dirty image is proportional to -∂χ²/∂c, the gradient of χ² with respect to the flux c_p of a point at each pixel. This function uses that gradient directly, so it works for any data virgil can fit: closure phases, kernel or DISCO phases, squared visibilities, or a mix. J e_p is the change in the whitened residuals per unit flux at pixel p. Each iteration picks the pixel whose Gauss–Newton step would lower χ² most, the largest g_p² / |J e_p|² with g_p < 0, and adds gain times that step, -g_p / (2 |J e_p|²). This is matching pursuit; for linear data, where |J e_p| is the same everywhere, it is exactly Högbom's CLEAN. The norms |J e_p| are computed once, at the start, with one Jacobian–vector product per pixel. Normalizing by them matters next to an analytic star, where flux in a pixel is nearly the same as the star's: the gradient there is small, but so is |J e_p|, and an unnormalized search would pile flux beside the star.

Components are added to a fixed base scene, usually an analytic star at flux 1; fit its parameters first. Their fluxes are relative to the base's, like a companion's in a System; for a System base they are siblings of its components, relative to its total. By default the components are grey, the same fraction of the base's flux at every wavelength; give a spectrum to make them follow one, as the image does in SPARCO (e.g. a star with its own spectral index and an environment with another). Without a base, the components alone make the image, starting from one at the centre of the grid (closure phases do not fix the position, so this is also the anchor).

Each iteration only adds flux, so an early step that overshoots cannot be undone by later ones. Every refit_every iterations a major cycle (as in Clark and Cotton–Schwab CLEAN) refits the fluxes of all the components at once, by non-negative least squares on the linearized residuals: flux can move between components, and components whose flux falls to zero are removed. The fluxes stay non-negative. A last major cycle runs when CLEAN stops at the target or at max_iterations, so that the fluxes it returns are refitted even when it stops early (when the data start close to the target, a major cycle may otherwise never run).

The base is fixed unless base_priors names some of its parameters, e.g. a companion's position or a star's diameter. Each major cycle then fits those together with the components' fluxes, with fit: fitting them in turn instead would let the components absorb an error in the base, or let the base absorb flux it does not model. A companion much closer to the star than λ/B has its flux and separation nearly degenerate, whatever fits it.

Multi-scale CLEAN (Cornwell 2008, IEEE JSTSP 2, 793): with several scales_mas, each component is a pixel convolved with a circular Gaussian of FWHM s (s = 0 is a point), and the image is Σ_s G_s * F_s. The search runs over every (scale, pixel), with the same score g² / |J e|², and the major cycles refit the components at all scales together. An extended source then takes a few broad components instead of many points. Each shape is cut to the support and renormalized there, so a component's flux is the flux it puts in the image. Cornwell multiplies the peak residual at each scale by a bias 1 - 0.6 s / s_max that favours small scales, because the peak of a smoothed residual grows with the scale's area. The score here is instead the χ² decrease of each component's Gauss–Newton step, which does not depend on its size, so no bias is needed; scale_bias applies one, 1 - scale_bias s / s_max, if the broad components grab flux that belongs to points.

Iteration stops when χ² per data point reaches target_chi2_red (the discrepancy principle; it relies on correct error bars), or after max_iterations. It also stops when χ² has stalled: when no pixel lowers it, or when it has fallen by less than a fraction stall_tolerance over the last stall_window iterations, and a major cycle does not help. That happens when noise puts the truth's own χ² per point above the target; chi2_red shows how close it came.

Parameters:

Returns:

A starting model for an image fit, sized from the data.

It fits a quick parametric model, an analytic star (if star) plus a circular Gaussian envelope with free flux and width, a low-order description of the visibilities, from a few starting widths. Then it chooses the image:

• a field of about six FWHMs of the envelope, but at least 500 mas, and never larger than the interferometric field of view (λ / B_min) or largest_mas;
• pixels oversample times finer than the Nyquist scale;
• an Image with the fitted flux, whose pixels are the fitted Gaussian (start="moments") or the positive part of the :func:dirty_image (start="dirty", with the star removed). A dirty start is better when the Fourier coverage is dense (it already has the right shape) and worse when it is sparse (its sidelobes dominate).

Fit it with fit(start, image_priors(start), data, regularizers), adding a prior on the image's flux, "env.flux" (or "flux" with star=False), to fit it too: with a star, this stops the fit from parking excess flux next to it.

Parameters:

Returns:

The dirty image: direct synthesis of the visibilities, unregularized.

I(x) = Σ_k w_k Re[V_k exp(+2πi u_k · x)] / Σ_k w_k, summing over the samples (each standing for itself and its conjugate) with uniform weights on the informed ones: the image convolved with the dirty beam, plus noise. It peaks at one for a lone point source. For AMIGO DISCO data the visibilities are the least-squares estimate from the modes; since the modes are blind to the total flux, the estimate is normalized to |V| = 1 on average on the shortest baselines.

Parameters:

Returns:

Priors on the log-brightness of every Image in a scene.

Returns a priors dict for fit. An Image with a plain log-brightness array gets a flat prior, {"env.log_brightness": ImproperUniform(...)}, so that its pixels are constrained only by the data and the regularizers. An Image with a GaussianField gets standard-normal priors on the field's latents, {"env.log_brightness.latent": Normal(0, 1)}: the Gaussian-process prior, which needs no regularizer. Add priors for any other free parameters (fluxes, offsets) to the dict.

The largest pixel scale (mas) that samples the data's finest fringes.

That is λ / (2 B) for the longest baseline B in units of wavelength, over all samples of data (an OIData or a sequence of them). Reconstructions normally use pixels 2–4 times smaller.

The default field of view for an image of data, in mas.

The smaller of largest_mas (500 by default) and the interferometric field of view, the largest λ / B over all samples with a non-zero baseline B (the shortest baseline at the longest wavelength, where they are observed together): structure larger than that is not measured, and on a uv lattice a larger field would alias. Use it with nyquist_pixel_scale to choose the image's size and pixels.

The resolution of a dataset: the FWHM ellipse of its beam's core.

The dirty beam is the image of a point source made by the data's uv sampling. Near its peak it falls off as 1 - 2π² xᵀ M x, where M is the weighted mean of u uᵀ over the samples, which matches a Gaussian of covariance M⁻¹ / 4π². Its FWHM ellipse is the standard "beam" drawn on interferometric images: roughly λ/B, and elongated where the coverage is. For complicated coverage the real beam can look quite different, but the ellipse still shows the size and shape of a resolution element.

Every sample that carries information is weighted equally ("uniform weighting"); in AMIGO DISCO data those are the uv points on which the modes have weight (at least 1e-3 of the largest), i.e. the splodges. Weighting by information instead would let the low-frequency central splodge dominate, and give a broader beam.

Parameters:

Returns:

An image convolved with a Gaussian beam: the image at the data's resolution.

Reconstructed images are super-resolved. A regularized image can put structure on scales finer than the beam, where the data constrain it only weakly, so it often looks clumpy or streaky. Convolved with the beam (the "restored" image of radio astronomy), it shows only what the data resolve. That is the fair way to compare two reconstructions, or a reconstruction with a model: convolve both.

The kernel is an elliptical Gaussian with the beam's FWHMs and position angle (North to East), normalized to unit sum. It is sampled on an odd grid centred on a pixel, so the convolution does not shift the image. Flux beyond the edge of the image is taken to be zero, and the result is cropped to the image, so the total flux is kept only for structure more than about a beam from the edge; structure nearer the edge is dimmed.

Parameters:

Returns:

Fit a model over a range of weights for one regularizer.

The weights are fitted from largest to smallest, each starting from the previous solution, which is faster and more stable than starting every fit afresh. Plot penalty against chi2 (both on log axes) to see the trade-off between fitting the data and regularizing the image, and compare LCurve.corner and LCurve.discrepancy with the images either side: there is usually a wide range of good weights.

Parameters:

Returns:

Laplace-approximated log evidence of a Gaussian-field image fit.

For an Image whose log-brightness is a GaussianField with standard-normal latents z, at the MAP model from fit,

log Z ≈ -½ χ² - L - ½ |z|² - ½ log det(I + JᵀJ),

up to a constant that depends only on the data (see below). J is the Jacobian of the whitened residuals with respect to z, so JᵀJ is the Gauss–Newton curvature of the likelihood. Other fitted parameters (fluxes, spectra) are held at their MAP values. Compare it across fits with different hyperparameters and choose the largest, as MacKay's evidence framework does; it is exact for a linear model, and assumes correct error bars. It uses the data's quoted errors, so it does not support fits with noise= terms, nor fits with one model per dataset.

Normalization. -½ χ² - L is model_loglike at the MAP less its data-only normalization, -Σ log σ - ½ n log 2π (in its von Mises form for unprojected phases, and with the correlated closure phases' determinant), which is the same for every model and hyperparameter on a given dataset and so cancels in differences of log Z, the only meaningful quantity. For plain data L = 0. L is the rest of the likelihood's normalizer, which can depend on the model and on nuisance widths: for data with calibration gains (OIData.with_gains) or closure offsets (OIData.with_closure_offsets), which the likelihood marginalizes, it is ½ log det of their covariance factor, and for extra observables (OIData.extras) the log ratio of their effective to quoted errors. It is evaluated at the MAP exactly as the likelihood evaluates it, so that evidences for data with different gain or offset widths can be compared; its curvature with respect to z is left out of J, as Gauss–Newton leaves out the residuals' own second derivatives.

Parameters:

Returns:

Draws from the Laplace posterior of a Gaussian-field image.

For an Image whose log-brightness is a GaussianField, the prior on its whitened latents z is N(0, I), so about the MAP z₀ the Gauss–Newton (Laplace) posterior is

z ~ N(z₀, (I + JᵀJ)⁻¹),

with J the Jacobian of the whitened residuals with respect to z: the same approximation as log_evidence. With the thin SVD J = U S Vᵀ, a draw is z₀ + V diag((1 + s²)^-½) ε₁ + (I - V Vᵀ) ε₂ with standard-normal ε₁ and ε₂, so directions the data do not constrain keep the prior's unit width. The draws show which features of the image the data support. Only the latents vary: every other fitted parameter (fluxes, the star, spectra) is held at its MAP value, so the spread is narrower than the full posterior's where they are correlated with the image. Like log_evidence, it uses the data's quoted errors.

Parameters:

Returns:

Re-estimate the scale of the error bars from a Gaussian-field fit.

What it does. It estimates the factor s by which every error bar should be multiplied for the data to be consistent with the fit. s < 1 means the error bars are too large; s > 1 that they are too small, or that the model is missing something.

The idea. In MacKay's evidence framework, the level of the noise is a hyperparameter, like the prior's sigma and length_mas. Write the noise precision as β = 1/s², so that the likelihood is exp(-β χ²/2), with χ² computed with the quoted errors. The Laplace-approximated evidence, as a function of β, is maximized when

N is the number of data. γ is the effective number of parameters the data measure: the λ_i are the eigenvalues of the Gauss–Newton curvature JᵀJ of the likelihood (with the quoted errors) in the field's whitened latents, in which the prior's curvature is the identity, so that βλ_i are those of the rescaled likelihood. A direction with βλ ≫ 1 is fixed by the data and counts as one parameter; one with βλ ≪ 1 is fixed by the prior and counts as none. Since γ depends on β, the equation is solved for β at the fitted model (by Newton's method, which converges monotonically from β = 0 because γ is concave in β).

Each measured parameter uses up one datum's worth of scatter. So of the N residuals, only N − γ are free to scatter, and an honest error bar gives χ² ≈ N − γ, not N. The ordinary "χ² per point" estimate, s² = χ²/N, is biased low for the same reason as the 1/N estimate of a sample variance; this is its Bayesian, nonlinear generalization. It is MacKay's re-estimation formula for β (MacKay 1992, eq. 4.10, with γ from eq. 4.9; Bishop 2006, eqs. 3.91–3.95).

How to use it. Fit at your chosen hyperparameters, call this, rescale the data with OIData.with_error_scale, and refit. The fit itself depends on the errors, so in principle fitting and rescaling is a fixed-point iteration; in practice one rescaling usually suffices. Error bars that are too large make the discrepancy principle, classic MaxEnt and the evidence all over-regularize, so rescale before choosing hyperparameters with any of them. The estimate assumes the model is adequate: if the data contain structure the model cannot fit, s absorbs it.

Only the field's latents are counted in γ. Each other fitted parameter (the image's flux, a star's position, a spectral index) that the data measure lowers N − γ by about one more, and so raises s by a fraction of about 1/(2N). That is negligible while such parameters are few compared with the data, as in every SPARCO fit.

One scale per observable. A single s assumes every error bar is off by the same factor. Often it is not: the V² and closure phases of one instrument are calibrated differently, and in a MATISSE N-band contest file the V² gave χ² per point 0.005 while the closure phases gave 0.49. One scale between the two leaves the V² overweighted and the phases underweighted. With by_observable=True each kind of observable b ("vis", "phi" and each kind in extras, pooled over the datasets) gets its own precision β_b on its own rows of the residual vector. The evidence is then maximized where, for every block,

with G_b = J_bᵀJ_b the curvature from block b's rows of the Jacobian: MacKay's re-estimation of β applied to each block (MacKay 1992, eqs. 4.9–4.10; Bishop 2006, §§3.5.2–3.5.3). Equivalently, γ_b is the sum over the block's rows of the diagonal of the hat matrix B½ J A⁻¹ Jᵀ B½ (B the diagonal matrix of each row's β_b): the share of the measured parameters that block b pays for. The γ_b add up to the total γ = tr(A⁻¹(A − I)). They are coupled, so the equations are solved together by iterating β_b ← N_b / (χ²_b + γ_b/β_b) (the same fixed point as MacKay's β_b ← (N_b − γ_b)/χ²_b, but monotone, so it cannot oscillate) from β_b = N_b/χ²_b to a relative change of 1e-10. With one block it reproduces the single scale. The blocks' covariances must be s_b² D: data with calibration gains, closure-phase offsets, marginalized flux scales or differential phases with a finite prior_width add nuisance covariance that does not scale with the quoted errors, and raise a ValueError (the single scale makes the same assumption, so treat it with care for such data). N_b counts independent data: closure phases from four or more telescopes count their independent combinations, while their periodic penalty residuals (see whitened_residuals) belong to the block's χ² and γ_b but are not counted as data.

Use it when the blocks' χ² per point differ markedly, e.g. V² and closure phases from different calibrations, or OI_FLUX spectra beside interferometry. Rescale with the dictionary it returns, data.with_error_scale(scales), and refit. Prefer the single scale when a block has few data: its s_b scatters by about 1/√(2(N_b − γ_b)). Like the single scale, it absorbs into s_b any structure in block b the model cannot fit.

It uses the data's quoted errors, so it does not support fits with noise= terms (which estimate the errors another way), nor fits with one model per dataset.

Parameters:

Returns:

Check a model and its data for common imaging pitfalls.

Nothing is printed or warned: print the returned Diagnosis to read it. Run it on the fitted model, with the regularizers used in the fit.

Parameters:

Returns:

Sky offsets (mas) of the pixel centres along one image axis.

Along columns this is dra, decreasing to the right (East left); along rows it is ddec, decreasing downwards (North up). Both are (centre - index) * pixel_scale_mas with the centre at index (npix - 1) / 2, as in :func:image_coordinates.