From 5591f22b4eb0950445a2d1cea45fba6004f99166 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Thi=C3=A9baut?= Date: Mon, 5 Feb 2024 16:49:02 +0100 Subject: [PATCH] Update doc. --- README.md | 104 +++++- notes.md | 575 ++++++++++++++++++++++++++++++++ old/beauty2016.jl | 793 +++++++++++++++++++++++++++++++++++++++++++++ old/metrics.jl | 380 ++++++++++++++++++++++ src/interpolate.jl | 110 +++++++ 5 files changed, 1961 insertions(+), 1 deletion(-) create mode 100644 notes.md create mode 100644 old/beauty2016.jl create mode 100644 old/metrics.jl create mode 100644 src/interpolate.jl diff --git a/README.md b/README.md index d98b21b..a5dc761 100644 --- a/README.md +++ b/README.md @@ -1 +1,103 @@ -# ImageMetrics [![Build Status](https://github.com/emmt/ImageMetrics.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/emmt/ImageMetrics.jl/actions/workflows/CI.yml?query=branch%3Amain) [![Coverage](https://codecov.io/gh/emmt/ImageMetrics.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/emmt/ImageMetrics.jl) +# Image metrics +[![Build Status](https://github.com/emmt/ImageMetrics.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/emmt/ImageMetrics.jl/actions/workflows/CI.yml?query=branch%3Amain) +[![Coverage](https://codecov.io/gh/emmt/ImageMetrics.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/emmt/ImageMetrics.jl) + +This repository collect notes and tools for image metrics that is means to +quantitatively compare images. + +## References + + +## Image Metrics Principles + +Assessing image quality for optical interferometry has been studied[^Gomnes+2016] +specific: + +* bright object on a dark background; +* extrapolation of the synthesized filed of view can be naturally done by zero-padding (i.e. assuming that missing pixels have a zero value). + + +The quality of an image has to be assessed by an objective quantitative +criterion. The literature shows that establishing a universal image quality +criterion is a very controversial subject. What is the best criterion also +largely depends on the context. Here we will assume that the *metric* +$\def\V#1{\boldsymbol{#1}}\def\M#1{\mathbf{#1}}\def\Dist{\mathcal{D}}\Dist(\V{x}, \V{y})$ is used to estimate the discrepancy between a reconstructed +image $\V{x}$ and a reference image $\V{y}$. To simplify the discussion, we also assume +that the lower $\Dist(\V{x}, \V{y})$ the better the agreement between $\V{x}$ and $\V{y}$. In +other words $\Dist(\V{x}, \V{y})$ can be thought as a measure of the distance between +$\V{x}$ and $\V{y}$. + + +When assessing image quality it is important that the result does not depend on +irrelevant changes. This however depends on the type of images and on the +context. For instance, for object detection or recognition, the image metric +should be insensitive to the background level, to a geometrical transform +(translation, rotation, magnification, etc.) or to a multiplication of the +brightness by some positive factor which does not affect the shape of the +object. In cases where image reconstruction has some degeneracies, these should +not have any incidence on the metric. The easier is then to minimize the metric +with respect to the degenerated parameters. These parameters include, but are +not limited to, brightness scaling, geometric transformation of coordinates, +etc. For optical interferometry and when only power-spectrum and closure phase +data are available, the images to be compared may have to be shifted for best +matching. + +When comparing a true image (with potentially an infinitely high resolution) to +a restored image, the effective resolution achievable by the instrument and the +image restoration process must be taken into account. Otherwise and because +image metrics are in general based on pixel-wise comparisons, the slightest +displacement of sharp features would lead to large loss of quality (according +to the metric) whereas the images may look very similar at a lower and more +realistic resolution. The easiest solution is then to define the reference +image to be the true image blurred by an effective point spread function (PSF) +whose shape corresponds to the effective resolution. The choice of the +effective resolution is then a parameter of the metric. + +To summarize and to be more specific, using the distance $\Dist(\V{x}, \V{y})$ +between the restored image $\V{x}$ and the reference image $\V{y}$, the discrepancy +between the restored image $\V{x}$ and the true image $\V{z}$ would be given by: + +``` math +\label{eq:general-discrepancy} +d(\V{x}, \V{z}) = \min_{\alpha,\beta,t} +\Dist\bigl(\alpha\,\V{x} + \beta, h_{\sigma,t} \ast \V{z} \bigr), +``` + +with $\alpha$ a brightness scale, $\beta$ a bias and $h_{\sigma,t}$ an +effective PSF of *width* $\sigma$ and centered at position $t$. The symbol +$\ast$ denotes the convolution. Assuming an effective PSF of Gaussian shape, +the parameter $\sigma > 0$ can be chosen to be the standard deviation of the +PSF; the full width at half maximum ($\def\FWHM{\operatorname{FWHM}}\FWHM$) is +then given by: + +``` math +\label{eq:gaussian-fwhm} +\FWHM = \sqrt{8\,\log2}\,\sigma \approx 2.355\,\sigma. +``` + +Note that the merit function could be minimized with respect to the width of +the PSF in order to estimate the effective resolution achieved by a given +restored image. Our choice to assigning the translation to the PSF is to avoid +relying on some particular method to perform sub-pixel interpolation (of $\V{x}$, +$\V{y}$ or $\V{z}$) for fine tuning the position. Not doing so would add another +ingredient to the metric. + +In the following, we first review the most common metrics found in the +literature and argue whether they are appropriate or not in the context of +optical interferometry. We then propose a family of suitable metrics. + +**Notations:** In this document, we denote scalars by Greek letters (e.g. +$\alpha$, $\beta$), collection of values, a.k.a. *vectors*, by italic bold +lower case Latin letters (e.g. $\V{x}$, $\V{y}$, $\V{z}$), and linear +operators, a.k.a. *matrices*, by upright bold upper case Latin letters (e.g., +$\M{W}$). + +[^Gomnes+2016]: N. Gomes, P. J. V. Garcia & É. Thiébaut, *Assessing the quality +of restored images in optical long-baseline interferometry* in Monthly Notices +of the Royal Astronomical Society, vol. **465**, pp. 3823-3839 (2016). + +[^Sanchez+2016]: J. Sanchez-Bermudez, É. Thiébaut, K.-H. Hofmann, M. Heininger, +D. Schertl, G. Weigelt, F. Millour, A. Schutz, A. Ferrari, M. Vannier, D. Mary, +J. Young & F. Malbet, F., *The 2016 interferometric imaging beauty contest* in +Optical and Infrared Interferometry and Imaging V, SPIE International +Conference, **9907**, 99071D (2016) [doi](https://doi.org/10.1117/12.2231982). diff --git a/notes.md b/notes.md new file mode 100644 index 0000000..c685194 --- /dev/null +++ b/notes.md @@ -0,0 +1,575 @@ + +\subsection{Merit functions} +\label{sec:merit_functions} + +\subsubsection{Quadratic metrics} +\label{sec:quadratic_metric} +The most general expression of a quadratic metric to measure the discrepancy +between two images $x$ and $y$ takes the form of a weighted (squared) +$\ell_2$-norm: +\begin{displaymath} + \wltwon(x, y; W) = \Vert x - y \Vert_{W}^2 +\end{displaymath} +where we denote by $\Vert q \Vert_{W}^2 = q\T\,W\,q$ the weighted squared +Euclidean norm with $W$ a positive (semi-)definite weighting operator. Using a +diagonal weighting operator $W = \Diag(w)$ yields: +\begin{equation} + \label{eq:quad} + %\wltwon(x, y; w) = \sum\nolimits_{i} w_i \, (x_i - y_i)^2, + \wltwon(x, y; w) = \sum_{i} w_i \, (x_i - y_i)^2, +\end{equation} +where the sum is carried out for all pixels of the images and where $w_i \ge 0$ +is the weight of the $i$-th pixel. + +By choosing specific weights, it is possible to mimic a number of commonly used +metrics. For instance, the metric of the \textit{Interferometric Imaging Beauty + Contest} \citep{Lawson2004} is: +\begin{align} + \label{eq:ibc} + \ibc(x, y) + &= \sqrt{\wltwon(x, y; w = y/{\textstyle \sum_{i} y_{i}})} \nonumber \\ + &= \left[\frac{\sum_{i} y_i \, (x_i - y_i)^2}{\sum_{i} y_i}\right]^{1/2} , +\end{align} +which amounts to taking the weights as being proportional to the reference +image: $w = y/\sum_{i} y_{i}$. One of the main drawback of this merit function +is that it discards pixels where the reference image $y$ is zero, which occurs +for many pixels for a compact astronomical source on a dark background. + +The most simple quadratic metric is the squared $\ell_2$-norm (also known as +the squared Euclidean norm) of the pixelwise differences between the images: +\begin{align} + \label{eq:l2-norm} + \ltwon (x, y) + &= \Vert x - y \Vert_2^2 \nonumber \\ + &= \sum_{i} (x_i - y_i)^2, +\end{align} +which is $\wltwon$ when $w = 1$. + +The Mean Squared Error (MSE) is directly derived from the Euclidean norm by +taking $w = 1/\Npix$ with $\Npix$ the number of pixels: +\begin{equation} + \label{eq:mse} + \mse(x, y) = \frac{1}{\Npix} \, \lVert x - y \rVert_2^2 \,. +\end{equation} +The MSE was used by \citet{Renard2011} to benchmark the effects of the +regularisation in the image reconstruction from interferometric data. For all +the metrics presented so far, the smaller the merit value, the more similar are +the images. \oops{$\Leftarrow$ this is important!} + +Some other commonly used metrics are also based on the Euclidean norm of the +differences. For instance, the Peak Signal to Noise Ratio (PSNR) is: +\begin{equation} + \label{eq:psnr} + \psnr(x, y) = 10 \times \log_{10} \! \left( + \frac{ \bigl( \max(y) - \min(y) \bigr)^2 }{ \mse(x, y) } + \right) \, . +\end{equation} +Here, $\min(y)$ and $\max(y)$ correspond respectively to the minimum and +maximum possible pixel value of the reference image $y$. The PSNR is given in +decibel (db) units and the higher the PSNR, the more similar are the images. + +Clearly, the MSE and the PSNR are the squared Euclidean norm of the pixel-wise +difference between the images ($\ltwon$) but expressed in different units. They +can be used interchangeably and we will only consider $\ltwon$ in what follows. +% \oops{(MSE or the $\ltwon$?)} + + +\subsubsection{Sum of absolute differences} +\label{sec:l1-norm} + +One of the drawbacks of quadratic metrics is that they strongly emphasize the +largest differences. To avoid this, an $\ell_p$-norm\footnote{possibly raised + to the power $p$ to simplify the equations and the computations} can be used +with an exponent $p < 2$. For instance, the \textit{sum of absolute + differences} or $\ell_1$-norm is given by: +\begin{align} + \label{eq:l1-norm} + \lonen (x, y) + &= \Vert x - y \Vert_1 \nonumber \\ + &= \sum_{i} \lvert x_i - y_i \rvert. +\end{align} + + +\subsubsection{Fidelity function} +\label{sec:fidelity} + +The fidelity function was introduced by \citet{Pety2001b} in the context of +image reconstruction for ALMA. It is defined as the ratio of the total flux of +the reference $y$ to the difference between the restored image $x$ and the +reference one: +\begin{equation} + \label{eq:fidelity} + \fid(x, y) = \frac{ + \sum_{i} y_i + }{ + \sum_{i} \max \{ \eta, \vert y_i - x_i \vert \} + } \, , +\end{equation} +where $\eta$ is some non-negative threshold. The higher the fidelity value, the +better the agreement. + +Choosing $\eta > 0$ avoids divisions by zero and \citet{Pety2001b} took +$\eta = 0.7 \, \RMS(x - y)$, where $\RMS$ is the root mean squared value: +\begin{displaymath} + \RMS(x) = \frac{1}{\Npix} \sum_{i} x_{i}^2 \, . +\end{displaymath} +We note that with $\eta > 0$, all differences smaller than $\eta$ have the same +incidence on the total cost and are therefore irrelevant. To avoid this, one +has to take $\eta = 0$, in which case the reciprocal of the fidelity function +is then just the $\ell_1$-norm defined in Eq.~(\ref{eq:l1-norm}) times some +constant factor which only depends on the reference $y$. As the fidelity +function would then yield the same results as the $\ell_1$-norm, we only +consider the latter in our study. + + +\subsubsection{\textsc{Universal Image Quality Index} and \textsc{Image Structural Similarity}} +\label{sec:quality_index} +\label{sec:ssim} + +The universal image quality index was proposed by \citet{Wang2002} to overcome +the MSE and PSNR, which were found to be very poor estimators of the image +quality for common brightness distortions and image corruptions (like +salt-and-pepper noise, lossy compression artefacts, \etc). The universal image +quality index is defined as: +\begin{equation} + \label{eq:quality_index} + \QualityIndex(x, y) = \frac{ + 4\,\AVG(x)\,\AVG(y)\,\Cov(x,y) + }{ + \bigl(\AVG(x)^2 + \AVG(y)^2\bigr)\, + \bigl(\Var(x) + \Var(y)\bigr) + } \, . +\end{equation} +where $\AVG(x)$, $\Var(x)$ and $\Cov(x,y)$ are the empirical average, +variance and covariance of $x$ and $y$ given by: +\begin{align*} + \AVG(x) &= \frac{1}{\Npix} \, \sum_{i} x_i \, , \\ + \Var(x) &= \Cov(x,x) \, , \\ + \Cov(x, y) &= \frac{1}{\Npix - 1} \, \sum_{i} + \bigl(x_i - \AVG(x)\bigr) \, \bigl(y_i - \AVG(y)\bigr) \, . +\end{align*} + +The universal image quality index takes values in the range $[-1,1]$. $Q(x,y)$ +is maximal for the best agreement which occurs when $y = \alpha\,x + \beta$ and +minimal when $y = -\alpha\,x + \beta$ for any $\alpha > 0$ and any $\beta$. +Although the universal image quality index was designed to cope with brightness +distortions such as mean shift or dynamic shrinkage, this indicator is not +exactly insensitive to any affine transform of the intensity as is (see the +demonstration in Section~\ref{sec:affine-correction}) the correlation +coefficient: +\begin{equation} + \label{eq:correlation} + \Corr(x, y) = \frac{\Cov(x,y)}{\sqrt{\Var(x) \, \Var(y)}} \, . +\end{equation} + +In order to improve over the universal image quality index, \citet{Wang2004} +introduced the image Structural SIMilarity (SSIM): +\begin{align} + \label{eq:ssim} + \SSIM(x, y) &= + \frac{% + 2 \, \AVG(x) \, \AVG(y) + \varepsilon_1 % + }{% + \AVG(x)^2 + \AVG(y)^2 + \varepsilon_1 + } \notag \\ + & \quad \times \frac{ + 2 \, \Cov(x, y) + \varepsilon_2% + }{ + \Var(x) + \Var(y) + \varepsilon_2% + } \, , +\end{align} +where $\varepsilon_1 > 0$ and $\varepsilon_2 > 0$ are small values introduced +to avoid divisions by zero. Note that with $\varepsilon_1 = 0$ and +$\varepsilon_2 = 0$, the SSIM is just the image quality index defined in +Eq.~(\ref{eq:quality_index}). The higher the SSIM, the better the agreement. In +principle SSIM and the quality index should be used *locally*, that is on +small regions of the images. + + +\subsubsection{Kullback-Leibler divergence} +\label{sec:kl} + +Being non-negative everywhere and normalised, the images can be thought as +distributions (over the pixels). The Kullback-Leibler divergence measures the +similarity between two distributions. When applied to our (normalised), images +it writes: +\begin{equation} + \label{eq:kl} + \KL(x, y) = \sum_{i} y_i \, \log(x_i / y_i) \, . +\end{equation} +A restriction for the Kullback-Leibler divergence is that $x$ and $y$ must be +strictly positive everywhere. It is however possible to account for +non-negative distributions by modifying the definition of the Kullback-Leibler +divergence as follows: +\begin{displaymath} + \KL(x, y) = \sum_{i} c_{\KL}(x_i, y_i) \, , +\end{displaymath} +where $c_{\KL}(q, r)$ extends $r\,\log(q/r)$ by continuity: +\begin{displaymath} + c_{\KL}(q, r) = \ + \begin{cases} + 0 & \text{if $q = r$ or $q > 0$ and $r = 0$,}\\ + r \, \log(q / r) & \text{otherwise.} + \end{cases} +\end{displaymath} + +Note that the Kullback-Leibler divergence is not symmetric, i.e.\ +$\KL(x,y) \not= \KL(y,x)$. The Kullback-Leibler divergence is less or equal to +zero, the lower the Kullback-Leibler divergence the worse is the agreement +between $x$ and $y$. The maximal value of the Kullback-Leibler divergence is +equal to zero and is achieved when $x = y$. + +Like the IBC metric, the Kullback-Leibler divergence disregards $x_i$ where +$y_i = 0$. In addition, any image $x$ with at least one pixel, say $i_0$, such +that $x_{i_0} = 0$ while $y_{i_0} > 0$ yields $\KL(x, y) = -\infty$ which +corresponds to the maximum possible discrepancy. These are serious drawbacks +for using the Kullback-Leibler divergence as an image metric because it could +not make a distinction between restored images such that $x_{i_0} = 0$ whatever +the values of the other pixels. + + + +\subsubsection{Minimizing the discrepancy with respect to the brightness distortion} +\label{sec:affine-correction} + +\oops{Make this part shorter (or move it in an appendix) and more didactic.} In +order to make a formal link between different metrics it is worth investigating +what happens when the minimization with respect to the brightness distortion +parameters $\alpha$ and $\beta$ is carried on. As we will show, this +minimization has a closed form solution with a quadratic metric: +\begin{displaymath} + \Vert \alpha \, x + \beta\,\One - y \Vert_{W}^2 \, , +\end{displaymath} +with $x$ and $y$ the images to compare, $\alpha \in \Reals$ a factor, +$\beta \in \Reals$ a bias and $\One$ an image where all pixels are equal to 1. + +Let us first consider the bias correction. Introducing $q = y - \alpha \, x$, +we want to minimize $\Vert q - \beta\,\One\Vert_{W}^2$ with respect to $\beta$. +Expanding the quadratic norm yields: +\begin{displaymath} + \Vert q - \beta\One\Vert_{W}^2 = \Vert q \Vert_{W}^2 + - 2\,(\One\T\,W\,q)\,\beta + (\One\T\,W\,\One)\,\beta^2 \, . +\end{displaymath} +This is a simple 2nd order polynomial in $\beta$ and the minimum is achieved +for the optimal bias: +\begin{equation} + \beta^\star = \frac{\One\T\,W\,q}{\One\T\,W\,\One} \, , +\end{equation} +which can be seen as a weighted averaging of $q$. Thus: +\begin{equation} + \min_{\beta} \Vert q - \beta\,\One\Vert_{W}^2 + = \Vert q - \beta^\star\,\One\Vert_{W}^2 + = \Vert C\,q\Vert_{W}^2 \, , +\end{equation} +where the linear operator $C$ is given by: +\begin{equation} + \label{eq:centering-operator} + C = I - \One\,\frac{\One\T\,W}{\One\T\,W\,\One} \, , +\end{equation} +where $I$ is the identity. The linear operator $C$ has the effect of removing +the weighted average of its argument. Replacing $q$ by $y - \alpha \, x$ +yields: +\begin{equation} + \min_{\beta} \Vert \alpha \, x + \beta\,\One - y \Vert_{W}^2 + = \Vert \alpha\,C\,x - C\,y\Vert_{W}^2 \, , +\end{equation} +which amounts to comparing the weighted average subtracted images. + +The expansion: +\begin{displaymath} + \Vert \alpha \, x - y \Vert_{W}^2 + = \Vert y \Vert_{W}^2 - 2\,(y\T\,W\,x)\,\alpha + \alpha^2 \, \Vert x \Vert_{W}^2 \, , +\end{displaymath} +readily shows that the optimal factor $\alpha$ is: +\begin{displaymath} + \argmin_{\alpha} \Vert \alpha \, x - y \Vert_{W}^2 + = \frac{y\T\,W\,x}{\Vert x \Vert_{W}^2} \, , +\end{displaymath} +and, after trivial simplifications, that: +\begin{displaymath} + \min_{\alpha} \Vert \alpha \, x - y \Vert_{W}^2 + = \Vert y \Vert_{W}^2 - \frac{(y\T\,W\,x)^2}{\Vert x \Vert_{W}^2} \, . +\end{displaymath} + +Putting all together we have shown that: +\begin{equation} + \min_{\alpha,\beta} \Vert \alpha \, x + \beta\,\One - y \Vert_{W}^2 + = \Vert C\,y \Vert_{W}^2 - \frac{(y\T\,C\T\,W\,C\,x)^2}{\Vert C\,x \Vert_{W}^2} \, , +\end{equation} +where the linear operator $C$ is given in Eq.~(\ref{eq:centering-operator}). If +no bias correction is wanted, it is sufficient to take $C = I$. The above +expression can be divided by $\Vert C\,y \Vert_{W}^2$ to obtain a symmetric +distance between $x$ and $y$ and which is independent of an affine transform of +the brightness of any of the two image: +\begin{equation} + d(x, y) = 1 - \Corr(x, y)^2 \, , +\end{equation} +with: +\begin{equation} + \label{eq:weighted-correlation} + \Corr(x, y) = \frac{y\T\,C\T\,W\,C\,x} + {\Vert C\,x \Vert_{W} \, \Vert C\,y \Vert_{W}} \, , +\end{equation} +the (weighted) correlation between the two images $x$ and $y$. If +$W \propto I$, then the usual definition of the correlation, given in +Eq.~(\ref{eq:correlation}), is retrieved. + +The distance $d(x, y)$ takes values in the range $[0,1]$, the smaller it is the +better is the agreement. Conversely, the better the agreement the larger the +absolute value of the (weighted) correlation. It is therefore clear now that +comparing images by mean of their (weighted) correlation coefficient is +equivalent to using a quadratic norm minimized with respect to an affine +transform of the image intensity. + + +\subsubsection{Human perception} +\label{sec:human_perception} + +\oops{Here the rationale and the results of the poll. But should be part of a + next section, say *Benchmarking the metrics*.} $\rightarrow$ \oops{See + section~\ref{sec:mu_sigma}.} + + +\subsubsection{Designing the metric} +\label{sec:designing_metric} + +We want to derive an image metric that is adapted to our particular case: we +consider images of compact objects (i.e.\ with finite size support) over a +constant background and which may be shifted by an arbitrary translation. + +We assume that $d(x,y,t)$ yields the discrepancy between the image $x$ and the +image $y$ shifted by a translation $t$. Quite naturally, we require that the +following properties hold: +\begin{itemize} +\item[(i)] the metric does not change if the images are extended with pixels + set with the background level; likewise, the metric does not change if + the images are truncated provided that the values of the removed pixels + equal the background level; +\item[(ii)] the metric is non-negative and equal to zero if the two images are + the same (for a given relative translation), in particular $d(x,x,0)=0$ + whatever the image $x$. +\item[(iii)] if translations are irrelevant, the metric must be + *stationary* in the sense that whatever the images $x$ and $y$ and + the translations $t$, $t'$ and $t''$: + \begin{equation} + d\bigl(s(x,t),s(y,t'),t''\bigr) = d(x,y,t'+t''-t) \, , + \end{equation} + where $s(x,t)$ yields image $x$ shifted by translation $t$: + \begin{displaymath} + s(x,t)_{i} = x_{i - t}, + \end{displaymath} + except that if $t$ is fractional then $s(x,t)$ implies some form of + *interpolation*. +\end{itemize} +A last requirement, although not strictly required, could be: +\begin{itemize} +\item[(iv)] the metric is *symmetric* in the sense that: + \begin{equation} + \label{eq:symmetry} + d(y,x,-t) = d(x,y,t) \, , + \end{equation} + whatever the images $x$ and $y$ and the translation $t$. +\end{itemize} + +To limit the number of possibilities, we consider that the metric is the sum of +a pixel-wise cost. Then, accounting for the Property~(i) yields: +\begin{equation} + \label{eq:additive_metric} + d(x,y,t) = \sum_{i \in \Integers^n} + c(\Extend{x}_{i}, \Extend{y}_{i - t}) +\end{equation} +where $n$ is the number of dimensions $x$ and $y$ (in the case of images, +$n=2$), $\Integers$ is the set of integers, $t \in \Integers^n$ is the +considered translation, $c(u,v)$ is the pixel-wise cost and $\Extend{x}$ +(resp.\ $\Extend{y}$) is the image $x$ (resp.\ $y$) infinitely extended with +the background level $\beta$: +\begin{equation} + \Extend{x}_{i} = + \begin{cases} + x_{i} & \text{if $i \in \Set{X}$;} \\ + \beta & \text{else;} \\ + \end{cases} +\end{equation} +where $\Set{X} \subset \Integers^n$ (resp.\ $\Set{Y} \subset \Integers^n$) is +the support of the image $x$ (resp.\ $y$). We note that the Property~(ii), +implies that $c(u,u) = 0$ whatever $u \in \Reals$ and also that the background +level must be the same for the two images. We also note that the Property~(iv) +implies that the pixel-wise cost be a symmetric function, i.e.\ +$c(u,v) = c(v,u)$ whatever $(u,v) \in \Reals^2$. Finally, the Property~(iii) +holds because the same pixel-wise cost is used whatever the index $i$. + +As $c(\beta,\beta) = 0$, the sum over the infinite set $\Integers^n$ in +Eq.~(\ref{eq:additive_metric}) simplifies to sums over three finite (and +possibly empty) subsets: +\begin{equation} + \label{eq:additive_metric_finite} + d(x,y,t) + = \sum_{\mathclap{i \in \Set{X} \cap \Set{Y}_{t}}} c(x_{i}, y_{i - t}) + + \sum_{\mathclap{i \in \Set{X} \backslash \Set{Y}_{t}}} c(x_{i}, \beta) + + \sum_{\mathclap{i \in \Set{Y} \backslash \Set{X}_{-t}}} c(y_{i}, \beta), +\end{equation} +where: +\begin{displaymath} + \Set{Y}_{t} = \{ i \in \Integers^n \given i - t \in \Set{Y} \}, +\end{displaymath} +is the support of $y$ shifted by $t$. An efficient implementation of the metric +may be achieved with: +\begin{equation} + \label{eq:additive_metric_fast} + d(x,y,t) = \gamma + \sum_{\mathclap{i \in \Set{X} \cap \Set{Y}_{t}}} + \bigl[c(x_{i}, y_{i - t}) - c(x_{i}, \beta) - c(y_{i - t}, \beta)\bigr] \, , +\end{equation} +where $c(x_{i}, \beta)$ (resp.\ $c(y_{i}, \beta)$) can be pre-computed for +all $i\in\Set{X}$ (resp.\ for all $i\in\Set{Y}$) and: +\begin{displaymath} + \gamma + = \sum_{i \in \Set{X}} c(x_{i}, \beta) + + \sum_{i \in \Set{Y}} c(y_{i}, \beta) \, . +\end{displaymath} + +\oops{It remains to choose the pixel-wise cost $c(u,v)$.} + + +\subsubsection{Choice of the candidates} +\label{sec:choice_candidates} + +\oops{We must justify this choice.} We mention all merit functions from +Section~\ref{sec:l1-norm} to Section~\ref{sec:accuracy} in order to be +exhaustive, but in our simulations we have used only Eqs.~\ref{eq:l1-norm}, +\ref{eq:ibc}, \ref{eq:l2-norm}, \ref{eq:fidelity}, and +\ref{eq:accuracy}.\footnote{For convenience, we took the reciprocal of the + fidelity function in our analysis, such that for all used merit functions, + the smaller the value, the better the agreement.} + + +\section*{Acknowledgements} +\label{sec:acknowledgments} + +\bibliography{/home/eric/work/tex/journals-short.bib, + /home/eric/work/articles/biblio.bib} +\bibliographystyle{plainnat} + +\end{document} + + +## 2016' Interferometric Imaging Beauty Contest + +For the 2016' Interferometric Imaging Beauty Contest held in Edinburgh, the +objective was to reconstruct multi-wavelength images. A metric was designed to +compare multi-channel images so as to be insensitive to irrelevant differences +due to: + +* **orientation**: image axes may be inverted, notably the E-W axis; +* **translation**: there may be an arbitrary shift between images (possibly a + different shift in every spectral channel); +* **pixel scale**: the size of the pixel used for the restored image is + arbitrary even though it should be sufficiently small to account for the + highest measured frequencies; +* **brightness**: not all entries are normalized according to the channel-wise + flux given in `OI_FLUX` data-block. + +, $x$, to a +reference image, $y$, + +In addition, the comparison must take into account that the measurements have +limited resolution. Thus the reference image $y$ is the true image $z$ +convolved with an effective PSF with a full width at half maximum (FWHM) chosen +to match the interferometric resolution. In the comparison, the restored images +are also convolved with a PSF whose FWHM is tuned to best match $y$. + +Let $R(\delta,t,\omega)$ be a linear operator which resamples an image at a +given common resolution and orientation. This operator depends on $\delta$ the +pixel size of the image to resample, $t$ the translation and $\omega$ the full +width at half maximum (FWHM) of the bur to apply to the image. All these +parameters are in the same units. Then the **reference image** is: + +``` math +y = R(\delta_{\mathrm{ref}},0,\omega_{\mathrm{ref}}) \cdot z, +``` + +where $z$ is the true image used to simulate the data, +$\delta_{\mathrm{ref}} = 3\,\mathrm{mas}/\mathrm{pixel}$ is the pixel size of the +image $z$ and $\omega_{\mathrm{ref}} \sim \lambda_{\mathrm{min}}/(2\,B_{\mathrm{max}})$ is the +FWHM of the objective resolution. We took $\omega_{\mathrm{ref}} = 6\,\mathrm{mas}$. + +The score for a given image $x$ is the sum of the squared differences between +the $\Gamma$-corrected images: + +``` math +\def\mybig#3{\bigl#1#2\bigr#3} +s(x) = \frac{ + \sum\limits_{\lambda} \min\limits_{\alpha_\lambda,t_\lambda,\omega_\lambda} + \sum\limits_{j} \left( + \Gamma\left( + \alpha_\lambda\, + \left[R(\delta,t_\lambda,\omega_\lambda)\,x_\lambda\right]_{j} + \right) + - \Gamma(y_{\lambda,j}) + \right)^2 + }{ + \sum\limits_{\lambda,j} \Gamma\bigl(y_{\lambda,j}\bigr)^2 + } +``` + +where $x_{\lambda}$ is the restored image in the spectral channel indexed by +$\lambda$, $j$ is the pixel index and $\Gamma: \Reals\to\Reals$ is a brightness +correction function to emphasizes the interesting parts of the images. We have +chosen: + +``` math +\def\Sign{\mathrm{sign}} +\Gamma(t) = \Sign(t)\,|t|^\gamma \, , +``` + +where $\Sign(t)$ is the sign of $t$ (\ie, $+1$ if $t > 0$, $-1$ if $t < 0$ and +$0$ if $t = 0$). We took $\gamma = 1$ and $\gamma = 0.7$. Note that all +restored images are nonnegative. The denominator is to normalize the score: 0 +is the best possible value and 1 is the score for a black image. The lower the +score the better. Because $\Gamma(\alpha\,t) = \Gamma(\alpha)\,\Gamma(t)$ +(whatever $t$ and $\alpha$), minimizing with respect to $\alpha_\lambda$ has a +closed form solution and the score simplifies to: + +``` math +s(x) = 1 - \frac{ + \sum\limits_{\lambda} \max\limits_{t_\lambda,\omega_\lambda} + c_\lambda(t_\lambda,\omega_\lambda) + }{ + \sum\limits_{\lambda,j} \Gamma\Paren[\big]{y_{\lambda,j}}^2 + } \, , +``` + +with: + +``` math + c_\lambda(t,\omega) = \frac{ + \Brack[\Big]{ + \sum_{j} + \Gamma\Paren[\big]{y_{\lambda,j}}\, + \Gamma\Paren[\big]{\Brack*{R(\delta,t,\omega)\,x_\lambda}_{j}} + } + }{ + \sum_{j} + \Gamma\Paren[\big]{\Brack*{R(\delta,t,\omega)\,x_\lambda}_{j}}^2 + } \, , +``` + +which is a normalized correlation between the $\Gamma$-corrected images. + + +To compensate for different pixel sizes, the image $x$ or $x$ which has the +larger pixel size is interpolated so that both images have the same (smallest) +pixel size. + +Separable linear interpolation with a triangle kernel is applied for magnifying +and fine shifting the images. + +The criterion is minimized for a translation between each images (for each +spectral channel). + +In every spectral channel, the brightness of the restored images is scaled +so that the total flux per channel is the same as in the reference image. + +Because the resolution of the reference image may change (see above) the score +is the ratio between the sum for all spectral channels of the scores between +the restored images and the reference images divided by the sum for all +spectral channels of the scores between a zero image and the reference images. diff --git a/old/beauty2016.jl b/old/beauty2016.jl new file mode 100644 index 0000000..7974acc --- /dev/null +++ b/old/beauty2016.jl @@ -0,0 +1,793 @@ + +if ! isdefined(:arcsec) + const arcsec = pi/(180*3600) +end + +if ! isdefined(:mas) + const mas = arcsec/1E3 +end + +if ! isdefined(:plt) + using PyCall + pygui(:gtk); # can be: :wx, :gtk, :qt + using LaTeXStrings; + import PyPlot + const plt = PyPlot; +end +fma() = (plt.cfl(); nothing) +fma(i::Integer) = (plt.figure(i); plt.cfl(); nothing) +function fig(i::Integer, clear::Bool=false) + plt.figure(i) + if clear + plt.clf() + end + nothing +end +function pli{T<:Real}(a::AbstractArray{T,2}; + cmin::Real=minimum(a), + cmax::Real=maximum(a), + interp::Symbol=:nearest) + plt.imshow(a, interpolation=interp, vmin=cmin, vmax=cmax, + origin=:lower) + nothing +end + + +module Test + +import Base: size + +include("metrics/metrics.jl") +import .Metrics: absdif, compare_images, resample +using TiPi +using OptimPack +using FITSIO +import PyPlot +const plt = PyPlot + + +const arcsec = pi/(180*3600) +const mas = arcsec/1E3 + +typealias Table Dict{Symbol,Any} + +type Cube + name::ASCIIString + pixelsize::Float64 + data::Array{Float64,3} + flux::Array{Float64,1} +end +size(c::Cube) = size(c.data) +size(c::Cube, i::Integer) = size(c.data, i) +pixelsize(c::Cube) = c.pixelsize +data(c::Cube) = c.data + + +#function waverange(first::Real, step::Real, number::Integer) +# +#end + +TrueObj2v2 = Table(:name => "truth", + :flipx => false, + :flipy => false, + :pixelsize => 3.0mas, + :dirname => "../dustydisk/", + :filenames => ("dustydisk-Nband-v2.fits",)) + +TrueObj2v3 = Table(:name => "truth", + :flipx => false, + :flipy => false, + :pixelsize => 3.0mas, + :dirname => "../dustydisk/", + :filenames => ("dustydisk-Nband-v3.fits",)) + +YoungObj2 = Table(:name => "Young", + :flipx => true, + :flipy => false, + :pixelsize => 2.0mas, + :dirname => "Young/", + :filenames => ("bsmem_Object2-v3.fits.gz",)) + +HofmannObj2 = Table(:name => "Hofmann", + :flipx => true, + :flipy => false, + :pixelsize => 3.125mas, + :dirname => "Hofmann/Object2.8h.E.1.3DCubus/", + :filenames => ("Object2_08.00mu.fits", + "Object2_08.38mu.fits", + "Object2_08.77mu.fits", + "Object2_09.15mu.fits", + "Object2_09.54mu.fits", + "Object2_09.92mu.fits", + "Object2_10.31mu.fits", + "Object2_10.69mu.fits", + "Object2_11.08mu.fits", + "Object2_11.46mu.fits", + "Object2_11.85mu.fits", + "Object2_12.23mu.fits", + "Object2_12.62mu.fits", + "Object2_13.00mu.fits")) + +SchutzObj2 = Table(:name => "Schutz et al.", + :flipx => false, + :flipy => false, + :pixelsize => 3.125mas, + :dirname => "Schutz/", + :filenames => ("Obj2_schutzetal.fits.gz",)) + +MillourObj2a = Table(:name => "Millour (a)", + :flipx => false, + :flipy => false, + :pixelsize => 6.0mas, + :dirname => "Millour/", + :filenames => ("Object2_combinImage_Median.fits.gz",)) + +MillourObj2b = Table(:name => "Millour (b)", + :flipx => false, + :flipy => false, + :pixelsize => 1000mas/255, # see text + :dirname => "Millour/", + :filenames => ("Object2_modelImagesCube.fits.gz",)) + +# Full width at half maximum of a cubic B-spline: +const FWHM = 4/3 * (1 + 2*cos((pi + atan(3*sqrt(7)))/3)) + +function resampleimage{T1<:Real,T2<:Real,T<:AbstractFloat}(x1::AbstractArray{T1,1}, + x2::AbstractArray{T2,1}, + A::Array{T,2}, + g1::Range, + g2::Range; + border::Union{Void,Real}=nothing, + fwhm::Real=3mas) + @assert length(g1) == size(A,1) + @assert length(g2) == size(A,2) + ker = TiPi.Kernels.CubicKernel(Float64) + @assert length(ker) == 4 + + @assert fwhm > 0 + const alpha = Float64(FWHM/fwhm) + + + # zero pad the image? + if border !== nothing + A = pad(T(border), A, (size(A,1) + 2, size(A,2) + 2)) + g1 = enlargerange(g1) + g2 = enlargerange(g2) + end + + # Resample the image. + op1 = Resampler(ker, alpha*x1, alpha*g1, quiet=true) + op2 = Resampler(ker, alpha*x2, alpha*g2, quiet=true) + return resample(op1, resample(op2, A, 2), 1) +end + +function correlateimages{T}(A0::AbstractArray{T,2}, A1::AbstractArray{T,2}, gamma::Real=1) + @assert size(A1) == size(A0) + gamma = Float64(gamma) + a::Float64 = 0.0 + b::Float64 = 0.0 + if gamma == 1.0 + @inbounds @simd for i in eachindex(A0,A1) + v0 = Float64(A0[i]) + v1 = Float64(A1[i]) + a += v1*v1 + b += v0*v1 + end + else + @inbounds @simd for i in eachindex(A0,A1) + v0 = Float64(A0[i])^gamma + v1 = Float64(A1[i])^gamma + a += v1*v1 + b += v0*v1 + end + end + return b*b/a +end + +function enlargerange(range::Range) + s = step(range) + linspace(first(range) - s, last(range) + s, length(range) + 2) +end + +function makerange(step::Real, i0::Integer, len::Integer) + linspace(step*(1 - i0), step*(len - i0), len) +end + +immutable Resampler{T<:AbstractFloat,M} + nrows::Int + ncols::Int + width::Int + idx::Array{Int,2} + wgt::Array{T,2} +end + +import TiPi.Kernels.Kernel + + + +doc""" + + Resampler(ker, pos, grd) + +yields a linear interpolator which is suitable to interpolate or resample +an array along a given dimension. `ker` is the interpolation kernel, `pos` +specifies the coordinates where to interpolate and `grd` specifies the +sampled coordinates. + +In 1-D (for the sake of clarity), the resampling writes: + + dst[i] = sum_j ker(grd[j] - pos[i])*src[j] + +where `dst` is the resulting array interpolated/resampled at coordinates +`pos` and `src` gives the values of a function sampled on the regular grid +`grd`. + +""" +function Resampler{T,T1<:Real}(ker::Kernel{T}, + pos::AbstractArray{T1,1}, grd::Range; + quiet::Bool=false) + + # All computations done in double precision. + const Flt = Float64 + + for x in pos + isfinite(x) || error("invalid position ($x)") + end + + # Determine the size of the kernel support in grid units. + d::Flt = step(grd) # grid increment + w::Flt = length(ker)/abs(d) # kernel width in grid index units + M::Int = floor(Int, w) # max. number of non-zero coefficients + if M < length(ker) + if ! quiet + warn("too few neighbors per interpolated values, will use Catmull-Rom cubic interpolation") + end + ker = TiPi.Kernels.CatmullRomKernel(Flt) + grd ./= d # = 1:length(grd) + pos ./= d + d = 1.0 + M = length(ker) + w = M + for x in pos + isfinite(x) || error("invalid position ($x)") + end + end + + # Number of rows and columns of the corresponding linear operator. + const nrows = length(pos) + const ncols = length(grd) + + # Coefficients of the reverse interpolation formula to convert + # coordinates in (fractional) grid index. + const j0 = 1 # first grid index + const j1 = length(grd) # last grid index + const t0 = Flt(first(grd)) # first grid coordinate + const t1 = Flt(last(grd)) # last grid coordinate + if (t1 - t0) == zero(Flt) + error("input grid must have two distinct endpoints") + end + const f0 = Flt(j0/(t1 - t0)) + const f1 = Flt(j1/(t1 - t0)) + + # Extrapolation is done by using the nearest position in the grid. In + # infinite precision, there would be no needs to take care of this, but + # when the fractional index `f` is outside the range `[fmin, fmax]`, we + # treat that specially to avoid `InexactError` exceptions. + const h = Flt(w/2) + const fmin = Flt(j0 - h) + const fmax = Flt(j1 + h) + const q = h + one(Flt) # offset to edge of support plus one + idx = Array(Int, (M, nrows)) + wgt = Array(T, (M, nrows)) + for i in 1:nrows + # Apply reverse formula to get + t::Flt = pos[i] # interpolation coordinate + f::Flt = f0*(t1 - t) + f1*(t - t0) # corresponding fractional grid index + if f <= fmin + f -= ceil(f - q) + for k in 1:M + idx[k,i] = j0 + wgt[k,i] = ker((f - Flt(k))*d) + end + elseif f >= fmax + f -= ceil(f - q) + for k in 1:M + idx[k,i] = j1 + wgt[k,i] = ker((f - Flt(k))*d) + end + else + l = try + ceil(Int, f - q) + catch + warn("f=$f, q=$q, t=$t, f0=$f0, f1=$f1, t1=$t1, t0=$t0") + unsafe_trunc(Int, f - q) + end + f -= l + for k in 1:M + idx[k,i] = min(max(1, l + k), ncols) + wgt[k,i] = ker((f - Flt(k))*d) + end + end + end + return Resampler{T,M}(nrows, ncols, M, idx, wgt) +end + +function check{T,M}(op::Resampler{T,M}) + nrows, ncols, idx, wgt = op.nrows, op.ncols, op.idx, op.wgt + if op.width != M + throw(ArgumentError("bad width")) + end + if nrows < 1 + throw(ArgumentError("bad output dimension")) + end + if ncols < 1 + throw(ArgumentError("bad input dimension")) + end + if size(idx) != (M, nrows) + throw(ArgumentError("bad dimensions of indexes")) + end + if size(wgt) != (M, nrows) + throw(ArgumentError("bad dimensions of weights")) + end + for j in idx + 1 ≤ j ≤ ncols || throw(BoundsError()) + end +end + +function resample{T,M,N}(op::Resampler{T,M}, + src::AbstractArray{T,N}, + dim::Int) + resample!(Array(T, ntuple(d -> (d == dim ? op.nrows : size(src, d)), N)), + op, src, dim) +end + +function resample!{T,M,N}(dst::AbstractArray{T,N}, + op::Resampler{T,M}, + src::AbstractArray{T,N}, + dim::Int) + check(op) + nrows, ncols, idx, wgt = op.nrows, op.ncols, op.idx, op.wgt + @assert 1 ≤ dim ≤ N + if size(src, dim) != ncols + error("bad source dimension") + end + for d in 1:N + if size(dst, d) != (d == dim ? nrows : size(src, d)) + error("bad destination dimension") + end + end + @inbounds begin + if N == 1 + for i in 1:nrows + s = zero(T) + @simd for k in 1:M + j = idx[k,i] + s += wgt[k,i]*src[j] + end + dst[i] = s + end + elseif dim == 1 + rng1 = CartesianRange(size(src)[dim+1:end]) + for i1 in rng1 + for i in 1:nrows + s = zero(T) + @simd for k in 1:M + j = idx[k,i] + s += wgt[k,i]*src[j,i1] + end + dst[i,i1] = s + end + end + elseif dim == N + rng0 = CartesianRange(size(src)[1:dim-1]) + for i in 1:nrows + for i0 in rng0 + s = zero(T) + @simd for k in 1:M + j = idx[k,i] + s += wgt[k,i]*src[i0,j] + end + dst[i0,i] = s + end + end + else + rng0 = CartesianRange(size(src)[1:dim-1]) + rng1 = CartesianRange(size(src)[dim+1:end]) + for i1 in rng1 + for i in 1:nrows + for i0 in rng0 + s = zero(T) + @simd for k in 1:M + j = idx[k,i] + s += wgt[k,i]*src[i0,j,i1] + end + dst[i0,i,i1] = s + end + end + end + end + end + return dst +end + +@noinline function sumalongdims!{T1,T2,N}(dst::AbstractArray{T1,N}, + src::AbstractArray{T2,N}) + # It's assumed that DST has size 1 along any dimension that we're summing + fill!(dst, 0) + srcmax = CartesianIndex(size(dst)) + @inbounds for i in CartesianRange(size(src)) + dst[min(srcmax,i)] += src[i] + end + dst +end + +function sumalongdims(A::AbstractArray, dims) + sz = [size(A)...] + sz[[dims...]] = 1 + B = Array(eltype(A), sz...) + sumalongdims!(B, A) +end + +function collapse{T,N}(A::AbstractArray{T,N}) + sz = Array(Int, N) + j = 0 + for i in 1:N + dim = size(A,i) + if dim > 1 + j += 1 + sz[j] = dim + end + end + j > 0 && return reshape(A, sz[1:j]...) + for i in eachindex(A) + return A[i] + end +end + +function spectrum{T<:Real}(arr::Array{T,3}) + range = CartesianRange(size(arr)[1:2]) + result = Array(Float64, size(arr,3)) + @inbounds begin + for j in 1:size(arr,3) + s::Float64 = 0.0 + @simd for i in range + s += arr[i,j] + end + result[j] = s + end + end + return result +end + +function loadcube(inp::Table; verb::Bool=false) + n = length(inp[:filenames]) + temp = Array(Any, n) + dim1 = -1 + dim2 = -1 + dim3 = -1 + data = Array(Float64,0,0,0) + for i in 1:n + verb && println("reading: ",inp[:dirname]*inp[:filenames][i]) + f = FITS(inp[:dirname]*inp[:filenames][i]) + a = read(f[1]) + rank = ndims(a) + verb && println("size: ", size(a)) + if i == 1 + @assert 2 <= rank <= 3 + dim1 = size(a,1) + dim2 = size(a,2) + dim3 = (rank >= 3 ? size(a,3) : n) + else + @assert rank == 2 + @assert dim1 == size(a,1) + @assert dim2 == size(a,2) + end + if eltype(a) != Float64 + a = map(Float64, a) + end + if i == 1 + # Allocate cube. + data = (rank == 3 ? a : Array(Float64, dim1, dim2, dim3)) + end + if rank == 2 + data[:,:,i] = a + end + end + if inp[:flipx] + data = data[end:-1:1,:,:] + end + if inp[:flipy] + data = data[:,end:-1:1,:] + end + return Cube(inp[:name], inp[:pixelsize], data, spectrum(data)) +end + +function resamplecube(c::Cube, mag::Real, align::Symbol=:center) + nchannels = size(c,3) + data = Array(Float64,0,0,0) + for i in 1:nchannels + a = resample(c.data[:,:,i], mag, align) + if i == 1 + data = Array(Float64, size(a)..., nchannels) + end + data[:,:,i] = a + end + return Cube(c.name, c.pixelsize/mag, data, spectrum(data)) +end + +function findmaximum{T<:Real,N}(a::AbstractArray{T,N}) + range = CartesianRange(size(a)) + imax = first(range) + amax = a[imax] + @inbounds for i in range + aval = a[i] + if aval > amax + amax = aval + imax = i + end + end + return (amax, imax) +end + +function findminimum{T<:Real,N}(a::AbstractArray{T,N}) + range = CartesianRange(size(a)) + imin = first(range) + amin = a[imin] + @inbounds for i in range + aval = a[i] + if aval < amin + amin = aval + imin = i + end + end + return (amin, imin) +end + +function findextrema{T<:Real,N}(a::AbstractArray{T,N}) + range = CartesianRange(size(a)) + imin = first(range) + amin = a[imin] + imax = imin + amax = amin + @inbounds for i in range + aval = a[i] + if aval < amin + amin = aval + imin = i + elseif aval > amax + amax = aval + imax = i + end + end + return (amin, imin, amax, imax) +end + +doc""" + + comparecube(c0, c1) + +return a map of the error for different translations. `c0` is the reference cube. +""" +function comparecube(cub0::Cube, cub1::Cube; + gamma::Real=0.5, # brightness correction + fwhm::Real=6mas, # reference blur + delta::Real=2mas, # sampling step for images + omega::Real=400mas) # width of field of view + + # FIXME: interpolate wavelength + @assert size(cub0,3) == size(cub1,3) + + const debugging = false + const plotting = true + + fwhm::Float64 = fwhm + delta::Float64 = delta + omega::Float64 = omega + gamma::Float64 = gamma + if gamma == 0.5 + Gamma(a::AbstractArray) = map(t -> sign(t)*sqrt(abs(t)), a) + elseif gamma == 1.0 + Gamma(a::AbstractArray) = a + elseif gamma > 0.0 + Gamma(a::AbstractArray) = map(t -> sign(t)*abs(t)^gamma, a) + else + error("invalid value for Gamma-correction") + end + + n = round(Int, omega/delta/2) + fov_x = makerange(delta, n+1, 2*n+1) + fov_y = makerange(delta, n+1, 2*n+1) + + # Start by aligning the maxima. + ref = Array(Float64, 0, 0, 0) + dat = Array(Float64, 0, 0, 0) + sumcorrel = 0.0 + nchannels = size(cub0,3) + for j in 1:nchannels + a0 = WorkImage(cub0.data[:,:,j], pixelsize(cub0)) + a1 = WorkImage(cub1.data[:,:,j], pixelsize(cub1)) + + # Reference image: + r0 = resample(a0, fov_x, fov_y, fwhm) + r0c = Gamma(r0) + if plotting + Main.fig(1,true) + Main.pli(r0c) + end + c0 = sum(r0c.*r0c) + + best_c::Float64 = -Inf + best_tx::Float64 = 0.0 + best_ty::Float64 = 0.0 + best_w::Float64 = fwhm + best_img::Array{Float64,2} = Array(Float64, 0, 0) + scale = [mas, mas, mas] + param = [best_tx/mas, best_ty/mas, best_w/mas] + function cost(a::Array{Float64,1}) + @assert length(a) == 3 + if ! isfinite(a[1]) || ! isfinite(a[2]) || ! isfinite(a[3]) + warn("tx = $(a[1]), ty = $(a[2]), w = $(a[3])") + return 0.0 + end + tx = a[1]*scale[1] + ty = a[2]*scale[2] + w = max(0.0, a[3])*scale[3] + r1 = resample(a1, fov_x - tx, fov_y - ty, w) + r1c = Gamma(r1) + c = correlateimages(r0c, r1c)/c0 + if c > best_c + best_c = c + best_tx = tx + best_ty = ty + best_w = w + best_img = r1 + if plotting + Main.fig(2,true) + Main.pli(r1c, cmin=0) + end + if debugging + @printf(" - score = %9.7f, t = (%6.3f,%6.3f) mas, FWHM = %6.3f mas\n", + 1.0 - best_c, (best_tx/mas), (best_ty/mas), (best_w/mas)) + end + end + return -c + end + + if debugging + const n = 101 + const i0 = div(n,2) + 1 + const s = 0.1*delta + w = fwhm + crt = Array(Float64,n,n) + c = -cost(param) + for i2 in 1:n + ty = (i2 - i0)*s + for i1 in 1:n + tx = (i1 - i0)*s + param = [tx, ty, w]./scale + crt[i1,i2] = -cost(param) + end + end + if plotting + Main.fig(3,true) + Main.pli(crt, cmin=0) + end + else + rslt = OptimPack.newuoa!(cost, param, 2*length(param)+1, 1.0, 1e-4; + verbose=0, maxeval=500) + @printf(" - spectral channel = %2d, score = %9.7f, t = (%6.3f,%6.3f) mas, FWHM = %6.3f mas\n", + j, 1.0 - best_c, (best_tx/mas), (best_ty/mas), (best_w/mas)) + end + sumcorrel += best_c + r1 = best_img + if j == 1 + dat = Array(Float64, size(r1)..., nchannels) + ref = Array(Float64, size(r0)..., nchannels) + end + dat[:,:,j] = r1 + ref[:,:,j] = r0 + end + return 1.0 - (sumcorrel/nchannels), dat, ref +end + +type WorkImage{T<:AbstractFloat} + data::AbstractArray{T,2} + x::AbstractArray{T,1} + y::AbstractArray{T,1} +end + +function WorkImage{T<:AbstractFloat}(a::AbstractArray{T,2}, pixelsize::Real) + amax, imax = findmaximum(a) + x = makerange(T(pixelsize), imax[1], size(a,1)) + y = makerange(T(pixelsize), imax[2], size(a,2)) + return WorkImage{T}(a, x, y) +end + +function resample(img::WorkImage, x, y, w) + resampleimage(x, y, img.data, img.x, img.y; border=0, fwhm=w) +end + +function runtests() + flag = true + for (name, filename, ref, inp) in + (("Florentin Millour (a)", "FlorentinMillourA.fits", trueobj2v3, millourobj2a), + ("Florentin Millour (b)", "FlorentinMillourB.fits", trueobj2v3, millourobj2b), + ("Karl-Heinz Hoffman", "KarlHeinzHoffman.fits", trueobj2v2, hofmannobj2), + ("John Young", "JohnYoung.fits", trueobj2v3, youngobj2), + ("Antony Schutz at al.", "AntonySchutz.fits", trueobj2v3, schutzobj2),) + + println("\n$(name):"); + score, dat, ref = comparecube(ref, inp, + gamma=0.7, # brightness correction + fwhm=6mas, # reference blur + delta=2mas, # sampling step for images + omega=400mas) # width of field of view + write(FITS(filename, "w"), map(Float32, dat)) + if flag + write(FITS("Truth.fits", "w"), map(Float32, ref)) + flag = false + end + @printf(" - final score: %10.8f\n", score) + write(FITS("/tmp/"*filename, "w"), map(Float32, dat)) + end +end + +doc""" + + comparecubeold(c0, c1) + +return a map of the error for different translations. `c0` is the reference cube. +""" +function comparecubeold(c0::Cube, c1::Cube) + # FIXME: interpolate wavelength + @assert size(c0,3) == size(c1,3) + if pixelsize(c0) > pixelsize(c1) + println("resampling C0") + c0 = resamplecube(c0, pixelsize(c0)/pixelsize(c1)) + elseif pixelsize(c0) < pixelsize(c1) + println("resampling C1") + c1 = resamplecube(c1, pixelsize(c1)/pixelsize(c0)) + end + + # Fix flux. + nchannels = size(c0,3) + result = Array(Float64,0,0,0) + temp = Array(Float64, size(c1.data)[1:2]) + for j in 1:nchannels + @assert isfinite(c0.flux[j]) && c0.flux[j] > 0.0 + @assert isfinite(c1.flux[j]) && c1.flux[j] > 0.0 + alpha = c0.flux[j]/c1.flux[j] + @assert isfinite(alpha) + @inbounds @simd for i in CartesianRange(size(temp)) + temp[i] = c1.data[i,j]*alpha + end + plt.clf() + plt.imshow(temp, vmin=0) + a = compare_images(absdif, c0.data[:,:,j], temp, 0.0) + if j == 1 + result = Array(Float64, size(a)..., nchannels) + end + range = CartesianRange(size(a)) + imin = first(range) + amin = a[imin] + for i in range + if a[i] < amin + imin = i + amin = a[imin] + end + end + @printf("%2d alpha = %.5e amin = %.5e imin = %s\n", j, alpha, amin, string(imin)) + result[:,:,j] = a + end + return result +end + +trueobj2v2 = loadcube(TrueObj2v2); +trueobj2v3 = loadcube(TrueObj2v3); +hofmannobj2 = loadcube(HofmannObj2); +youngobj2 = loadcube(YoungObj2); +schutzobj2 = loadcube(SchutzObj2); +millourobj2a = loadcube(MillourObj2a); +millourobj2b = loadcube(MillourObj2b); + +end # module diff --git a/old/metrics.jl b/old/metrics.jl new file mode 100644 index 0000000..a47761a --- /dev/null +++ b/old/metrics.jl @@ -0,0 +1,380 @@ +module Metrics + +using TiPi + +export interpolate + +# FIXME: @assert n >= 2 +immutable LinearInterpolationCoefficients{T<:AbstractFloat} + i0::Int + i1::Int + w0::T + w1::T +end + +function LinearInterpolationCoefficients(vals::AbstractArray{Float64,1}, dim::Int) + c = Array(LinearInterpolationCoefficients{Float64}, length(vals)) + k = 0 + @inbounds for x in vals + k += 1 + c[k] = LinearInterpolationCoefficients{Float64}(interpcoefs(x, dim)...) + end + return c +end + + +# immutable SeparableInterpolator{T,N} <: LinearOperator{AbstractArray{T,N}, +# AbstractArray{T,N}} +# ops::Vector{Interpolator{T,N}} +# inpdims::NTuple{N,Int} +# outdims::NTuple{N,Int} +# end +# +# eltype{T,N}(op::SeparableInterpolator{T,N}) = T +# input_size{T,N}(op::SeparableInterpolator{T,N}) = op.inpdims +# output_size{T,N}(op::SeparableInterpolator{T,N}) = op.outdims +# +# function apply_direct{T,N}(A::SeparableInterpolator{T,N}, +# src::AbstractArray{T,N}) +# dst = Array(T, output_size(op)) +# apply_direct!(dst, A, src) +# return dst +# end +# +# function apply_direct!{T,N}(dst::AbstractArray{T,N}, +# A::SeparableInterpolator{T,N}, +# src::AbstractArray{T,N}) +# @assert size(src) == input_size(A) +# @assert size(dst) == output_size(A) +# +# for i in 1:N +# end +# end +# +# function apply_adjoint{T,N}(A::SeparableInterpolator{T,N}, +# src::AbstractArray{T,N}) +# dst = Array(T, input_size(A)) +# apply_adjoint!(dst, A, src) +# return dst +# end +# +# function apply_adjoint!{T,N}(dst::AbstractArray{T,N}, +# op::Interpolator{T,N}, +# src::AbstractArray{T,N}) +# @assert size(src) == output_size(A) +# @assert size(dst) == input_size(A) +# end + + +function interpcoefs{T<:AbstractFloat}(x::T, n::Int) + if x <= one(T) + #return LinearInterpolationCoefficients{T}(1, 1, zero(T), one(T)) + return (1, 1, zero(T), one(T)) + elseif x >= T(n) + #return LinearInterpolationCoefficients{T}(n, n, one(T), zero(T)) + return (n, n, one(T), zero(T)) + else + f = floor(x) + i = Int(f) + u = x - f + #return LinearInterpolationCoefficients{T}(i, i + 1, one(T) - u, u) + return (i, i + 1, one(T) - u, u) + end +end + +# FIXME: make this work for complex-valued arrays +# FIXME: the 2D interpolation is still ~40% slower than Yorick +# FIXME: multi-dimensionalseparable interpolation is faster (and more generic) +# when done one dimension at a time (ordering of dimensions is important) + +function interpolate{T<:Real}(src::AbstractArray{T,1}, + xvals::AbstractArray{Float64,1}) + const dim1 = size(src, 1) + dst = Array(Float64, length(xvals)) + k = 0 + for x in xvals + (i0, i1, u0, u1) = interpcoefs(x, dim1) + k += 1 + dst[k] = src[i0]*u0 + src[i1]*u1 + end + return dst +end + +function interpolate{T<:Real}(src::AbstractArray{T,2}, + xvals::AbstractArray{Float64,1}, + yvals::AbstractArray{Float64,1}; + slow::Bool=false) + const dim1 = size(src, 1) + const dim2 = size(src, 2) + dst = Array(Float64, (length(xvals), length(yvals))) + if slow + k = 0 + for y in yvals + (j0, j1, v0, v1) = interpcoefs(y, dim2) + for x in xvals + (i0, i1, u0, u1) = interpcoefs(x, dim1) + k += 1 + dst[k] = ((src[i0,j0]*u0 + src[i1,j0]*u1)*v0 + + (src[i0,j1]*u0 + src[i1,j1]*u1)*v1) + end + end + else + c1 = LinearInterpolationCoefficients(xvals, dim1) + k = 0 + for y in yvals + (j0, j1, v0, v1) = interpcoefs(y, dim2) + for i in 1:length(c1) + i0 = c1[i].i0 + i1 = c1[i].i1 + u0 = c1[i].w0 + u1 = c1[i].w1 + k += 1 + dst[k] = ((src[i0,j0]*u0 + src[i1,j0]*u1)*v0 + + (src[i0,j1]*u0 + src[i1,j1]*u1)*v1) + end + end + end + return dst +end + +function interpolate{T<:Real}(src::AbstractArray{T,3}, + xvals::AbstractArray{Float64,1}, + yvals::AbstractArray{Float64,1}, + zvals::AbstractArray{Float64,1}; + slow::Bool=false) + const dim1 = size(src, 1) + const dim2 = size(src, 2) + const dim3 = size(src, 3) + dst = Array(Float64, (length(xvals), length(yvals), length(zvals))) + if slow + k = 0 + for z in zvals + (k0, k1, w0, w1) = interpcoefs(z, dim3) + for y in yvals + (j0, j1, v0, v1) = interpcoefs(y, dim2) + for x in xvals + (i0, i1, u0, u1) = interpcoefs(x, dim1) + k += 1 + dst[k] = ((((src[i0,j0,k0]*u0 + src[i1,j0,k0]*u1)*v0 + + (src[i0,j1,k0]*u0 + src[i1,j1,k0]*u1)*v1))*w0 + + (((src[i0,j0,k1]*u0 + src[i1,j0,k1]*u1)*v0 + + (src[i0,j1,k1]*u0 + src[i1,j1,k1]*u1)*v1))*w1) + end + end + end + else + c1 = LinearInterpolationCoefficients(xvals, dim1) + c2 = LinearInterpolationCoefficients(yvals, dim2) + k = 0 + for z in zvals + (k0, k1, w0, w1) = interpcoefs(z, dim3) + for i in 1:length(c2) + j0 = c2[i].i0 + j1 = c2[i].i1 + v0 = c2[i].w0 + v1 = c2[i].w1 + for i in 1:length(c1) + i0 = c1[i].i0 + i1 = c1[i].i1 + u0 = c1[i].w0 + u1 = c1[i].w1 + k += 1 + dst[k] = ((((src[i0,j0,k0]*u0 + src[i1,j0,k0]*u1)*v0 + + (src[i0,j1,k0]*u0 + src[i1,j1,k0]*u1)*v1))*w0 + + (((src[i0,j0,k1]*u0 + src[i1,j0,k1]*u1)*v0 + + (src[i0,j1,k1]*u0 + src[i1,j1,k1]*u1)*v1))*w1) + end + end + end + end + return dst +end + +type DiscrepancyMap{N} + d::Array{Float64,N} + r::NTuple{N,UnitRange{Int}} +end + +center(len::Integer) = div(len,2) + 1 + +function resample_range(mag::Real, dim::Integer, align::Symbol) + mag >= 1 || throw(ArgumentError("magnification should be at least one")) + newdim = ceil(Int, mag*dim) + + if align == :center + # x' = alpha*i + beta with alpha = 1/magn + # alpha*center(n') + beta = center(n) + # beta = center(n) - alpha*center(n') + alpha = one(Float64)/Float64(mag) + beta = center(dim) - alpha*center(newdim) + a = alpha + beta + b = alpha*newdim + beta + elseif align == :middle + margin = ((newdim - 1)/mag - (dim - 1))/2 + a = Float64(1 - margin) + b = Float64(dim + margin) + elseif align == :first + a = one(Float64) + b = Float64(newdim - 1)/Float64(mag) + a + elseif align == :last + b = Float64(dim) + a = b - Float64(newdim - 1)/Float64(mag) + else + throw(ArgumentError("invalid alignment")) + end + linspace(a, b, newdim) +end + +function resample{T<:Real}(src::AbstractArray{T,1}, mag::Real, + align::Symbol=:center) + interpolate(src, + resample_range(mag, size(src, 1), align)) +end + +function resample{T<:Real}(src::AbstractArray{T,2}, mag::Real, + align::Symbol=:center) + interpolate(src, + resample_range(mag, size(src, 1), align), + resample_range(mag, size(src, 2), align)) +end + +function resample{T<:Real}(src::AbstractArray{T,3}, mag::Real, + align::Symbol=:center) + interpolate(src, + resample_range(mag, size(src, 1), align), + resample_range(mag, size(src, 2), align), + resample_range(mag, size(src, 3), align)) +end + +""" + resample(A, mag) + +yields array `A` resampled with magnification `mag` using linear interpolation. +The result has dimensions `ceil(mag*size(A))`, so there may be some extra +margin to fill. The alignment can be specified as follows: + + resample(A, mag, :first) + resample(A, mag, :last) + resample(A, mag, :middle) + resample(A, mag, :center) + +to align the result at the first element of each dimension, at the last element +of each dimension or to geometrically center the result exactly at the middle +of the interval or at the "center" as defined by the conventions of the +`fftshift` function. For a dimension of length `dim` this is given by the +index: + + icen = div(dim,2) + 1 + +The default alignment is `:center`. + + """ resample + +function compare_images{T}(c::Function, + a::AbstractArray{T,2}, + b::AbstractArray{T,2}, + beta) + xrange = -size(b, 1) : size(a, 1) + yrange = -size(b, 2) : size(a, 2) + compare_images(c, a, b, T(beta), xrange, yrange) +end + +function compare_images{T<:Real}(c::Function, + a::AbstractArray{T,2}, + b::AbstractArray{T,2}, + beta::Real, + xrange::UnitRange{Integer}, + yrange::UnitRange{Integer}) + compare_images(c, a, b, T(beta), + UnitRange{Int}(xrange), UnitRange{Int}(yrange)) +end + +function compare_images{T}(c::Function, + a::AbstractArray{T,2}, + b::AbstractArray{T,2}, + beta::T, + xrange::UnitRange{Int}, + yrange::UnitRange{Int}) + # Build the discrepancy maps with respect to the background level. + (ca, sa) = c(a, beta) + (cb, sb) = c(b, beta) + gamma = sa + sb # maximal discrepancy + + # Compute the discrepancy map for all requested translations. + result = Array(Float64, length(xrange), length(yrange)) + iy = 0 + for ty in yrange + iy += 1 + ix = 0 + for tx in xrange + ix += 1 + result[ix,iy] = gamma + c(a, ca, b, cb, tx, ty) + end + end + return result +end + +absdif(a::Float64, b::Float64) = abs(a - b) + +abs2dif(a::Float64, b::Float64) = (d = a - b; d*d) + +# Metric used for the former Interferometric Imaging Beauty Contests. +# A is the reconstructed image, B is the reference image (which should be +# normalized but this is not essential). +iibc(a::Float64, b::Float64) = (d = a - b; b*d*d) + +kullback_leibler(a::Float64, b::Float64) = + (a == b || (a > zero(Float64) && b == zero(Float64)) + ? zero(Float64) + : b*log(a/b)) + +for cost in (:absdif, :abs2dif, :iibc, :kullback_leibler) + + @eval begin + + $cost(a::Real, b::Real) = $cost(Float64(a), Float64(b)) + + # Compute discrepancy map with respect to the background level. + function $cost{T,N}(a::AbstractArray{T,N}, beta::T) + c = Array(Float64, size(a)) + s = 0.0 + j = 0 + @inbounds begin + @simd for i in eachindex(a) + t = $cost(a[i], beta) + j += 1 + c[j] = t + s += t + end + end + return (c, s) + end + + function $cost{T,N}(a::AbstractArray{T,N}, beta) + $cost(a, T(beta)) + end + + function $cost{T}(a::AbstractArray{T,2}, ca::Array{Float64,2}, + b::AbstractArray{T,2}, cb::Array{Float64,2}, + tx::Int, ty::Int) + @assert size(ca) == size(a) + @assert size(cb) == size(b) + s = 0.0 + xrange = max(1, 1 + tx) : min(size(a, 1), size(b, 1) + tx) + length(xrange) > 0 || return s + yrange = max(1, 1 + ty) : min(size(a, 2), size(b, 2) + ty) + length(yrange) > 0 || return s + @inbounds begin + for y in yrange + for x in xrange + s += $cost(a[x,y], b[x-tx,y-ty]) - ca[x,y] - cb[x-tx,y-ty] + end + end + end + return s + end + end + +end + +end # module diff --git a/src/interpolate.jl b/src/interpolate.jl new file mode 100644 index 0000000..0c37d06 --- /dev/null +++ b/src/interpolate.jl @@ -0,0 +1,110 @@ + +""" + interpolate(A, i...) -> B + +interpolates array `A` along its dimensions at fractional indices `i...`. + +""" +function interpolate(A::AbstractArray{T,N}, + i::Vararg{N,AbstractVector{<:Real}}) where {T,N} + return interpolate(A, Val(1), i...) +end + +function interpolate(A::AbstractArray{T,N}, + ::Val{d}, + xd::AbstractVector{<:Real}, + xs::AbstractVector{<:Real}...) where {T,N,d} + B = Array{float(T),N}(undef, ) + unsafe_interpolate!(B, A, xd, Val(d)) + return interpolate(A, Val(d+1), xs...) +end +function interpolate(A::AbstractArray{T,N}, + ::Val{N}, + xd::AbstractVector{<:Real}, + xs::AbstractVector{<:Real}...) where {T,N,d} + B = Array{float(T),N}(undef, ) + unsafe_interpolate!(B, A, xd, d) + return interpolate(A, Val(d+1), xs...) +end + + +function interpolate!(dst::AbstractArray{T,N}, + src::AbstractArray{T,N}, + grd::AbstractVector{<:Real}, + dim::Integer) where {T,N} + @assert 1 ≤ dim ≤ N + I_dst = axes(dst) + I_src = axes(src) + I_head = I_src[1:dim-1] + I_tail = I_src[dim+1:N] + I_run_src = I_src[dim] + I_run_dst = I_dst[dim] + + I_dst[1:dim-1] == I_head || throw(DimensionMismatch( + "leading axes are not the same")) + I_dst[dim+1:N] == I_tail || throw(DimensionMismatch( + "trailing axes are not the same")) + axes(grd) == (I_run_src,) || throw(DimensionMismatch( + "vector of positions has incompatible indices")) + + unsafe_foo!(dst, CartesianIndices(I_head), I_run_dst, CartesianIndices(I_tail), + src, I_run_src, grd) + return dst + , + const dim1 = size(src, 1) + dst = Array(Float64, length(xvals)) + k = 0 + for x in xvals + (i0, i1, u0, u1) = interpcoefs(x, dim1) + k += 1 + dst[k] = src[i0]*u0 + src[i1]*u1 + end + return dst +end + +@Inline function interp_coefs(x::T, I::AbstractUnitRange{Int}) where {T<:AbstractFloat} + if x < first(I) + i1 = first(I) + i2 = i2 + w1 = zero(T) + w2 = one(T) + elseif x ≥ last(I) + i1 = last(I) + i2 = i1 + w1 = one(T) + w0 = zero(T) + else + i1 = floor(Int, x) + i2 = i1 + 1 + w2 = x - floor(x) + w1 = one(T) - w2 + end + return i1, w1, i2, w2 +end + +function unsafe_interpolate(dst::AbstractArray{<:Any,N}, + A::AbstractArray{<:Any,N}, + x::AbstractVector{<:Real}, + ::Val{d}) where {N,d} + inds = axes(A) + I = CartesianIndices(inds[1:d-1]) + J = axes(dst)[d] + K = CartesianIndices(inds[d+1:N]) + Jp = inds[d] + + zerofill!(dst) + + # Assume A and dst are in column-major order. + for k ∈ K + for j ∈ J + j1, w1, j2, w2 = interp_coefs(x[j], Jp) + for i ∈ I + dst[i,j,k] = w1*A[i,j1,k] + w2*A[i,j2,k] + end + end + end + return dst +end + + +zerofill!(A::AbstractArray) = fill!(A, zero(eltype(A)))