Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix off-by-one Error and related WCS changes #276

Merged
merged 47 commits into from
Oct 3, 2023
Merged

Conversation

teutoburg
Copy link
Contributor

@teutoburg teutoburg commented Sep 22, 2023

Two important bugs were fixed

  • FITS headers are 1-indexed, meaning the first pixel is [1, 1], unlike the 0-indexed nature of Python [0, 0]. Previously, there were several occurrences of this not being taken into consideration, which is fixed by this PR. Pixel coordinates which are written to a fits.Header must always follow the 1-based convention. This also ensures headers are correctly interpreted by other software.
  • Related to this, the center of an image with an even number of pixels along at least one axis cannot be at an integer pixel number, but rather must fall on a half-pixel. This was previously also ignored in many cases.

Other changes for everything to work as expected

  • In some places, it was necessary to convert float to int in a different way. I.e., int(), np.ceil() and round() are all used, but great care must be taken in deciding which is the correct function for which case.
  • Some of the code depends on the simplification that the conversion between pixel and world coordinates (and vice versa) is a linear operation. This is the case for e.g. detector coordinate systems, but not for sky coordinates, which usually use TAN projection. Previously, sky coordinate information in headers was almost always given as TAN, but the implementation of the conversion functions just ignored that and converted everything linearly. Changing this behavior however breaks a number of things, because the linear outcome is expected. A more permanent (and geometrically accurate) solution is needed here, but for the time being, two changes were introduced as a stop-gap solution:
    • FITS headers that are created dynamically are initialized with LINEAR wcs axes, even if they are dealing with sky coordinates.
    • A new wcs_suffix X is introduced, which uses sky coordinates, but with LINEAR projection. This is currently only used for mock objects in testing.

As a consequence of these changes, many test cases and mock objects had to be adapted, as they were based on erroneous assumptions and sometimes "bent" to pass with an actually wrong result.

Additional changes

  • The output format of optics.image_plane_utils.calc_footprint() was changed to match that of astropy.wcs.WCS.calc_footprint(), i.e. a single array of shape (4, 2) containing the corner points in clockwise order starting from the bottom left. A new function optics.image_plane_utils.calc_table_footprint() is introduced to do the same for astropy.table.Table objects.
  • In a number of functions (including the one just mentioned), astropy.wcs.WCS is now used internally, instead of custom implementations. In terms of performance, it seems that once the WCS instance is created, the actual operations are faster then the self-implemented ones. The only drawback is the performance hit that results from creating the WCS objects in the first place. This can be mitigated somewhat by keeping them around as much as possible, instead of creating them from scratch every time, in cases where this is possible. Sometimes, it may be sufficient to pass the WCS instead of a whole header.
  • Several operations were parameterized, or "array-ized", which results in a number of simplifications.
  • Many small changes and improvements to auxiliary functions, mostly but not exclusively in optics.image_plane_utils.

What is still to do

  • More consistency for coordinate-based operations is needed in general.
  • Tolerances in various tests should be re-evaluated.
  • Some tests use random numbers, think about using a fixed seed for reproducability.
  • Both deg and arcsec are now used internally, which actually caused the most headache in this PR. The checks and conversions need to be more robust!
  • The mentioned change to the output format of calc_footprint() necessitated a number of adaptions in the code. Those were done very crudely and could be improved, in particular many tests use this as an optional feature for plotting result footprints, all of this could be factored out.
  • astropy.wcs.WCS.calc_footprint() has an optional argument center, which controls if the footprint is considered to go until the center of the corner pixels (where the coordinates are defined) or to the very edge, adding one full CDELT to the whole image. We need to decide which version to use in which cases.
  • As mentioned above, WCS instances should be passed along as much as possible for performance.
  • We are using "canvas" headers and HDUs a lot. Maybe think about a Canvas class.
  • Make full use of the features available through astropy.wcs.WCS.
  • As discussed in other places, use less units internally.
  • The function optics.fov_utils.is_field_in_fov (and that whole module) could use some refactoring.
  • In this PR, the auxiliary function optics.image_plane_utils._fix_360 was introduced to deal with "overflowing" values of deg, i.e. e.g. 359.995 instead of -0.005. This needs some more attention, unit tests, and also think about how to deal with deg vs. arcsec in this context.

Before, this was just a multiplication by the zoom factor,
which can produce wrong results even in simple cases.

The solution in place now will probably only work for LINEAR WCSs.
…ntral

pixel, which gets all the flux. However, this was a wrong assumption
based on an effect of the off-by-one error. Given an even number of pixels
along both axes of the canvas, a source perfectly in the center should be
spread out over the 4 central pixels, which now works. With this change,
the test still tests the exact same thing, just in a more general and
correct way.
Change API of calc_footprint to return corner in the same format as
astropy.wcs.WCS.calc_footprint, which is now used internally.
The affected parts could be improved, but for now just get it to work.
@teutoburg teutoburg marked this pull request as ready for review September 26, 2023 14:12
@teutoburg teutoburg changed the title [DRAFT] Fix off-by-one Error and related WCS changes Fix off-by-one Error and related WCS changes Sep 26, 2023
@teutoburg
Copy link
Contributor Author

This should fix #209

Revert changes to extract_area_from_table and change mock data instead
Copy link
Collaborator

@hugobuddel hugobuddel left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

https://github.com/AstarVienna/ScopeSim/blame/a54514e4ab737cf836b96dcfc65530820a43da96/scopesim/effects/psfs.py#L627C18-L627C18

                kernel_pixel_scale = hdr["CDELT1"] * unit_factor

now fails, because CDELT1 is not there, but CDELT1D is. I think you made some helper function for this; I'll have a look

Sorry for not being very constructive, I'll get something better.

Co-authored-by: Hugo Buddelmeijer <hugo@buddelmeijer.nl>
@hugobuddel hugobuddel self-requested a review September 27, 2023 09:59
@teutoburg teutoburg merged commit 7da0dd3 into dev_master Oct 3, 2023
19 checks passed
@teutoburg teutoburg deleted the fh/off-by-one branch October 3, 2023 16:28
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
Archived in project
Development

Successfully merging this pull request may close these issues.

WCS in Scopesim, off-by-one error, specifically image_plane_utils.header_from_list_of_xy
2 participants