diff options
| -rw-r--r-- | lib/stwcs/wcsutil/hstwcs.py | 354 |
1 files changed, 271 insertions, 83 deletions
diff --git a/lib/stwcs/wcsutil/hstwcs.py b/lib/stwcs/wcsutil/hstwcs.py index e6c812d..8f66283 100644 --- a/lib/stwcs/wcsutil/hstwcs.py +++ b/lib/stwcs/wcsutil/hstwcs.py @@ -1,6 +1,6 @@ from __future__ import division # confidence high -import os.path +import os from pywcs import WCS import pyfits import instruments @@ -79,12 +79,12 @@ class NoConvergence(Exception): solution appears to be divergent. If the solution does not diverge, `divergent` will be set to `None`. - nonconvergent : None, numpy.array + failed2converge : None, numpy.array Indices of the points in :py:attr:`best_solution` array for which the solution failed to converge within the specified maximum number of iterations. If there are no non-converging poits (i.e., if the required accuracy has been achieved for all points) then - `nonconvergent` will be set to `None`. + `failed2converge` will be set to `None`. """ def __init__(self, *args, **kwargs): @@ -94,7 +94,7 @@ class NoConvergence(Exception): self.accuracy = kwargs.pop('accuracy', None) self.niter = kwargs.pop('niter', None) self.divergent = kwargs.pop('divergent', None) - self.nonconvergent = kwargs.pop('nonconvergent', None) + self.failed2converge= kwargs.pop('failed2converge', None) # @@ -454,7 +454,9 @@ class HSTWCS(WCS): Parameters ---------- accuracy : float, optional (Default = 1.0e-4) - Required accuracy of the solution. + Required accuracy of the solution. Iteration terminates when the + correction to the solution found during the previous iteration + is smaller (in the sence of the L2 norm) than `tolerance`\ . maxiter : int, optional (Default = 20) Maximum number of iterations allowed to reach the solution. @@ -462,25 +464,94 @@ class HSTWCS(WCS): adaptive : bool, optional (Default = False) Specifies whether to adaptively select only points that did not converge to a solution whithin the required accuracy for the - next iteration. Default is recommended for HST instruments. + next iteration. Default is recommended for HST as well as most + other instruments. .. note:: - The :py:meth:`all_sky2pix` uses a vectorized implementation of - the method of consecutive approximations (see `Notes` section - below) in which it iterates over *all* input poits *regardless* - untill the required accuracy has been reached for *all* input - points. In some cases it may be possible that *almost all* points - have reached the required accuracy but there are only a few - of input data points for which additional iterations may be - needed. In this situation it may be advantageous to set - `adaptive` = `True`\ in which case :py:meth:`all_sky2pix` will - continue iterating *only* over the points that have not yet - converged to the required accuracy. However, for the HST's - ACS/WFC detector, which has the strongest distortions of all - HST instruments, testing has shown that enabling this option - would lead to a 10-30\% penalty in computational time. - Therefore, for HST instruments, it is recommended to set - `adaptive` = `False`\ . + The :py:meth:`all_sky2pix` uses a vectorized implementation + of the method of consecutive approximations (see `Notes` + section below) in which it iterates over *all* input poits + *regardless* until the required accuracy has been reached for + *all* input points. In some cases it may be possible that + *almost all* points have reached the required accuracy but + there are only a few of input data points left for which + additional iterations may be needed (this depends mostly on the + characteristics of the geometric distortions for a given + instrument). In this situation it may be + advantageous to set `adaptive` = `True`\ in which + case :py:meth:`all_sky2pix` will continue iterating *only* over + the points that have not yet converged to the required + accuracy. However, for the HST's ACS/WFC detector, which has + the strongest distortions of all HST instruments, testing has + shown that enabling this option would lead to a about 10-30\% + penalty in computational time (depending on specifics of the + image, geometric distortions, and number of input points to be + converted). Therefore, for HST instruments, + it is recommended to set `adaptive` = `False`\ . The only + danger in getting this setting wrong will be a performance + penalty. + + .. note:: + When `detect_divergence` is `True`\ , :py:meth:`all_sky2pix` \ + will automatically switch to the adaptive algorithm once + divergence has been detected. + + detect_divergence : bool, optional (Default = True) + Specifies whether to perform a more detailed analysis of the + convergence to a solution. Normally :py:meth:`all_sky2pix` + may not achieve the required accuracy + if either the `tolerance` or `maxiter` arguments are too low. + However, it may happen that for some geometric distortions + the conditions of convergence for the the method of consecutive + approximations used by :py:meth:`all_sky2pix` may not be + satisfied, in which case consecutive approximations to the + solution will diverge regardless of the `tolerance` or `maxiter` + settings. + + When `detect_divergence` is `False`\ , these divergent points + will be detected as not having achieved the required accuracy + (without further details). In addition, if `adaptive` is `False` + then the algorithm will not know that the solution (for specific + points) is diverging and will continue iterating and trying to + "improve" diverging solutions. This may result in NaN or Inf + values in the return results (in addition to a performance + penalties). Even when `detect_divergence` is + `False`\ , :py:meth:`all_sky2pix`\ , at the end of the iterative + process, will identify invalid results (NaN or Inf) as "diverging" + solutions and will raise :py:class:`NoConvergence` unless + the `quiet` parameter is set to `True`\ . + + When `detect_divergence` is `True`\ , :py:meth:`all_sky2pix` will + detect points for + which current correction to the coordinates is larger than + the correction applied during the previous iteration **if** the + requested accuracy **has not yet been achieved**\ . In this case, + if `adaptive` is `True`, these points will be excluded from + further iterations and if `adaptive` + is `False`\ , :py:meth:`all_sky2pix` will automatically + switch to the adaptive algorithm. + + .. note:: + When accuracy has been achieved, small increases in + current corrections may be possible due to rounding errors + (when `adaptive` is `False`\ ) and such increases + will be ignored. + + .. note:: + Setting `detect_divergence` to `True` will incurr about 5-10\% + performance penalty (in our testing on ACS/WFC images). + Because the benefits of enabling this feature outweigh + the small performance penalty, it is recommended to set + `detect_divergence` to `True`\ , unless extensive testing + of the distortion models for images from specific + instruments show a good stability of the numerical method + for a wide range of coordinates (even outside the image + itself). + + .. note:: + Indices of the diverging inverse solutions will be reported + in the `divergent` attribute of the + raised :py:class:`NoConvergence` object. quiet : bool, optional (Default = False) Do not throw :py:class:`NoConvergence` exceptions when the method @@ -497,9 +568,9 @@ class HSTWCS(WCS): Notes ----- - Inputs can either be (RA, Dec, origin) or (RADec, origin) where RA and Dec - are 1-D arrays/lists of coordinates and RADec is an array/list of pairs - of coordinates. + Inputs can either be (RA, Dec, origin) or (RADec, origin) where RA + and Dec are 1-D arrays/lists of coordinates and RADec is an + array/list of pairs of coordinates. Using the method of consecutive approximations we iterate starting with the initial approximation, which is computed using the @@ -510,8 +581,8 @@ class HSTWCS(WCS): highly efficient (>30x) when *all* data points that need to be converted from sky coordinates to image coordinates are passed at *once*\ . Therefore, it is advisable, whenever possible, to pass - as input a long array of all points that need to be converted to - :py:meth:`all_sky2pix` instead of calling :py:meth:`all_sky2pix` + as input a long array of all points that need to be converted + to :py:meth:`all_sky2pix` instead of calling :py:meth:`all_sky2pix` for each data point. Also see the note to the `adaptive` parameter. Examples @@ -520,6 +591,7 @@ class HSTWCS(WCS): >>> hdulist = pyfits.open('j94f05bgq_flt.fits') >>> w = stwcs.wcsutil.HSTWCS(hdulist, ext=('sci',1)) >>> hdulist.close() + >>> ra, dec = w.all_pix2sky([1,2,3],[1,1,1],1); print(ra); print(dec) [ 5.52645241 5.52649277 5.52653313] [-72.05171776 -72.05171295 -72.05170814] @@ -537,10 +609,71 @@ class HSTWCS(WCS): [[ 1.00000233 0.99999997] [ 2.00000232 0.99999997] [ 3.00000233 0.99999998]] - >>> xy = w.all_sky2pix(radec,1, maxiter=3, accuracy=1.0e-10, quiet=False) + >>> xy = w.all_sky2pix(radec,1, maxiter=3, accuracy=1.0e-10, \ +quiet=False) NoConvergence: 'HSTWCS.all_sky2pix' failed to converge to requested \ accuracy after 3 iterations. + >>> + Now try to use some diverging data: + >>> divradec = w.all_pix2sky([[1.0,1.0],[10000.0,50000.0],\ +[3.0,1.0]],1); print(divradec) + [[ 5.52645241 -72.05171776] + [ 7.15979392 -70.81405561] + [ 5.52653313 -72.05170814]] + + >>> try: + >>> xy = w.all_sky2pix(divradec,1, maxiter=20, accuracy=1.0e-4, \ +adaptive=False, detect_divergence=True, quiet=False) + >>> except stwcs.wcsutil.hstwcs.NoConvergence as e: + >>> print("Indices of diverging points: {}".format(e.divergent)) + >>> print("Indices of poorly converging points: {}".format(e.failed2converge)) + >>> print("Best solution: {}".format(e.best_solution)) + >>> print("Achieved accuracy: {}".format(e.accuracy)) + >>> raise e + Indices of diverging points: + [1] + Indices of poorly converging points: + None + Best solution: + [[ 1.00006219e+00 9.99999288e-01] + [ -1.99440907e+06 1.44308548e+06] + [ 3.00006257e+00 9.99999316e-01]] + Achieved accuracy: + [[ 5.98554253e-05 6.79918148e-07] + [ 8.59514088e+11 6.61703754e+11] + [ 6.02334592e-05 6.59713067e-07]] + Traceback (innermost last): + File "<console>", line 8, in <module> + NoConvergence: 'HSTWCS.all_sky2pix' failed to converge to the requested accuracy. + After 5 iterations, the solution is diverging at least for one input point. + + >>> try: + >>> xy = w.all_sky2pix(divradec,1, maxiter=20, accuracy=1.0e-4, \ +adaptive=False, detect_divergence=False, quiet=False) + >>> except stwcs.wcsutil.hstwcs.NoConvergence as e: + >>> print("Indices of diverging points: {}".format(e.divergent)) + >>> print("Indices of poorly converging points: {}".format(e.failed2converge)) + >>> print("Best solution: {}".format(e.best_solution)) + >>> print("Achieved accuracy: {}".format(e.accuracy)) + >>> raise e + Indices of diverging points: + [1] + Indices of poorly converging points: + None + Best solution: + [[ 1. 1.] + [ nan nan] + [ 3. 1.]] + Achieved accuracy: + [[ 0. 0.] + [ nan nan] + [ 0. 0.]] + Traceback (innermost last): + File "<console>", line 8, in <module> + NoConvergence: 'HSTWCS.all_sky2pix' failed to converge to the requested accuracy. + After 20 iterations, the solution is diverging at least for one input point. + """ nargs = len(args) @@ -571,13 +704,15 @@ accuracy after 3 iterations. .format(nargs)) # process optional arguments: - accuracy = kwargs.pop('accuracy', 1.0e-4) - maxiter = kwargs.pop('maxiter', 20) - quiet = kwargs.pop('quiet', False) - adaptive = kwargs.pop('adaptive', False) + accuracy = kwargs.pop('accuracy', 1.0e-4) + maxiter = kwargs.pop('maxiter', 20) + adaptive = kwargs.pop('adaptive', False) + detect_divergence = kwargs.pop('detect_divergence', True) + quiet = kwargs.pop('quiet', False) # initialize iterative process: - x0, y0 = self.wcs_sky2pix(ra, dec, origin) # initial approximation (WCS based only) + x0, y0 = self.wcs_sky2pix(ra, dec, origin) # <-- initial approximation + # (WCS based only) # see if iterative solution is required (when any of the # non-CD-matrix corrections are present). If not required @@ -585,7 +720,8 @@ accuracy after 3 iterations. if self.sip is None and \ self.cpdis1 is None and self.cpdis2 is None and \ self.det2im1 is None and self.det2im2 is None: - # no non-WCS corrections are detected - return initial approximation + # no non-WCS corrections are detected - return + # initial approximation if vect1D: return [x0, y0] else: @@ -608,8 +744,8 @@ accuracy after 3 iterations. y -= dy # norn (L2) squared of the correction: - dn2prev = dx**2+dy**2 - dn2 = dn2prev + dn2prev = dx**2+dy**2 + dn2 = dn2prev # process all coordinates simultaneously: iterlist = range(1, maxiter+1) @@ -617,7 +753,12 @@ accuracy after 3 iterations. ind = None inddiv = None - divergent = False + npts = x.shape[0] + + # turn off numpy runtime warnings for 'invalid' and 'over': + old_invalid = np.geterr()['invalid'] + old_over = np.geterr()['over'] + np.seterr(invalid = 'ignore', over = 'ignore') if not adaptive: for k in iterlist: @@ -625,66 +766,113 @@ accuracy after 3 iterations. if np.max(dn2) < accuracy2: break - # check for divergence: - inddiv, = np.where((dn2 > dn2prev) & (dn2 >= accuracy2)) - if inddiv.shape[0] > 0: - divergent = True - break - # find correction to the previous solution: dx, dy = self.pix2foc(x, y, origin) # If pix2foc does not apply all the required distortion # corrections then replace the above line with: #r0, d0 = self.all_pix2sky(x, y, origin) #dx, dy = self.wcs_sky2pix(r0, d0, origin ) - dx -= x0 dy -= y0 + # update norn (L2) squared of the correction: + dn2 = dx**2+dy**2 + + # check for divergence (we do this in two stages + # to optimize performance for the most common + # scenario when succesive approximations converge): + if detect_divergence: + ind, = np.where(dn2 <= dn2prev) + if ind.shape[0] < npts: + inddiv, = np.where( + np.logical_and(dn2 > dn2prev, dn2 >= accuracy2)) + if inddiv.shape[0] > 0: + # apply correction only to the converging points: + x[ind] -= dx[ind] + y[ind] -= dy[ind] + # switch to adaptive iterations: + ind, = np.where((dn2 >= accuracy2) & \ + (dn2 <= dn2prev) & np.isfinite(dn2)) + iterlist = iterlist[k:] + adaptive = True + break + #dn2prev[ind] = dn2[ind] + dn2prev = dn2 + # apply correction: x -= dx y -= dy - # update norn (L2) squared of the correction: - dn2prev = dn2.copy() - dn2 = dx**2+dy**2 + if adaptive: - else: - ind, = np.where(dn2 >= accuracy2) + if ind is None: + ind = np.asarray(range(npts), dtype=np.int64) for k in iterlist: # check convergence: if ind.shape[0] == 0: break - # check for divergence: - inddiv = ind[np.where(dn2[ind] > dn2prev[ind])] - if inddiv.shape[0] > 0: - divergent = True - break - # find correction to the previous solution: dx[ind], dy[ind] = self.pix2foc(x[ind], y[ind], origin) # If pix2foc does not apply all the required distortion # corrections then replace the above line with: #r0[ind], d0[ind] = self.all_pix2sky(x[ind], y[ind], origin) - #dx[ind], dy[ind] = self.wcs_sky2pix(r0[ind], d0[ind], origin ) + #dx[ind], dy[ind] = self.wcs_sky2pix(r0[ind], d0[ind], origin) dx[ind] -= x0[ind] dy[ind] -= y0[ind] + # update norn (L2) squared of the correction: + dn2 = dx**2+dy**2 + + # update indices of elements that still need correction: + if detect_divergence: + ind, = np.where((dn2 >= accuracy2) & (dn2 <= dn2prev)) + dn2prev[ind] = dn2[ind] + else: + ind, = np.where(dn2 >= accuracy2) + #ind = ind[np.where(dx[ind]**2+dy[ind]**2 >= accuracy2)] + # apply correction: x[ind] -= dx[ind] y[ind] -= dy[ind] - # update norn (L2) squared of the correction: - dn2prev = dn2.copy() - dn2 = dx**2+dy**2 + #FINAL DETECTION OF INVALID, DIVERGING, AND FAILED-TO-CONVERGE POINTS: + invalid = (((~np.isfinite(y)) | (~np.isfinite(x)) | \ + (~np.isfinite(dn2))) & \ + (np.isfinite(ra)) & (np.isfinite(dec))) + if detect_divergence: + inddiv, = np.where(((dn2 >= accuracy2) & (dn2 > dn2prev)) | \ + invalid) + # identify diverging and/or invalid points: + if inddiv.shape[0] == 0: + inddiv = None + # identify points that did not converge within + # 'maxiter' iterations: + if k >= maxiter: + ind, = np.where((dn2 >= accuracy2) & (dn2 <= dn2prev) & \ + (~invalid)) + if ind.shape[0] == 0: + ind = None + else: + ind = None - # update indices of elements that still need correction: - ind, = np.where(dn2 >= accuracy2) - #ind = ind[np.where(dx[ind]**2+dy[ind]**2 >= accuracy2)] + else: + # dn2prev is not available in this case. + inddiv, = np.where(invalid) + # identify diverging and/or invalid points: + if inddiv.shape[0] == 0: + inddiv = None + # identify points that did not converge within + # 'maxiter' iterations: + if k >= maxiter: + ind, = np.where((dn2 >= accuracy2) & (~invalid)) + if ind.shape[0] == 0: + ind = None + else: + ind = None - if k >= maxiter and not quiet: + if (ind is not None or inddiv is not None) and not quiet: if vect1D: sol = [x, y] err = [np.abs(dx), np.abs(dy)] @@ -692,31 +880,31 @@ accuracy after 3 iterations. sol = np.dstack( [x, y] )[0] err = np.dstack( [np.abs(dx), np.abs(dy)] )[0] - if ind is None: - ind, = np.where(dx**2+dy**2 >= accuracy2) - - if inddiv is None: - inddiv, = np.where(dn2[ind] > dn2prev[ind]) - - if ind.shape[0] == 0: - ind = None - inddiv = None - - elif inddiv.shape[0] == 0: - inddiv = None - - assert(ind is not None or inddiv is not None) # <-- sanity check + # restore previous numpy error settings: + np.seterr(invalid = old_invalid, over = old_over) + + if inddiv is None: + raise NoConvergence("'HSTWCS.all_sky2pix' failed to " \ + "converge to the requested accuracy after {:d} " \ + "iterations.".format(k), best_solution = sol, \ + accuracy = err, niter = k, failed2converge = ind, \ + divergent = None) + else: + raise NoConvergence("'HSTWCS.all_sky2pix' failed to " \ + "converge to the requested accuracy.{0:s}" \ + "After {1:d} iterations, the solution is diverging " \ + "at least for one input point." \ + .format(os.linesep, k), best_solution = sol, \ + accuracy = err, niter = k, failed2converge = ind, \ + divergent = inddiv) - raise NoConvergence("'HSTWCS.all_sky2pix' failed to converge to "\ - "requested accuracy after {:d} iterations." \ - .format(k), best_solution = sol, \ - accuracy = err, niter = k, \ - nonconvergent = ind, divergent = inddiv) + # restore previous numpy error settings: + np.seterr(invalid = old_invalid, over = old_over) if vect1D: return [x, y] else: - return np.dstack([x,y])[0] + return np.dstack( [x, y] )[0] def _updatehdr(self, ext_hdr): |