Deblurring galaxy images with Tikhonov regularization on magnitude domain

Natural images such as astronomical images are known to have smooth structures. In other words, values of adjacent pixels are similar to each other and norms of the image gradient should be small. This knowledge can be a priori information and often adopted as regularization terms, e.g. TV and Tikhonov regularization. The Tikhonov regularization is a quadratic norm with a linear operator;

where we adopt L as a differential operator for both horizontal and vertical directions, and $\|\cdot\|$ indicates the Euclidean norm.

Although the Tikhonov regularization successfully improves the image quality in astronomy (e.g. Kuramochi et al. \yearcite2018ApJ…858…56K¹¹1They call this term total squared variation.), it is not necessarily optimal for galaxy images. This is because the Tikhonov regularization linearly calculates the image gradients despite the fact that galaxies have a super linear profile such as a Sérsic profile. As a natural consequence, the image gradients in logarithmic space can be more appropriate for the regularization. Hence, we suggest a Tikhonov regularization term on magnitude-image domain. One problem in the use of magnitude is a divergence when pixel values are zero or close to zero. To avoid this problem, we apply asinh magnitudes suggested in Lupton et al. (1999).

where $a$ = $-2.5\log e$ = 1.08574, and b is an arbitrary parameter to determine the linearity behaviour. While $\mu(x)$ is linear in $x$ for $|x|\ll b$, it shows logarithmic behaviour for $|x|\gg b$. Although it is natural to choose $b$ to be an flux uncertainty as discussed in Lupton et al. (1999), we discuss the optimal choice of $b$ in section 4.3.

In summary, we propose the following formulation for a galaxy-image deblurring problem,

The difficulty in solving the problem (7) is the handling of the second term. It has three functions to be calculated; the linear operator $\mathbf{L}$, Euclidean norm, and asinh magnitude. To address this issue, we take a strategy that divides the variable $\mathbf{X}$ into two terms, $\mathbf{X_{flux}}$ and $\mathbf{X_{mag}}$, and updates them separately. Furthermore, the updates are conducted via a primal-dual splitting method (Condat, 2013), with which the calculations of linear operator and norm can be separated. This method has also another advantage that it can avoid an inner loop caused by a sub-problem, which is often appeared in the alternating direction method of multipliers (Boyd et al., 2011). In the following, we describe our strategy in detail.

Dividing the variable $\mathbf{X}$ into flux and magnitude terms, we modify the problem (7) as the following.

In the modified problem, the non-linear function (6), is removed from the objective function.

We solve the problem (2.2) via the primal-dual splitting method. To adopt the primal-dual splitting to our problem, here we define some variables and functions. We define the primal variable ${\bf\mathcal{X}}$, which have flux and magnitude parts,

To match the dimension, we also modify the observed image, dual variable, PSF matrix, and linear differential operator as, ${\bf\mathcal{Y}}=({\bf Y}\;{\bf 0})^{\top}$, ${\bf\mathcal{V}}=({\bf 0}\;{\bf V})^{\top}$, ${\bf\mathcal{A}}=({\bf A}\;{\bf 0})$, and ${\bf\mathcal{L}}=({\bf 0}\;{\bf L})$. Using these definitions, we also introduce data fidelity term F, regularization term H, and constraints G.

With the above definitions, the problem(2.2) can be expressed as the following.

This minimization problem is equivalent to a problem finding the saddle point of a Lagrangian,

where $\mathbf{\mathcal{V}}$ is the dual variable of $\mathbf{\mathcal{X}}$, and $H^{*}$ is the convex conjugate of $H$. As can be seen, the linear operator $\mathcal{L}$ and the function $\mathcal{H}$ (corresponds to Euclidean norm) are separated in equation (20), and they can be easily calculated. This is the greatest advantage of the primal-dual splitting method.

The primal and dual variables can be updated as the following.

where $\sigma$, $\tau$, and $\rho_{n}$ are the update parameters described later, and prox indicates the proximal operator calculated as follows. The proximal operator with $G_{1}$ calculates the positivity constraint,

The proximal operator with $G_{2}$ calculates the consistency between $\mathbf{X_{flux}}$ and $\mathbf{X_{mag}}$,

where equation (22) indicates the weighted average of flux and magnitude with the weight vector $\mathbf{w_{mag}}$, calculated as follows.

