Confirm whether ePSF residual mapping differing from Anderson & King (2000) is correct with regards to how far samplings are drawn per grid point for sigma clipping

[Onoddil]

At present, in EPSFBuilder, when calculating star residual samplings and sigma clipping multiple samplings to find the robust residual for a given grid point -- used to update the interpolation spline points with more accurate PSF samples -- the current implementation uses a nearest neighbour fitting (_py2intround), but the algorithm as detailed in Anderson & King (2000) uses all star samplings within an oversampled grid spacing (i.e., a quarter of a pixel for oversampling=4) to compute the sigma clipped average residual value. Note that nearest neighbour fitting is essentially all data points within half a grid spacing, instead of a whole grid spacing, reducing the number of samples used to compute each grid point residual by a factor 4, and removing some implicit data smoothing without the "rolling average".

There were some other, undocumented changes to EPSFBuilder relative to those in the AK00 paper, so is this another case of a fix made for unknown but valid reasons, or should this be changed to align with the algorithm's intent and swap to pulling samplings from within a whole grid spacing?