Variational Image Feature Extraction for the Event Horizon Telescope

Imaging algorithms form powerful analysis tools for very long baseline interferometry (VLBI) data analysis. However, these tools cannot measure certain image features (e.g., ring diameter) by their nonparametric nature. This is unfortunate since these image features are often related to astrophysically relevant quantities such as black hole mass. This paper details a new general image feature-extraction technique that applies to a wide variety of VLBI image reconstructions called variational image domain analysis. Unlike previous tools, variational image domain analysis can be applied to any image reconstruction regardless of its structure. To demonstrate its flexibility, we analyze thousands of reconstructions from previous Event Horizon Telescope synthetic data sets and recover image features such as diameter, orientation, and ellipticity. By measuring these features, our technique can help extract astrophysically relevant quantities such as the mass and orientation of the central black hole in M87.


INTRODUCTION
Generating quantitative measurements about intrinsic radio images from very long baseline interferometry (VLBI) is a computationally and theoretically difficult task.Even in principle, it is impossible to associate a single unique image with a given set of visibility data.In practice, the small number of stations participating in the Event Horizon Telescope (EHT), and thus sparse coverage in the u-v plane, results in variety of potential image structures.As a result, the process of reaching quantitative conclusions about image features requires significant additional analysis.For the EHT analysis of the 2017 M 87 data, this has taken two, complementary forms (Event Horizon Telescope Collaboration et al. 2019a).
The first is a traditional Bayesian parametric modeling approach (Event Horizon Telescope Collaboration et al. 2019e, hereafter Paper VI). Simple geometric models are fit to the visibility data, and direct quantitative inferences about the image properties encoded within Corresponding author: Paul Tiede paul.tiede@cfa.harvard.eduthe model, e.g., diameter, are possible.Higher model fidelity improves the accuracy of the derived measurements, and thus this requires an input assumption that the key model features capture the "truth".
The second is a nonparametric approach, usually referred to as "imaging" (Event Horizon Telescope Collaboration et al. 2019d, hereafter Paper IV).This category of methods is broadly defined, and includes deconvolution algorithms like CLEAN (Högbom 1974;Schwarz 1978;Clark 1980;Schwab 1984) and forward modeling approaches like "maximum entropy" (Frieden 1972;Gull & Daniell 1978;Narayan & Nityananda 1986), regularized maximum likelihood (RML) analyses (Chael et al. 2016(Chael et al. , 2018;;Akiyama et al. 2017a,b), and Bayesian imaging (Broderick et al. 2020).The output of imaging by the EHT has been an ensemble of image reconstructions that reproduce the observed visibility data Event Horizon Telescope Collaboration et al. 2019d.The imaging methods have the significant advantage that they are extremely flexible, and therefore are reasonably expected to cover the "truth".However, unlike parametric modeling, imaging methods do not give direct quantitative measurements of the image features of interest, e.g., ring diameter, width, orientation, etc.Therefore, reaching quantitative conclusions about the image properties requires an additional processing step, which we call "feature extraction".It is this final step that is the subject of this paper.
Feature extraction is similar to the Bayesian parametric modeling applied in the image, rather than visibility, domain.However, it differs in one important respect: by virtue of being performed after imaging, the class of applicable "models", i.e., image features to be measured, is already well known.For example, in Paper IV and Paper VI, quantitative image features were extracted by the algorithm "ring extractor" REx (Chael 2019).However, REx is only applicable to images that have a dominant ring-like feature.In general, images from the EHT can have a complex structure and are dependent on the intrinsic source.The active galactic nuclei 3C 279 is displays a jet morphology which while poorly described by a ring, can be described by a set of Gaussians (Kim et al. 2020).
One possible approach to feature extraction is to "template" relevant image features using a transformation.For example, the Hough transform (Hough 1964;Duda & Hart 1972), is used to extract rings and other shapes from images using template matching.A related method is to approximate the complicated image reconstruction with parametric templates that describe the features of interest.This idea is more akin to the parametric/geometric modeling approach of visibility data in Paper VI and is the approach taken in this paper.
Any comparison requires a suitable quality metric, i.e., objective function.Because the total flux is typically an arbitrary rescaling and the image brightness is positive definite, there is a natural identification between the flux-normalized image and probability distribution.This motivates the use of "divergences" as an extremely flexible class of objective functions for comparing images; for this reason divergences have been used extensively in image processing (e.g., Goudail et al. 2004;Choi & Lee 2003;Aherne et al. 1998).
We adopt the method similar to variational inference (Blei et al. 2017), in which complicated distributions are approximated by simple parametric forms, with parameters estimated via the minimization of an appropriate divergence.In this paper, we develop a number of appropriate parametric forms and explore the performance of a variety of divergence for application to image feature extraction.Therefore, we call this method variational image domain analysis, or VIDA.The VIDA algorithm has been implemented in the open source package VIDA.jl 1 written in Julia (Bezanson et al. 2017).
There are two main reasons that imaging reconstruction followed by feature extraction is advantageous to directly fitting simple geometric models to the data.First, choosing the correct geometric model is difficult before imaging, meaning imaging is the first step in both methods.Second, fitting simple geometric models can "underfit" the data leading to biased results.For instance, in Paper VI additional "nuisance" Gaussians were required to obtain a reasonable reduced chisquare.On the other hand, Bayesian imaging techniques (e.g., (Broderick et al. 2020)) make minimal assumptions about source structure.Combining the Bayesian imaging posterior with VIDA then provides a deterministic map from image to feature posteriors.VIDA's feature posteriors can then be compared to the geometric modeling results testing whether the features are robust across methods.Therefore, VIDA fills a gap in the EHT modeling pipeline and is generic, unlike current EHT tools.
The layout of the paper is as follows: In Section 2, we present the details of VIDA.Different types of templates implemented, and the objective function used to find the best approximation to the true image are detailed.Section 3 applies VIDA to a variety of ring-like image reconstructions from the test set of Paper IV and compares the results to REx.This is an empirical demonstration that we can recover the optimal template, even though the objective function is non-convex.In Section 4, we demonstrate VIDA's flexibility by applying it to non-ring images from the test set of Paper IV. Finally, the conclusions are detailed in Section 5. Additional material related to the model and results for CLEAN is presented in Appendices A -C.

VARIATIONAL IMAGE DOMAIN ANALYSIS
The critical insight behind VIDA is that images (sans polarization and modulo total flux) and probability densities are in one-to-one correspondence.Namely, images are point-wise positive and integrable.Probability divergences are a natural class of objective functions used to compare two probability distributions.Furthermore, they have been used in image feature extraction and similarity measures before.Divergences thus form a natural objective function between the image reconstruction and the template.VIDA consists of three ingredients: 1. Image I(α, β) whose features we want to extract The choice of template used will depend on the structure of the image.For example, the images of M 87 from Paper IV are ring-like, while the reconstructions of 3C 279 from Kim et al. (2020) can be described by several Gaussian brightness distributions.In this section, we present the various templates that are implemented in VIDA.jl.

Gaussian template
To model a source of compact flux we include an asymmetric Gaussian template.The parameters of the Gaussian template are: 1.The size, σ = √ σ a σ b , where σ2 a,b are the variances in the principal directions of the Gaussian.
3. ξ, rotation angle (relative to the Gaussian center) of the principal axes measured east of north.
4. x 0 , y 0 , the center of the Gaussian.

Disk template
In Section 3 we will test VIDA.jl on a number of synthetic data tests.One of the test images is a disk.To approximate disks we use the template: ) where r 0 is the radius of the flat disk, α controls the smoothness of fall off and N ensures the template is normalized.The radial distance, r, is relative to the center (x 0 , y 0 ).When r 0 = 0 this template reduces to a symmetric Gaussian with standard deviation α.

