Add a function to estimate covariance matrix of output image

[astrofrog]

It would be nice to add a function that estimates a covariance matrix of output pixels for any of the reprojection functions. We could make this a 'brute force' way which uses Monte-Carlo simulations to reproject a number of noise images, and empirically measure the covariance between output pixels. This would not be the most efficient way to do it but could work for all of the reprojection algorithms.