Uploaded image for project: 'Data Management'
  1. Data Management
  2. DM-15431

Add Jacobian BoundedField accessor to SkyWcs

    Details

    • Type: Story
    • Status: Done
    • Resolution: Done
    • Fix Version/s: None
    • Component/s: afw
    • Labels:
      None

      Description

      We've agreed that it should not be the responsibility of the PhotoCalib to include a Jacobian derived from the WCS, because we don't always want to apply that.  Instead, SkyWcs should have the ability (given a pixel bounding box) to return a BoundedField that evaluates to the determinant of the Jacobian of the pixelToSky transform (aka the pixel area in angular units) at every point in that bounding box efficiently.  This probably involves a new BoundedField class (though it's possible there's a way to do this with TransformBoundedField).

      A reference implementation can be found in meas_mosaic's calculateJacobian functions, though a different approach to vectorization may be necessary to get good performance from AST.

        Attachments

          Issue Links

            Activity

            Hide
            Parejkoj John Parejko added a comment -

            We may be able to use skyWcs.getTransform().getJacobian() to accomplish this. Whether it is performant is a different question.

            Show
            Parejkoj John Parejko added a comment - We may be able to use skyWcs.getTransform().getJacobian() to accomplish this. Whether it is performant is a different question.
            Hide
            Parejkoj John Parejko added a comment -

            Summary of what I've learned so far regarding the determinant of the jacobian of the WCS (hereafter det(J)):

            meas_mosaic.calculateJacobian(point) uses the square of the ratio of getPixelScale(point)/getPixelScale(origin). Conceptually, this is the value that we're interested in: how much a pixel is distorted at some pixel location vs. at the "center" of the WCS.

            Alternate methods for computing this same value for a SkyWcs object (one would again take the ratio at the point vs. the pixel origin):

            1.
            np.linalg.det(wcs.getTransform().getJacobian(point))}}
             
            2.
            linearMap = wcs.getFrameDict().getMapping("pixels", "iwc")
            t = lsst.afw.geom.TransformPoint2ToPoint2(linearMap)
            detJ = np.linalg.det(t.getJacobian(point))
             
            3.
            linear = wcs.linearizePixelToSky(point, lsst.geom.radians)
            detJ = np.linalg.det(linear.getMatrix())
            

            When compared with the pixel scale method, 1. gives results that can differ at the percent-ish level, 2. gives results that are closer, but slowly start to diverge as the point of interest is further from the origin, and 3. seems to stay within floating point precision even for points far from the origin. For example, using a skyWcs loaded from the testdata_jointcal, these are the ratios of the three methods above to the pixel scale calculation:

            In [37]: print(compare(skyWcs, point3, rawJacobianFull)) 
                ...: print(compare(skyWcs, point3, rawJacobianIwc)) 
                ...: print(compare(skyWcs, point3, rawJacobianLinearized))                                     
            1.0027808211807387
            1.0000102561363005
            0.9999999954468474
            

            I have been trying to come up with a comparable calculation using the gtransfo object in jointcal itself, but have not yet been able to produce any similar numbers. If I were to implement det(J) as an explicit non-fitted term in the photometric calibration, I would have done it similar to the "IWC" calculation above: take the partial derivatives of the AstrometryModel, which goes from pixel->tangent plane (where the tangent plane is the IWC in the resulting SkyWcs object). But this seems to produce a very different det(J).

            Show
            Parejkoj John Parejko added a comment - Summary of what I've learned so far regarding the determinant of the jacobian of the WCS (hereafter det(J)): meas_mosaic.calculateJacobian(point) uses the square of the ratio of getPixelScale(point)/getPixelScale(origin) . Conceptually, this is the value that we're interested in: how much a pixel is distorted at some pixel location vs. at the "center" of the WCS. Alternate methods for computing this same value for a SkyWcs object (one would again take the ratio at the point vs. the pixel origin): 1. np.linalg.det(wcs.getTransform().getJacobian(point))}}   2. linearMap = wcs.getFrameDict().getMapping("pixels", "iwc") t = lsst.afw.geom.TransformPoint2ToPoint2(linearMap) detJ = np.linalg.det(t.getJacobian(point))   3. linear = wcs.linearizePixelToSky(point, lsst.geom.radians) detJ = np.linalg.det(linear.getMatrix()) When compared with the pixel scale method, 1. gives results that can differ at the percent-ish level, 2. gives results that are closer, but slowly start to diverge as the point of interest is further from the origin, and 3. seems to stay within floating point precision even for points far from the origin. For example, using a skyWcs loaded from the testdata_jointcal, these are the ratios of the three methods above to the pixel scale calculation: In [37]: print(compare(skyWcs, point3, rawJacobianFull)) ...: print(compare(skyWcs, point3, rawJacobianIwc)) ...: print(compare(skyWcs, point3, rawJacobianLinearized)) 1.0027808211807387 1.0000102561363005 0.9999999954468474 I have been trying to come up with a comparable calculation using the gtransfo object in jointcal itself, but have not yet been able to produce any similar numbers. If I were to implement det(J) as an explicit non-fitted term in the photometric calibration, I would have done it similar to the "IWC" calculation above: take the partial derivatives of the AstrometryModel, which goes from pixel->tangent plane (where the tangent plane is the IWC in the resulting SkyWcs object). But this seems to produce a very different det(J).
            Hide
            erykoff Eli Rykoff added a comment -

            What is the relative compute time of these four methods? 

            And maybe I'm missing something in why you can't just use method 3 inside jointcal.  Doesn't this only have to be done once to compute these corrections for each star observation, once you've computed the astrometric solution?

            Show
            erykoff Eli Rykoff added a comment - What is the relative compute time of these four methods?  And maybe I'm missing something in why you can't just use method 3 inside jointcal.  Doesn't this only have to be done once to compute these corrections for each star observation, once you've computed the astrometric solution?
            Hide
            rowen Russell Owen added a comment -

            Note that SkyWcs already offers getPixelScale(pixel) which returns the square root of the pixel area (as an Angle so you can pick your favorite units for the area). Whether or not it is presently efficient enough is a valid question, but if it's not, it can be the basis of an interpolation of some kind.

            Show
            rowen Russell Owen added a comment - Note that SkyWcs already offers getPixelScale(pixel) which returns the square root of the pixel area (as an Angle so you can pick your favorite units for the area). Whether or not it is presently efficient enough is a valid question, but if it's not, it can be the basis of an interpolation of some kind.
            Hide
            Parejkoj John Parejko added a comment -

            Jim Bosch: do you think I need to implement the persistence parts of this BoundedField? I feel like it will only be used in-memory.

            Show
            Parejkoj John Parejko added a comment - Jim Bosch : do you think I need to implement the persistence parts of this BoundedField? I feel like it will only be used in-memory.
            Hide
            jbosch Jim Bosch added a comment -

            I don't think it needs to be done immediately, and whether we need it in the end depends on:

            • whether we decide to make the Jointcal photometric fit include it explicitly as a fixed model component;
            • how we want to persist that composite model.

            So I'm happy at least deferring until we know the answers to those questions.

             

            Show
            jbosch Jim Bosch added a comment - I don't think it needs to be done immediately, and whether we need it in the end depends on: whether we decide to make the Jointcal photometric fit include it explicitly as a fixed model component; how we want to persist that composite model. So I'm happy at least deferring until we know the answers to those questions.  
            Show
            Parejkoj John Parejko added a comment - - edited Jenkins run: https://ci.lsst.codes/blue/organizations/jenkins/stack-os-matrix/detail/stack-os-matrix/28985/pipeline
            Hide
            Parejkoj John Parejko added a comment -

            Paul Price: do you mind taking on this medium sized review (~350 new lines)?

            The handful of new lines in jointcal were from my initial testing, but I figured they were useful enough to be added to the python layer.

            If Jira doesn't pick up the afw PR, it's here: https://github.com/lsst/afw/pull/411

            Show
            Parejkoj John Parejko added a comment - Paul Price : do you mind taking on this medium sized review (~350 new lines)? The handful of new lines in jointcal were from my initial testing, but I figured they were useful enough to be added to the python layer. If Jira doesn't pick up the afw PR, it's here: https://github.com/lsst/afw/pull/411
            Hide
            price Paul Price added a comment -

            Reviews completed on the GitHub PRs. Minor comments in afw that should be addressed before merging. I don't understand how the changes in jointcal relate, since they don't use the new functionality in afw, but you say it's useful so I won't stand in your way; nevertheless, the commit messages should be fixed.

            Show
            price Paul Price added a comment - Reviews completed on the GitHub PRs. Minor comments in afw that should be addressed before merging. I don't understand how the changes in jointcal relate, since they don't use the new functionality in afw, but you say it's useful so I won't stand in your way; nevertheless, the commit messages should be fixed.
            Hide
            Parejkoj John Parejko added a comment -

            Thanks for the review.

            Merged and done.

            Show
            Parejkoj John Parejko added a comment - Thanks for the review. Merged and done.
            Hide
            Parejkoj John Parejko added a comment -

            While reviewing fgcm, I discovered SingleFrameJacobianPlugin in meas_base plugins.py. It uses linearizePixelToSky().determinant (my version 2 above), so it should be very close to the pixel area calculation implemented here. However, it does result in the question of whether we need both of these approaches (the plugin itself doesn't appear to be used anywhere).

            Show
            Parejkoj John Parejko added a comment - While reviewing fgcm, I discovered SingleFrameJacobianPlugin in meas_base plugins.py . It uses linearizePixelToSky().determinant (my version 2 above), so it should be very close to the pixel area calculation implemented here. However, it does result in the question of whether we need both of these approaches (the plugin itself doesn't appear to be used anywhere).
            Hide
            jbosch Jim Bosch added a comment -

            I suspect the plugin as it exists now would be better off using something based on getPixelScale() (and perhaps renaming it to make that clear).   For completely different reasons (computing shapes and centroids in sky coordinates from pixel-coordinate measurements) we're going to also want to start saving all of the matrix elements returned from linearizePixelToSky independently, and I imagine we'll do that in the same plugin.

            Show
            jbosch Jim Bosch added a comment - I suspect the plugin as it exists now would be better off using something based on getPixelScale() (and perhaps renaming it to make that clear).   For completely different reasons (computing shapes and centroids in sky coordinates from pixel-coordinate measurements) we're going to also want to start saving all of the matrix elements returned from linearizePixelToSky independently, and I imagine we'll do that in the same plugin.

              People

              • Assignee:
                Parejkoj John Parejko
                Reporter:
                jbosch Jim Bosch
                Reviewers:
                Paul Price
                Watchers:
                Eli Rykoff, Jim Bosch, John Parejko, John Swinbank, Paul Price, Russell Owen
              • Votes:
                0 Vote for this issue
                Watchers:
                6 Start watching this issue

                Dates

                • Created:
                  Updated:
                  Resolved:

                  Summary Panel