Ring templates
One of the principal quantities of interest in images of M 87 is the ring diameter, d 0 , since it is related to the mass of the central black hole.Additionally, the ring is expected to have some thickness, w, due to the emitting material around the black hole.The simplest template would be a circular Gaussian ring with some thickness.Doppler boosting however can cause the emission to appear asymmetric.To model this a slash can be added to the ring template to include the brightness asymmetry.Additionally, the ring itself does not have to be circular.Ring ellipticity could occur from e.g., the emitting material not being azimuthally symmetric around the black hole.As well, due to the sparse coverage of the EHT array and imaging algorithms, ellipticity may be introduced into the reconstructions.Consolidating each of these features into a template parameter, we get the following: • d 0 : the geometric mean of the semi-major, a, and semi-minor, b, axis d 0 = 2 √ ab which is related to the area of the ellipse, π(d 0 /2) 2 .
• τ : the ellipticity of the ellipse, τ = 1 − b/a • ξ τ : the position angle of the semi-major axis measured north of east.
• w: the width of the Gaussian ring, defined to be the full width half max (FWHM) of the Gaussian, i.e. w = 2 √ 2 log 2σ, where σ is the standard deviation.
• s: the strength of the slash described in equation ( 4).
• ξ s : the position angle of the slash measured east of north.
The functional form of the template is given by: where S(x, y; s, ξ) is the slash function and d θ (x, y) is the minimum distance between the ellipse with parameters θ = d 0 , τ, x 0 , y 0 and the point x, y.
However, for an ellipse there is no analytical equation and instead once has to numerically minimize the function, subject to the constraint that e x , e y are points on the ellipse with parameters d 0 , τ, ξ τ .For the slash function S, we use a first order cosine expansion in azimuthal angle ϕ around the center x 0 , y 0 : where N 0 is a normalization factor to ensure the template is unit normalized.We restrict s ∈ [0, 1] to prevent image brightness from becoming negative.In the VIDA.jlpackage, this template is called GeneralGaussianRing (GGR), and an example reconstruction using said template is shown in Figure 1.
Additionally, VIDA.jl has a number of other ring-like templates currently implemented: • GaussianRing: Symmetric Gaussian ring with constant azimuthal intensity (i.e.GGR with τ, s = 0) • SlashedGaussianRing: Symmetric Gaussian ring with azimuthal slash described by Equation 4(i.e.GGR with τ = 0) • EllipticalGaussianRing: Elliptical Gaussian ring with constant azimuthal brightness (i.e.GGR with s = 0) • TIDAGaussianRing: GGR template where the slash and ellipticity position angle are either aligned or anti-aligned.
We also include a more general version of the GGR called the CosineRing{N,M}.This template is similar to the GGR template but where the width, σ, and slash function ( 4) are replaced by a higher order cosine expansion in azimuthal angle ϕ: where s, σ, ξ (.) ,, are vectors with the cosine expansion coefficients of the slash, standard deviation, and angular offset3 .We can reproduce the GeneralGaussianRing template by setting M = 1 and N = 0 in ( 5) and ( 6) respectively.This template can be used if the image has a ring-like feature that has a bumpy azimuthal profile.

