%%% $Id: dualconv.tex,v 1.1 2007-12-14 02:22:06 price Exp $
\documentclass[panstarrs]{panstarrs}

% basic document variables
\title{Dual Convolution Image Subtraction}
\subtitle{Formulation}
\shorttitle{Dual Convolution}
\author{Paul Price}
\audience{Pan-STARRS PMO}
\group{Pan-STARRS Algorithm Group}
\project{Pan-STARRS Image Processing Pipeline}
\organization{Institute for Astronomy}
\version{DR}
\docnumber{PSDC-430-???}

\setlength{\topsep}{-2pt}
  
\begin{document}
\maketitle
\sloppy

% -- Revision History --
% provide explicit values for the old versions
% use '\theversion' for the current version (set above)
% use \hline between each table row
\RevisionsStart
% version  Date            Description
DR & 2007 Dec 13 & Draft \\ \hline
\RevisionsEnd

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\DocumentsInternal
\DocumentsExternal
\DocumentsEnd

\tableofcontents
\pagebreak 
\pagenumbering{arabic}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}

The current state-of-the-art in image subtraction uses the algorithm
of Alard (\code{2000A&AS..144..363A}), which builds on the original
algorithm of Alard \& Lupton (\code{1998ApJ...503..325A}).  However,
these algorithms cannot deal with the case where the PSFs are
misaligned.  That is, if the PSFs are both elongated, and the
direction of the elongation is not identical, then the Alard (2000)
algorithm cannot provide a good subtraction.

We might attempt to get around this case by using a dual convolution
--- convolving both images to a common PSF.

\section{Formulation}

Given images $I_1(x,y)$ and $I_2(x,y)$, we wish to solve for the
convolution kernels $k_1(u,v)$ and $k_2(u,v)$, that produce:
\begin{equation}
k_1(u,v) \otimes I_1(x,y) + \Delta bg = k_2(u,v) \otimes I_2(x,y) + I_2(x,y)
\end{equation}

We have explicitly included an unconvolved $I_2(x,y)$ on the RHS
because our intention is that the convolution kernel $k_2$ has a sum
of zero (i.e., its purpose is merely to broaden the PSF in $I_2$).
Convolution kernel $k_1$ will have a non-zero sum, and allow the
matching of the normalisation between the images.  

We proceed by minimising:
\begin{equation}
\chi^2 = \sum_{x,y} \left[ \frac{k_1(u,v) \otimes I_1(x,y) + \Delta bg - \sum_{x,y} k_2(u,v) \otimes I_2(x,y) - I_2(x,y)}{\sigma(x,y)} \right]^2
\label{eqn:chi2}
\end{equation}

If we express the kernels as a linear combination of basis functions,
this becomes a linear problem which can be solved (relatively)
quickly:
\begin{eqnarray}
k_1(u,v) & = & \sum_i a_i f_i(u,v) \\
k_2(u,v) & = & \sum_i b_i g_i(u,v)
\end{eqnarray}

In order to get the normalisation behaviour we desire, we actually
use:
\begin{eqnarray}
k_1(u,v) & = & a_1 \delta_{u0} \delta_{v0} + \sum_{i>1} a_i f_i' \\
k_2(u,v) & = & \sum_{i>1} b_i g_i'
\end{eqnarray}
where $\delta_ij$ is the Kronecker delta, and:
\begin{equation}
f_i' = \begin{cases}
  f_i / \sum_{u,v} f_i - \delta_{u0} \delta_{v0} & \text{if $f_i$ is an even function} \\
  f_i & \text{if $f_i$ is an odd function}
  \end{cases} \\
\end{equation}
and similarly for $g_i$.

Then we take the derivative of equation \ref{eqn:chi2} w.r.t. $a_j$
and $b_j$ to yield the following two equations:
\begin{eqnarray}
\sum_{x,y} [ \sum_i a_i A_i(x,y) - \sum_i b_i B_i(x,y) - I_2(x,y) ] A_j(x,y) & = & 0 \\
\sum_{x,y} [ \sum_i a_i A_i(x,y) - \sum_i b_i B_i(x,y) - I_2(x,y) ] B_j(x,y) & = & 0
\end{eqnarray}
where we have written:
\begin{equation}
A_i(x,y) = \begin{cases}
  1 & \text{if $i = 0$ (background matching);} \\
  I_1(x,y) & \text{if $i = 1$ (normalisation matching);} \\
  f_i(u,v) \otimes I_1(x,y) & \text{for all the rest.}
