---
abstract: 'We present a Bayesian Voronoi image reconstruction technique (VIR) for interferometric data. Bayesian analysis applied to the inverse problem allows us to derive the *a-posteriori* probability of a novel parameterization of interferometric images. We use a variable Voronoi diagram as our model in place of the usual fixed pixel grid. A quantization of the intensity field allows us to calculate the likelihood function and *a-priori* probabilities. The Voronoi image is optimized including the number of polygons as free parameters. We apply our algorithm to deconvolve simulated interferometric data. Residuals, restored images and $\chi^2$ values are used to compare our reconstructions with fixed grid models. VIR has the advantage of modeling the image with few parameters, obtaining a better image from a Bayesian point of view.'
author:
- 'G. F. Cabrera, S. Casassus'
- 'N. Hitschfeld'
title: |
    Bayesian Image Reconstruction\
    Based on Voronoi Diagrams
---

Introduction
============

Astronomical interferometric data result from the addition of instrumental noise to the convolution of the sky image and the instrumental response. Because of incomplete sampling in the $(u, v)$ plane, obtaining sky images from interferometric data is an instance of the inverse problem, and involves reconstruction algorithms.

The CLEAN method consists of modeling the side-lobe disturbances and subtract them iteratively from the dirty map [@CLEAN]. The CLEAN method works well for low noise and simple sources. But if the source has many complex features, or if the data is too noisy, CLEAN will do only a few iterations returning a noisy image [@CLEAN]. Another shortcoming is that CLEAN involves some ad-hoc parameters (the loop gain, stopping criteria, clean beam) that bias the final reconstruction, in the sense that CLEAN can give many different reconstructions for the same dataset. The maximum entropy method (MEM) finds the image that simultaneously best fits the data, within the noise level, and maximizes the entropy $S$. This is done by minimizing $$L_\mathrm{MEM} = \chi^2 - \lambda S,\label{eq:LMEM}$$ where, for the case of interferometric data, $\chi^2$ can be calculated as $$\chi^2 = \sum_{k=1}^{N_\mathrm{Vis}}\frac{||V_k^\mathrm{obs} -
    V_k^\mathrm{mod}||^2}{\sigma_k^2},$$ where the sum runs over all the $N_\mathrm{Vis}$ visibilities, the symbol $||z||$ stands for the modulus of the complex number $z$ and $\sigma_k$ is the root mean square (rms) noise of the corresponding visibility. $\lambda$ is a control parameter and the entropy $S$ varies for different implementations . The entropy is used as a regularizing term in a degenerate inverse problem, when there are more free parameters than data. Different formulations for $S$ appear in the literature. Some examples are $\sum_i\ln(I_i)$, $\sum_iI_i\ln(I_i)$, $\sum_i\ln(p_i)$, $\sum_ip_i\ln(p_i)$, where $I_i$ is the specific intensity value at pixel $i$ and $p_i = I_i/\sum_iI_i$ .

used MEM in the AIPS VM task. Their method makes some approximations that diagonalize the Hessian matrix required to optimize their merit function. They used an entropy of the form $S =
-\sum_iI_i\log{(I_i/ m_i)}$, where the sum extends over all the pixels $i$, $\{I_i\}_{i=1}^n$ is the model image and $\{m_i\}_{i=1}^n$ is a prior image. However, the neglect of the side-lobe contribution to the Hessian may lead the optimization to local minima that still bear instrumental artifacts. [@Casassus]implemented a MEM algorithm based on the conjugate gradient method, without the use of the Cornwell and Evans approximation. They used an entropy of the form $S = -\sum_iI_i\log{(I_i/ M)}$, where $\{I_i\}_{i=1}^N$ is the model image and $M$ is a small intensity value, i.e they start with a blank image prior, and $M$ is an intensity value much smaller than the noise. Bayesian analysis is a powerful tool for image reconstruction techniques. In this application, our goal is to find the most probable image by maximizing its *a-posteriori* probability. For Bayesian methods, the *a-priori* and likelihood distributions are needed. To derive the *a-priori* probability the definition of an intensity quantum is needed. This quantum represents the minimum measurable intensity unit. The intensity in each pixel can be interpreted as a number of quanta $I_i = \sigma_\mathrm{q}N_i$, where $I_i$ is the intensity in pixel $i$, $\sigma_\mathrm{q}$ is the quantum size and $N_i$ the number of quanta in pixel $i$.

used Bayesian analysis in the Pixon algorithm. They use a variable model and maximize $P(I,M|D)$, that is, the probability of the image $I$ and model $M$ given the data $D$. In their approach the model used to parameterize the image is a set of Gaussians which are used to average a pseudo-image. The pseudo-image starts as a maximum residual likelihood reconstruction and a local Gaussian pixon is assigned to each of its pixels. The number of pixons, and hence the number of free parameters, is reduced in each iteration.

have used Bayesian analysis for interferometric data, but using a fixed pixel grid to parameterize the model image. They use Gibbs sampling to determine the posterior density distribution.

The most typical model used in astronomy to represent the sky brightness distribution consists of a pixel grid. A big disadvantage of this grid is that the number of pixels remains fixed as well as their size. Often, uniform pixel grids involve more free parameters than really needed to fit the data.

The purpose of this paper is to explore Bayesian reconstruction with image models based on Voronoi tessellations in place of the usual pixelated image. We call this new deconvolution method “Voronoi image reconstruction” (VIR, hereafter). The advantage of using Voronoi models is that it is possible to use a smaller number of free parameters, as required by Bayesian theory. Our purpose is not optimal CPU efficiency; we search for the optimal image and model from a Bayesian point of view.

We used the Cosmic Background Imager [CBI, @pad02] to illustrate our method. The CBI is a planar interferometer array with 13 antennas, each 0.9 m in diameter, mounted on a 6 m tracking platform. An example of CBI baselines is shown in Figure \[fig:baselines\]. The radius of the hole at the center of the $(u,
v)$ plane is the reciprocal of the minimum distance between two antennas, measured in wavelengths. The side-lobes of the CBI are caused mainly by this central hole in the $(u, v)$ baselines.

We briefly summarize the elements of Bayesian theory that determine the probability distributions concerning our problem (Section \[sec:Bayes\]). The new model based on Voronoi tessellations is described (Section \[sec:Voronoi\]), as well as optimization issues involved in our problem (Section \[sec:Optimization\]). We discuss implementation details such as the optimal quantum size and number of Voronoi polygons (Section \[sec:Implementation\]), compare reconstructions made with MEM and VIR (Section \[sec:Results\]) and finally summarize our results (Section \[sec:Conclusions\]).

