Speed up adaptive reprojection

[svank]

This PR applies on top of PR #274 and is the second of three I'll be submitting.

These commits provide very significant speed-ups to the adaptive resampling algorithm. In my test case (which arose naturally from real-world use but was unintentionally a bit pathological [1]), runtime decreased from 972 seconds to 3.3 seconds (comparing well to the 'exact' algorithm, at 15 seconds). The final commit also benefits the interpolation algorithm, reducing its runtime for my test case from 7s to 2.6s.

Here's some quick discussion on each change:

• The first commit here is "Improve handling of zero-weight pixels". The algorithm computes an averaging weight for each input sample. It's frequently the case that the computed weight is zero (since sampling must be done across a region larger that the "non-zero weight" region to be sure to include everything). This commit adds an explicit and cheap check for that condition, and the weighting function returns 0 directly rather than computing out the zero weight the long way. This alone took my test case runtime from 972s to 160s.

• The next commit trades a very small amount of accuracy in computing the Jacobian matrix in order to reduce the amount of coordinate transforms done by two thirds. This is described in the documentation additions in this commit. I've made this speedier option the default, but left the option to use the slower, more accurate calculations. For real-world situations, where the coordinate transformation doesn't change much between neighboring output pixels, the loss in accuracy from the faster version should be negligible. For my PSP/WISPR test case, the output image changes by less than 0.002% per pixel, and the runtime drops from 160s to 145s. For the implemented test case with an AIA image, the largest few per-pixel changes are a few percent where the transformation is most extreme, while 99.6% of non-NaN pixels change by less than 1%, and 97% by less than 0.1%.

• The next commit reduces the number of input pixels sampled by cropping out many zero-weight pixels from the sampling region. The region where the filter function produces non-zero weights is a square in the output plane, so I transform the corners of that square to the input plane with the Jacobian, and sample through a bounding box containing the transformed corners. This is the commit that depends on the previous PR, as it makes use of that Jacobian. This cuts down on the number of zero-weight pixels that are iterated over, and reduces my runtime from 145s to 49s.

• The next commit clamps the sampling region in the input plane to the bounds of the input image. This prevents attempts to sample non-existent input pixels---those attempts didn't cause any problems, but they did waste time. An extra benefit is that in cases where the output image has many pixels that map to regions outside the input image, zero sampling has to be done for those output pixels. This was significant in my test case, bringing runtime from 49s to 7.7s.

• The last commit adds a flag to disable roundtrip checking of the coordinate transforms. By default, for every output pixel coordinate that's transformed to input-image coordinates, that coordinate is immediately transformed back to the output plane to verify that the transform is defined in both directions. In my case, that check isn't necessary, and using this flag to disable the check cuts my runtime from 7.7s to 3.3s. I've set the default to keep doing the check, to ensure that situations that need the check get correct behavior by default. The interpolation algorithm, reproject_interp, does this same check, and so I've added the same flag to that function, where it takes my test case from 7.0s to 2.6s.

The commit speeding up the Jacobian calculation is the only one that changes the outputs of the test cases, since the speedier behavior is on by default. The other commits have no effect on the output images, and do not change the reference images for the tests.

---

[1] I'm using Parker Solar Probe's WISPR instrument, which includes two cameras, providing an inner (more Sun-ward) and outer field of view. Combining images from the two to form a composite field of view requires reprojecting the images into a common projection. My test case is reprojecting the inner FOV image into an array the size of the composite FOV, which leaves about half of the output array empty. (The reprojected outer-FOV image would then be pasted into this space.) The coordinate transform is well-behaved for the portion of the output image that maps to the inner FOV, but for the pixels on the far side of the output array, the coordinate transform produces some extreme warping and causes each output pixel to attempt to sample thousands of (non-existent) pixels from the input image. My runtimes could have been reduced quite a bit by reprojecting the image into an output array sized for just the inner-FOV image, but after this PR, that sort of care isn't required from the user.

[astrofrog]

@svank - could you rebase this, which will remove the first commit since that is now merged?

[svank]

@astrofrog Done

[astrofrog]

Looks great! Just a minor comment relating to docs syntax below

[astrofrog]

Can you write True and False in this paragraph with backticks to get proper syntax highlighting/linking?

[svank]

Good call, I've fixed and force-pushed. (The other PRs need to be rebased now, but I'll wait until after this PR is merged.)

[astrofrog]

As above

[astrofrog]

Thanks!