The proximal operation of Euclidean norm can be calculated as follows.

where we used Moreau’s identity (e.g. Bauschke & Combettes (2011)) in the first equality and defined $\alpha{\rm=}2\lambda/(2\lambda+\sigma)$ in the last equality.

Initial estimation of $\mathbf{X_{flux,0}}$ is arbitrary but we choose the observed image $\mathbf{Y}$. We confirmed that images with zero or random values lead to the same restoration results although more iterations are required for convergence. It indicates that our method is not sensitive to the initial estimation.

There are three tuning parameters in the primal-dual splitting, $\tau$, $\sigma$, and $\rho_{n}$. All of these parameters are mainly related to convergence speed, and not sensitive to the resultant image quality. The parameters $\tau$ and $\sigma$ are related to updates of the primal and dual variables, and should be $0<\tau<2/\beta$ and $\sigma\leq\left(\frac{1}{\tau}-\frac{\beta}{2}\right)/\|L\|^{2}$, where $\beta$ is the Lipschitz constant of the equation (10),

Although the choice of $\tau$ and $\sigma$ is arbitrary, we define them as the following.

With the above definition, these parameters always satisfy their constraint. The parameter $\rho_{n}$ is an update rate and can be modified during the iteration. The optimal choice of $\rho_{n}$ may accelerate the convergence, making use of the previous iteration. Nonetheless, we simply adopt $\rho_{n}$ of one, that is, we do not use the information from the previous update.

With the above setting, we can solve the problem (20). Nonetheless, the convergence might be very slow, depending on the parameters $\sigma$ and $\lambda$. This is because, if $\lambda\gg\sigma$, $\alpha$ in equation (29) converges to 1 and not sensitive to the regularization parameter $\lambda$. This means that many iterations are required to reflect the regularization into the restored image. To avoid this situation, we scale the image and the variance beforehand. The scaling parameter is as follows.

This scaling leads to $\sigma=\lambda$ and $\alpha=0.66$. Although $\alpha$ is independent of $\lambda$, the primal variable strongly depends on $\lambda$ because the scaling factor is proportional to $1/\sqrt{\lambda}$. This scaling as a pre-processing significantly accelerates the convergence of our problem.

The main objective of this paper is to investigate the feasibility of Tikhonov regularization on magnitude domain. For this purpose, we compare our method with a widely applied RL method and a method with conventional Tikhonov regularization on flux domain. We explain these methods in the following.

Simulation studies are suitable for a principle demonstration of our method because we have ground-truth images and perfect PSFs. The former can be directly compared with the reconstructed images, and owing to the latter we can also directly evaluate our method without uncertainty in PSFs.

We retrieved one spiral and one elliptical galaxies from the Subaru HSC-SSP DR2 (Aihara et al., 2019) observed with the i band. Hereafter we call them Spiral-1 and Elliptical-1, respectively. The use of different morphologies is suitable for performance verification. We visually chose the galaxies using the HSC map²²2https://hsc-release.mtk.nao.ac.jp/doc/index.php/tools-2/ and show their coordinates and redshifts in the upper side of Table 3.3. The galaxies have moderate redshifts, so that we can ignore the PSF blurring and avoid saturation. We resampled and shrank the images to images with 64$\times$64 pixels. The magnification factors of the pixel scales were 4.2 and 2.1 for Spiral-1 and Elliptical-1, respectively³³3It corresponds to pixel scales of 0.7 and 0.35 arcsec/pixel while the original is 0.168 arcsec/pixel.. Figures 1a and 2a show their resampled images.

The images were blurred with a realistic PSF. While Gaussian functions are often applied as a PSF, they are clearly distinct from realistic ones. We used a more realistic PSF. We combined the PSFs for Spiral-1 and Elliptical-1, which were provided by the HSC project⁴⁴4https://hsc-release.mtk.nao.ac.jp/psf/pdr2/. The FWHM of the combined PSF was $\sim$3 pixel. Because the images were shrank as mentioned in the previous paragraph, the convolution with the combined PSFs significantly blurred the images.

