DM-20371: Create new shift/rot only WCS fitter #127

morriscb · 2019-07-08T22:35:02Z

Initial try at a new, stripped down fitter for use with @parejkoj 's new distortion model code. Not meant for full review but for @parejkoj 's information if he would like to take up the ticket while I am on vacation.

Debug and adjust fitting algorithm. Initial try at new affine fitter with comments.

Remove failing tests. Modify and add tests.

Debug last tests. Fix linting.

TallJimbo

Only big-picture comment is that you've got a lot of snake_case variables and methods, which aren't allowed by our coding conventions.

TallJimbo · 2019-07-24T17:34:28Z

python/lsst/meas/astrom/fitAffineWcs.py

+
+class FitAffineWcsTask(pipeBase.Task):
+    """Fit a TAN-SIP WCS given a list of reference object/source matches.
+    """


I'd mention somewhere here that it fits on top of the cameraGeom distortion model, maybe just by moving some of the docs from fitWcs here.

TallJimbo · 2019-07-24T17:36:30Z

python/lsst/meas/astrom/fitAffineWcs.py

+        # Grab reference coordinates and source centroids. Compute the average
+        # direction and separation between the reference and the sources.
+        # I'm not sure if bearingTo should be computed from the src to ref
+        # or ref to source.


Has this question been resolved? I'd have thought your testing wouldn't go so well if you had it backwards, but maybe I'm misunderstanding what the impact of that would be.

It's resolved. The code is fine as is. There are a few comments round the code that were for @parejkoj to look into if he had time over the 2 weeks I was gone.

Please remove those comments if you haven't already.

Already is.

TallJimbo · 2019-07-24T17:39:50Z

python/lsst/meas/astrom/fitAffineWcs.py

+        # minimize). Fits all current test cases with a scatter of a most 0.15
+        # arcseconds. exits early because of the xTol value which cannot be
+        # disabled in scipy1.2.1. Matrix starting values are non-zero as this
+        # results in better fit off-diagonal terms.


That 0.15 arcseconds strikes me still pretty big (almost a pixel for most cameras we care about), though maybe it's as good as we can get with the camera geometry we have. Are you saying you think it could be better without the xTol/scipy1.2.1 issue, or is that unrelated?

This was from the unittests copied from fitTanSip code which as 1 arcsecond pixels. I think things would be a bit better with the newer versions and removed tolerances as the code tends to exit quickly. With the unittests written in this ticket, the scatter is much much smaller and passes the test criteria without modification.

TallJimbo · 2019-07-24T17:40:51Z

python/lsst/meas/astrom/fitAffineWcs.py

+
+    output_separations = []
+    # Fit both sky to pixel and pixel to sky to avoid any non-invertible
+    # affine matrixes.


👍 good idea!

TallJimbo · 2019-07-24T17:45:03Z

python/lsst/meas/astrom/fitAffineWcs.py

+    """Function to minimize to fit the shift and rotation in the WCS.
+
+    Parameters
+    ----------


Should document x here; seems like a particularly good place to state which parameters are at each index in the array.

TallJimbo · 2019-07-24T19:46:27Z

python/lsst/meas/astrom/fitAffineWcs.py

+
+    def makeWcs(self, crval_offset, rot_matrix):
+        """Apply a shift and rotation to the WCS internal to this class,
+        a new Wcs with these transforms applied.


Missing "return" in this sentence?

TallJimbo · 2019-07-24T19:47:11Z

python/lsst/meas/astrom/fitAffineWcs.py

+        rot_matrix : 'numpy.ndarray', (3, 3)
+            Rotation matrix to apply to the mapping/transform to add to the
+            WCS. Rotation matrix need not be unitary and can therefor
+            stretch/squish as well as rotate.


If it's not unitary I think we probably shouldn't call it a rotation matrix.

TallJimbo · 2019-07-24T19:52:36Z

