16.4 Image reconstruction process

As the experimental frequency list is finite, and in addition the experimental data blurred, the object representation that can be obtained from these data is of course incomplete. This simple remark shows that the inverse problems of Fourier synthesis must be regularized: the high-frequency components of the image to be reconstructed must be negligible.

The central problem is to specify in which conditions it is possible to extrapolate or interpolate, in some region of the Fourier domain, the Fourier transform of a function $\phi$ whose support is contained in some finite region of $H_o$. It is now well established that extrapolation is forbidden, and interpolation allowed to a certain extent. The corresponding regularization principle is then intimately related to the concept of resolution: the interpolation is performed in the frequency gaps of the frequency coverage to be synthesized.

16.4.1 Synthesized aperture

Let  $\mathcal{H}$ be the Fourier domain: $\mathcal{H}=(-\Delta u/2, \Delta u/2)^2$. In Fourier synthesis, the frequency coverage to be synthesized is a centro-symmetric region  $\mathcal{H}_s \subset \mathcal{H}$ (see Fig. 16.4).

CLEAN and WIPE share a common objective, that of the image to be reconstructed. This image, $\Phi_{\! s}$, is defined so that its Fourier transform is quadratically negligible outside  $\mathcal{H}_s$. More explicitly, $\Phi_{\! s}$ is defined by the convolution relation:

$$\Phi_{\! s} = \Theta_s \star \Phi_{\! o}.$$ (16.7)

The ``synthetic beam'' $\Theta _s$ is a function resulting from the choice of  $\mathcal{H}_s$: the well-known clean beam in CLEAN, the neat beam in WIPE.

16.4.2 Synthetic beam

The neat beam can be regarded as a sort of optimal clean beam: the optimal apodized point-spread function that can be designed within the limits of the Heisenberg principle. More precisely, the neat beam $\Theta _s$ is a centro-symmetric function lying in the object space $H_o$, and satisfying the following properties:

• The energy of  $\widehat{\Theta}_s$ is concentrated in  $\mathcal{H}_s$. In other words, $\widehat{\Theta}_s$ has to be small outside  $\mathcal{H}_s$ in the mean-square sense: we impose the fraction $\chi^2$ of this energy in  $\mathcal{H}_s$ to be close to $1$ (say  $\chi^2=0.98$).
• The effective support  $\mathcal{D}_s$ of $\Theta _s$ in $H_o$ is as small as possible with respect to the choice of  $\mathcal{H}_s$ and $\chi^2$. The idea is of course to have the best possible resolution.

This apodized point-spread function is thus computed on the grounds of a trade-off between resolution and efficiency, with the aid of the power method.

Figure: Experimental frequency coverage and frequency coverage to be synthesized  $\mathcal{H}_s$ (left hand). The experimental frequency list $\mathcal{L}_e$ includes $N_e=2862$ frequency points. The frequency coverage to be synthesized  $\mathcal{H}_s$ is centred in the Fourier grid  $\mathbb{G} \delta u$, where $\delta u = \Delta u/N$ with $N=128$ (here, the diameter of the circle is equal to  $40 \delta u$). The neat beam $\Theta _s$ (right hand) represented here corresponds to the frequency coverage to be synthesized $\mathcal{H}_s$ for a given value of  $\chi ^2=0.97$. It is centred in the object grid $\mathbb{G} \delta x$ where  $\delta x=1/\Delta u$ (here, the full width of $\Theta _s$ at half maximum is equal to  $5 \delta x$).

16.4.3 Regularization frequency list

As extrapolation is forbidden, and interpolation only allowed to a certain extent in the frequency gaps of the frequency coverage to be synthesized, the experimental frequency list  $\mathcal{L}_e$ should be completed by high-frequency points. These points, located outside the frequency coverage to be synthesized  $\mathcal{H}_s$, are those for which the high-frequency components of the image to be reconstructed are practically negligible.