We added noises to the blurred images. To represent photon noises from both the background and objects, we applied Gaussian and Poisson distributions, respectively. We note that the sky background had been subtracted from the images and were assumed to be uniform across the images, and the sky noise was able to be expressed as a Gaussian distribution. We set the standard deviation of the Gaussian noise to 1.0% of the peak signal, whereas we produced two Poisson noises: the standard deviations of 1.0% and 5.0% of the peak signal. Hereafter the noise levels are referred to as 1.0% and 5.0% despite the additional 1.0% Gaussian noise. The procedure of the noise production is the same for both Spiral-1 and Elliptical-1. Panels (b) and (f) in figures 1 and 2 show the blurred Spiral-1 and Elliptical-1 images with 1.0% and 5.0% noises, respectively. The S/N at the centres are about 70 and 20 for both morphologies, respectively.

We adopted our deblurring method to the images. We set the hyper parameter $\lambda$ to 2.0 for Spiral-1 while 16.0 and 4.0 for Elliptical-1 with 1.0% and 5.0% noises, respectively. The other parameter b was set to 1.0$\bar{\sigma}_{y}$, where $\bar{\sigma}_{y}$ is the median flux uncertainty in the blurred image. The number of iterations was 1000 unless otherwise stated. We discuss the appropriate parameter values in section 4.3.

To evaluate our method with more realistic situations, we applied it on real observation images. In contrast to the simulation, observation images could have unknown and/or systematic errors so that PSF deconvolution becomes a more difficult problem. We chose one spiral and one elliptical galaxies from the HSC-SSP DR2 as listed in the lower side of Table 3.3. Hereafter we call them Spiral-2 and Elliptical-2. They have an appropriate size to fit to 64$\times$64 pixel images as shown in panels (b) and (f) of figures 3 and 4. We obtained the images, the variances, and the PSFs of each galaxy from the public data release site⁵⁵5https://hsc-release.mtk.nao.ac.jp/doc/index.php/tools-2/. We retrieved the images from both HSC-DUD and HSC-Wide surveys to obtain deep and shallow images, respectively. The S/N at the centres of Spiral-2 and Elliptical-2 were 409 and 178 for HSC-DUD, and 105 and 53 for HSC-Wide, respectively. We used only i-band images since other band images have no reference HST images as described later. The typical seeing at the i-band was $\sim$ 0.6 arcsec (Aihara et al., 2019). We did not use mask information because the images were produced with a sufficient number of observations and no noticeable artefacts were found.

These observed images were deblurred with our method. We set $\lambda$ to 2.0 and b to 1.0$\bar{\sigma}_{y}$ for both Spiral-2 and Elliptical-2. We discuss the appropriate parameter values in section 4.3.

Since we did not have ground-truth images, we instead used images taken with the HST. The HST has the F814W filter, whose wavelength coverage is similar to the i-band, so that the images can be used as a proxy of the ground-truth. The pixel scale was 0.05 arcsec, which were resampled to 0.168 arcsec to match the pixel scale of the Subaru/HSC. To avoid uncertainties in the flux calibration and sky subtraction, we scaled the images and added offsets to them to match the restored HSC i-band images. The HST images were obtained from the Hubble Legacy Archive (Jenkner et al., 2006; Whitmore et al., 2008, 2016).

For quantitative image evaluation, we applied two widely applied metrics, the peak signal to noise ratio (PSNR) and structural similarity (SSIM; Wang et al. (2004)). The PSNR of image $\boldsymbol{X}$ is defined as

where $N$ is the number of elements, and the subscript $ref$ indicates the ground-truth image. By definition, images similar to the ground-truth have a high PSNR.

The definition of SSIM is the following.

where $\mu$ is the average value, $\sigma_{x}$ is the standard deviation, $\sigma_{x,ref}$ is the covariance, $L$ is the dynamic range of the image, and $K_{1}$ and $K_{2}$ are arbitrary small values, where we adopted 0.001, following Wang et al. (2004). The terms $l$, $c$, and $s$ correspond to comparisons between the image and reference in terms of the luminance, contrast, and structure, respectively. When the image and reference are similar, these terms should be close to one, and the SSIM should also be close to one.

Since our simulation has a reference image for each galaxy, these metrics can be appropriate for our evaluation. On the other hand, the HST images are not perfect ground-truth ones for the observational evaluation. For example, when the restored HSC images have better image quality, the HST images do not work as the ground-truth. Hence, a care must be taken for the evaluation results.

