Improve stability of SIP fitting. Fix constant term - CRPIX - in SIP #427
Conversation
Codecov ReportBase: 86.84% // Head: 86.93% // Increases project coverage by
Additional details and impacted files@@ Coverage Diff @@
## master #427 +/- ##
==========================================
+ Coverage 86.84% 86.93% +0.09%
==========================================
Files 24 25 +1
Lines 3769 3888 +119
==========================================
+ Hits 3273 3380 +107
- Misses 496 508 +12
Flags with carried forward coverage won't be shown. Click here to find out more.
Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here. ☔ View full report at Codecov. |
|
Here I will attach an ASDF file with two test WCS: NIRCAM and MIRI. After unzipping, do the following: import asdf
with asdf.open('test_wcs.asdf') as af:
w_nircam = af["nircam_wcs"]
w_miri = af["miri_wcs"]
hdr = w_nircam.to_fits_sip(max_pix_error=0.001, degree=None, max_inv_pix_error=None, npoints=128, crpix=[500, 500]) |
mcara
force-pushed
the
sip-use-lu-solve-fix-crpix
branch
3 times, most recently
from
89baa8c
to
5c37e28
Compare
gwcs/wcs.py
Outdated
| # 0 < i + j <= degree. | ||
| pseudo_vander = np.empty((x.size, nterms), dtype=float) | ||
|
|
||
| def crd_pwr(p, q): |
There was a problem hiding this comment.
Is all this conditional logic really necessary? Why not just use the general expression?
| a[i, j] = np.sum(coord_pq, dtype=flt_type) | ||
| a[j, i] = a[i, j] | ||
|
|
||
| with warnings.catch_warnings(record=True): |
There was a problem hiding this comment.
Should we add this fitter to astropy modeling?
gwcs/wcs.py
Outdated
| @@ -2752,67 +2756,197 @@ def _calc_approx_inv(self, max_inv_pix_error=5, inv_degree=None, npoints=16): | |||
| (Shift(crpix[0]) & Shift(crpix[1]))) | |||
|
|
|||
|
|
|||
| def _fit_2D_poly(ntransform, npoints, degree, max_error, plate_scale, | |||
| def _poly_fit_lu(xin, yin, xout, yout, degree): | |||
There was a problem hiding this comment.
This deserves a significant docstring.
mcara
force-pushed
the
sip-use-lu-solve-fix-crpix
branch
from
34dbdc8
to
bfa163c
Compare
| # normal system proved to be significantly more accurate, efficient, | ||
| # and stable than SVD. | ||
| # | ||
| # coord_pow - a dictionary used to store powers of coordinate arrays |
There was a problem hiding this comment.
I guess that this is not a significant memory issue for large arrays, but what is the highest power we envision?
There was a problem hiding this comment.
doesn't matter. memory caching cannot store/use more than 2x memory used for the pseudo-vandermonde matrix. Just an estimate, for 9th (highest) order SIP it will store 54*npoints**2 values in addition to the same number computed anyway. For 1024 points, it would require about 50-100MB extra depending on longdouble size.
There was a problem hiding this comment.
For an expected typical usage with npoints=256 and longdouble being 128 bits and SIP degree of 5, caching would require 20MB
There was a problem hiding this comment.
Right. Somehow I was thinking it was for every pixel.
mcara
force-pushed
the
sip-use-lu-solve-fix-crpix
branch
from
cd12f32
to
b76d061
Compare
mcara
force-pushed
the
sip-use-lu-solve-fix-crpix
branch
from
b76d061
to
a2960c3
Compare
mcara
force-pushed
the
sip-use-lu-solve-fix-crpix
branch
from
a2960c3
to
9f57c17
Compare
|
|
Closes #426
Closes #425
This PR addresses most of the issues in #425 and replaces #426. It has simplified the logic on when to raise errors or warnings, fixes (constrains) the first (constant) term in the SIP polynomials being fitted to distortions, and most importantly, it uses a different method for solving equations for the polynomial fit that in my testing seem to significantly outperform
astropy.modeling.fitting.LinearLSQFitterin accuracy and stability(and possibly timing - not tested)and time/performance.