python/lsst/meas/astrom/fitAffineWcs.py

+        self.frame_max = np.max(self.frame_idxs)
+
+        # Find frame just before the final mapping to sky and store those
+        # indices and mappings for later.


I was hoping @parejkoj's work would let us rely on some named frame we could expect to be present in the initial WCS, rather than assuming the mapping to sky is implemented by some fixed number of actual AST mapping instances. If that's not the case, is it an oversight or is there something that makes it hard?

I mean that could be, however for now to use with the generic stack produced WCS, this is what I came up with to catch most cases. Most of the code in the stack seems to flatten the transform as pixels->IWC->Sky. If this changes and Parejko's new code will have more Frames, the code above will still work.

👍 Working with generic stack-produced WCSs is probably more important than named frames vs indices. But I'd also like to get @parejkoj's thoughts on this, if only to wonder whether to bring this up again after this ticket is merged: will we have a named frame that can consistently be used for this purpose in all WCSs we produced (including, but not limited to the new initial ones)?

I dislike this approach, too. makeSkyWcsFrameDict (which is used under the hood to construct the new "initial" WCS) loses the names of the intermediate frames; we would normally want to put the new Transform right after FIELD_ANGLE, but that becomes IWC (which is the only "intermediate" frame that FITS WCS knows about). I filed DM-20315 about this issue.

At minimum, we should file a ticket triggered on DM-20315 to cleanup this code to insert the transformations between FIELD_ANGLE and IWC explicitly. That would make it incompatible with SkyWcs objects created in other ways, but I think that's fine: we're already well past the point of being "true FITS WCS" compatible with this work anyway.

I think we probably could make the names explicit here by explicitly inserting the new Transform between IWC and Sky, and correcting that to be between FIELD_ANGLE and IWC once 20315 is done.

TallJimbo · 2019-07-24T19:59:23Z

python/lsst/meas/astrom/fitAffineWcs.py

+                               crval_offset[1] * arcseconds),
+            np.array([[1., 0.], [0., 1.]]))
+        iwc_to_sky_map = iwcs_to_sky_wcs.getFrameDict().getMapping("PIXELS",
+                                                                   "SKY")


The fact that you need to make a bogus WCS just to extract this mapping from it suggests we should have a direct way to construct that mapping available in afw. Would you mind creating a ticket to add that and then update this code accordingly?

I think this requires a deeper change the C++ API of SkyWcs and the Sky frame. Sky frame is pretty special it seems and the code you see here is taken from unittests others wrote to create exactly this case. I believe the one I took here was written by @parejkoj and used in jointcal?

Yeah, I'm not advocating addressing this on this ticket. But if we've got the same workaround appearing in jointcal as well, that'd make it even more important to factor it out into its own routine.

Russell suggested to me that the easiest way to create the correct IWC->SKY Mapping is to create a blank SkyWcs centered at the CRVAL you want and extract the Mapping and Frame from that. I don't know if we need another way to go about it or not: this is only two lines, and is a pretty specialized use.

This is the relevant code in jointcal:

https://github.com/lsst/jointcal/blob/master/src/ConstrainedAstrometryModel.cc#L234

TallJimbo · 2019-07-24T20:00:21Z

tests/test_fitAffineWcsTask.py

@@ -0,0 +1,274 @@
+# LSST Data Management System
+# Copyright 2008, 2009, 2010 LSST Corporation.


Should copy this copyright blurb from a newer file (like one of the new modules you added).

morriscb · 2019-07-25T01:05:00Z

Take a look at the changes. I've changed all the method and class variables to camel case. Is that enough?

I have one last improvement I thought of as the initial guess for the angle and distance for the fitter is computed naively and wrong. I'll have to do some testing tomorrow on this.

TallJimbo · 2019-07-25T14:06:06Z

I've changed all the method and class variables to camel case. Is that enough?

Not up to me. I don't think that as a reviewer I have the authority to grant exceptions to the style guide for no particular reason, especially in new code.