The elements of the regularization frequency list  $\mathcal{L}_r$ are the frequency points  $\mathbf{u}_r$ located outside the frequency coverage to be synthesized  $\mathcal{H}_s$ at the nodes of the Fourier grid  $\mathbb{G} \delta u$:

$$\mathcal{L}_r = \{ \mathbf{u}_r=\mathbf{q} \delta u, \mathbf{q}\in\mathbb{G}: \mathbf{q} \delta u \not\in \mathcal{H}_s \}.$$ (16.8)

The global frequency list  $\mathcal{L}$ is then the concatenation of  $\mathcal{L}_e$ with  $\mathcal{L}_r$.

16.4.4 Data space

According to the definition of the image to be reconstructed, the Fourier data corresponding to  $\Phi_{\! s}$ are defined by the relationship:

$$\Psi_{\! s}(\mathbf{u}) = \widehat{\Theta}_s(\mathbf{u}) \Psi_{\! e}(\mathbf{u}) \quad\forall \mathbf{u} \in \mathcal{L}_e.$$ (16.9)

Clearly, $\Psi_{\! s}$ lies in the experimental data space $K_e$.

Let us now introduce the data vector:

$$\Psi_{\! d}(\mathbf{u}) = \begin{cases}\Psi_{\! s}(\mathbf{u}) & \text{on $\mathcal{L}_e$;}\ 0 & \text{on $\mathcal{L}_r$.} \end{cases}$$ (16.10)

This vector lies in the data space $K_d$, the real Euclidian space underlying the space of complex-valued functions $\psi$ on  $\mathcal{L}$, such that  $\psi(-\mathbf{u}) = \bar{\psi}(\mathbf{u})$. This space is equipped with the scalar product:

$$(\psi_1 \mid \psi_2)_d = \sum_{\mathbf{u}\in\mathcal{L}_e} \bar{\... ...{u}\in\mathcal{L}_r} \bar{\psi_1}(\mathbf{u}) \psi_2(\mathbf{u}) (\delta u )^2;$$ (16.11)

$W(\mathbf{u})$ is a given weighting function that takes into account the reliability of the data via the standard deviation  $\sigma_e(\mathbf{u})$ of  $\Psi_{\! e}(\mathbf{u})$, as well as the local redundancy  $\rho(\mathbf{u})$ of  $\mathbf{u}$ up to the sampling interval $\delta u$.

The Fourier sampling operator $A$ is the operator from the object space $H_o$ into the data space $K_d$:

$$A: H_o \longrightarrow K_d, \qquad (A\phi)(\mathbf{u}) = \begin... ...idehat{\phi}(\mathbf{u}) & \text{on $\mathcal{L}_r$.} \end{cases}$$ (16.12)

As the experimental data  $\Psi_{\! e}(\mathbf{u})$ are blurred values of  $\widehat{\Phi}_o(\mathbf{u})$ on  $\mathcal{L}_e$, this operator will play a key role in the image reconstruction process. The definition of this Fourier sampling operator suggests that the action of $A$ should be decomposed into two components: $A_e$ on the experimental frequency list  $\mathcal{L}_e$, and $A_r$ on the regularization frequency list  $\mathcal{L}_r$.

16.4.5 Object representation space

The reconstructed image is defined as the function $\Phi _E$ of the object space $H_o$ minimizing some objective functional. The definition of this functional takes into account the nature of the data, as well as other constraints. For example, the image to be reconstructed may be confined to a subspace, or more generally to a convex set, of the object space $H_o$: this convex set is the object representation space $E$. It may be defined from the outset (in an interactive manner, for example), or step by step throughout the image reconstruction procedure (this is the case of the current implementation of WIPE). In both cases, the projection operator onto this space, the projector $P_E$, will play an essential role in the image reconstruction process.

REMARK 1: positivity constraint.

In most cases encountered in practice, the scalar components of $\Phi _E$ in the interpolation basis of $H_o$ must be non-negative (cf. Eq.). In the current implementation of WIPE this constraint is taken into account. The object representation space $E$ is then built, step by step, accordingly.

16.4.6 Objective functional

