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

Replace griddata with RegularGridInterpolator for distortions #429

Merged
merged 7 commits into from
Mar 29, 2021

Conversation

JarronL
Copy link
Collaborator

@JarronL JarronL commented Mar 23, 2021

Previously, the add_distortion function used scipy's griddata to transform from an irregularly gridded 'sci' coordinate array to a regularly gridded 'sci' coordinate array. However, the input 'sci' coordinates are derived from a regularly gridded 'idl' array. Since, it would be much faster to use RegularGridInterpolator (RGI), we could instead start with these 'idl' coordinates and interpolate onto some set of irregular 'idl' points, which map to our desired 'sci' coord values.

I've written a generalized distortion function that uses this concept to transform an image between 'idl' coordinates to some other pysiaf frame ('sci', 'tel', 'det'). This is useful for other purposes, like distorting model images (e.g, disks) into the detector frame before convolving with PSFs.

Benchmarks on my laptop for fov_pixels=512 and oversample=4 for a single PSF:

  1. without distortion: ~11 sec
  2. with original distortion (griddata): ~75 sec
  3. with updated distortion (RGI): ~13 sec

While the original implementation added a factor of ~6x to the PSF calculation time, the RGI only takes a 10-20% overhead hit. The resulting distortions calculations are also very similar:

Above is the difference of the corner PSFs relative to the PSF at the center of the detector (for A5, F405N). I've set include_si_wfe=False so that the only change would be from geometric distortions rather than WFE differences.

So, I think solves the excessive distortion overhead in #426.

Remove constraints on detector_pixel that it should be an integer from 0 to npix. Since these values correspond to SIAF aperture 'sci' coords and only get converted to 'tel' (V2/V3), there's no reason to constrain to the detector focal plane, especially since it is possible to simulate PSFs outside this region via SI WFE extrapolation. Also, fractional values are necessary to track small shifts (e.g., coronagrahic masks) and provide sub-pixel coordinate transformations.
Explicitly set x,y=map() in order to force map to return the values, otherwise it gets saved as an object, and we are unable to subscript the values as expected.
@JarronL JarronL requested a review from mperrin March 23, 2021 21:50
@codecov
Copy link

codecov bot commented Mar 23, 2021

Codecov Report

Merging #429 (a0a69ab) into develop (ce518be) will decrease coverage by 0.12%.
The diff coverage is 60.00%.

Impacted file tree graph

@@             Coverage Diff             @@
##           develop     #429      +/-   ##
===========================================
- Coverage    48.39%   48.26%   -0.13%     
===========================================
  Files           14       14              
  Lines         5837     5865      +28     
===========================================
+ Hits          2825     2831       +6     
- Misses        3012     3034      +22     
Impacted Files Coverage Δ
webbpsf/distortion.py 81.81% <60.00%> (-13.42%) ⬇️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update ce518be...a0a69ab. Read the comment docs.

@JarronL JarronL added enhancement New feature or request JWST Affects JWST models in WebbPSF labels Mar 23, 2021
@mperrin mperrin assigned JarronL and unassigned mperrin Mar 25, 2021
Copy link
Collaborator

@mperrin mperrin left a comment

Choose a reason for hiding this comment

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

This looks great! I have a few very minor change requests; it looks like there's a couple extra unused imports, and one parameter missing in the doc string.

Thanks very much for doing this.

webbpsf/distortion.py Outdated Show resolved Hide resolved
webbpsf/distortion.py Outdated Show resolved Hide resolved
Comment on lines +76 to +91
xnew_coords : None or ndarray
Array of x-values in new coordinate frame to interpolate onto.
Can be a 1-dimensional array of unique values, in which case
the final image will be of size (ny_new, nx_new). Or a 2d array
that corresponds to full regular grid and has same shape as
`ynew_coords` (ny_new, nx_new). If set to None, then final image
is same size as input image, and coordinate grid spans the min
and max values of siaf_ap.convert(xidl,yidl,'idl',to_frame).
ynew_coords : None or ndarray
Array of y-values in new coordinate frame to interpolate onto.
Can be a 1-dimensional array of unique values, in which case
the final image will be of size (ny_new, nx_new). Or a 2d array
that corresponds to full regular grid and has same shape as
`xnew_coords` (ny_new, nx_new). If set to None, then final image
is same size as input image, and coordinate grid spans the min
and max values of siaf_ap.convert(xidl,yidl,'idl',to_frame).
Copy link
Collaborator

Choose a reason for hiding this comment

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

What are the use cases for passing in coordinate arrays here? I can see this adds flexibility, but also complexity. I'm curious if you have particular needs or scenarios in mind that will use this.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Useful for plotting purposes onto custom grid spacing. I am currently using this for visualizing model images onto already existing V2/V3 coordinate grids. If you want, I can remove these options here, since I've actually already made a slightly different version of this function for webbpsf_ext.

Copy link
Collaborator

Choose a reason for hiding this comment

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

No, definitely don't remove. I could tell you had put in the option for some reason, and wanted to understand what it was.

webbpsf/distortion.py Show resolved Hide resolved
Copy link
Collaborator

@mperrin mperrin left a comment

Choose a reason for hiding this comment

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

Looks good!

I checked the timing and results of a 3x3 PSF grid calculation using NIRISS.

  • With current develop, took 46.6 s.
  • With this PR, takes 25.5 s.
    The resulting PSFs are visually indistinguishable in the view from display_psf_grid.

The test calculation was a 3x3 grid with default properties per each PSF (using 21 wavelengths and a 101x101 PSF 4x oversampled). Hence the distortion part is only part of the total runtime. In any case, this suffices to verify a nice speedup from this PR. Thanks again @JarronL

@mperrin mperrin merged commit ad9c2de into develop Mar 29, 2021
@mperrin mperrin deleted the distortion_grid branch November 18, 2021 18:47
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request JWST Affects JWST models in WebbPSF
Projects
None yet
Development

Successfully merging this pull request may close these issues.

WebbPSF's distortion model is really slow to compute - we should improve it
2 participants