Merge pull request #48 from keflavich/issue47

Fix for issue 47

image_registration/chi2_shifts.py

@@ -29,7 +29,7 @@ def chi2_shift(im1, im2, err=None, upsample_factor='auto', boundary='wrap',
           & = & \Sigma_{ij} \left[ X_{ij}^2/\sigma_{ij}^2 - 2X_{ij}Y_{ij}(dx,dy)/\sigma_{ij}^2 + Y_{ij}(dx,dy)^2/\sigma_{ij}^2 \\right]  \\\\
 
 
-    Equation 2-4: blahha
+    Equation 2-4:
 
     .. math::
        :nowrap:
@@ -63,6 +63,14 @@ def chi2_shift(im1, im2, err=None, upsample_factor='auto', boundary='wrap',
     since term 1 is the same for all shifts, and the quantity of interest is
     :math:`\Delta \chi^2` when determining the best-fit shift and error.
 
+    The resulting shifts are limited to an accuracy of +/-0.5 pixels in the
+    upsampled image frame.  That is not a Gaussian uncertainty but a quantized
+    limit: the best solution will always be +/-0.5 pixels offset because we're
+    zooming in on an even grid, so the "best" position is required to be a
+    discrete pixel center.  If you're looking at an image with exactly zero
+    shift, it will have exactly +/- 1/usfac/2 systematic error in the resulting
+    solution.
+
 
     Parameters
     ----------
@@ -128,26 +136,29 @@ def chi2_shift(im1, im2, err=None, upsample_factor='auto', boundary='wrap',
     >>> dx2,dy2,edx2,edy2 = chi2_shift(image, shifted2, upsample_factor='auto')
 
     """
-    chi2,term1,term2,term3 = chi2n_map(im1, im2, err, boundary=boundary,
-                                       nthreads=nthreads, zeromean=zeromean,
-                                       use_numpy_fft=use_numpy_fft,
-                                       return_all=True, reduced=False)
+    chi2, term1, term2, term3 = chi2n_map(im1, im2, err, boundary=boundary,
+                                          nthreads=nthreads, zeromean=zeromean,
+                                          use_numpy_fft=use_numpy_fft,
+                                          return_all=True, reduced=False)
     ymax, xmax = np.unravel_index(chi2.argmin(), chi2.shape)
 
     # needed for ffts
     im1 = np.nan_to_num(im1)
     im2 = np.nan_to_num(im2)
 
-    ylen,xlen = im1.shape
-    xcen = xlen/2-(1-xlen%2)
-    ycen = ylen/2-(1-ylen%2)
+    ylen, xlen = im1.shape
+
+    # this is the center pixel - it's an integer pixel ID (not the center
+    # coordinate)
+    xcen = xlen/2 - (1 if xlen % 2 == 0 else 0.5)
+    ycen = ylen/2 - (1 if ylen % 2 == 0 else 0.5)
 
     # original shift calculation
-    yshift = ymax-ycen # shift im2 by these numbers to get im1
-    xshift = xmax-xcen
+    yshift = ymax - ycen # shift im2 by these numbers to get im1
+    xshift = xmax - xcen
 
     if verbose:
-        print("Coarse xmax/ymax = %i,%i, for offset %f,%f" % (xmax,ymax,xshift,yshift))
+        print("Coarse xmax/ymax = %i,%i, for offset %f,%f" % (xmax, ymax, xshift, yshift))
 
     # below is sub-pixel zoom-in stuff
 
@@ -167,7 +178,7 @@ def chi2_shift(im1, im2, err=None, upsample_factor='auto', boundary='wrap',
         m_auto = 2.6088233328527037 # slightly >1 sigma
 
     # biggest scale = where chi^2/n ~ 9 or 11.8 for M=2?
-    if upsample_factor=='auto':
+    if upsample_factor == 'auto':
         # deltachi2 is not reduced deltachi2
         deltachi2_lowres = (chi2 - chi2.min())
         if verbose:
@@ -199,9 +210,11 @@ def chi2_shift(im1, im2, err=None, upsample_factor='auto', boundary='wrap',
             s1 = zoom_factor*upsample_factor
             s2 = zoom_factor*upsample_factor
 
-    (yshifts_corrections,xshifts_corrections),chi2_ups = zoom.zoomnd(chi2,
-            usfac=upsample_factor, outshape=[s1,s2], offsets=[yshift,xshift],
-            return_xouts=True)
+    (yshifts_corrections, xshifts_corrections), chi2_ups = zoom.zoomnd(chi2,
+                                                                       usfac=upsample_factor,
+                                                                       outshape=[s1, s2],
+                                                                       offsets=[yshift, xshift],
+                                                                       return_xouts=True)
 
     # deltachi2 is not reduced deltachi2
     deltachi2_ups = (chi2_ups - chi2_ups.min())
@@ -217,7 +230,7 @@ def chi2_shift(im1, im2, err=None, upsample_factor='auto', boundary='wrap',
 
     # BELOW IS TO COMPUTE THE ERROR
 
-    errx_low,errx_high,erry_low,erry_high = chi2map_to_errors(chi2_ups, upsample_factor)
+    errx_low, errx_high, erry_low, erry_high = chi2map_to_errors(chi2_ups, upsample_factor)
 
     yshift_corr = yshifts_corrections.flat[chi2_ups.argmin()]-ycen
     xshift_corr = xshifts_corrections.flat[chi2_ups.argmin()]-xcen

image_registration/fft_tools/zoom.py

@@ -52,7 +52,7 @@ def zoom1d(inp, usfac=1, outsize=None, offset=0, nthreads=1,
         return result
 
 def zoom_on_pixel(inp, coordinates, usfac=1, outshape=None, nthreads=1,
-        use_numpy_fft=False, return_real=True, return_xouts=False):
+                  use_numpy_fft=False, return_real=True, return_xouts=False):
     """
     Zoom in on a 1D or 2D array using Fourier upsampling
     (in principle, should work on N-dimensions, but does not at present!)
@@ -83,27 +83,30 @@ def zoom_on_pixel(inp, coordinates, usfac=1, outshape=None, nthreads=1,
     """
 
     inshape = inp.shape
-    if outshape is None: outshape=inshape
+    if outshape is None:
+        outshape = inshape
 
-    outarr = np.zeros((inp.ndim,)+tuple(outshape),dtype='float')
-    for ii,(insize, outsize, target) in enumerate(zip(inshape,outshape,coordinates)):
+    outarr = np.zeros((inp.ndim,)+tuple(outshape), dtype='float')
+    for ii, (insize, outsize, target) in enumerate(zip(inshape, outshape, coordinates)):
         # output array should cover 1/usfac *  the range of the input
         # it should go from 1/2.-1/usfac to 1/2+1/usfac
         # plus whatever offset is specified
         # outsize is always 1+(highest index of input)
         outarr_d = np.linspace(target - (outsize-1.)/usfac/2.,
                                target + (outsize-1.)/usfac/2.,
                                outsize)
-        
+
         # slice(None) = ":" or "get everything"
         # [None] = newaxis = add a blank axis on this dim
         dims = [None]*ii + [slice(None)] + [None]*(inp.ndim-1-ii)
         outarr[ii] = outarr_d[tuple(dims)]
 
     # temporary hack
     if inp.ndim == 1:
-        result = scale.fourier_interp1d(inp, outarr.squeeze(), nthreads=nthreads,
-                use_numpy_fft=use_numpy_fft, return_real=return_real)
+        result = scale.fourier_interp1d(inp, outarr.squeeze(),
+                                        nthreads=nthreads,
+                                        use_numpy_fft=use_numpy_fft,
+                                        return_real=return_real)
     elif inp.ndim == 2:
         result = scale.fourier_interp2d(inp, outarr, nthreads=nthreads,
                                         use_numpy_fft=use_numpy_fft,
@@ -112,7 +115,7 @@ def zoom_on_pixel(inp, coordinates, usfac=1, outshape=None, nthreads=1,
         raise NotImplementedError("Can't do more than 2D yet")
 
     if return_xouts:
-        return outarr,result
+        return outarr, result
     else:
         return result
 
@@ -137,7 +140,7 @@ def zoomnd(inp, offsets=(), middle_convention=float, **kwargs):
     outshape : int
         Number of pixels in output array
         (passed to :func:`zoom_on_pixel`)
-    
+
     Other Parameters
     ----------------
     return_xouts : bool
@@ -165,7 +168,7 @@ def zoomnd(inp, offsets=(), middle_convention=float, **kwargs):
     # outsize is always 1+(highest index of input)
 
     middlepix = [middle_convention((insize-1)/2.) + off
-                 for insize,off in zip(inp.shape,offsets)]
+                 for insize, off in zip(inp.shape, offsets)]
 
     return zoom_on_pixel(inp, middlepix, **kwargs)