# Libhires libhires is a C++ image reconstruction library for Time Ordered Information (TOI) data from the Planck satellite. # Theory The basis of all of the algorithms used in libhires is minimap. This simply bins the TOI's into the output pixels and averages the results per pixel. It does not attempt a true deconvolution, which would require taking into account the PSF of the detectors at different frequencies. The more sophisticated methods take into account the detector specific PSF. They all start by computing a minimap image. Then they deconvolve that minimap image. At the boundaries, the PSF can depend on pixels outside of the data. To reduce the effect from the edges, all of the methods reflect the boundaries. So U(-x) = U(x) U(x_boundary + x) = U(x_boundary - x) The basic equation that the deconvolution methods start with is the imaging equation b=Ax+N where 'b' is the measured data, 'x' is the true signal, 'A' maps the signal to the measurement using the PSF, and 'N' is the noise. All of the deconvolution methods implemented in libhires can be understood as trying the maximize the Bayesian a posteriori probability P(x|b). This probability can be expressed in terms of the probability of the measurement given the signal P(b|x) and the prior probability of the signal P(x). P(x|b) ~ P(b|x) P(x) The probability of the measurement P(b) is unity here. The differences between the deconvolution methods arise mainly because they choose different P(b|x) and P(x). ## MCM MCM stands for the [Maximum Correlation Method](http://labs.adsabs.harvard.edu/adsabs/abs/1990AJ.....99.1674A/). This is also known as the HIRES method for IRAS. MCM assumes that the noise follows a Poisson distribution. The prior probability for the signal P(x) is the improper constant distribution over the positive numbers. P(x) ~ x >0 constant x<=0 0 These two assumptions require that both the measured flux 'b' and the solution 'x' are positive. This implies that every pixel of the minimap image must be positive. For data from Planck's high frequency detectors, this is almost always true. This is less true for the lower frequencies. MCM then iteratively smoothes an initial guess to find a solution. The iteration is constructed so that flux is always conserved. One defect of this method is that there is no unique solution. Repeated iterations can lead to divergent results. This is a natural consequence of P(x) being a constant. For many cases, the solution does not actually wander off to infinity because flux is conserved and can never be negative. That means that the recovered signal always depends on the resolution of the recovered image. Also, the original MCM did not have binning. We have added binning in order to eliminate a bias whereby the recovered signal in regions with more samples would be brighter than regions with fewer samples. With that said, MCM is relatively fast, easy to compute, and widely used. It also conserves flux. It is not clear if that is a good thing, since it conserves the flux of the noise. ## Elastic Net [Elastic Net](https://en.wikipedia.org/wiki/Elastic_net_regularization) starts by assuming the noise follows a Gaussian distribution. Using only this would give the traditional least-squares solution. However, this is not sufficient, because this would lead to overfitting. The solution would become wildly varying in trying to fit the noise. To alleviate this, Elastic Net preferentially chooses smaller signals by setting the prior probability of the signal to a combination of Laplace and Gaussian distributions. P(x) ~ exp((x^2)/(2*g^2)) exp(-|x/l|) The Gausssian prior makes the recovered solution less prone to large variations. The Laplacian prior tends to select some regions for signal while setting the rest of the image to zero. Using the variance of the measurements v=Var(b) their maximum size m=max(|b|) and the ratio of the number of pixels in the measurement 'N_b' and recovered signal 'N_x' r=N_b/N_x we set the constants 'g' and 'l' as g=r*v/(m*m) l=r*v/m This results in a matrix that is, in principle, invertible. In practice, we use mlpack's [LARS](http://www.mlpack.org/doxygen.php?doc=classmlpack_1_1regression_1_1LARS.html) class to compute a solution. Elastic Net has a well defined solution with relatively few artifacts. However, it is much, much slower than MCM, making it impractical for large image reconstructions. # Implementation Libhires is written in C++ and depends on * [mlpack](http://www.mlpack.org/) * [armadillo](http://arma.sourceforge.net/). * [boost](http://www.boost.org/) * [cfitsio](http://heasarc.gsfc.nasa.gov/fitsio/fitsio.html) * [CCfits](http://heasarc.gsfc.nasa.gov/fitsio/CCfits/) * [libxml2](http://xmlsoft.org/) Libhires uses C++11 features, and has been tested with gcc versions 4.7 and 4.9. # Testing There is a simple test program in test/deconvolve.cxx that gets built into build/deconvolve Running this compares a deconvolution of a gaussian elliptical measurement against expected results in test/expected/ # License Copyright © 2012-2015, California Institute of Technology Authors: Walter Landry, Peregrine McGehee, and Angela Zhang Based on Government Sponsored Research NAS7-03001 and NNN12AA01C. This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.