Bayesian Theory {#sec:Bayes}
===============

An image model is required to parameterize the sky brightness distribution. The most typical model used in astronomy is a rectangular grid of uniform pixels. That configuration of pixels is the model $M$, and the distribution of brightness in the model is called an image $I$. We search for the image that represents as accurately as possible the visibility data $D$. The Bayesian image reconstruction approach, using a fixed model, tries to find the image that maximizes the probability $P(I | D, M)$, i.e. find the most probable image given the data and the model.

Using the Bayes theorem, we obtain $$P(I | D, M) = \frac {P (D| I, M)P(I | M)}{P (D | M)}.$$ Since the data is fixed, $P(D | M)$ is a constant in the problem when the model is not considered as a variable. Thus, the fixed image model optimization problem reduces to $$\max_I P(I | D, M) = \max_I P (D| I, M)P(I | M). \label{eq:P(I|D|M)}$$ The first term, $P (D|I, M)$ is called the likelihood, and measures how well our data represents our image. The second term, $P (I | M)$ is called the image prior, and gives the *a-priori* probability of the image given the model, i.e. how probable is the image given only the model.

In the case of having a variable model, what we would like to find is the image and model that maximize $P(I, M| D)$, i.e. find the most probable image and model given the data. In this case we find $$\begin{aligned}
   P(I, M| D) & = & P(I|D,M)P(M|D)\nonumber\\
   & = & \frac {P (D| I, M)P(I| M)P(M|D)}{P (D | M)}\nonumber\\
   & = & \frac {P (D| I, M)P(I| M)P(M)}{P (D)}.
 \end{aligned}$$ Since the data is fixed, $P(D)$ is constant in our problem. As we cannot privilege one model over another in the absence of image and data, $P(M)$ is the same for all models, so it is not important for our analysis. This way, our optimization problem reduces to $$\max_{I, M} P(I, M| D) = \max_{I, M} P (D| I, M)P(I | M).$$

Probability Distributions
-------------------------

Our data is a set of $N_\mathrm{Vis}$ observed visibilities $\{V_1^\mathrm{obs}, V_2^\mathrm{obs}, \cdots,
V_{N_\mathrm{Vis}}^\mathrm{obs}\}$. If we have a certain model $M$ and image $I$, we obtain model visibilities $\{V_k^\mathrm{mod}\}$ by simulating the interferometric observations over our image: $$V^\mathrm{mod}_k = V^\mathrm{mod}(u_k,v_k) = \int_{-\infty}^{+\infty} A(x,y)
I(x,y)\exp\left[2\pi i (u_kx+v_ky)\right]
\frac{dx\,dy}{\sqrt{1-x^2-y^2}} ~, \label{eq:vmodel}$$ where $\{u_k, v_k\}$ are the coordinates of baseline $k$ in the $(u,
v)$ plane and $A$ is the primary beam. We thus have a set of $N_\mathrm{Vis}$ model visibilities. Assuming that each visibility is independent from the others and Gaussian noise, the likelihood is $$\begin{aligned}
  P (D|I, M) & = &
  P(\{V_k^\mathrm{obs}\}_{k=1}^{N_\mathrm{Vis}}|\{V_k^\mathrm{mod}(I, M)\}_{k=1}^{N_\mathrm{Vis}})
  = \prod_{k=1}^{N_\mathrm{Vis}} P (V_k^\mathrm{obs}|V_k^\mathrm{mod}) \nonumber\\
  & = &
  \prod_{k=1}^{N_\mathrm{Vis}}\frac{1}{2\pi\sigma_k^2}e^{-||V_k^\mathrm{obs}
  - V_k^\mathrm{mod}||^2/2\sigma_k^2}.\end{aligned}$$ To obtain the image prior, $P (I | M)$, we calculate the statistical weight of a given distribution of counts . Consider a model consisting of $n$ cells. In the case of a traditional image, each pixel would be a cell. There is a number of $N$ quanta falling into these cells. These are intensity quanta of some size $\sigma_\mathrm{q}$. In the case of a pixelated image, the intensity in each pixel $i$ would be $I_i = \sigma_\mathrm{q} N_i$, where $I_i$ is the intensity in cell $i$. Each quantum could fall into any of the $n$ cells, so the total number of possible configuration for the $N$ quanta will be $n^N$. The probability of the image given the model is the probability of a certain state $\{N_1, N_2, \cdots, N_n\}$ that represents that image, where $N_i$ is the number of quanta in cell $i$. Consider a given image configuration defined by a particular distribution $\{N_i\}$. The image distribution is not changed in the $N!$ possible redistributions of counts between cells, provided each $N_i$ is constant. The $\prod_i N_i!$ swaps of counts within each cell keep the same image configuration. The model $M$ consists of the Voronoi diagram and the total number of quanta (i.e. n, the position of the generators and N), thus the *a-priori* probability is $$P (I | M) = P (\{N_i\}|n, N) = \frac{N!}{n^N\prod_i N_i!}.
\label{eq:apriori}$$ As explained above, $\sigma_\mathrm{q}$ is an intensity quantum. It is also possible to describe the number of quanta per cell using a flux quantum $\sigma_i^\mathrm{F}$, where $i$ is the index of the cell to which we associate the quantum. This flux quantum can be expressed in terms of the intensity quantum as $\sigma_i^\mathrm{F} =
\sigma_\mathrm{q}A_i$, where $A_i$ is the area of cell $i$. In this case, the number of quanta per cell is $N_i =
F_i/\sigma_i^\mathrm{F}$, where $F_i = I_i A_i$ is the flux of cell $i$. This leads to $N_i = I_i/\sigma_\mathrm{q}$, which is the same expression for $N_i$ obtained using the intensity quantum $\sigma_\mathrm{q}$. Using these cell-dependent flux quanta, the probability of a quantum falling into each cell will be $\frac{1}{n}$ for every cell, leaving the *a-priori* probability the same as Eq. \[eq:apriori\].

MEM and Natural Entropy
-----------------------

In Bayesian theory, for a fixed model, the image $I$ can be found by optimizing the *a-posteriori* probability: $$\begin{aligned}
  \max_I P (I | D, M) & = & \min_I(-\ln{P (D| I, M)P(I | M)})\nonumber\\
  & = & \min_I\sum_{k=1}^{N_\mathrm{Vis}}\frac{||V_k^\mathrm{obs} -
  V_k^\mathrm{mod}||^2}{2\sigma_k^2} - \ln\bigg(\frac{N!}{n^N\prod_i
  N_i!}\bigg)\nonumber\\
  & = & \min_I\frac{1}{2}\chi^2 - S,\label{eq:funcL}\end{aligned}$$ where we have defined the natural entropy $S =
\ln\bigg(\frac{N!}{n^N\prod_i N_i!}\bigg)$. call the term $\ln{(N!/\prod_i N_i!)}$ the multiplicity prior. In the limit of large $N_i$, $$\begin{aligned}
S 
& \simeq & N\ln\frac{N}{n} - \sum_i N_i\ln{N_i},\label{eq:SStirling}\end{aligned}$$ and it can be seen that the Bayesian method is very similar to MEM in the sense that we are adjusting the image to the data while maximizing an entropy of the form of Eq. \[eq:SStirling\]. VIR uses the natural entropy as a regularizing term.

A New Image Model based on Voronoi Diagrams {#sec:Voronoi}
===========================================

A Voronoi diagram is a division of the Euclidian plane into $n$ regions $\mathcal{V}_i$ defined by $n$ points $\vec{x_i}$ (called sites or generators) such that every coordinate $\vec{x}$ in the space belongs to $\mathcal{V}_i$ if and only if $||\vec{x} - \vec{x_i}|| <
||\vec{x} - \vec{x_j}||\ \forall\ j\neq i$. The result of the above definition is a set of polygons defined by the generators. Figure \[fig:Voronoi\] shows an example of a Voronoi diagram. For further details on Voronoi diagrams see [@Voronoi].

We propose a 2D Voronoi diagram in place of the usual pixelated, uniform grid, image as our model. We associate an intensity $I_i$ to each of these polygons. The advantage of using a Voronoi diagram is that we can use just as many cells (i.e. free parameters) as the data requires. Our optimization parameters will be the position of each generator $\vec{x_i} = (x_i, y_i)$, and the intensity at each cell, $I_i$.

With our new model $M$ consisting of $n$ generators ($3\times n$ parameters, $x_i, y_i$ and $I_i$ for each generator), we can vary the number of free parameters as required by the optimization problem. We can see in equation (\[eq:funcL\]) that the entropy $S$ increases as the number of cell $n$ decreases.

Optimization {#sec:Optimization}
============

The optimization problem can be seen as a maximization of the *a-posteriori* probability $\max_{I, M} P (I, M | D)$, or equivalently as a minimization of the more convenient merit function $L = \frac{1}{2}\chi^2 - S$. The conjugate gradient method (CG) is often used for minimization problems where derivatives can be easily calculated. Though it is usually fast in convergence, CG has the problem of converging on local minima depending on the initial condition. The use of other optimization algorithms is postponed to future work.

The CG method searches parameters space using the gradient of the function to be minimized. The derivatives of this function are $$\begin{aligned}
  \frac{\partial L}{\partial x} & = & \frac{1}{2}\frac{\partial \chi^2
  }{\partial x} - \frac{\partial S}{\partial x},\\ 
  \frac{\partial \chi^2}{\partial x} & = & 2\sum_{k=1}^{N_\mathrm{Vis}}
  \frac{1}{\sigma_k^2} \mathrm{Re}\left((V_k^{\mathrm{mod}}
  - V_k^{\mathrm{obs}})^* \frac{\partial
  V_k^{\mathrm{mod}}}{\partial x}\right) \label{dLdx},\end{aligned}$$ where $x$ is any of the optimization parameters ($x_i$, $y_i$ or $I_i$). The derivatives of the visibilities with respect to the position $\vec{x}_i = (x_i, y_i)$ of the $i$ generator are $$\begin{aligned}
  \frac{\partial V_k^{\mathrm{mod}}}{\partial x_i} & = & 
  \sum_{j\in J_i} \bigg[(I_i - I_j)
    \sum_{l|\mathrm{pixel~}l\epsilon a_{ij}}A_l\Delta t_l
    (M_xt_l + b_x)e^{(t_lc_2 + s_0c_1)}\bigg],\\ 
  \frac{\partial V_k^{\mathrm{mod}}}{\partial y_i} & = & 
  \sum_{j\in J_i} \bigg[(I_i - I_j)
    \sum_{l|\mathrm{pixel~} l\epsilon a_{ij}}A_l\Delta t_l
    (M_yt_l + b_y)e^{(t_lc_2 + s_0c_1)}\bigg],\end{aligned}$$ where $I_i$ is the intensity in cell $i$, $J_i$ is a set of the indices of the polygons adjacent to $\mathcal{V}_i$, $a_{ij}$ is the edge which divides polygons $\mathcal{V}_i$ and $\mathcal{V}_j$, $l$ sums over the pixels which intersect $a_{ij}$, $A$ is the CBI primary beam. For further details see Sec. \[ap:derivatives\].

The derivative of the visibilities with respect to the intensity of each cell $I_i$ is $$\frac{\partial V_k^{\mathrm{mod}}}{\partial I_i} =
  \frac{\sin{(\pi u_k\Delta x)}\sin{(\pi v_k\Delta y)}}{\pi^2u_kv_k}
  \sum_{\textrm{\scriptsize{pixels }}l\epsilon \mathcal{V}_i}A_l
  e^{2\pi i(u_kx_l+v_ky_l)} \label{eq:dVdI},$$ where $\vec{k}_k = (u_k, v_k)$ is the baseline corresponding to the pair of antennas $k$, $\Delta x$ and $\Delta y$ are the pixel width and height, and the sum extends over all the pixels inside $\mathcal{V}_i$.

The entropy only depends of the intensities $I_i$, so $\frac{\partial
S}{\partial x_i} = \frac{\partial S}{\partial y_i} = 0$, then (see Sec. \[ap:derivativesdS\]) $$\frac{\partial S}{\partial I_i} = \frac{1}{\sigma_\mathrm{q}}
  (\sum_{k=1}^{N}\frac{1}{k} - \ln{n} - \sum_{k=1}^{N_i}\frac{1}{k}).$$

VIR Design and Implementation {#sec:Implementation}
=============================

We have designed, and implemented in c++, VIR with 6 modules which include algorithms for:

-   the generation of the Voronoi diagram

-   calculation of model visibilities

-   calculation of the merit function $L$ to be optimized as well as its derivatives

-   fitting a Voronoi diagram to an image

-   the CG method

-   the optimization of the number of polygons

VIR uses the CG method from [@NumericalRecipes] and searches for the position and intensities of the Voronoi polygons, ${x_i, y_i,
I_i}$, that minimize our merit function $L$. The CG method modifies the intensities and also moves the positions of the Voronoi generators. This causes the shape of the Voronoi polygons to change as well. A general problem with CG is that it usually converges on local minima. For VIR in particular, though Voronoi polygons intensities adjust quite fine, the positions of the generators are difficult to modify substantially. The VIR parameter space is smooth enough in intensity space to converge to a good solution. But the parameter space in cell generator positions is very structured, and CG is quickly stuck on local minima.

Due to the fact that CG easily falls into local minima, we needed a good approximation for the initial Voronoi diagram. For this purpose we used a pixelated version of the Bayesian algorithm, where the model was a uniform grid. We decided to do a pure $\chi^2$ (maximum likelihood, ML) reconstruction and use the fifth CG iteration as our starting image. We chose this particular iteration because on inspection the modeled images were still smooth. Pure $\chi^2$ reaches convergence with noisy images, where the true image is unrecognizable. We then fitted a Voronoi diagram to the image (see Sec. \[ap:fitting\]) and ran CG using the positions and intensities of the generators as our free parameters, which led to our final reconstruction. Truncation to a level of $10^{-5}$ quanta was used to enforce positivity.

An important issue to consider is the size of the quantum $\sigma_\mathrm{q}$. treat $\sigma_\mathrm{q}$ as a free parameter. But, as we now explain, $\sigma_\mathrm{q}$ was held constant in this implementation of VIR. We treat the number of quanta per cell as a continuous variable in order to use the CG method. Entropy is maximized at $\sigma_\mathrm{q} = \infty$, where, for a given configuration of intensities $\{I_i\}$, $N = 0$ and $S = 0$. For every other value of $N$, the entropy will be negative. This means that even for large $\sigma_\mathrm{q}$, the intensities $I_i =
\sigma_\mathrm{q}N_i$ can have reasonable values (using small $N_i$). Figure \[fig:SvsN\] shows $S$ as a function of $N$ for $51$ Voronoi generators and 3 different intensity distributions using the model tessellation of Figure \[fig:reconstruccion\]c. We considered: 1- the VIR intensities of Figure \[fig:reconstruccion\]c, 2- a uniform intensity distribution image ($N_i = \frac{N}{n}$ $\forall$ $i$), 3- a spike where all $N$ are only in one cell ($N_i = N$, $N_j
=0$ $\forall$ $j\neq i$). The curves of Figure \[fig:SvsN\] are obtained by keeping the intensities fixed and modifying $\sigma_\mathrm{q}$ in order to obtain different $N$. It can be seen on Figure \[fig:SvsN\] that the entropy is maximized at $N = 0$, independently of the intensities $\{I_i\}$ of the model, where the optimal value of $\sigma_\mathrm{q} = \infty$ is achieved if the number of quanta per cell is treated as a continuous variable. If the number of quanta per cell were discrete variables, as in , the choice of a big $\sigma_\mathrm{q}$ would admit only zero values for every cell. Otherwise, if one or more quanta fell in a given cell, the intensity of that cell would diverge as $\sigma_\mathrm{q}$ for arbitrarily large $\sigma_\mathrm{q}$, causing a big $\chi^2$ value. Therefore, in our continuous optimization the intensity quantum must be determined a-priori.

In the Bayesian description of the entropy we count events that fall in each cell. It seems reasonable to take the noise level as the minimum value of intensity we can distinguish. So, $\sigma_\mathrm{q}$ should approximate the estimated thermal noise in the naturally weighted dirty map. The definition of the weighted dirty map [e.g. @Briggs] is $$\begin{aligned}
  I^\mathrm{D} (x, y) \equiv
  \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}W(u, v)V(u, v)e^{-2\pi
  i(ux + vy)}dudv,\\
  W(u, v) = \frac{1}{\sum_kw_k} \sum_kw_k\delta(u-u_k, v-v_k),\end{aligned}$$ where the sums extend over all visibilities, $w_k$ are the weights given to visibility $k$ and $\delta$ is the two-dimensional Dirac delta function. Propagating the thermal noise, we get for the standard deviation of the dirty map $$\sigma_\mathrm{rms}^\mathrm{D} =
\sqrt{\frac{\sum_kw_k^2\sigma_k^2}{(\sum_kw_k)^2}}, \label{eq:sigmaD}$$ where $\sigma_k$ are the visibilities standard deviations. To take into account model pixels correlated by the interferometer beam, we should multiply the previous expression by $\sqrt{N_\mathrm{beam}}$, where $N_\mathrm{beam}$ is the number of pixels inside a beam pattern. This leads to $$\sigma_\mathrm{rms} = \sqrt{\frac{\sum_kw_k^2\sigma_k^2}
  {(\sum_kw_k)^2}}\sqrt{N_\mathrm{beam}}.$$ For natural weights, $\sigma_k^2 = \frac{1}{w_k}$, $$\sigma_\mathrm{rms} = \sqrt{\frac{N_\mathrm{beam}}{\sum_kw_k}} =
  \sqrt{\frac{N_\mathrm{beam}}{\sum_k\frac{1}{\sigma_k^2}}}.$$ We calculated the noise with natural weighting, $w_k =
\frac{1}{\sigma_k^2}$, because this is the weight we give to each individual visibility data in the optimization of the merit function.

Once we have the value of $\sigma_\mathrm{q}$ we search for the optimal number of cells $n$. In Figure \[fig:Lvsn\] we plot the optimal merit function for different $n$ and $\sigma_\mathrm{q}$. These reconstructions were made over a simulation of CBI observations on a mock sky image (Figure \[fig:reconstruccion\]a). We averaged over 100 reconstructions with different realizations of Gaussian noise. The average curves shown in Figure \[fig:Lvsn\], start with $n = 10$ and end with $n = 100$ for even $n$. One single reconstruction for all $n$ took about two hours using an AMD Athlon64 XP3000 processor with 1GB of DDR RAM at 333 MHz, so the 300 reconstructions took about $25$ CPU days, but we distributed the work in six computers, so it took about $5$ real days in total. It can be seen that for a signal to noise ratio (SNR) of $\sim 52$, on average, the optimal number of polygons $n$ is between 50 and 55. When $\sigma_\mathrm{q}$ is diminished to $\frac{1}{10}\sigma_\mathrm{q}$, on average, the optimal merit function is found at $n$ close to $30$. For $\sigma_\mathrm{q} = 10\sigma_\mathrm{q}$, the optimal $n$ is found between $80$ and $90$. It can be seen that as we increase the value of $\sigma_\mathrm{q}$ we reach lower values for our function, as discussed above. Furthermore, the optimal number of polygons increases.

Example Reconstruction {#sec:Results}
======================

Mock Dataset
------------

The mock sky image we used for simulations is a $256\times 256$ image consisting of three Gaussians and a rectangle. Figure \[fig:reconstruccion\]a shows this image on a $128\times128$ pixel field. Pixels are $0.75'\times0.75'$, while the CBI’s primary beam is of $45'$ FWHM (60 pixels), so most of the emission lies under the beam. We simulated a CBI observation of 3620 visibilities over this image and added Gaussian noise to the visibilities in order to reach a SNR of $\sim52$. This SNR was calculated by taking the maximum intensity from the dirty map using natural weights, and using the noise $\sigma_\mathrm{rms}^\mathrm{D}$ (see Eq. \[eq:sigmaD\]). Simulation of the CBI observations is performed with the MockCBI program (Pearson 2000, private communication), which calculates the visibilities $V(u,v)$ on the input images $I(x,y)$ with the same $uv$ sampling as a reference visibility dataset (Eq. \[eq:vmodel\]). Thus MockCBI creates the visibility dataset that would have been obtained had the sky emission followed the true image. Figure \[fig:reconstruccion\]b shows the dirty map calculated over these simulated visibilities using the DIFMAP package [@she97]. The CBI’s primary beam is drawn as a dashed circle. The secondary side-lobes due to the central discontinuity in $u$-$v$ coverage can be distinguished in Figure \[fig:reconstruccion\]b at a level comparable to the true emission.

MEM Reconstruction
------------------

The VIR method was compared with the MEM algorithm described in [@Casassus]. To fit the model image to the observed visibilities, MEM calculates the model visibilities required by its merit function $L_\mathrm{MEM}$. The model visibilities are those obtained by a simulation of CBI observations had the sky followed the model image . The free-parameters of our MEM model are the pixels in the model $64\times64$ image. The model functional we minimize is $L_\mathrm{MEM} = \chi^2 - \lambda S$, with the entropy $S= - \sum_i
I_i \log I_i/M$, where $M$ is a default pixel value well below the noise level, and $\{I_i\}_{i=1}^{N}$ is the model image. We started with the fifth iteration of a pure $\chi^2$ reconstruction ($\lambda =
0$) as initial condition for the CG minimization. This is the same ML initial condition used in our VIR method. Figure \[fig:reconstruccion\]g shows the reconstructed image using $\lambda
= \frac{100}{\sigma_\mathrm{rms}}$ and $M =
10^{-2}\sigma_\mathrm{rms}$ inset on a larger $128\times128$ image [^1].

VIR Reconstruction
------------------

The MEM algorithm described above requires the prior assignment of the $\lambda$ and $M$ parameters as well as the entropy formula. In contrast, our VIR algorithm is free from such arbitrary parameters (provided the optimal $\sigma_\mathrm{q}$ is indeed equal to $\sigma_\mathrm{rms}$). For our VIR method, we only need to find the number of polygons to be used. In order to find the optimal number of polygons we reconstructed with different numbers of generators in a range covering each natural number from $n = 6$ to $n = 100$. We found a minimum at $n = 51$. Figure \[fig:Lvsn\] summarizes this search. The whole search for a particular realization of noise took about 10 hours on the AMD Athlon64 XP3000 processor with 1GB of DDR RAM at 333 MHz. The VIR reconstruction using 51 polygons is shown in Figure \[fig:reconstruccion\]c, where the Voronoi cells have also been drawn. Figure \[fig:reconstruccion\]d shows the same model but without drawing the Voronoi mesh.

Results
-------

The quality of each reconstruction can be assessed by visual inspection, comparing the VIR and MEM model images with the true image. The MEM model looks similar to the true image but is noisy. The density of Voronoi generators in the VIR model is greater where there is more emission in the true image, approximating the true image with only a few polygons. We calculated $\chi^2_\mathrm{im} =
\sum_i(I_i^\mathrm{mod} - I_i^\mathrm{true})^2$, where $I_i^\mathrm{mod}$ is the intensity at pixel $i$ of the model image (MEM or VIR), $I_i^\mathrm{true}$ is the intensity at pixel $i$ of the true image, and the sum extends over all pixels in the images. $\chi^2_\mathrm{im}$ gives a measure of how well the model fits the true image. It can be seen in Table \[table:funciones\] that the VIR reconstruction has a better $\chi^2_\mathrm{im}$ than MEM, showing that the VIR model is closer to the true image than the MEM model.

[crrrrr]{} MEM & 7354.85 & 1.016 & 12192.6 & 0.001608\
VIR & 7221.04 & 0.997 & 3753.28 & 0.001396\

Figures \[fig:reconstruccion\]e and \[fig:reconstruccion\]h show the VIR and MEM models residuals. Residual images are the dirty map of the residuals of the visibilities, calculated over the optimal model visibilities. It can be noted on Figure \[fig:reconstruccion\]e that the VIR residuals are very good, showing only noise. On the other hand, in the MEM residuals (Figure \[fig:reconstruccion\]h) the object shape can clearly be distinguished as well as the CBI’s side-lobes. The object seems to be more compact in the model than in its MEM residuals; as expected these residuals are convolved with the synthetic beam.

Restored images are shown in Figures \[fig:reconstruccion\]f and \[fig:reconstruccion\]i. These images are obtained by convolving the models with a Gaussian point spread function (PSF) given by DIFMAP and adding the dirty map of the residuals visibilities. On Figures \[fig:reconstruccion\]f and \[fig:reconstruccion\]i it can be assessed that VIR produces improved restored images relative to MEM. The VIR restored image is similar to that expected given the instrumental noise: it approximates the true image convolved with a Gaussian PSF plus a uniform noise level. In the MEM restored image, on the other hand, the CBI side-lobes can still be distinguished.

The number of optimization parameters in MEM are $64\times64 = 4096$, while the VIR method has only $51$ triplets (cell’s $(x, y)$ position and intensity) i.e. $153$ free parameters. This smaller number of parameters causes the Bayesian entropy to be greater than the pixelated version, obtaining a smaller value for our merit function $L$ to be minimized.

Table \[table:funciones\] also shows $\frac{\chi^2}{n_\mathrm{data}}$ values, where $n_\mathrm{data}$ is the number of data points ($3620\times 2$ in our case). A good reconstruction should have a $\frac{\chi^2}{n_\mathrm{data}}$ value close to $1$. It can be seen that the VIR model gives a value of $\frac{\chi^2}{n_\mathrm{data}}$ closer to $1$ than the MEM reconstruction.

Conclusions {#sec:Conclusions}
===========

We have introduced a Bayesian Voronoi image reconstruction (VIR) technique for interferometric data where the image is represented by a Voronoi tessellation in place of the usual pixelated image. The advantage of Voronoi models is that we can use a smaller number of free parameters, as required by the Bayesian analysis of a discretized intensity field. Our purpose is not optimal CPU efficiency; we search for the optimal image and model from a Bayesian point of view. The free parameters of our model are the Voronoi generators positions $(x_i, y_i)$ and intensities $I_i$. The following points summarize our work:

-   We discretized the intensity field in order to calculate *a priori* probabilities. We defined a quantum intensity value $\sigma_\mathrm{q}$ such that $I_i = \sigma_\mathrm{q} N_i$, where $I_i$ is the intensity at cell $i$ and $N_i$ the number of quanta in cell $i$.

-   We calculated the analytical derivatives required by the conjugate gradient and cross checked them by finite differences. Because the parameter space in cell generators positions is very structured, the positions of the Voronoi generators are difficult to change. As initial condition we took a Voronoi tessellation adjusted to an interrupted maximum likelihood reconstruction.

-   We simulated a CBI observation over a true image and reconstructed sky images from this mock visibility dataset using MEM and VIR.

-   We defined the value of $\sigma_\mathrm{q}$ as the estimated noise of the dirty map and searched for the optimal number of Voronoi polygons for our example dataset.

-   We finally compared the MEM and VIR models, residuals and restored images. The VIR model is closer to our true image than the MEM model. Residuals and restored images are also better in VIR than in MEM. We found that VIR model visibilities give a better fit to the data than MEM, in the sense that $\chi^2$ is closer to its expected value.

We are grateful to Tim Pearson for advice on FFTs and the use of MOCKCBI. G.F.C. and S.C. acknowledge support from FONDECYT grant 1060827, and from the Chilean Center for Astrophysics FONDAP 15010003.

Derivatives {#ap:derivatives}
===========

Our merit function for minimization is $$\begin{aligned}
L & = &
\frac{1}{2}\sum_{j=1}^{N_\mathrm{Vis}}\frac{||V_j^{\mathrm{mod}}
  - V_j^{\mathrm{obs}}||^2}{\sigma_j^2}
-\ln\left(\frac{N!}{n^N \prod_{i=1}^{n}N_{i}!}\right)\\ & = &
\frac{1}{2}\chi^2 - S.\end{aligned}$$ So, the derivative of $L$ with respect to any variable $x$ is $$\frac{\partial L}{\partial x} = \frac{1}{2}\frac{\partial
\chi^2}{\partial x} - \frac{\partial S}{\partial x}$$

Calculation of the Derivatives of $\chi^2$
------------------------------------------

$\chi^2$ derivatives with respect to any variable $x$ can be obtain as follows $$\begin{aligned}
\frac{\partial}{\partial x} \frac{1}{2}\chi^2
  & = & \frac{\partial}{\partial x}
    \left(\frac{1}{2}\sum_{k=1}^{N_\mathrm{Vis}}\frac{||V_k^{\mathrm{mod}} -
    V_k^{\mathrm{obs}}||^2}{\sigma_k^2}\right)\nonumber\\
  & = & \sum_{k=1}^{N_\mathrm{Vis}} \frac{1}{\sigma_k^2}
    \left(\mathrm{Re}(V_k^{\mathrm{mod}} - V_k^{\mathrm{obs}})
    \mathrm{Re}\left(\frac{\partial V_k^{\mathrm{mod}}}{\partial x}\right)
        + \mathrm{Im}(V_k^{\mathrm{mod}} - V_k^{\mathrm{obs}})
    \mathrm{Im}\left(\frac{\partial V_k^{\mathrm{mod}}}{\partial x}\right)
    \right),\nonumber\\
  && \label{eq:dLdI}\end{aligned}$$ where its necessary to calculate the model visibilities derivatives with respect to $x$.

### Calculation of $\frac{\partial V_k^{\mathrm{mod}}}{\partial I_i}$

In our Voronoi tessellation representation of the sky image $$\label{eq:VisDiscreta}
V(\vec{k}) = \sum_i^{N_\mathrm{V}} I_i\int_{\mathcal{V}_i} A(\vec{x}) e^{2\pi
i\vec{k}\vec{x}}d\vec{x},$$ where $N_\mathrm{V}$ is the number of polygons, $\mathcal{V}_i$ is polygon $i$ and $I_i$ its intensity. We neglected the $\sqrt{1 - x^2 -
y^2}$ term which is close to $1$, but it can easily be included in $A(\vec{x})$. After derivation and defining $f_k(\vec{x})\equiv
A(\vec{x}) e^{2\pi i\vec{k_k}\vec{x}}$ we obtain $$\begin{aligned}
\frac{\partial V_k^{\mathrm{mod}}}{\partial I_i} & = &
\int\int_{\mathcal{V}_i} f_k(\vec{x})d^2x, \\
& = & \frac{\sin{(\pi u_k\Delta x)}\sin{(\pi v_k\Delta y)}}{\pi^2u_kv_k}
  \sum_{\textrm{\scriptsize{pixels }}l\epsilon \mathcal{V}_i}A_l
  e^{2\pi i(u_kx_l+v_ky_l)}, \\
& \simeq & \Delta x\Delta y
  \sum_{\textrm{\scriptsize{pixels }}l\epsilon \mathcal{V}_i}A_l
  e^{2\pi i(u_kx_l+v_ky_l)}\end{aligned}$$ for small $\Delta x$ and $\Delta y$.

### Calculation of $\frac{\partial V_k^{\mathrm{mod}}}{\partial x_i}$ and $\frac{\partial V_k^{\mathrm{mod}}}{\partial y_i}$

To evaluate $\frac{\partial V_k}{\partial x_i}$ we move the generator $\vec{x}_i$ an infinitesimal quantity $\delta_x$ parallel to the $\hat{x}$ axis as in Figure \[fig:VoronoiDelta\]. We will calculate $$\frac{\partial V_k}{\partial x_i} =
\lim_{\delta_x\to 0}\frac{\Delta V}{\delta_x},\label{eq:DerivadaLimite}$$ where $\Delta V_k = V_k(\vec{x}_1, \cdots,
\vec{x}_i + \vec{\delta_x}, \cdots, \vec{x}_{N_V}) -
V_k(\vec{x}_1, \cdots, \vec{x}_i, \cdots, \vec{x}_{N_V})$.

It can be seen in Figure \[fig:VoronoiDelta\] that when moving the generator $\vec{x}_i$, the only polygons modified are $\mathcal{V}_i$ and its neighbors. Using this, Eq. \[eq:VisDiscreta\] leads to $$\begin{aligned}
\Delta V_k & = & I_i'\int_{\mathcal{V}_i'}f_k(\vec{x})d\vec{x} -
I_i\int_{\mathcal{V}_i}f_k(\vec{x})d\vec{x}\nonumber\\
& & + \sum_{j\in J_i}\bigg(I_j'\int_{\mathcal{V}_j'}f_k(\vec{x})d\vec{x} -
I_j\int_{\mathcal{V}_j}f_k(\vec{x})d\vec{x}\bigg),\label{eq:deltaVis}\end{aligned}$$ where $\mathcal{V}_i$ is the polygon generated by $\vec{x}_i$ before moving, $\mathcal{V}_i'$ is the same polygon after moving $\vec{x}_i$, $J_i$ is the set of indices of the polygons that are neighbors to $\mathcal{V}_i$ and $J_i'$ is the set of indices of the polygons that are neighbors to $\mathcal{V}_i'$.

It can be seen in Figure \[fig:VoronoiDelta\] that $$\begin{aligned}
\mathcal{V}_i = (\mathcal{V}_i\cap \mathcal{V}_i')\cup(\mathcal{V}_i\setminus \mathcal{V}_i\cap \mathcal{V}_i'), & \mathcal{V}_i' =
(\mathcal{V}_i\cap \mathcal{V}_i')\cup(\mathcal{V}_i'\setminus \mathcal{V}_i\cap \mathcal{V}_i'),\\
\mathcal{V}_j = (\mathcal{V}_j\cap \mathcal{V}_j')\cup(\mathcal{V}_i'\cap \mathcal{V}_j), & \mathcal{V}_j' = (\mathcal{V}_j\cap
\mathcal{V}_j')\cup(\mathcal{V}_i\cap \mathcal{V}_j'),\end{aligned}$$ so, Eq. \[eq:deltaVis\] is $$\begin{aligned}
\Delta V_k & = & (I_i'-I_i)\int_{\mathcal{V}_i'\cap
\mathcal{V}_i}f_k(\vec{x})d\vec{x}\nonumber\\
&& + \sum_{j\in J_i}\bigg[(I_i' - I_j)\int_{\mathcal{V}_i'\cap
\mathcal{V}_j}f_k(\vec{x})d\vec{x} + (I_j' - I_i)\int_{\mathcal{V}_i\cap
\mathcal{V}_j'}f_k(\vec{x})d\vec{x}\bigg].\end{aligned}$$ In our case the cells’ intensities don’t depend of the position of the generators, so we obtain $$\label{eq:DeltaVis}
\Delta V_k = \sum_{j\in J_i}\bigg[(I_i-I_j)\bigg(\int_{\mathcal{V}_i'\cap
\mathcal{V}_j}f_k(\vec{x})d\vec{x} - \int_{\mathcal{V}_i\cap
\mathcal{V}_j'}f_k(\vec{x})d\vec{x}\bigg)\bigg].$$

It can be seen in Figure \[fig:VoronoiDelta\] that to obtain $\Delta
V_k$ we must integrate only over the shaded regions. For this purpose, for each region between $\vec{x}_i$ and $\vec{x}_j$ we will define a coordinate system $$\begin{aligned}
\hat{s} = -\cos{\alpha_j}\hat{x} + \sin{\alpha_j}\hat{y}. & \hat{t}
= \sin{\alpha_j}\hat{x} + \cos{\alpha_j}\hat{y},\label{eq:sistema_coordenadas}\end{aligned}$$ where $\alpha_j$ is the angle formed by the $-\hat{x}$ axis and the edge $a_{ij}$ between $\vec{x}_i$ and $\vec{x}_j$ (see Figure \[fig:VoronoiCoordenadas\]). Using this change of coordinates, the integral over the region of interest is $$\int_{\mathcal{V}_i'\cap \mathcal{V}_j}f_k(x, y)dxdy = \int_{\mathcal{V}_i'\cap
\mathcal{V}_j}f_k(s,t)dsdt.\label{eq:int}$$

Let $\vec{x_i} = (x_i, y_i)$ be the position of the $i$ cell’s generator, $\vec{x_j} = (x_j, y_j)$ one of its neighbor, and $\vec{x_i}' = (x_i + \delta_x, y_i)$ the site’s position after moving it a quantity $\delta_x$. We define $\vec{x_0} \equiv (x_0, y_0)$ as the point in the intersection of the segment formed by $\vec{x_i}$ and $\vec{x_j}$ and its respective edge $a_{ij}$. The same way, we define $\vec{x_0}' = (x_0', y_0')$ as the point in the intersection of the segment formed by $\vec{x_i}'$ and $\vec{x_j}$ and its respective edge $a_{ij}'$. It can be seen on Figure \[fig:VoronoiCoordenadas\] that $x_0 = \frac{x_i + x_j}{2}$ , $x_0' = x_0 + \frac{\delta}{2}$ and $y_0'
= y_0 = \frac{y_i + y_j}{2}$.

The edge $a_{ij}$ is defined in the new coordinate system by $$s = s_0 = -x_0\cos{\alpha_j} + y_0\sin{\alpha_j}.$$ In the same way, the edge $a_{ij}'$ is defined in the original coordinate system by $$y = m(x - x_0') + y_0,$$ where $$m \equiv \frac{x_i + \delta_x - x_j}{y_j - y_i}.$$ We can define the same line in our new coordinate system as $$s = m't+b',$$ where $$\begin{aligned}
m' & \equiv & -\frac{\cos{\alpha_j} + m\sin{\alpha_j}}{\sin{\alpha_j}
  - m\cos{\alpha_j}},\\ 
b' & \equiv & \frac{-mx_0' +
  y_0}{\sin{\alpha_j} - m\cos{\alpha_j}}.\end{aligned}$$ This can be approximated to first order in $\delta_x$ as $$\begin{aligned}
m' & \simeq & \delta_xM_x,\\
b' & \simeq &s_0 + \delta_xB_x,\end{aligned}$$ where $$\begin{aligned}
M_x & \equiv &\frac{\sin^2\alpha_j}{y_j - y_i} =
\frac{\sin{\alpha_j}\cos{\alpha_j}}{x_i - x_j},\\ 
B_x & \equiv &
\frac{\sin{\alpha_j}}{y_j-y_i}(s_0\cos{\alpha_j} + x_i) =
\frac{\cos{\alpha_j}}{x_i-x_j}(s_0\cos{\alpha_j} + x_i).\end{aligned}$$

The integral in Eq. \[eq:DeltaVis\] using our new coordinate system will be $$\begin{aligned}
\mathcal{I} & = & \int_{\mathcal{V}_i'\cap \mathcal{V}_j}f_k(\vec{x})d\vec{x} - \int_{\mathcal{V}_i\cap
\mathcal{V}_j'}f_k(\vec{x})d\vec{x}\\
& = & \int\int_{a_{ij}}^{a_{ij}'}A(\vec{x})e^{2\pi i(ux + vy)}dxdy.\label{eq:Ixy}\end{aligned}$$ If we use $A(\vec{x})$ in the $(s, t)$ coordinate system as a pixelated image, Eq. \[eq:Ixy\] will be $$\mathcal{I} = \sum_{l\ \epsilon \mathrm{\ pixeles\ de\ }
a_i}A_l\int_{t_{ijl}^1}^{t_{ijl}^2}\int_{s_0}^{m't+b'} e^{2\pi i(ux(s,
t) + vy(s, t))}dsdt,$$ where $t_{ijl}^1$ and $t_{ijl}^2$ are the $t$ coordinate of the beginning and end of the portion of the edge $a_{ij}$ that intersects pixel $l$. Developing the previous expression, $$\begin{aligned}
\mathcal{I} & = & \sum_{l}A_l\int_{t_{ijl}^1}^{t_{ijl}^2}\int_{s_0}^{m't+b'}
e^{2\pi i(u(-s\cos{\alpha_j} + t \sin{\alpha_j}) + v(s\sin{\alpha_j} +
t \cos{\alpha_j}))}dsdt,\label{eq:integralPixel}\\ 
& \simeq & 
\sum_{l}\frac{A_l}{\pi c_2}e^{2\pi i(s_0c_1 + \bar{t}_{ijl}c_2)}\kappa_{ijl}\delta_x,
\label{eq:IExacto}\end{aligned}$$ where we defined $$\begin{aligned}
  c_1 & \equiv & -u\cos{\alpha_j} + v\sin{\alpha_j},\\ 
  c_2 & \equiv & u\sin{\alpha_j} + v\cos{\alpha_j},\\ 
  \kappa_{ijl} & \equiv &
  (M_x\bar{t}_{ijl} + B_x)\sin{(\pi c_2\Delta t_{ijl})} \\\nonumber 
  &&+ i\frac{M_x}{2} \bigg(\frac{\sin(\pi c_2\Delta t_{ijl})}{\pi c_2} -
  \Delta t_{ijl}\cos{(\pi c_2\Delta t_{ijl})}\bigg),\\ 
  \bar{t}_{ijl} & \equiv & \frac{t_{ijl}^1 + t_{ijl}^2}{2},\\ 
  \Delta t_{ijl} & \equiv & \frac{t_{ijl}^2 - t_{ijl}^1}{2}.\end{aligned}$$

In the calculation above we integrated over the fraction of the edge that falls inside pixel $l$ and then summed these integrals over the whole edge $a_i$. It is also possible to approximate the integral of Eq. \[eq:integralPixel\] as $\int_{t_{ijl}^1}^{t_{ijl}^2}g(t) dt =
g(\bar{t}_{ijl})\Delta t_{ijl}$, which is equivalent to taking the limit over the integral $\mathcal{I}$ of Eq. \[eq:IExacto\], $\lim_{\Delta t_{ijl}\rightarrow 0} \mathcal{I}$, obtaining $$\mathcal{I} = \sum_{l} A_l\Delta t_{ijl} (M_x\bar{t}_{ijl} + B_x)e^{2\pi
    i(\bar{t}_{ijl}c_2 + s_0c_1)}\delta_x.\label{eq:IAprox}$$

We found by direct evaluation that the difference between Eq. \[eq:IAprox\] and Eq. \[eq:IExacto\] is negligible, so, for simplicity, we will use Eq. \[eq:IAprox\]. Introducing Eq. \[eq:IAprox\] in Eq. \[eq:DeltaVis\], we obtain $$\Delta V_k = \delta_x\sum_{j\in
    J_i}\bigg[(I_i-I_j)\sum_{l} A_l\Delta t_{ijl} (M_x\bar{t}_{ijl} +
    B_x)e^{2\pi i(\bar{t}_{ijl}c_2 + s_0c_1)}\bigg],$$ so, according to Eq. \[eq:DerivadaLimite\], the derivative of the $k$ visibility with respect to the position $x$ of polygon $i$ is $$\begin{aligned}
  \frac{\partial V_k}{\partial x_i} & = &
  \lim_{\delta_x\to 0}\frac{\Delta V}{\delta_x},\\
  & = & \sum_{j\in
    J_i}\bigg[(I_i-I_j)\sum_{l} A_l\Delta t_{ijl} (M_x\bar{t}_{ijl} +
    B_x)e^{2\pi i(\bar{t}_{ijl}c_2 + s_0c_1)}\bigg].\end{aligned}$$ Similarly, for the derivative with respect to the position $y$ of the $i$ polygon we obtain $$\frac{\partial V_k}{\partial y_i} = \sum_{j\in
    J_i}\bigg[(I_i-I_j)\sum_{l} A_l\Delta t_{ijl} (M_y\bar{t}_{ijl} +
    B_y)e^{2\pi i(\bar{t}_{ijl}c_2 + s_0c_1)}\bigg],$$ where $$\begin{aligned}
M_y & \equiv &\frac{\cos^2\alpha_j}{x_i - x_j} =
\frac{\sin{\alpha_j}\cos{\alpha_j}}{y_j - y_i},\\ 
B_y & \equiv & \frac{\sin{\alpha_j}}{y_j-y_i}(s_0\sin{\alpha_j} - y_i)
= \frac{\cos{\alpha_j}}{x_i-x_j}(s_0\sin{\alpha_j} - y_i).\end{aligned}$$

Calculation of the Derivatives of $S$ {#ap:derivativesdS}
-------------------------------------

We defined our entropy as $$\begin{aligned}
S & = & \ln\left(\frac{N!}{n^N
  \prod_{i=1}^{n}N_{i}!}\right) \\
  & = & \ln(N!) - N\ln(n) - \sum_{i=1}^{n}\ln(N_i!)\\
  & = & \ln\Big(\Gamma(N + 1)\Big) - N\ln(n) - \sum_{i=1}^{n}\ln\Big(\Gamma(N_i +
  1)\Big),\end{aligned}$$ where $N_i = \frac{I_i}{\sigma_\mathrm{q}}$ is the number of quanta in cell $i$, $N = \sum_iN_i$ and $\Gamma$ is the Gamma function. It can be seen that this function does not depend on the position of the Voronoi generators, so $$\frac{\partial S}{\partial x_i} = \frac{\partial S}{\partial y_i} = 0.$$

Using Weierstrass’ definition of the Gamma function $$\Gamma(z) = z^{-1}e^{-\gamma z}\prod_{n=1}^\infty
  \left[\left(1 + \frac{z}{n}\right)^{-1} e^{z/n}\right],$$ where $\gamma$ is Euler’s constant, we can obtain $$\frac{\partial \ln \Big(\Gamma(z + 1)\Big)}{\partial z} = -\gamma +
  \sum_{n=1}^z\frac{1}{n}$$ so, the derivative of $S$ with respect to $I_i$ is $$\frac{\partial S}{\partial I_i} = \frac{1}{\sigma_\mathrm{q}}
  (\sum_{k=1}^{N}\frac{1}{k} - \ln{n} - \sum_{k=1}^{N_i}\frac{1}{k}).$$

Finite Difference Cross Check on the Derivatives
------------------------------------------------

Numerical calculation of the derivatives by finite differences is not very accurate, in particular for the position of the generators. Finite difference derivatives are calculated as $\frac{\partial L}{\partial x} = \frac{ L(x + \delta) -
L(x)}{\delta}$, where $\delta$ is a small displacement of $x$. In the case of the positions of the generators, if $\delta$ is too small, the pixelization of the Voronoi diagram (needed to obtain the model visibilities) will not change after the displacement $\delta$. On the other hand, if $\delta$ is too big, the generator displacement may cause the function to change abruptly, as explained below. That is why we calculated the analytical expression for the derivatives.

To verify that our derivatives are correctly calculated and programmed, we compared our analytical result with a numerical calculation. We created a Voronoi tessellation of $50$ polygons with random positions and intensities and calculated the analytical and numerical derivatives using these parameters $\{x_i, y_i, I_i\}$. For the case of $\frac{\partial L}{\partial x_i}$ and $\frac{\partial
L}{\partial y_i}$ this numerical cross check consists of moving each Voronoi generator a quantity $\delta$ from -0.1 to 0.1 with an interval of $10^{-3}$ in units of the total size of the square image. We evaluate the merit function $L$ at each position intervals, thus obtaining two sequences $\{L_i\}_{i=1}^{2\times 10^2}$. We then fitted a polynomial of order four to the curve defined by each sequence $\{L_i\}$ and calculated the derivative of the polynomial at $\delta = 0$. For the case of $\frac{\partial L}{\partial I_i}$ we varied the intensity of cell $i$ from $-\sigma_\mathrm{q}$ to $\sigma_\mathrm{q}$ and did the same approximation to a polynomial of order four and calculated its derivative. Figure \[fig:dLdx\] shows this cross check for $\frac{\partial L}{\partial x_i}$ and $\frac{\partial L}{\partial I_i}$. Although the derivatives are similar, they are not exactly the same for $\frac{\partial L}{\partial
x_i}$. This is caused by the polynomial coarseness fit, as explained below.

Figure \[fig:ajuste\] shows the curve fit for $\frac{\partial
L}{\partial x_i}$ for three different generators (generator number 37, 36 and 18 respectively). It can be seen in Figure \[fig:ajuste\] that the polynomial fit adjusts quite well to the function values for polygon number 37, so on Figure \[fig:dLdx\] both derivatives are the same. On the contrary, for polygons number 36 and 18, the fitted polynomial does not resemble the function $L$ at $\delta = 0$, causing a slight difference in their derivatives on Figure \[fig:dLdx\]. For polygon number 18 the polynomial does not fit the curve at all. This is the main problem of using a numerical approximation for the derivatives of $\{\vec{x_i}\}$: when two polygons are closer than $\delta$, the generator displacement causes the function $L$ to change abruptly (see Figure \[fig:VorProblem\]).

It can be seen that the analytical and numerical derivatives on Figure \[fig:dLdx\] are almost the same. As explained above, differences are produced because there are cases where the polynomials do not fit well to the variations of the merit function $L$ (for example, when two generators are too close). In an accurate calculation it is necessary to use the analytical derivatives.

Fitting a Voronoi Tessellation to an Image {#ap:fitting}
==========================================

Once we have a reasonable reconstruction for a pixelated image, we would like to fit a Voronoi tessellation to it in order to have a good initial starting point for the CG. This is done in an incremental way.

We start with a mesh consisting in only one polygon. We calculate the error per polygon as $$e_i^2 = \sum_l
  (I_i - I^\mathrm{im}_l)^2,$$ where the sum runs over all the pixels that fall inside polygon $i$, $I_i$ is the intensity of that polygon and $I^\mathrm{im}_l$ is the intensity of pixel $l$ in the image to be fitted. In each iteration we add a new polygon inside the one with the greatest error. The new generator is inserted in the position of the pixel that has the most different intensity value with respect to the mesh intensity.

Briggs, D. S., Schwab, F. R. & Sramek R. A. 1999, ASP Conf. Ser., 180, 127 Casassus, S., Cabrera, G. F., Förster, F, Pearson, T. J., Readhead, A. C. S., Dickinson, C. 2006, , 639, 951 Cornwell, T. J. & Evans, K. F. 1985, A&A, 143, 77 Högbom, J. A. 1974, A&AS, 15, 417 Narayan, Ramesh & Nityananda, Rajaram 1986, , 24, 127 Okabe, A., Boots, B. & Sugihara, K. 1992, Spacial Tessellations Concepts and Applications of Voronoi Diagrams, John Wiley & Sons Padin, S., et al, 2002, , 114, 83 Piña, R. K. & Puetter, R. C. 1993, , 105, 630 Press, W. H., Flannery, B. P., Teukilsky, S. A., Vettering, W. Y. 1992, Numerical Recipes in C, C. Cambridge University Press Shepherd, M.C., 1997, in ASP Conf. Ser. 125, Astronomical Data Analysis Software and Systems VI, ed. G. Hunt & H.E. Payne (San Francisco: ASP), 77 Sutton, E. C. & Wandelt, B. D. 2006, , 162, 401

[^1]: We choose to display the sky images in a larger field than the domain of free parameters; larger fields are required to highlight secondary side-lobes