Constant Template
We found that many image reconstructions had a diffuse intensity throughout the image due to poor dynamic range from sparse coverage and regularization effects.To model the background, we added a constant intensity template.This template is typically required to be included in any analysis to ensure reliable feature extraction.

Composite Templates
A general image reconstruction from a VLBI observation may have multiple image features.As such, VIDA.jl allows the user to combine multiple features into composite templates, where each individual component is given a relative flux4 .An example of this is shown in Section 4.1 where three Gaussian templates are used to model the reconstructed image.

Probability Divergences
As mentioned above, VIDA.jl uses an analogy between images and 2-D probability distributions to motivate the use of divergences as objective functions.Divergences form a general way to measure the similarity between two distributions.A divergence can be thought of as a functional F q [p] = D(p||q), comparing the distribution p (template) to a reference q (image), and is required to be non-negative, D(p||q) ≥ 0, and non-degenerate, D(p||q) = 0 if and only if p = q.Note that this definition is more general than a metric.Namely, a divergence does not have to be symmetric, i.e.D(p||q) ̸ = D(q||p) or satisfy the triangle equality.
One of the most well-known divergences is the Kullback-Leiber (KL) divergence (Kullback & Leibler 1951) or relative entropy, One issue with the KL divergence is its definition when the support of q and p differ, i.e. if q(x) = 0.In this case we set the contribution to the integral to be zero.
In addition to the KL divergence VIDA.jl also includes the Bhattacharyya divergence (Bh) (Bhattacharyya 1943), The Bh divergence is related to a well known metric on probability spaces, the Square Hellinger Distance : (9) Therefore, minimizing the Bh divergence is simply a least squares fit in the space of the square root of distributions.
In this paper we will present the results from optimizing the Bh divergence for two reasons.First, we found that while the KL and Bh divergence produce near-identical results, the Bh divergence required ∼ 25% less evaluations to converge.Second, the Bh divergence has preferable theoretical properties compared to the KL divergence.Namely, it is well defined when the image pixels have zero intensity and is symmetric.