Target galaxies used for the evaluation with simulation and observation. Name R.A. Dec. z Simulation Spiral-1 150.746077 2.342989 0.044 Elliptical-1 150.8170172 3.4772409 0.106 Observation Spiral-2 150.1728195 2.6090672 0.265 Elliptical-2 150.687837 2.61864322 0.860^a

• a

  photometric redshift from Tanaka et al. (2018).

Figures 1 and 2 show the image-restoration results for Spiral-1 and Elliptical-1, respectively. While all the three methods (RL, Tikhonov, and this work) successfully restore the blurred and noisy images (panels b and f), our method shows the best performance. Remarkable distinction can be seen especially in noisier images (bottom panels). We show the quantitative comparison in table 1. In the following, we discuss the image qualities in detail.

Reconstructed images with the RL method (panels c and g) show higher noise despite the higher spatial resolution than the blurred ones (panels b and f). This is because the iteration of the RL method amplifies the noise while improving the resolution. Owing to the trade-off, the number of iterations should be optimized. We determined the number to maximize the PSNR. While the best numbers of iterations were slightly different for PSNR and SSIM, we show the maximum values of them in table 1. Whereas we can determine the best number of iterations in a simulation study, we could not in general, so that the use of the best numbers leads to unfair comparisons. Nonetheless, the restored images with our method still show higher image qualities than those from the RL method. For example, for the Spiral-1 and Elliptical-1 images with 1% noises, the RL method provides the PSNRs of 40.39 and 43.98 dB whereas ours show 43.07 and 57.11 dB, respectively. These results strongly suggest that our method outperforms the RL method. We have to note that the image quality is affected by the number of iterations even in our method, which we discuss in subsection 4.2. Since the increase of the number of iterations after exceeding 1000 affects only a few central pixels, which could not be seen visually, we show only images with the 1000 iterations in the figures 1 and 2.

Panels (d) and (h) in figures 1 and 2 show the Spiral-1 and Elliptical-1 images that have 1 and 5% noises and are restored with the Tikhonov regularization. Compared with the results from the RL method, high frequency components caused by noises are reduced owing to the regularization, but the images are blurred. This results in lower PSNRs and SSIMs as shown in table 1. The low image qualities resulted from the regularization suggest that the assumed prior is not appropriate. Since both spiral and elliptical galaxies tend to have intensity profiles that exponentially change with the position, minimizing the gradients in linear space leads to blurring detailed structures, as mentioned in the introduction. In contrast, images restored with our method (panels e and i) are appropriately smoothed but remarkably reproduce the galaxy structures. The quantitative comparison in table 1 also shows that our method leads to the best qualities. These results suggest that the proposed regularization, a quadratic norm of the differential image in magnitude space (i.e. log space), is appropriate.

Since the elliptical galaxies are compact and the intensity profiles rapidly change, the difference of the restored images may not be obvious. To see the difference more clearly, we plot in figure 5a the horizontal profiles at the centre of the Elliptical-1 images, where the added noise level is 5.0%. At the first look, any de-blurring method produces sharper profiles than the blurred one (grey). Among them, our method (bright cyan) shows the most similar profile to the original (black).

Figures 3 and 4 show the results on real images (Spiral-2 and Elliptical-2, respectively). Similar to the simulation results, our method best restores the detailed structures in the galaxies. The quantitative evaluation in table 2 also shows that our method basically produces the highest PSNRs and SSIMs. The only case our method leads to lower PSNR than the RL is for Elliptical-2 in HSC-DUD. This is because the S/N exceeds 170 as mentioned in section 3.2. Although the RL method is sensitive to a noise, it can work for high S/N images. On the other hand, the lower S/N case, i.e. in the HSC-Wide, our method produces a higher PSNR than the RL. It indicates that our method can work even for low S/N images.

Here we discuss the optimal number of iterations for our method. It is generally difficult to determine the number of iterations. While increasing the number of iterations decreases the objective function, it directly leads to a high computational cost and does not necessarily improve the quality of restored images. Hence, the number of iterations should be optimized by balancing the computational costs with the resulting image quality.

The left panels of figure 6 show objective functions with the components against the number of iterations for the four galaxies. In all cases, the total cost and its components seem to converge at about 100 iterations. On the other hand, the PSNR for our method continues to increase until the number of iterations reaches to 1000$\sim$10000, as shown in the black lines in the right panels of figure 6.

