Better algorithm to get list of HiPS tiles

[cdeil]

With the one example for make_sky_image that we have in the docs we get this iamge: https://github.com/hipspy/hips/pull/30#issuecomment-312290956

There's two tiles missing at the edge, because the algorithm we have is to do a disk query with center at the center of the image, and radius such that all pixels in the sky image are covered by the disk, and then it can happen that tile centers are outside that disk: http://hips.readthedocs.io/en/latest/api/hips.utils.compute_healpix_pixel_indices.html

The current algorithm is following the description of: http://hips.readthedocs.io/en/latest/drawing_algo.html#naive-algorithm

@tboch - Do you have a suggestion for a better algorithm?

[tboch]

Right, it seems that query_disc method sometimes misses some pixels, especially for small orders (large field of view). In Aladin Lite, the workaround I used was to extend the radius by a factor 1.5: radius *= 1.5 This will return additional unneeded tiles, but they can be filtered out easily using method described in step 5 of http://hips.readthedocs.io/en/latest/drawing_algo.html#naive-algorithm

[cdeil]

@tboch - I was thinking about this a bit: extending the radius by 1.5 can still result in missing tiles if the sky image is small, no?

Here's an alternative algorithm: change compute_healpix_pixel_indices and call http://healpy.readthedocs.io/en/latest/generated/healpy.pixelfunc.ang2pix.html for the coordinates of all the pixels in the image (we already have those as pixel_coords = wcs_geometry.pixel_skycoords in the function currently), followed by making making a unique list.

That should always give the exact list of tiles needed, and also means we don't need the tile list trimming step that's needed with the disk query. When drawing it could also be useful to have this "image" of ipix knowing which tile will / should be used to draw on a given pixel in the sky image.

Possible drawbacks of that method are that it could be slower (ang2pix and computing unique) or memory intensive (extra image array of ipix).

I think my suggestion would be to implement the algorithm with ang2pix first, because it should be just <~ 5 lines of code, and always yield correct results.

@tboch - Thoughts?

[tboch]

Yes indeed, this should always return the exact correct list of needed tiles. Please @adl1995 go ahead with @cdeil 's suggestion, we'll see later how we cope if computation is too slow.

[adl1995]

@cdeil Should I overwrite the current compute_healpix_pixel_indices function, or should there be a separate function for this?

[cdeil]

Should I overwrite the current compute_healpix_pixel_indices function, or should there be a separate function for this?

I think either changing the existing function or adding a new one would be reasonable. My slight preference would maybe be to just replace the existing version to keep our codebase simple. @tboch - Do you have a preference?

Concerning implementation of the new function, I think that maybe calling https://docs.scipy.org/doc/numpy/reference/generated/numpy.bincount.html to compute the list of unique ipix values would be best. It should be very fast because it doesn't have to sort the array.

[tboch]

I agree with Christoph, just replace it @adl1995

[cdeil]

This was implemented by @adl1995 in #50 .