Index: trunk/doc/misc/dualconv.tex
===================================================================
--- trunk/doc/misc/dualconv.tex	(revision 15822)
+++ trunk/doc/misc/dualconv.tex	(revision 15822)
@@ -0,0 +1,200 @@
+%%% $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}
+