The result that the PSNRs continues to increase despite the convergence of the objective function can be interpreted as the PSNR is sensitive to a small difference in the restored image. This is because the PSNR is expressed with a logarithmic of the MSE as shown in equation (35). In other words, the PSNR is sensitive to the detailed structures after the objective function almost converges and global structures are reproduced. It may lead to an idea of including a logarithmic function into the objective function to accelerate the convergence speed. However, it makes the calculation of the Lipschitz constant difficult and may also make the convergence unstable. Fortunately, our analysis shows that even the moderate number of iterations (i.e. $\sim$1000) produces a good image quality as shown in figure 6 and tables 1–2. Hence, the conclusion that our regularization method is effective for deblurring galaxy images is robust.

We apply the RL method with the optimal number of iterations that leads to the highest PSNR and SSIM. As mentioned earlier, it is usually not possible to know the optimal number, and one has to stop the iteration at some time. However, the optimal number of iterations for the RL method significantly depends on cases. As shown in the red lines in the right panels of figure 6, it was 23, 34, 127, and 35 for Spiral-1 with 5.0% noise, Elliptical-1 with 1.0% noise, Spiral-2 from the HSC-DUD, and Elliptical-2 from the HSC-Wide, respectively. The serious problem is an over-optimization, which leads to a significantly low PSNR and degrades the image as shown in figure 7a and 7c.

In contrast, our method is stable during the iterations. As shown in figure 7b and 7d, the large number of iteration does not degrade the restored images. We note that the PSNR of our method decreases with the number of iterations for real observation images as shown in figure 6g and 6h. It suggests that the best number of iterations is 572 for Elliptical-2 from HSC-Wide, and the adopted number of 1000 is too large. Although PSNRs from both RL and our methods decrease with increasing the iterations, the causes are different as indicated from figure 7. To explore the cause of the decrease in PSNR, we show the horizontal profile of Elliptical-2 from the HSC-Wide in figure 5b. It shows that, in our method, increasing the number of iterations leads to a sharper profile than the HST image. One possible interpretation is that the PSF is slightly fatter than the reality because this behaviour is only seen in case of real images. If this is the case, the optimization leads to too sharp images and it might be better to adopt a moderate number of iterations. Nonetheless, as indicated in figure 7b and 7d, too many iterations do not lead to noise amplification. Hence, we conclude that our method can produce stable solutions.

The qualities of restored images depend on hyper parameters and our method have four parameters: $b$, $\lambda$, $\tau$ and $\sigma$. The latter two parameters, $\tau$ and $\sigma$, are the step lengths of the primal-dual splitting method and, in this study, they are automatically determined by equations (31) and (32). In the following, we discuss the general characteristics of the former parameters, $b$ and $\lambda$, and investigate the dependence of our method on these parameters.

The parameter $b$ is used in the definition of asinh magnitude which is proportional to $x$ and $\log{x}$ when $|x|\gg b$ and $|x|\ll b$, respectively. The concept of using asinh magnitude is to appropriately treat pixels having zero or negative flux caused by a noise. In this context, it is natural to choose $b$ to be related to flux uncertainty as suggested in Lupton et al. (1999). In the regularization point of view, the parameter $b$ should be less than the brightness which fine structures to be restored would have. As discussed in section 2, galaxies have much steeper profiles than linear and this is why we apply asinh magnitude in our regularization. Nonetheless, galaxy structures fainter than $b$ is calculated in linear space and they can be smoothed out if $b$ is inappropriately large. In the optimization point of view, it is favourable that $b$ is proportional to the flux uncertainty. In this case, both the data-fidelity and regularization terms in our objective function (7) include the square of flux uncertainty in the denominator and balancing these terms becomes less sensitive to S/N of the measurements.

The parameter $\lambda$ is used to balance the data fidelity and regularization terms in equation (7). In general, a higher $\lambda$ leads to smoothing out detailed structures while reducing noises, and this parameter should be chosen to balance these two effects. As mentioned in the previous paragraph, both terms in our cost function include the square of flux uncertainty in the denominator. In this case, the balance between the data fidelity and regularization terms is less sensitive to $\lambda$. Exception is the case where some structures have a remarkably high S/N. The asinh magnitude in our regularization shows a logarithmic behaviour at higher S/N regions, leading to less weight to the regularization than the data fidelity. That is, high S/N images may favour high $\lambda$, contrary to intuition.

