From: fmasci@ipac.caltech.edu Subject: Re: attempt to back out scales using CD breakdown Date: January 30, 2010 11:21:54 PM PST To: astrocat1@gmail.com, jwf@ipac.caltech.edu, tim@ipac.caltech.edu, hlm@ipac.caltech.edu, jarrett@ipac.caltech.edu, carl@ipac.caltech.edu, roc@ipac.caltech.edu Hi all, Apologies if the below comes across as obvious, or has been discussed in previous email threads. I sense that concepts are not yet fully grasped. Here are some facts as well as suggestions on how to proceed with minimal code changes, assuming of course we don't yet have a clear path. (1) Following the formalism on the distortion webpage, the absolute pixel scales (s1, s2) cannot be explicitly extracted from the f(u,v) and g(u,v) polynomial fits. However, one can *refine* prior s1,s2 values using the linear terms in these polynomials, see below. These polynomials only pertain to pixel space. The CD matrix only gets us to the sky and vice-versa, on a per-frame basis, and that's it. (2) By definition, the pixel scales s1,s2 pertain to the Coordinate Reference PIXels only (CRPIX1,2), usually taken as the center of the FOV. Deviations in scale towards the edges can arise from (i) map projections (e.g., gnomonic, orthographic, etc..) and are generally handled by WCS s/w libraries. This is negligible for small fields, including WISE. (ii) real focal plane distortion will change effective scales over the FOV. But remember, there is only one s1,s2 in the SIP representation. Parameters describing position dependent scale changes are handled by the distortion representation and the (negligible) map projection geometry. (3) To recap, the linear terms in the distortion polynomials are A01, A10, B01, B10 and are extracted from: Delta_u = f(u,v) = A01*v + A10*u Delta_v = g(u,v) = B01*v + B10*u where u = x_distorted - CRPIX1 and v = y_distorted - CRPIX2. Furthermore, I've assumed for simplicity that skew=0, so that the corrections Delta_u, Delta_v are precisely equal to f, g respectively. I am also ignoring the offset terms (A00, B00) since they represent residuals in positioning of the center reference pixel, and should vanish following absolute pointing reconstruction/refinement. The constant and linear coefficients are included when fitting for f and g so they can absorb any dependencies that cannot be explained by the higher order terms. Including them generally improves the higher order estimates. (4) In the SIP convention, the A10, B01 coefficients in (3) are *implicitly* represented by s1,s2 when they are *precisely zero*. Of course, right. Also, the linear cross terms A01, B10 are taken care of by the skew, precisely when these are zero. If any of these coefficients significantly differ from zero, then they must be included as part of the distortion polynomial in FITS headers, otherwise the distortion representation is incomplete. The usual working assumption is that these linear terms are ~zero and can be dropped, but they're important if they're not zero because it means the s1,s2 alone cannot explain the whole FOV. This is particularly true if you're dealing with a large FOV like WISE. Therefore, an important sanity check (as I mentioned last week) is to check that the linear terms are ~0 to warrant having the linear part described solely by s1,s2 (and skew) in the CD matrix. (5) Given estimates of s1,s2 (e.g., from priors or SFPRex), one can *refine* these using the A10,B01 coefficients from a distortion fit as follows. For a skewless FOV (an ok assumption for now, otherwise I'll have to hit the Maple code), the linear change in pixels caused by distortion across the FOV along orthogonal x and y directions from the center is: Delta_u = A10*u, Delta_v = B01*v, which really means: x_true - x_distorted = A10*(x_distorted - CRPIX1), y_true - y_distorted = B01*(y_distorted - CRPIX2) respectively. Again, if s1,s2 were ~correct, then Delta_u ~ 0 and Delta_v ~ 0. If we call the old (or prior) pixel scales s1_old, s2_old and the new (refined) values s1_new, s2_new, then the following relationships are also true: x_true - x_distorted = (s1_new - s1_old)*(x_distorted - CRPIX1) y_true - y_distorted = (s2_new - s2_old)*(y_distorted - CRPIX2) Combining with the above, This means that: s1_new = [1 + A10 ]*s1_old; s2_new = [1 + B01]*s2_old; A subsequent distortion fit with these new scales should give ~zero linear coefficients in the poly fit and *better* estimates for the higher order terms. If not, then it's another round of refinement as described above. For relatively large FOV's (like here), there's greater potential for any distortion derivation to be uncertain. That's because the higher order distortion terms depend sensitively on what was assumed for scales in the first place (i.e., a chicken-and-egg type of problem). An uncertain scale gives progressively larger (albeit systematic) uncertainties in the mapped 2MASS positions towards a frame edge. (6) Now to SFPRex. This must assume some prior distortion, scales, and twist (for a given frame) to map 2MASS references to pixel space. The 5 param fit corrects scales, twist, and absolute pointing with the assumption that the prior distortion has been characterized as best it can over the full FOV. A bad distortion assumption (in particular, for the higher-order terms) will corrupt all parameter estimates from SFPRex. The least corrupted will be corrections to the absolute pointing. To get better estimates of s1,s2, my suggestion is to run SFPRex in 5 parameter, single-band mode but only using the central 400-500 hundred pixels when matching WISE-to-2MASS to minimize the effects effects of the "uncertain" (non-linear) distortion. The polynomials look ~linear in that region for the most part. The caveat is that there'll be fewer matches per frame, but many frames will need to be processed and trended to get some workable scales. Remember, the goal is to get better estimates of s1,s2 to kick things off. These don't have to be super accurate at first, because as outlined in (5) above, they can be refined by factoring out the linear coefficients of the polynomial fits. The pixel-usage window for SFPRex can then be opened up of course when scales and distortion are known to sufficient accuracy. How well will be determined from the residuals of polynomial fits. It's going to be difficult, but I'm confident we'll get there. Again, these are large(ish) FOVs, and the scale only needs to be in error by ~0.4% (or ~ +/-0.01 arcsec) to give a residual of ~2 pixels at the edge of frame (considering only the linear part of the distortion). There's not much margin for error. Regards, Frank ----------------------------------------------------------------------------- From: fmasci@ipac.caltech.edu Subject: Re: attempt to back out scales using CD breakdown Date: January 31, 2010 10:34:53 AM PST To: jwf@ipac.caltech.edu Cc: astrocat1@gmail.com, tim@ipac.caltech.edu, hlm@ipac.caltech.edu, jarrett@ipac.caltech.edu, carl@ipac.caltech.edu, roc@ipac.caltech.edu The constant terms should be ~zero if there were no global bias in the WISE-to-reference matches over the FOV, or the match separation vectors were ~globally symmetric about CRPIX1,2. But it's okay if these are not zero since they will absorb any dependence that cannot be explained by the linear+higher order terms. The assumption is that this bias will be taken out by future SFPRex runs, i.e., when the correct RA,Dec is assigned to the CRPIX1,2, and the distortion has converged to a state to enable unbiased matches over the full FOV in the first place. So, you can see that it's iterative. You can't ensure bias free solutions from SFPRex without first nailing the distortion, everywhere! But don't forget to check those constant + linear terms in your distortion fits. As these get smaller, you know were getting closer and closer, and SFPRex is being fed the right fodder. Regards, Frank I'm not sure how many iterations it will take so that A10, B01 are low enough to be confident that the final s1,s2 will not mess up the non-linear part of the distortion in subsequent fits. Should we be concerned at all that the constant term for X (it would be A_0_0 in SIP notation if such a term existed) keeps coming out around 0.6 to 0.7 pixel? The constant term for Y tends to run around 0.1 pixel. Shouldn't these converge toward zero? Regards, John ----------------------------------------------------------------------------- From: fmasci@ipac.caltech.edu Subject: Re: back to the frame level Date: February 7, 2010 9:30:40 AM PST To: jwf@ipac.caltech.edu, tim@ipac.caltech.edu, jarrett@ipac.caltech.edu, roc@ipac.caltech.edu, astrocat1@gmail.com, hlm@ipac.caltech.edu, carl@ipac.caltech.edu John is correct. The linear cross terms soak up any slack in a rotation residual from a distortion fit, just like the constant terms soak up any global bias in reference vs. wise x,y position differences from imperfect WCS solutions. The goal is to obviously get these as small as possible after repeated SFPrex -> gnDSTR runs. I have derived some relations to refine scales using _all_ the linear terms. These nicely reduce to the simple refinement relations emailed last week. The main relations are now: s1new = s1old + sqrt[ (s1old*A10)**2 + (s2old*B10)**2 ]; s2new = s2old + sqrt[ (s2old*B01)**2 + (s1old*A01)**2 ]; where A01, B10 are the cross terms in zero-index based notation (of course, why use any other notation!). One can also estimate the rotation error from the linear terms (see attached). If you're interested in the derivation, here it is. Apologies for the font. http://wise2.ipac.caltech.edu/staff/fmasci/figures_codeVdist/refined_scales.pdf Regards, Frank Plus the cross terms are sizeable, which I don't believe we know how to interpret. I have found that to relate to a rotation error. The cross-coefficients are basically the cosine and sine of the rotation, at least when there's no scale-factor error. If there are both, then they get blended. Regards, John ------------------------------------------------------------------------------- From: fmasci@ipac.caltech.edu Subject: Re: back to the frame level Date: February 8, 2010 2:46:29 PM PST To: astrocat1@gmail.com Cc: hlm@ipac.caltech.edu Right there's an ambiguity in using those relations. The sqrt hides the +/- direction. Those relations just reproduce the case when cross terms are exactly zero, giving: s1new = s1old*(1 + A10); s2new = s2old*(1 + B01); and the sign of the A10 and B10 determines the direction of the refinement. However, here are some new relations that use all the linear coefficients: s1new = s1old*(1 + A10) - s2old*B10*Tan(theta); s2new = s2old*(1 + B01) + s1old*A01*Tan(theta); Tan(theta) = -x / [1 + sqrt(1 + x**2)], where: x = (s2old*B10)/(s1old*A10); what I'm not sure about now is the sign of sqrt in the Tan(theta) formula. I bet its direction can be calibrated using a simple test case. I'm going to bet it's a "+". More simply, try both +/- for the sqrt in the Tan(theta) formula and pick the one that gives you closer to your expectation. When this is done, there should be no ambiguity in using these new relations. Thanks for mentioning that. I will update my documentation. Frank Frank- Haven't been through your derivation as yet, but it doesn't appear to me that there's a mechanism in this formulation for scales to adjust downward. That is, the adjustment term can never be negative. Am I reading it wrong? Regards, Howard ----------------------------------------------------------------------------- From: astrocat1@gmail.com Subject: backing scales out of the gnDSTR fit Date: February 12, 2010 4:06:49 PM PST To: fmasci@ipac.caltech.edu Cc: hlm@ipac.caltech.edu, jwf@ipac.caltech.edu, tim@ipac.caltech.edu, jarrett@ipac.caltech.edu, carl@ipac.caltech.edu, roc@ipac.caltech.edu Frank- Finally got a chance to go through your approach to backing the scales out of the linear terms in the gnDSTR fit. http://wise2.ipac.caltech.edu/staff/fmasci/figures_codeVdist/refined_scales.pdf Thanks again for taking the time to do it. Nice derivation. It's great that the expressions for scale change in equations 13 & 14 come out so clean. My concern, as I mentioned in our discussion a few days back, is whether, when we get down close to the final scales, the uncertainties in the coefficients need to be accounted for in the solution. As a test, I went to a case which came up back on January 31 when we were trying to nail down the W1 scales: A12= 0.00064466 A21= -0.0010697729 B12= -0.0006919832 B21= 0.0003678687 Note that the cross terms are about half the magnitude of the primary terms, so it should serve as a good test case. The TRSrvb solution for this case was: dsx = -0.0016084889; dsy = -.0006006777 These are fractional scale changes. In the application of your equations which follows I assumed the apriori scales could be taken to be 1.0, making the computed scale changes fractional. I used your equation 13 to solve for dsx & 14 to get dsy. Eqn 13&14 => dsx= 0.001131 or -0.001131; dsy= 0.000946 or -0.000946 These are not terribly different from the TRSrvb results. Once the scales are known, any one of the equations 9, 10, 11 or 12 can be used to solve for the rotation angle. I did so separately for each equation to see how consistent the answers would be. Since the signs on dsx & dsy are ambiguous, I did it both ways. For each equation below I give the derived tangents, followed by the rotation angles: Eqn 9 => tan= -5.983193 or 0.167135; the= 99.488419 or -170.511581 Eqn 10 => tan= -1.115574 or 0.490754; the= 131.873036 or 26.139697 Eqn 11 => tan= 0.167135 or -5.983193; the= 9.488419 or 99.488419 Eqn 12 => tan= -0.540183 or 1.227934; the= 151.622844 or -129.158549 The angles are not only inconsistent, but huge. Assuming I haven't made a mistake, something's going on here. I'm thinking we might need to go to a chi-square fit which factors uncertainties on the coefficients into the solution. Or possibly as you mentioned, unintended skew is somehow slipping in. Regards, Howard -----------------------------------------------------------------------------