From 75ff3078c9cd63717793617cebbd02016d4942d8 Mon Sep 17 00:00:00 2001 From: Jarron Leisenring Date: Fri, 12 May 2023 13:21:28 -0700 Subject: [PATCH 1/7] non-standard pixel sizes for distortion Previously, distortion code assumed that webbpsf and pysiaf apertures had similar pixel sizes. However, it is possible to request any non-standard pixel size when generated a PSF. This fix correctly applies distortion for any arbitrary pixel size. --- webbpsf/distortion.py | 38 +++++++++++++++++++++++++------------- 1 file changed, 25 insertions(+), 13 deletions(-) diff --git a/webbpsf/distortion.py b/webbpsf/distortion.py index 6beb805d..7b41297c 100644 --- a/webbpsf/distortion.py +++ b/webbpsf/distortion.py @@ -63,7 +63,7 @@ def distort_image(hdulist_or_filename, ext=0, to_frame='sci', fill_value=0, Value used to fill in any blank space by the skewed PSF. Default = 0. If set to None, values outside the domain are extrapolated. to_frame : str - Type of input coordinates. + Requested of output coordinate frame. * 'tel': arcsecs V2,V3 * 'sci': pixels, in conventional DMS axes orientation @@ -116,7 +116,6 @@ def distort_image(hdulist_or_filename, ext=0, to_frame='sci', fill_value=0, # Pixel scale information ny, nx = hdu_list[ext].shape pixelscale = hdu_list[ext].header["PIXELSCL"] # the pixel scale carries the over-sample value - oversamp = hdu_list[ext].header["DET_SAMP"] # PSF oversampling relative to detector # Get 'sci' reference location where PSF is observed xsci_cen = hdu_list[ext].header["DET_X"] # center x location in pixels ('sci') @@ -129,16 +128,15 @@ def distort_image(hdulist_or_filename, ext=0, to_frame='sci', fill_value=0, ylin = np.linspace(-1*ny_half, ny_half, ny) xarr, yarr = np.meshgrid(xlin, ylin) - # Convert the PSF center point from pixels to arcseconds using pysiaf - xidl_cen, yidl_cen = aper.sci_to_idl(xsci_cen, ysci_cen) - - # Get 'idl' coords - xidl = xarr * pixelscale + xidl_cen - yidl = yarr * pixelscale + yidl_cen - # ############################################### # Create an array of indices (in pixels) that the final data will be interpolated onto - xnew_cen, ynew_cen = aper.convert(xsci_cen, ysci_cen, 'sci', to_frame) + if to_frame=='sci': + xnew_cen, ynew_cen = (xsci_cen, ysci_cen) + else: + # Choose conversion function + func_convert = getattr(aper, f'sci_to_{to_frame}') + xnew_cen, ynew_cen = func_convert(xsci_cen, ysci_cen) + # If new x and y values are specified, create a meshgrid if (xnew_coords is not None) and (ynew_coords is not None): if len(xnew_coords.shape)==1 and len(ynew_coords.shape)==1: @@ -147,9 +145,18 @@ def distort_image(hdulist_or_filename, ext=0, to_frame='sci', fill_value=0, assert xnew_coords.shape==ynew_coords.shape, "If new x and y inputs are a grid, must be same shapes" xnew, ynew = xnew_coords, ynew_coords elif to_frame=='sci': - xnew = xarr / oversamp + xnew_cen - ynew = yarr / oversamp + ynew_cen + osamp_x = aper.XSciScale / pixelscale + osamp_y = aper.YSciScale / pixelscale + xnew = xarr / osamp_x + xnew_cen + ynew = yarr / osamp_y + ynew_cen else: + # Convert the PSF center point from pixels to arcseconds using pysiaf + xidl_cen, yidl_cen = aper.sci_to_idl(xsci_cen, ysci_cen) + + # Get 'idl' coords + xidl = xarr * pixelscale + xidl_cen + yidl = yarr * pixelscale + yidl_cen + xv, yv = aper.convert(xidl, yidl, 'idl', to_frame) xmin, xmax = (xv.min(), xv.max()) ymin, ymax = (yv.min(), yv.max()) @@ -171,7 +178,12 @@ def distort_image(hdulist_or_filename, ext=0, to_frame='sci', fill_value=0, ynew += ynew_cen - np.median(ynew) # Convert requested coordinates to 'idl' coordinates - xnew_idl, ynew_idl = aper.convert(xnew, ynew, to_frame, 'idl') + if to_frame=='idl': + xnew_idl, ynew_idl = (xnew, ynew) + else: + # Choose conversion function + func_convert = getattr(aper, f'{to_frame}_to_idl') + xnew_idl, ynew_idl = func_convert(xnew, ynew) # ############################################### # Interpolate using Regular Grid Interpolator From 64348ab1aa9a9f6bf9c6af220cffb6994e62a9ec Mon Sep 17 00:00:00 2001 From: Jarron Leisenring Date: Fri, 12 May 2023 13:59:52 -0700 Subject: [PATCH 2/7] xidl_cen and yidl_cen need to be defined earlier --- webbpsf/distortion.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/webbpsf/distortion.py b/webbpsf/distortion.py index 7b41297c..435e1112 100644 --- a/webbpsf/distortion.py +++ b/webbpsf/distortion.py @@ -121,6 +121,9 @@ def distort_image(hdulist_or_filename, ext=0, to_frame='sci', fill_value=0, xsci_cen = hdu_list[ext].header["DET_X"] # center x location in pixels ('sci') ysci_cen = hdu_list[ext].header["DET_Y"] # center y location in pixels ('sci') + # Convert the PSF center point from pixels to arcseconds using pysiaf + xidl_cen, yidl_cen = aper.sci_to_idl(xsci_cen, ysci_cen) + # ############################################### # Create an array of indices (in pixels) for where the PSF is located on the detector nx_half, ny_half = ( (nx-1)/2., (ny-1)/2. ) @@ -150,9 +153,6 @@ def distort_image(hdulist_or_filename, ext=0, to_frame='sci', fill_value=0, xnew = xarr / osamp_x + xnew_cen ynew = yarr / osamp_y + ynew_cen else: - # Convert the PSF center point from pixels to arcseconds using pysiaf - xidl_cen, yidl_cen = aper.sci_to_idl(xsci_cen, ysci_cen) - # Get 'idl' coords xidl = xarr * pixelscale + xidl_cen yidl = yarr * pixelscale + yidl_cen From 005a60fa1d742e193169a76adeaa4860805ea45e Mon Sep 17 00:00:00 2001 From: Jarron Leisenring Date: Fri, 12 May 2023 14:42:25 -0700 Subject: [PATCH 3/7] fix test_apply_distortion_pixel_scale() failure Need to apply the same pixel sampling along the x- and y- axis, otherwise test_apply_distortion_pixel_scale() fails. May need to update test? --- webbpsf/distortion.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/webbpsf/distortion.py b/webbpsf/distortion.py index 435e1112..a481ba44 100644 --- a/webbpsf/distortion.py +++ b/webbpsf/distortion.py @@ -150,8 +150,9 @@ def distort_image(hdulist_or_filename, ext=0, to_frame='sci', fill_value=0, elif to_frame=='sci': osamp_x = aper.XSciScale / pixelscale osamp_y = aper.YSciScale / pixelscale - xnew = xarr / osamp_x + xnew_cen - ynew = yarr / osamp_y + ynew_cen + osamp = (osamp_x + osamp_y) / 2 + xnew = xarr / osamp + xnew_cen + ynew = yarr / osamp + ynew_cen else: # Get 'idl' coords xidl = xarr * pixelscale + xidl_cen From 31926e759cfc9d9ccc70b47a6d9388a8ddd10a89 Mon Sep 17 00:00:00 2001 From: Jarron Leisenring Date: Tue, 16 May 2023 11:00:54 -0700 Subject: [PATCH 4/7] revert to aper.convert --- webbpsf/distortion.py | 16 +++------------- 1 file changed, 3 insertions(+), 13 deletions(-) diff --git a/webbpsf/distortion.py b/webbpsf/distortion.py index a481ba44..8bdf5730 100644 --- a/webbpsf/distortion.py +++ b/webbpsf/distortion.py @@ -63,7 +63,7 @@ def distort_image(hdulist_or_filename, ext=0, to_frame='sci', fill_value=0, Value used to fill in any blank space by the skewed PSF. Default = 0. If set to None, values outside the domain are extrapolated. to_frame : str - Requested of output coordinate frame. + Requested type of output coordinate frame. * 'tel': arcsecs V2,V3 * 'sci': pixels, in conventional DMS axes orientation @@ -133,12 +133,7 @@ def distort_image(hdulist_or_filename, ext=0, to_frame='sci', fill_value=0, # ############################################### # Create an array of indices (in pixels) that the final data will be interpolated onto - if to_frame=='sci': - xnew_cen, ynew_cen = (xsci_cen, ysci_cen) - else: - # Choose conversion function - func_convert = getattr(aper, f'sci_to_{to_frame}') - xnew_cen, ynew_cen = func_convert(xsci_cen, ysci_cen) + xnew_cen, ynew_cen = aper.convert(xsci_cen, ysci_cen, 'sci', to_frame) # If new x and y values are specified, create a meshgrid if (xnew_coords is not None) and (ynew_coords is not None): @@ -179,12 +174,7 @@ def distort_image(hdulist_or_filename, ext=0, to_frame='sci', fill_value=0, ynew += ynew_cen - np.median(ynew) # Convert requested coordinates to 'idl' coordinates - if to_frame=='idl': - xnew_idl, ynew_idl = (xnew, ynew) - else: - # Choose conversion function - func_convert = getattr(aper, f'{to_frame}_to_idl') - xnew_idl, ynew_idl = func_convert(xnew, ynew) + xnew_idl, ynew_idl = aper.convert(xnew, ynew, to_frame, 'idl') # ############################################### # Interpolate using Regular Grid Interpolator From e5110734cd27688c707abdab4a325fd37fb98b16 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Wed, 17 May 2023 09:39:05 -0400 Subject: [PATCH 5/7] add test to verify fix of distortion with nonstandard pixelscale --- webbpsf/tests/test_distortion.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/webbpsf/tests/test_distortion.py b/webbpsf/tests/test_distortion.py index 3b3d44a8..b0bed7d6 100644 --- a/webbpsf/tests/test_distortion.py +++ b/webbpsf/tests/test_distortion.py @@ -281,3 +281,19 @@ def test_miri_conservation_energy(): assert pytest.approx(psf_sum, 0.005) == psf_cross_sum, "The energy conversation of the PSF before/after the " \ "scattering is added is greater than the tolerance of " \ "0.005" + +def test_distortion_with_custom_pixscale(): + """ Verifies the distortion model works properly even if the pixel scale is changed to + a nonstandard value for the calculation. This tests/verifies the fix in PR 669: + https://github.com/spacetelescope/webbpsf/pull/669 + """ + + miri = webbpsf_core.MIRI() + miri.pixelscale = 0.061 + psf = miri.calc_psf(fov_arcsec=2) + + # A symptom of the prior bug was the total sum of a distorted PSF would be very + # discrepant from the sum of the undistorted PSF. So verif that symptom is not the case: + + assert np.isclose(psf[0].data.sum(), psf[3].data.sum(), rtol=0.001) + assert np.isclose(psf[1].data.sum(), psf[3].data.sum(), rtol=0.001) From 454f1c4c561b84c038285608e391fa17c65dafe6 Mon Sep 17 00:00:00 2001 From: Jarron Leisenring Date: Wed, 26 Jul 2023 15:42:34 -0700 Subject: [PATCH 6/7] updated test_apply_distortion_pixel_scale test_apply_distortion_pixel_scale() now follows what its description Also updated test_apply_distortion_skew() so it actually creates a downsampled version of the data (previously was just creating a copy of the oversampled data). --- webbpsf/tests/test_distortion.py | 43 +++++++++++++++++--------------- 1 file changed, 23 insertions(+), 20 deletions(-) diff --git a/webbpsf/tests/test_distortion.py b/webbpsf/tests/test_distortion.py index b0bed7d6..bc1338de 100644 --- a/webbpsf/tests/test_distortion.py +++ b/webbpsf/tests/test_distortion.py @@ -36,7 +36,7 @@ def test_apply_distortion_skew(): # Rebin data to get 3rd extension fgs.options["output_mode"] = "Both extensions" - fgs.options["detector_oversample"] = 1 + fgs.options["detector_oversample"] = psf[0].header['DET_SAMP'] webbpsf_core.SpaceTelescopeInstrument._calc_psf_format_output(fgs, result=psf_siaf, options=fgs.options) # Test the slope of the rectangle @@ -78,6 +78,15 @@ def test_apply_distortion_pixel_scale(): fgs.options["output_mode"] = "Oversampled image" psf = fgs.calc_psf(add_distortion=False) + # Replace data with a fake image of row values equal to the row number + data = np.zeros_like(psf[0].data) + ny, nx = data.shape + for i in np.arange(ny): + data[i, :] = i + + # Replace the data in the PSF with the fake image + psf[0].data = data + # Set up new extensions (from webbpsf_core.JWInstrument._calc_psf_format_output) n_exts = len(psf) for ext in np.arange(n_exts): @@ -89,9 +98,9 @@ def test_apply_distortion_pixel_scale(): # Run data through the distortion function psf_siaf = distortion.apply_distortion(psf) - # Rebin data to get 3rd extension + # Rebin data to get 3rd extension (DET_DIST) fgs.options["output_mode"] = "Both extensions" - fgs.options["detector_oversample"] = 1 + fgs.options["detector_oversample"] = psf[0].header['DET_SAMP'] webbpsf_core.SpaceTelescopeInstrument._calc_psf_format_output(fgs, result=psf_siaf, options=fgs.options) # Test that the change caused by the pixel distortion is approximately constant along the row @@ -99,31 +108,25 @@ def test_apply_distortion_pixel_scale(): i = 20 ext = 3 - psf_arr = psf_siaf[ext].data - size = psf_siaf[ext].data.shape[0] - inds = np.arange(len(psf_arr)) + # Crop off the edges due to skew / rotation that brings in 0s from beyond edge of detector + psf_arr = psf_siaf[ext].data[5:-5, 5:-5] + ncol = psf_arr.shape[1] + inds = np.arange(ncol) # Model the skew with a basic linear function - yN = psf_arr[i, -1] - y0 = psf_arr[i, 0] - slope = (yN - y0) / len(psf_arr) - linear = (slope * inds) + y0 + slope, intercept = np.polyfit(inds, psf_arr[i, :], 1) + linear = (slope * inds) + intercept # Create a new 1D array that's your 20th row with the linear skew subtracted out - # Add y0 because we want the values compared to 0th value, not subtracted down to 0 (just preference) - final = psf_arr[i, :] - linear + y0 + final = psf_arr[i, :] - linear # Check the difference between adjacent values is the same to 1 decimal place - for i in range(size - 1): - a = final[i] - b = final[i + 1] - - # This is the same as assert round(a - b, 1) == 0 - assert pytest.approx(a, abs=0.1) == b, \ - "FGS PSF does not have expected pixel scale distortion for adjacent pixels" + diff = final[:-1] - final[1:] + assert pytest.approx(diff, abs=0.1) == 0, \ + "FGS PSF does not have expected pixel scale distortion for adjacent pixels" # Check that the difference between the first and last value is also the same to 1 decimal - assert pytest.approx(final[-1], 0.1) == final[0], "FGS PSF does not have expected pixel scale distortion in the " \ + assert pytest.approx(final[-1], abs=0.1) == final[0], "FGS PSF does not have expected pixel scale distortion in the " \ "entire row" From b4fdee2bffd31e05aacbc16406e0f9d99739ead7 Mon Sep 17 00:00:00 2001 From: Jarron Leisenring Date: Wed, 26 Jul 2023 15:43:17 -0700 Subject: [PATCH 7/7] use per axis pixel sampling during distortion --- webbpsf/distortion.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/webbpsf/distortion.py b/webbpsf/distortion.py index 8bdf5730..f73fd262 100644 --- a/webbpsf/distortion.py +++ b/webbpsf/distortion.py @@ -145,9 +145,8 @@ def distort_image(hdulist_or_filename, ext=0, to_frame='sci', fill_value=0, elif to_frame=='sci': osamp_x = aper.XSciScale / pixelscale osamp_y = aper.YSciScale / pixelscale - osamp = (osamp_x + osamp_y) / 2 - xnew = xarr / osamp + xnew_cen - ynew = yarr / osamp + ynew_cen + xnew = xarr / osamp_x + xnew_cen + ynew = yarr / osamp_y + ynew_cen else: # Get 'idl' coords xidl = xarr * pixelscale + xidl_cen