We investigated the dependence of restored image quality on the two hyper parameters, $b$ and $\lambda$, via simulations. We used Spiral-1 and Elliptical-1 with noises of 1% and 5%, the same as in section 4.1. The images were restored with our method with $\lambda=0.125-32$ and $b=0.05-3.0\,\bar{\sigma}_{y}$. The number of iterations was fixed to 10000. Figure 8 shows the PSNRs at various values of $\lambda$ and $b$ for the above galaxies.

It can be seen that the PSNR is not sensitive to b at $b<1.0\,\bar{\sigma}_{y}$ while a higher b leads to lower PSNRs. The latter indicates that fainter structures than b are smoothed out owing to unsuitable b as they can be treated in linear space as previously mentioned. From these results, we chose $b=1.0\,\bar{\sigma}_{y}$ throughout this paper.

Next, we discuss the dependence on $\lambda$, which is not as simple as $b$. The figure 8 shows that the best $\lambda$ is typically $2\sim 4$, i.e. the lower and higher $\lambda$ lead to noisy and blurred images, respectively. On the other hand, for Elliptical-1 with 1% noises, $\lambda=16$ leads to the highest PSNR. It can be interpreted as this galaxy has a steep profile and high S/N regions, which leads to less weight in the regularization term owing to the asinh magnitude. In this case, despite using flux uncertainty in the regularization term, the best regularization parameter still depends on the S/N.

On the other hand, the use of $\chi^{2}_{\nu}$ could mitigate the problem. Figure 9 shows the PSNR against $\chi^{2}_{\nu}$. A higher $\chi^{2}_{\nu}$ corresponds to lower weights in the data-fidelity term due to higher $\lambda$. We can see that the PSNR achieves the highest value at $\chi^{2}_{\nu}$ of slightly lower than 1.0, regardless of the galaxy type and the noise level. Hence, $\lambda$ can be chosen to obtain $\chi^{2}_{\nu}$ of slightly less than 1.

In section 4.1, we demonstrated that the proposed regularization, which minimizes image gradients in magnitude domain, leads to the best restored image quality. Nonetheless, our method still has some limitations, which we discuss in the following.

The most considerable drawback of our method is the slow convergence; our method often requires a large number of iterations to obtain the best image quality, as discussed in section 4.2. Although images with a high quality can still be obtained with a moderate number of iterations, N=100-1000, more iterations are required to obtain the highest PSNR (figure 6 e-h). Calculation time under KM’s computer setup (Intel Core i5, quad-core 1.1 GHz) was 15s and 180s for 1000 and 10000 iterations, respectively. One thing to be considered is the rapid decrease of the cost function (figure 6 a-d). It indicates that the global structures are rapidly reproduced while the local detailed structures are not because it is less sensitive to the objective function, as discussed in section 4.2. Modifying the cost function to be sensitive to the detailed structures may accelerate the convergence, which could be a future work.

As in many other deblurring methods, our method provides only maximum-point estimation of the posterior probability distribution and no covariance matrix, although astronomical analysis usually requires uncertainty in pixel values. Besides, the variance could not be estimated from the restored image, since the pixel values are not independent of each other. Hence, some artefacts might be regarded as a “significant detection”. To address this problem, a Monte-Carlo technique or a bootstrap method with a different set of raw images could be a solution. A variational Bayes method could also be useful (e.g. Babacan et al. (2008), Ruiz et al. (2015), Sonogashira et al. (2017), \yearciteSonogashira18). In this method, not only the best point in each pixel value, but also the probability distribution can be simultaneously estimated. Although both the cost function and optimization method should be changed and the computational cost increases, this method is worth devising in a future work.

Finally, simultaneous restoration of multi-band images can be considered to produce resolution-matched images. To produce such data set, sharper images are usually blurred to match with the worst case, which leads to a loss of information and is not appropriate. However, deblurring all the images does not guarantee spatial resolution matching, especially in case of a different PSF size and/or S/N. To address this problem, the fact that galaxy structures are very similar in different band images can be used as a priori information. Of course the structures are not identical between different band images, but the colour gradients (in magnitude domain) should be smooth. If this is the case, our method can naturally be extended towards the colour space, which could also be a future work.