Made the new 'all_sky2pix()' method in HSTWCS dramatically faster, by performing the fit in the linear pixel space rather than in decimal degrees on the sky. This allows the method to use a simple (limited) for loop for the iterations instead of a potentially infinite while loop, while generally requiring much fewer than 10 iterations as opposed to over 350 for the fit on the sky. The added complication, though, is that this method now relies on a linear WCS derived using the 'utils.output_wcs()' function.

git-svn-id: http://svn.stsci.edu/svn/ssb/stsci_python/stsci_python/trunk/stwcs@12004 fe389314-cf27-0410-b35b-8c050e845b92

1 files changed, 20 insertions, 15 deletions

diff --git a/wcsutil/hstwcs.py b/wcsutil/hstwcs.py
index 9de9020..1d38f42 100644
--- a/wcsutil/hstwcs.py
+++ b/wcsutil/hstwcs.py
@@ -4,7 +4,7 @@ import os.path
 from pywcs import WCS
 import pyfits
 import instruments
-from stwcs.distortion import models, coeff_converter
+from stwcs.distortion import models, coeff_converter, utils
 import altwcs
 import numpy as np
 from pytools import fileutil
@@ -335,24 +335,31 @@ class HSTWCS(WCS):
         y = np.zeros(3,dtype=np.float64)
     
         # define delta for comparison
-        err = 0.00001
+        err = 0.0001
+        
+        # Use linear WCS as frame in which to perform fit
+        # rather than on the sky
+        wcslin = utils.output_wcs([self])
         
         # We can only transform 1 position at a time
         for r,d,n in zip(ra,dec,xrange(len(ra))):
             # First guess for output 
             x[0],y[0] = self.wcs_sky2pix(r,d,origin)
+            # also convert RA,Dec into undistorted linear pixel positions
+            lx,ly = wcslin.wcs_sky2pix(r,d,origin)
 
-            # Loop around until we are close enough
-            still_iterating = True
+            # Loop around until we are close enough (max 20 iterations)
             ev_old = None
-            while still_iterating:
+            for i in xrange(20):
                 x[1] = x[0] + 1.0
                 y[1] = y[0] 
                 x[2] = x[0]
                 y[2] = y[0] + 1.0
                 # Perform full transformation on pixel position
-                xo,yo = self.all_pix2sky(x,y,1)
-                
+                rao,deco = self.all_pix2sky(x,y,origin)
+                # convert RA,Dec into linear pixel positions for fitting
+                xo,yo = wcslin.wcs_sky2pix(rao,deco,origin)
+
                 # Compute deltas between output and initial guess as 
                 # an affine transform then invert that transformation
                 dxymat = np.array([[xo[1]-xo[0],yo[1]-yo[0]],
@@ -360,18 +367,16 @@ class HSTWCS(WCS):
                          
                 invmat = np.linalg.inv(dxymat)
                 # compute error in output position
-                delta_ra = r - xo[0]
-                delta_dec = d - yo[0]
+                dx = lx - xo[0]
+                dy = ly - yo[0]
 
                 # record the old position
                 xold = x[0]
                 yold = y[0]
 
                 # Update the initial guess position using the transform
-                delta_nx = delta_ra*dxymat[0][0] + delta_dec*dxymat[1][0]
-                delta_ny = delta_ra*dxymat[0][1] + delta_dec*dxymat[1][1]
-                x[0] = xold + delta_nx*3600*3600/self.pscale
-                y[0] = yold + delta_ny*3600*3600/self.pscale
+                x[0] = xold + dx*dxymat[0][0] + dy*dxymat[1][0]
+                y[0] = yold + dx*dxymat[0][1] + dy*dxymat[1][1]
 
                 # Work out the error vector length
                 ev = np.sqrt((x[0] - xold)**2 + (y[0] - yold)**2)
@@ -379,11 +384,11 @@ class HSTWCS(WCS):
                 # initialize record of previous error measurement during 1st iteration
                 if ev_old is None:
                     ev_old = ev
+                    
                 # Check to see whether we have reached the limit or 
                 # the new error is greater than error from previous iteration
                 if ev < err or (np.abs(ev) > np.abs(ev_old)): 
-                    still_iterating = False
-                    
+                    break
                 # remember error measurement from previous iteration
                 ev_old = ev