Improve draw algorithm to handle tile edges in a good way

[cdeil]

As discussed previously, our current method of combining projected tiles isn't good. We simply make full sky images for each projected tile and sum them up. This often leads to "bright lines" at tile edges, because pixels there are drawn and summed twice, i.e. the intensity is a factor 2 too high.

Example: https://gist.github.com/cdeil/796875efe4a3620987b9f62eb278cfbc
Output:

---

Using the RGB example from the high-level docs, one can see that it can also happen that we get "dark lines" or holes between tiles when drawing with the current method:

---

And the current SimpleTilePainter is super slow, presumably because it applies the projective transform for every tile for the whole sky image (even though the tile only contributes to a small subset of pixels in the output image).

---

So we need to implement a better method in the SimpleTilePainter, to avoid these artifacts and get reasonable performance. I have some ideas, @adl1995 - we can discuss in the call this afternoon.

[cdeil]

I'm moving this to the 0.2 milestone. I didn't have time to implement this before going on vacation and want to release v0.1 now.

[cdeil]

Coming back to this question how to draw tiles.

The core of the current implementation is here:

hips/hips/draw/paint.py

Lines 199 to 203 in 8d6078f

	for tile in tiles:
	tile_image = self.warp_image(tile)
	# TODO: put better algorithm here instead of summing pixels
	# this can lead to pixels that are painted twice and become to bright
	image += tile_image

So the algorithm is to project each tile into a full sky image, and to sum up those projected sky images.
That is very inefficient in terms of storage space and CPU time.
And using + is also incorrect, it leads to the "bright edges" where pixel values are double of what they should be, as shown in the example above.

I think a much better algorithm would be to first compute a "plan" for a given sky image, where for every pixel in the sky image a decision is made which tile to use to compute it's value. This will give you a label image which you can integer-encode e.g. like this:

0 0 1 1
0 1 1 1
0 0 1 2

where 0 means use tile 0 and 1 means use tile 1 and 2 means use tile 2.

Then do something like this:

for tile_idx in range(0, label_image.max()):
    pix_idx = np.nonzero(image == tile_idx)
    # project and fill just these pixels for that tile

This should be much faster, memory efficient, and I think avoid the bright edge issue.

[adl1995]

@cdeil - I'm trying to get back to this issue, and would like to work on it during the weekend. Regarding the computation of tile plan:

I think a much better algorithm would be to first compute a "plan" for a given sky image, where for every pixel in the sky image a decision is made which tile to use to compute it's value.

So, currently we apply projective transformation to project the tile corners to the sky image. In this case, do we just add a preprocessing step which creates the label image? Could you please elaborate how the label image should be constructed? Do we find the area spanned by a particular tile (using tile corners) and assign a unique label to that part of the image?

[keflavich]

It's been a long time, but this is still an open issue - anyone considering tackling it?

[adl1995]

Hi @keflavich, I'm afraid none of the original developers are actively maintaining this project any longer. If you want send us an PR, we will be happy to review it (or add you to the list of maintainers).