Index: /trunk/ppSub/TO_DO
===================================================================
--- /trunk/ppSub/TO_DO	(revision 14508)
+++ /trunk/ppSub/TO_DO	(revision 14508)
@@ -0,0 +1,119 @@
+Some ideas for improving the image subtraction
+==============================================
+
+Paul Price
+15 August 2007
+
+
+1. Dual convolution
+
+The Alard & Lupton (1998) and Alard (2000) algorithms are unable to
+produce useful subtractions when the PSFs of the input and reference
+images have misaligned position angles (e.g., / vs \), since this
+would involve a deconvolution along one of the axes, which does not
+work well.  One way to solve this problem is through a dual
+convolution.  Instead of solving
+
+	I(x,y) = k(u,v) * R(x,y)
+
+for k(u,v), we can solve
+
+	k1(u,v) * I(x,y) = k2(u,v) * R(x,y)
+
+for k1(u,v) and k2(u,v).  If we write
+
+      k1(u,v) = sum_i a_i f_i(u,v)
+
+and
+
+      k2(u,v) = sum_i b_i g_i(u,v)
+
+where f_i(u,v) and g_i(u,v) are sets of basis functions, then the
+solution boils down to two equations:
+
+      sum_j a_j A_j A_i = sum_j b_j B_j A_i
+
+and
+
+      sum_j a_j A_j B_i = sum_j b_j B_j B_i
+
+where
+
+      A_i = sum_x,y f_i(u,v) * I(x,y) / sigma(x,y)
+
+and
+
+      B_i = sum_x,y g_i(u,v) * R(x,y) / sigma(x,y)
+
+Noting that the matrix terms of the RHS of the first equation and the
+LHS of the second equation (A_i B_j) are the same leads us to attempt
+to solve this system by iteration:
+
+  1. Set a_i = 1 if i = 0, otherwise a_i = 0; f_0 = delta(u,v); and
+     solve for b in the first equation.  This corresponds to doing the
+     usual Alard & Lupton solution.
+
+  2. Use b in the second equation, and solve for a.
+
+  3. Use a in the first equation, and solve for b.
+
+  4. Proceed in this manner until the change in b between iterations
+     is small.
+
+This is basically doing A&L to get the reference image to match the
+input image in the usual manner, then doing A&L on the input image to
+match the convolved reference image, then doing A&L on the reference
+image to match the convolved input image, etc., until the system
+settles down to the solution.  A more refined method of solving the
+equations may exist, but this should work.
+
+There is a need to normalise one of the kernels: the solution is
+currently non-unique because a and b can be scaled by some arbitrary
+value.  If we write
+
+      k1(u,v) = delta(u,v) + sum_i a_i [k_i - delta(u,v)]
+
+and the k_i are normalised to a sum of unity, then
+
+      sum_u,v k1(u,v) = 1.
+
+So k1 has a sum of unity, always.  This means that k2 supplies the
+scaling to match the photometry (as happens in Alard & Lupton), and k1
+is merely used to broaden the input image as required for the best
+subtraction.
+
+Note that the expense for this method over the usual Alard & Lupton
+roughly amounts to the extra image convolution, since that is
+typically the dominant factor.  Accumulation of the sums is not much
+more than Alard & Lupton, and the iteration should be fairly fast.
+
+
+2. Data-based kernel selection
+
+The quality of the subtraction is highly sensitive to the choice of
+basis functions.  In the case of ISIS kernels (which seem to be the
+most useful of those experimented with so far, because it takes a
+small number of parameters to generate a large kernel), the choice of
+the Gaussian widths is very important.  Some recipes exist for
+choosing these widths, but these are generally motivated by experience
+(through much trial and error) rather than directly from the data.  It
+would be nice to be able to throw down a large number of widths and
+allow the least-squares solution to choose the best for the data at
+hand (i.e., the most dominant contributors), but this is prohibitively
+expensive --- at least for a full solution.
+
+What might be possible is to do a quick and dirty solution by reducing
+the dimensionality.  Instead of working with the full two-dimensional
+kernel, k(u,v), let's work in one dimension, k(u).  For each of our
+stamps, let's take a cut through them in a consistent direction (e.g.,
+the x direction), and solve
+
+	I(x) = k(u) * R(x)
+
+This is not nearly as computationally expensive as the full solution,
+so we can pack k(u) with multiple Gaussian widths, solve the
+least-squares problem, and identify the most important kernel
+contributions which we will use (in suitable two-dimensional versions)
+in solving the full problem.  This method could even be applied with
+multiple direction cuts to ensure the full range of required Gaussian
+widths is obtained.