Optimizing the Divergence
A problem when using probability divergences is that they give non-convex, non-linear optimization problems.Furthermore, the nature of the problem will change if the template changes, making an analytic analysis difficult.Therefore, to extract the globally optimal template, we turned to heuristic global optimizers, such as Genetic/Evolutionary strategies (see Das & Suganthan 2011, for a review) and simulated annealing (Goffe 1996).For this paper we used the Julia package BlackBoxOptim.jl5 6 .BlackBoxOptim.jl uses naturalevolution, and differential evolution strategies to perform a stochastic search of the parameter space.In the next section we will validate that our chosen optimizer is able to recover the optimal template reliably.

VALIDATING VIDA
To validate VIDA we need to analyze two related quantities.First, we need to verify that the objective function, i.e. the Bh divergence Equation 8, is robust to the artifacts that occur in image reconstructions.Namely whether the recovered parameter distribution contains the true value.Additionally, given the complex nature of the optimization problem we need to ensure that the chosen optimizer can recover the global minimum of the divergence.Our validation procedure will consist of: 1. Selecting an applicable ground truth image I truth (see the left column of Figure 2 2. For each truth image create a simulated EHT observation matching the observation characteristic of the EHT M 87 2017 observations.
3. Create an ensemble of image reconstructions of the simulated observations using the same procedure as Paper IV.
The infinitely thin ring is then convolved with a circular Gaussian with FWHM 10 µas.Two orientations (measured east-of-north) ξ = 180 • and ξ = 150 • are considered in this paper, and are shown in the top and left panels of Figure 2.For both orientations we took r 0 = 22 µas, s = 0.46, and F 0 = 0.6 Jy.After blurring the ring, it is important to note that the effective radius (the intensity peak) is smaller than the original radius of the ring (see Paper IV).The amount the diameter is biased inwards is given approximately by: where d true is the diameter of the non-convolved ring (44 µas), and α is the FWHM of the Gaussian kernel (α = 10 µas).Using (11) with d true = 44 µas gives d blur ≈ 43 µas.If we also consider the finite resolution of the EHT array (∼ 20µas) this is further decreased to ≈ 42 µas.Therefore, we expect both VIDA and REx to recover a diameter of 42 µas.Additionally, the slash strength is also modified by the convolution.Fitting the crescent with the GGR template we find s = 0.32, which is the value we will take as the ground truth below.

Step 2: Creating Simulated EHT Observations
While VIDA could be applied to the ground-truth images shown in the left column of Figure 2, this isn't applicable to what the EHT observes.The EHT is a verylong-baseline inferometer and instead observes complex visibilties, V (u, v) which are related to the on sky image through the van Cittert-Zernike theorem (Thompson 2017): In addition, atmospheric and telescope effects can further corrupt the signal.To model these corruption effects we use the eht-imaging package (Chael et al. 2016(Chael et al. , 2018) ) to generate realistic simulated data.

Step 3: Generating an Ensemble of Reconstructed Images from Simulated VLBI Data
To validate VIDA we used image reconstructions from the forward modeling or "regularized maximum likelihood methods" (RML), e.g.Honma et al. ( 2014 2018) and more specifically the eht-imaging package (Chael et al. 2016(Chael et al. , 2018)).The goal of RML methods is to find the image, I, that minimizes the objective function Following Paper IV, each χ 2 d is defined solely from the data products from the EHT telescope, e.g., complex visibilities.The second term encapsulates our additional assumptions, or regularizers, that are placed on the image.The α d β r , are the "hyperparameters" that control the relative weighting of the regularizers and data products.For the list of the regularizers used, see Paper IV.
In an attempt to model the uncertainty in the image reconstructions, we used the same set of imaging hyperparameters as Paper IV.The resulting set of image reconstructions is called the "topset" and results in 1572 reconstructions per dataset.
For completeness we also apply VIDA to the CLEAN M 87 reconstructions for the Crescent 180, GRMHD, double and Disk in Appendix C. Note that for M 87 there are only 30 reconstructions for each simulated dataset, making the comparison more uncertain.

3.4.
Step 4: Applying VIDA to the Ring-Like Image Reconstruction Ensembles VIDA was run on each set of image reconstruction ensembles using the GGR template with a constant brightness background whose relative flux was a free parameter, giving 9 parameters in total.Some example reconstructions and corresponding optimal templates are shown in Figure 2. The results are shown in Figure 3.For the crescent models, we were able to recover the expected diameter d, width w, and azimuthal orientation ξ s (black dotted lines in Figure 3.In addition, we compared the VIDA results to the REx method used in Paper VI. REx assumes that a single ring-like feature dominates the image reconstruction, and then finds the ring by finding the image location that leads to a ring with minimal radial dispersion.REx characterizes (see Appendix A and Paper IV for definitions) the ring through a diameter d, width w, brightness moments s and orientation ξ s and a fractional dispersion of the diameter f d .The diameter and width and brightness profile of the REx measurement are similar to VIDA's measurement with the GGR template.However, the fractional dispersion is not directly measured.Instead, VIDA measures the ellipticity of the ring.In Appendix A we demonstrate how f d and τ are related if the dominant source of radial dispersion in the ring is due to ellipticity.
The agreement between REx and VIDA is excellent for both the crescent images.The peak and overall width of the distribution for each parameter in Figure 3 are consistent between REx and VIDA.The ground truth values (black vertical lines) for the diameter, width, and brightness orientation ξ s are also consistent with the REx and VIDA results.We also analyzed the pairwise linearcorrelations between all the parameters and found that only the diameter and width were correlated (see Figure 9).This correlation is expected from Equation 11and occurs due to the finite resolution of the EHT array.For the crescent models we found good agreement between VIDA (blue) and REx (orange) for the ring diameter, width, slash, and location of the azimuthal brightness.Furthermore, the crescent distributions for the diameter (after accounting for Equation 11), width, slash strength s and position angle ξs are consistent with the truth image (black dashed line).For the ellipticity, τ , we found VIDA's value was systematically lower than REx's as expected since REx measures the fractional dispersion of the ring radius and random fluctuations can then create a ellipticity floor around 0.05 (see Appendix A for a detailed discussion).The ellipticity orientation, ξτ , is only recovered by VIDA, so there is no REx comparison.For the GRMHD simulations, the agreement between REx and VIDA is weaker.Namely, the distributions for the parameters are sometimes significantly different, although the portions of the distributions do overlap.Similar results were found using the simulated EHT coverage from April 5, 6, and 10.
To compare VIDA and REx's measurement of ellipticity, we first note that an additional processing step is needed since the two definitions differ.REx doesn't directly measure τ but instead measures the fractional diameter dispersion of the ring f d (see Equation A5 for a definition).If we assume that the ring's ellipticity dominates f d , then f d and τ are related by an invertible map.For more information about this conversion, see Appendix A. In Figure 3 we show REx's results after converting from f d to τ .Comparing the two measurements of τ , we see that REx's measurement is consis-tently greater than VIDA's.This bias is not unexpected given that when τ is small, the conversion described in Appendix A no longer applies.Instead, the fractional dispersion is dominated by random fluctuations in the ring diameter, creating a floor in f d .If we then naively apply the previous conversion, as was done in Figure 3, we will overestimate τ (see Figure 8).
VIDA also recovers the orientation of the ring ellipticity ξ τ .Interestingly, in all instances, we measure a similar distribution for ξ τ irrespective of the intrinsic image.This distribution is biased so that the semi-major axis of the ellipticity is in the north-south direction.By visually inspecting the image reconstruction ensembles (see Figure 2 for a typical example reconstruction) we confirmed this feature was present in the reconstructions and was not a bias from VIDA.We also find a non-linear correlation between the measured ring ellipticity and orientation.Namely, τ is significantly larger when ξ τ ≈ 0. The origin of this bias and implications for M 87 will be discussed in Tiede et al. (in prep.).

APPLYING VIDA TO ADDITIONAL SIMULATED IMAGE RECONSTRUCTIONS
In the previous section, we saw that VIDA and REx gave remarkably similar answers to problems that demonstrated similar ring-like structures.While REx is limited to ring extraction, VIDA can be applied to any image given a suitable template function.This section will explore VIDA's capabilities of extracting features from a broader range of potential sources.To accomplish this, we will consider the other non-ring test images from Paper IV: the symmetric disk and double Gaussian (see Figure 4).We will follow the same steps in the previous section to evaluate VIDA's performance.

Double Gaussian
Here we consider a source composed of a compact double with two circular Gaussian components.Each Gaussian has an FWHM of 20 µas.One of the Gaussians is placed at the origin and has a flux of 0.27Jy, which we will call the NW component.The other Gaussian is at ∆RA = 30µas and ∆DEC = −12µas and has a flux of 0.33Jy and will be called the SE component.This type of source could arise when looking at AGN using VLBI images, such as the recent 3C 279 results (Kim et al. 2020).
To extract the reconstruction's compact components, we used a template with three asymmetric Gaussian components and constant background.Two of the Gaussian components were allowed to be arbitrary, while the third Gaussian component was forced to be large (r 0 > 15 µas).The reason for the third Gaussian was there tended to be a region of additional emission around the two dominant Gaussian components in the image reconstructions.This diffuse emission can be seen in the top middle panel of Figure 4.If we didn't include this third component, we found that the Gaussian components tended to be quite large to soak up the extra emission.
VIDA's results for the double Gaussian are shown in Figure 5. Overall the on-sky size of each Gaussian, their separation, and flux ratio are contained in the support of the distribution.However, there do exist biases in the results, such as the FWHM of the Gaussian being biased low7 .Additionally, the ellipticity τ , is larger than zero; however, this is to be expected as asymmetry can only be added to the images.
For the simulated EHT sampling from April 6, we see that the ellipticity appears to be bimodal, and the parameter uncertainties are greater than the other days.This uncertainty was unexpected given the relatively good EHT coverage on April 6.After analyzing the image reconstructions we deemed this not to be due to .VIDA results of the diameter for the disk topset.The diameter is given by ( 14).Ignoring the simulated EHT April 10 sampling, we consistently find that the diameter is 4µas smaller than the original image.The origin of this discrepancy is discussed in Figure 7.For the simulated data using the April 10 sampling, which has poor coverage compared to the other days, the imaging gives two modes.One mode is similar to the other days, while the second fails to show a coherent disk structure giving the second peak in diameter at 60µas.
VIDA.Instead, a subset (10 − 15%) of image reconstructions on April 6 that exhibit a third bright Gaussian feature.If we remove these reconstructions, we find that the results for the simulated April 6 coverage are consistent with the other days.

Disk Image
The intrinsic image is a symmetric flat disk with a diameter of 70 µas, which is then convolved with a Gaussian with an FWHM of 10 µas.The true image, an example reconstruction, and optimal template for that reconstruction are shown in the lower left, middle panel, and right panel of Figure 4.
To encode the diameter of the disk we use FWHM of the disk template: where r 0 and α are described in Equation 1.We fit the disk template to the ground truth image to calibrate the diameter definition to the disk's true diameter.We found that the optimal template for the true image had a d disk ≈ 69 µas.If we convolved the image by an additional 20 µas to take into account the finite resolution of the EHT array, we found d disk ≈ 68 µas.This is the value we use as the ground truth diameter in all comparisons below.Figure 6 displays the results for VIDA applied to each day.Ignoring the simulated EHT sampling of April 10, which has poor coverage compared to the other days, we find that the results are very consistent between days.For the simulated April 5 EHT sampling, we find the median diameter d disk = 65.7 +0.42 −0.76 µas where the range are the 68% interval about the median.Similarly, for the simulated EHT sampling of April 6th, we find d disk = 65.6 +0.40 −0.47 µas, and for April 11th d disk = 65.5 +0.45 −0.80 µas.This demonstrates that VIDA is robust to the slight dif-ference in image reconstructions from different baseline coverage.
For the simulated EHT sampling of April 10, we had a different result finding a bimodal diameter.Analyzing the reason for this, we found that images with d disk ≈ 60 µas had a markedly different structure than the rest of the images.Given the distinct non-disk structure of the image, it is no surprise that VIDA struggles at recovering the correct diameter.
Comparing our result for the diameter to the true value, 68µas, we find that our result has a consistent bias of ≈ 2.4µas for the simulated April 5, 6, and 11 EHT sampling.Again, this appears to be an artifact of the imaging process.In Figure 7, the radial profiles of the truth (dotted lines), averaged image reconstructions8 , and optimal templates are compared.VIDA does an excellent job of recovering the size of the images, which are similarly biased toward smaller radii.This suggests that the diameter bias is intrinsic to the topset used for this disk.Furthermore, as shown in Appendix C the CLEAN reconstructions of the disk do not suffer from the same bias.Given this, it is unlikely that observed bias is due to the finite resolution of the EHT array and is intrinsic to the top set used.

SUMMARY AND CONCLUSIONS
We present VIDA, a new image feature extraction technique appropriate for use by the EHT.VIDA adopts a forward modeling approach to extract quantitative image properties by approximating the image with a parameterized family of functions that encode the desired image properties.
A key feature of VIDA is its flexibility.Multiple image components have already been implemented, from which composite models of significant complexity can be constructed.These include ring-like templates of particular relevance to EHT images, Gaussians, and constants.The ability of VIDA has been demonstrated for several sources, each with over a thousand reconstructions.These include image reconstructions from simulated data produced from double Gaussians, slashed rings, and GRMHD simulations.In all cases, key quantitative features were accurately recovered where they appeared in the underlying image reconstructions.These include separations, orientations, ring diameters, widths, brightness profiles, and multiple measures of ellipticity.Application of these to the EHT observations of M 87 will be explored in future work.
The applicability of VIDA extends beyond EHT observations of M 87.The ability to create composite models with multiple components is naturally relevant to VLBI reconstruction of AGN, such as 3C 279, that is composed of multiple compact features (e.g., Kim et al. 2020).
It should be noted that image feature extraction methods, like VIDA, are generally most useful when strong pri-ors may be placed on the image structure itself.That is, VIDA is primarily a method for quantifying what is already qualitatively apparent.Poorly chosen models can lead to significant parameter biases, as seen in Section 4.1, where an extra Gaussian blob was required to achieve acceptable results.However, because VIDA is an image characterization tool, not an imaging tool in and of itself, this presents only a very modest limitation on its utility.
Software: BlackBoxOptim.jl, eht-imaging, GR (Heinen et al. 1985(Heinen et al. -2019)), Julia (Bezanson et al. 2017), matplotlib 3.3 (Hunter 2007) In this case we see that CLEAN is able to recover the diameter of the disk unlike eht-imaging, which had underestimated the diameter by a few µas (see Figure 10).For the double we see that VIDA and the CLEAN reconstructions are able to recover the separation and flux ratio.The true ellipticity is also recovered.Note that the larger τ and fwhm for the Gaussian's follows from incorporating the contributions of the CLEAN beam (gray dotted line).

Figure 1 .
Figure 1.An example of VIDA run.Left: image reconstruction of a GRMHD simulation from Paper V and Paper IV. Middle:VIDA reconstruction using the GeneralGaussiantemplate and Bh divergence.Right: vertical (red) and horizontal (blue) chords through through the center of light of the truth image.The dashed lines are for the VIDA optimal template, and solid for the reconstruction.This plot can be made in the VIDA.jlpackage using the triptic function and is used to assess the quality of the template approximation.

Figure 2 .
Figure 2. Images used for the imaging validation from Paper IV.We considered 3 models: top row crescent with position angle ξ = 180 • north of east, middle row crescent with position angle ξ = 150 • and bottom row GRMHD simulation.The left column shows the truth image, the middle columns an example reconstruction, and right the optimal VIDA template applied to the reconstruction.hydrodynamical (GRMHD) simulation from Event Horizon Telescope Collaboration et al. (2019f).The geometric crescent model is described by: ); Bouman et al. (2016); Akiyama et al. (2017a,b); Ikeda et al. (2016); Kuramochi et al. (

Figure 3 .
Figure3.Results from VIDA (blue) and REx (orange) being applied to the Paper IV reconstructions using the simulated April 11 EHT sampling of Figure2.The top and middle rows show the geometric crescent results and the bottom the GRMHD.For the crescent models we found good agreement between VIDA (blue) and REx (orange) for the ring diameter, width, slash, and location of the azimuthal brightness.Furthermore, the crescent distributions for the diameter (after accounting for Equation11), width, slash strength s and position angle ξs are consistent with the truth image (black dashed line).For the ellipticity, τ , we found VIDA's value was systematically lower than REx's as expected since REx measures the fractional dispersion of the ring radius and random fluctuations can then create a ellipticity floor around 0.05 (see Appendix A for a detailed discussion).The ellipticity orientation, ξτ , is only recovered by VIDA, so there is no REx comparison.For the GRMHD simulations, the agreement between REx and VIDA is weaker.Namely, the distributions for the parameters are sometimes significantly different, although the portions of the distributions do overlap.Similar results were found using the simulated EHT coverage from April 5, 6, and 10.

Figure 4 .
Figure 4. Results from applying VIDA to the non-ring test images from Paper IV where the left column shows the ground truth image, the middle an example reconstruction from the topset, and right the optimal template from applying VIDA to the reconstruction.Top: Shows the results for the double image.Bottom: Shows the results for the disk.Overall VIDA can recover the intrinsic structure of both images.

Figure 5 .
Figure5.VIDA results for the two compact Gaussian components (blue for SE component and orange for the NW) in the double Gaussian test image.Each row corresponds to a different EHT sampling from the 2017 M 87 observations.The green curves are for parameters that are a combination of the SE and NW components.On all days the true values are included in the support of the parameter distributions found by VIDA.Note, that the broad distribution found using the simulated EHT April 6 sampling is due to an imaging artifact as discussed in the paper.
Figure6.VIDA results of the diameter for the disk topset.The diameter is given by (14).Ignoring the simulated EHT April 10 sampling, we consistently find that the diameter is 4µas smaller than the original image.The origin of this discrepancy is discussed in Figure7.For the simulated data using the April 10 sampling, which has poor coverage compared to the other days, the imaging gives two modes.One mode is similar to the other days, while the second fails to show a coherent disk structure giving the second peak in diameter at 60µas.

Figure 9 .Figure 12 .
Figure9.Joint density plot for the crescent image reconstruction ensemble using the April 11 EHT simulated coverage.For VIDA (blue, lower), we find that the truth values (black lines) are within the recovered distributions.For REx (orange, upper) we find that for all parameter except the ellipticity τ are recovered.REx's ellipticity has a floor around τ = 0.05 which is expected since REx measures fractional dispersion as is explained in Appendix A. Additionally, no ξτ is measured for REx so those columns are not shown.
Results when applying VIDA to the average reconstruction from the topset using the simulated April 11 EHT sampling.The images were normalized to unit flux and centered before averaging.The average radial profile is shown in solid blue.Comparing this to the optimal template (orange dashed line) and the true profile (dotted lines), we see that the optimal template matches the reconstruction's radial profile but underestimates the ground truth image.