\end{cases}
\end{equation}
and
\begin{equation}
B_i(x,y) = g_i(u,v) \otimes I_2(x,y)
\end{equation}

More general background matching might be accomplished by generalising
the above expression for $A_0$ (replacing with multiple polynomial
terms), but we will content ourselves with a constant background
difference here.

We therefore obtain the following equations:
\begin{eqnarray}
\sum_{x,y} \sum_i a_i A_i(x,y) A_j(x,y) & = & \sum_{x,y} \sum_i b_i B_i(x,y) A_j(x,y) + \sum_{x,y} I_2(x,y) A_j(x,y) \\
\sum_{x,y} \sum_i a_i A_i(x,y) B_j(x,y) & = & \sum_{x,y} \sum_i b_i B_i(x,y) B_j(x,y) + \sum_{x,y} I_2(x,y) B_j(x,y)
\end{eqnarray}
which may be written as a pair of coupled matrix equations:
\begin{eqnarray}
A a & = & C^T b + d \\
C a & = & B b + e
\end{eqnarray}
where $a$ and $b$ are column vectors containing the coefficients of interest, and:
\begin{eqnarray}
A_{ij} & = & \sum_{x,y} A_i(x,y) A_j(x,y) \\
B_{ij} & = & \sum_{x,y} B_i(x,y) B_j(x,y) \\
C_{ij} & = & \sum_{x,y} A_i(x,y) B_j(x,y) \\
d_i    & = & \sum_{x,y} A_i(x,y) I_2(x,y) \\
e_i    & = & \sum_{x,y} B_i(x,y) I_2(x,y)
\end{eqnarray}

Note that matrix $C$ is not necessarily square, even if we are using
the same basis functions in $k_2$ as in $k_1$, because of the
additional elements for background and normalisation matching.

The solution to the matrix equations are:
\begin{eqnarray}
a & = & (C^T B^{-1} C - A)^{-1} (C^T B^{-1} e - d) \\
b & = & B^{-1}(C a - e)
\end{eqnarray}

\section{Additional considerations}

In practise, the above sums over $x,y$ are performed over selected
sub-images (``stamps''), rather than the entire image.  Once a
solution is obtained, it is applied to the stamps, which may be
rejected and replaced in favour of another (e.g., the original stamp
was a cosmic ray, or unmasked saturated star); the process is repeated
until a satisfactory solution is found.

Spatial variation of the kernels may be achieved by writing
\begin{eqnarray}
a_i(x,y) & = & \sum_j \alpha_j P_j(x,y) \\
b_i(x,y) & = & \sum_j \beta_j P_j(x,y)
\end{eqnarray}
where $P_j$ is a polynomial term in $x$ and $y$.  Then we do the same
calculations as before, except we add an additional iteration over the
polynomial terms, and end up solving for the $\alpha_i$ and $\beta_i$.
As pointed out in Alard (2000), this doesn't add much in terms of
computation time, because the main operations are the generation of
the $A_i$ and $B_j$, and the summing over the products of these to
generate the least-squares matrices and vectors.

If the basis functions for the kernels $k_1$ and $k_2$ are the same,
then the solution might be close to degenerate because a solution can
be found by making the $k_1$ and $k_2$ as wide as possible using the
same elements of the basis function.  In this case, we can add an
additional term to equation \ref{eqn:chi2} to force the coefficients
of the same terms to be different:
\begin{equation}
\chi^2_{\rm force} = \chi^2 - w \sum_i (a_i^2 - b_i^2)
\end{equation}
where $w$ is a weighting parameter (read: fudge factor) to make the
forcing term to be of the same magnitude as $\chi^2$ (so its effect is
neither diluted nor emphasised).  Note that this forumulation means
that the sign of $a_i$ and $b_i$ doesn't matter, but what matters is
that $b_i$ is small when $a_i$ is large.  This change to $\chi^2$ has
the effect of adding an additional term to the diagonals of the square
matrices:
\begin{eqnarray}
A_{ij} & = & \sum_{x,y} A_i(x,y) A_j(x,y) - w \delta_{ij} \\
B_{ij} & = & \sum_{x,y} B_i(x,y) B_j(x,y) + w \delta_{ij}
\end{eqnarray}


\bibliographystyle{plain}
\bibliography{panstarrs}

\end{document}