Sorry - I hate playing style police and I'm not a huge fan of camelCase either, and I personally very much agree that making the public interfaces conform is what matters most. If you want to RFC allowing more snake_case, I think it'd probably get a lot of support (and depending on how we handle migrating existing code (or not) I'd probably get behind it, too).

Add statement that this fitter is to be used in conjunction with a cameraGeom distortion model. Remove debugging comments meant for earlier inspection. Add description of x to chi_func Remove inherit from object. Remove mentions of shift/rot matrix. Update copyright. Change all class and api variables to camel case. convert variables to camelCase.

Make offsetDir calculation safe. Debug direction calculation. Remove print statements.

morriscb · 2019-07-25T22:23:00Z

Just finished up. Take another look if you want and I'll merge and close the ticket.

morriscb force-pushed the tickets/DM-20371 branch from 562905a to 40aa9e3 Compare July 22, 2019 18:06

morriscb added 3 commits July 23, 2019 11:30

Add new FitAffineWcs task.

839654e

Debug and adjust fitting algorithm. Initial try at new affine fitter with comments.

Add log statements and fix matchList update.

fcdd8ff

Add specific affine test.

9ebe71f

Remove failing tests. Modify and add tests.

morriscb force-pushed the tickets/DM-20371 branch from 95aeb4c to 1d146ed Compare July 23, 2019 18:30

Add reverse transform unittests.

a808119

Debug last tests. Fix linting.

morriscb force-pushed the tickets/DM-20371 branch from eda1b74 to a808119 Compare July 23, 2019 18:39

TallJimbo approved these changes Jul 24, 2019

View reviewed changes

morriscb force-pushed the tickets/DM-20371 branch 2 times, most recently from 2997d0e to deb3b05 Compare July 25, 2019 22:12

Modify initial offset calculation.

e9a4a67

Make offsetDir calculation safe. Debug direction calculation. Remove print statements.

morriscb force-pushed the tickets/DM-20371 branch from f9e1c28 to e9a4a67 Compare July 25, 2019 22:19

TallJimbo approved these changes Jul 26, 2019

View reviewed changes

Add TODO for changing which FRAME the trans is appended to.

38b992f

morriscb merged commit 38b992f into master Jul 29, 2019

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

DM-20371: Create new shift/rot only WCS fitter #127

DM-20371: Create new shift/rot only WCS fitter #127

morriscb commented Jul 8, 2019

TallJimbo left a comment

TallJimbo Jul 24, 2019

TallJimbo Jul 24, 2019

morriscb Jul 25, 2019

TallJimbo Jul 25, 2019

morriscb Jul 25, 2019

TallJimbo Jul 24, 2019

morriscb Jul 25, 2019

TallJimbo Jul 24, 2019

TallJimbo Jul 24, 2019

TallJimbo Jul 24, 2019

TallJimbo Jul 24, 2019

TallJimbo Jul 24, 2019

morriscb Jul 25, 2019

TallJimbo Jul 25, 2019

parejkoj Jul 26, 2019

TallJimbo Jul 24, 2019

morriscb Jul 25, 2019

TallJimbo Jul 25, 2019

parejkoj Jul 26, 2019

parejkoj Jul 26, 2019

TallJimbo Jul 24, 2019

morriscb commented Jul 25, 2019

TallJimbo commented Jul 25, 2019

morriscb commented Jul 25, 2019

		@@ -0,0 +1,274 @@
		# LSST Data Management System
		# Copyright 2008, 2009, 2010 LSST Corporation.

DM-20371: Create new shift/rot only WCS fitter #127

DM-20371: Create new shift/rot only WCS fitter #127

Conversation

morriscb commented Jul 8, 2019

TallJimbo left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

morriscb commented Jul 25, 2019

TallJimbo commented Jul 25, 2019

morriscb commented Jul 25, 2019