The reconstructed image is defined as the function $\Phi _E$ minimizing on $E$ the objective functional:

$$q(\phi) = \Vert \Psi_{\! d} - A \phi\Vert^2_d.$$ (16.13)

According to the definition of the data vector  $\Psi_{\! d}$ and to that of the Fourier sampling operator $A$, this quantity can be written in the form:

$$q(\phi) = q_e(\phi) + q_r(\phi) \begin{cases}q_e(\phi) = \displaystyle \!\!\sum_{\mathbf{u}\in\... ...rt \widehat{\phi}(\mathbf{u}) \vert^2 (\delta u)^2. & \end{cases}$$ (16.14)

The experimental criterion $q_e$ constraints the object model $\phi$ to be consistent with the damped Fourier data  $\Psi_{\! s}$, while the regularization criterion $q_r$ penalizes the high-frequency components of $\phi$.

Let now $F$ be the image of $E$ by $A$ (the space of the $A\phi$'s, $\phi$ spanning $E$), $A_E$ be the operator from $E$ into $F$ induced by $A$, and $\Psi _F$ the projection of  $\Psi_{\! d}$ onto $F$ (see Fig. 16.7). The vectors $\phi$ minimizing $q$ on $E$, the solutions of the problem, are such that  $A_E\phi = \Psi_F$. They are identical up to a vector lying in the kernel of $A_E$ (by definition, the kernel of $A_E$ is the space of vectors $\phi$ such that  $A_E\phi = 0$).

As  $\Psi_{\! d} - \Psi_F$ is orthogonal to $F$, the solutions $\phi$ of the problem are characterized by the property: $\forall \varphi \in E, (A\varphi \mid \Psi_{\! d} - A\phi)_d = 0$. On denoting by $A^*$ the adjoint of $A$, this property can also be written in the form:

$$\forall \varphi \in E, \quad (\varphi \mid r)_o = 0, r = A^*(\Psi_{\! d} - A\phi).$$ (16.15)

where $r$ is regarded as a residue. This condition is of course equivalent to $P_E r = 0$, where $P_E$ is the projector onto the object representation space $E$. The solutions of the problem are therefore the solutions of the normal equation on $E$:

$$A^*_E A^{\vphantom{*}}_E \phi = A^*_E \Psi_{\! d},$$ (16.16)

where  $A^*_E = P_E A^*$.

Many different techniques can be used for solving the normal equation (or minimizing $q$ on $E$). Some of these are certainly more efficient than others, but this is not a crucial choice.

REMARK 2: beams and maps.

The action of $A^*A$ involved in  $A^*_E A^{\vphantom{*}}_E$ is that of a convolutor. As the two lists  $\mathcal{L}_e$ and  $\mathcal{L}_r$ are disjoints, we have: $A^*A = A_e^*A_e + A_r^*A_r$. Thus, the corresponding point-spread function, called the dusty beam, has two components: the traditional dirty beam $\Theta_d$ and the regularization beam. The latter corresponds to the action of $A_r^*A_r$, the former to that of $A_e^*A_e$ (see Fig. 16.5). Likewise, according to the definition of the data vector, $A^*\Psi_{\! d} = A_e^*\Psi_{\! s}$ is called the dusty map (as opposed to the traditional dirty map  $A_e^*\Psi_{\! d}$ because it is damped by the neat beam).

REMARK 3: construction of the object representation space.

With regard to the construction of the object representation space $E$, CLEAN and WIPE are very similar: it is defined through the choice of the (discrete) object support. It is important to note that this space may be constructed, in a global manner or step by step, interactively or automatically. In the last version of WIPE implemented at IRAM, the image reconstruction process is initialized with a few iterations of CLEAN. The support selected by CLEAN is refined throughout the iterations of WIPE by conducting a matching pursuit process at the level of the components of $r$ in the interpolation basis of $H_o$: the current support is extended by adding the nodes of the object grid  $\mathbb{G} \delta x$ for which these coefficients are the largest above a given threshold (half of the maximum value, for example). The objective functional is then minimized on that new support, and the global residue $r$ updated accordingly. The object representation space of the reconstructed image is thus obtained step by step in a natural manner.

The simulation presented on Fig.16.5-16.6 corresponds to the conditions of Fig. 16.4. The Fourier data  $\Psi_{\! e}$ were blurred by adding a Gaussian noise: for all  $\mathbf{u}\in\mathcal{L}_e$, the standard deviation of  $\Psi_{\! e}(\mathbf{u})$ was set equal to $5 \%$ of the total flux of the object ( $\widehat{\Phi}_o(\mathbf{0)}/20$). The image reconstruction process was initialized with a few iterations of CLEAN, and the construction of the final support of the reconstructed image was made as indicated in Remark 3. At the end of the reconstruction process, a final smoothing of the current object support was performed. In this classical operation of mathematical morphology, the effective support of $\Theta _s$, $\mathcal{D}_s$, is of course used as a structuring element. The boundaries of the effective support of the reconstructed neat map are thus defined at the appropriate resolution. In particular, the connected entities of size smaller than that of  $\mathcal{D}_s$ are eliminated.

Figure: Dirty beam (left hand) corresponding to the experimental frequency list  $\mathcal{L}_e$ of Fig. 16.4, and dusty map (right hand) of a simulated data set (the simulated Fourier data  $\Psi_{\! e}$ were blurred by adding a Gaussian noise with a standard deviation $\sigma _e$ equal to $5 \%$ of the total flux of the object  $\Phi_{\! o}$).
Figure: Image to be reconstructed  $\Phi_{\! s}$ (left hand) at the resolution level defined in Fig. 16.4, and reconstructed neat map $\Phi _E$ (right hand) at the same resolution: the final condition number $\kappa _E$ is equal to $2.46$ (cf. Eq. 16.17 and 16.18).

16.4.7 Uniqueness and robustness

When the problem is well-posed, $A_E$ is a one-to-one map ( $\ker A_E = \{0\}$) from $E$ onto $F$; the solution is then unique: there exists only one vector $\phi\in E$ such that  $A_E\phi = \Psi_F$. This vector, $\Phi _E$, is said to be the least-squares solution of the equation $A_E\phi$ ``=''  $\Psi_{\! d}$.

In this case, let  $\delta\Psi_F$ be a variation of $\Psi _F$ in $F$, and  $\delta\Phi_E$ be the corresponding variation of $\Phi _E$ in $E$ (see Fig. 16.7). It is easy to show that the robustness of the reconstruction process is governed by the inequality:

$$\frac{\Vert\delta\Phi_E\Vert_o}{\Vert\Phi_E\Vert_o} \leq \kappa_E \frac{\Vert\delta\Psi_F\Vert_d}{\Vert\Psi_F\Vert_d}.$$ (16.17)

The error amplifier factor $\kappa _E$ is the condition number of $A_E$:

$$\kappa_E = \frac{\sqrt{\lambda'}}{\sqrt{\lambda}};$$ (16.18)

here $\lambda$ and $\lambda'$ respectively denote the smallest and the largest eigenvalues of  $A^*_E A^{\vphantom{*}}_E$. The closer to $1$ is the condition number, the easier and the more robust is the reconstruction process (see Fig. 16.8 and 16.9).

Figure 16.7: Uniqueness of the solution and robustness of the reconstruction process. Operator $A$ is an operator from the object space $H_o$ into the data space $K_d$. The object representation space $E$ is a particular subspace of $H_o$. The image of $E$ by $A$, the range of $A_E$, is denoted by $F$. In this representation, $\Psi _F$ is the projection of the data vector  $\Psi_{\! d}$ onto $F$. The inverse problem must be stated so that $A_E$ is a one-to-one map from $E$ onto $F$, the condition number $\kappa _E$ having a reasonable value.

The part played by inequality 16.17 in the development of the corresponding error analysis shows that a good reconstruction procedure must also provide, in particular, the condition number $\kappa _E$. This is the case of the current implementation of WIPE which uses the conjugate gradient method for solving the normal equation 16.16.

To conduct the final error analysis, one is led to consider the eigenvalue decomposition of  $A^*_E A^{\vphantom{*}}_E$. This is done, once again, with the aid of the conjugate gradient method associated with the QR algorithm. At the cost of some memory overhead (that of the $M$ successive residues), the latter also yields approximations of the eigenvalues $\lambda_k$ of  $A^*_E A^{\vphantom{*}}_E$. It is thus possible to obtain the scalar components of the associated eigenmodes $\Phi_k$ in the interpolation basis of $H_o$. The purpose of this analysis is to check whether some of them (in particular those corresponding to the smallest eigenvalues) are excited or not in $\Phi _E$. If so, the corresponding details may be artefacts of the reconstruction.

Figure: Reconstructed neat map $\Phi _E$ (left hand) and eigenmode $\Phi_1$ (right hand) corresponding to the smallest eigenvalue  $\lambda_1=0.165$ of $A^*_E A^{\vphantom{*}}_E$. The conditions of the simulations are those of Fig. 16.4 and 16.5: in particular, the diameter of  $\mathcal{H}_s$ is equal to  $40 \delta u$. The final condition number is  $\kappa_E=2.46$ (the eigenvalues of  $A^*_E A^{\vphantom{*}}_E$ are plotted on the bar code below). This eigenmode is not excited in $\Phi _E$: the separation angle $\theta_1$ between $\Phi _E$ and $\Phi_1$ is greater than $89^\circ$. In other situations, when the final condition number is greater, this mode may be at the origin of some artefacts in the neat map (see Fig. 16.9).

The reconstructed map is then decomposed in the form:

$$\Phi_E = \sum_{k=1}^M w_k \Phi_k, \quad w_k = ( \Phi_k \mid \Phi_E ).$$ (16.19)

The separation angle $\theta_k$ between $\Phi _E$ and $\Phi_k$ is explicitly given by the relationship:

$$\cos\theta_k = \frac{w_k}{\sqrt{\sum_{k=1}^M w_k^2}} \quad (0 \le \theta_k \le \pi / 2).$$ (16.20)

The closer to $\pi/2$ is $\theta_k$, the less excited is the corresponding eigenmode $\Phi_k$ in the reconstructed neat map $\Phi _E$.

To illustrate in a concrete manner the interest of equations 16.19 and 16.20, let us consider the simulations presented in Fig. 16.4 and 16.9. Whatever the value of the final condition number is, the error analysis allows the astronomer to check if there exists a certain similitude between some details in the neat map and some features of the critical eigenmodes. This information is very attractive, in particular when the resolution of the reconstruction process is greater than a reasonable value (the larger is the aperture to be synthesized  $\mathcal{H}_s$, the smaller is the full width at half-maximum of $\Theta _s$). In such situations of ``super resolution,'' the error analysis will suggest the astronomer to redefine the problem at a lower level of resolution, or to keep in mind that some details in the reconstructed neat map may be artefacts of the reconstruction process.

Figure: Reconstructed neat map $\Phi _E$ (left hand) and related critical eigenmode $\Phi_1$ (right hand). The latter corresponds to the smallest eigenvalue $\lambda_1=0.057$ of $A^*_E A^{\vphantom{*}}_E$. The conditions of the simulations are those of Fig. 16.5, but here the diameter of  $\mathcal{H}_s$ is taken equal to  $48 \delta u$: the final condition number is  $\kappa_E=4.19$ (the eigenvalues of  $A^*_E A^{\vphantom{*}}_E$ are plotted on the bar code below). The critical eigenmode $\Phi_1$ is at the origin of the oscillations along the main structuring entity of $\Phi _E$. This mode is slightly excited (the separation angle $\theta_1$ between $\Phi _E$ and $\Phi_1$ is less than $86^\circ$), thus the corresponding details may be artefacts. In this case of ``super-resolution'' the error analysis provided by WIPE suggests that the procedure should be restarted at a lower level of resolution (see Fig. 16.8), so that the final solution be more stable and reliable.