Add classes to build RGB images based on Lupton et al. (2004)

[parejkoj]

This code implements the Lupton et al. (2004) algorithm for RGB image stacks, to make optimally colored images from 3 separate bands. Not all functionality is implemented (some LSST C++ code for better interpolation over saturated pixels will need to be ported eventually), but the code has been used by a few different people to make color JPEGs and PNGs from 1-3 FITS images.

Yet to do:

• I still need to convert the code to use astropy's new interval.ZScaleInterval in place of ZScaleMapping.
• Many of the tests just run the code and check that a file was generated, they don't actually check the contents of the generated image. We could do so, but there's not an independent set of images to check against, so the best I can think is to generate some images and save those as a baseline to test against. Whether that's necessary for incorporation into astropy, I don't know.
• Convert replaceSaturatedPixels from lsst::afw::display. This should probably be a future feature: it's going to be a bit of work to convert that C++ into relatively fast python.

This code was pulled from lsst.afw.display.

[parejkoj]

I don’t know if astropy prefers that feature commits be rebase flattened, but I can do so if you’d prefer.

[MSeifert04]

Thank you for the contribution! 👍

I didn't review this in detail but just some general comments:

• user facing functions expecting numpy.ndarray or astropy.units.Quantity should convert the input at the beginning with np.asarray(in) or np.asanyarray(in).
• assertions shouldn't be used to control the inputs, just raise custom TypeError or ValueErrors. If I remember correctly some pytest versions have problems with assert outside of test functions.
• importing the __future__ functions absolute_import, division, print_function, unicode_literals at the beginning of the files save some float-casts.
• doubly nested np.where are really hard to understand. Is there a way to make this more readable?
• docstrings should have a Returns section if they return something, even if it's just to state what type and how many returns there are.
• classes should have either a docstring or their __init__ should have a docstring. They are probably combined in the sphinx docs but I think it's better to have one docstring instead of 2.
• Whitespaces around operators like * or / enhance readability.

This looks really nice overall but I think someone else is probably more fit to review the actual code.

[saimn]

Very nice ! I was looking at this algorithm a few days ago, and wanted to test the code you sent to the mailing list a few months ago, so this is perfect 😄
Zscale can indeed be replaced by the Astropy one, but I think you could also use Astropy's stretching functions for the linear/asinh mappings (with ImageNormalize).

[bsipocz]

@MSeifert04 - classes should have either a docstring or their init should have a docstring. They are probably combined in the sphinx docs but I think it's better to have one docstring instead of 2.

It should be fine to have at both places, last time I've checked (when fixed a bug for it) sphinx combined the two.

[cdeil]

@parejkoj - Thank you!
I think this would be appropriate for Astropy core, given that it's a published method and is widely used.

There's already comments by yourself and @MSeifert04 with quite a few suggestions. Maybe implement those first, and then we do more rounds of review?

I only very quickly browsed the code. There's stuff like def writeRGB(fileName, rgbImage) which is just two lines of saving a matplotlib figure. I don't think that's appropriate for Astropy, I see Astropy mostly as a library that should implement the building blocks (functions / classes), not tons of small convenience functions. How to make and save a PNG can be explained in the docs.

In my experience API discussions and review work best when approached from the user perspective. I.e. I'll wait for high-level docs to be here before reviewing this PR.
(clearly the code "works", that's not really in question here, right?)

[bsipocz]

TODO list that I volunteer to finish:

• - clean up tests, transform unittest to pure pytest
• use astropy zscale
• making mpl optional dependency for tests, etc
• remove writeRGB and displayRGB as @cdeil suggested

Changes are in a PR to this PR: parejkoj#2, tests are running on travis on my fork: https://travis-ci.org/bsipocz/astropy/builds/

Edit: tests are passing now

[bsipocz]

Some more commits to address some of @MSeifert04's comments are added to my PR.

[larrybradley]

Can LinearMapping and AsinhMapping be replaced with astropy.visualization.LinearStretch and astropy.visualization.AsinhStretch?

http://docs.astropy.org/en/stable/api/astropy.visualization.LinearStretch.html#astropy.visualization.LinearStretch
http://docs.astropy.org/en/stable/api/astropy.visualization.AsinhStretch.html#astropy.visualization.AsinhStretch?

[parejkoj]

Ok, I've pushed some cleanups, and merged a big set of changes from @bsipocz .

Responding to a few points:

• I'd rather not touch the multiply-wrapped np.where calls. They're not the clearest, but unwrapping them would be a pain and might not be much clearer. If you insist, I can try to unravel them.
• LinearMapping/AsinhMapping: there detailed behavior of the methods I'm adding here are different from those in astropy.stretch, so I'm not sure those can be readily replaced. The ZScale replacement that @bsipocz did made sense because they really were doing the same thing. I suppose we could lift these mapping methods up into astropy.stretch, but there's an API difference we'd have to contend with.
• Does astropy prefer spaces around all operators? LSST has spaces around everything except * and /, to make things look like their order-of-operations would.

I've got an example usage at the top of the file, but I'm not sure what more would be needed for "high-level docs".

[bsipocz]

@parejkoj - High level docs would mean to add it to somewhere here: http://docs.astropy.org/en/latest/visualization/index.html#astropy-visualization, possibly as a sub page, with a paragraph and at least one example.

[parejkoj]

Does anyone have a triple of publicly-shareable >8bit images I could use to make the example?

[bsipocz]

yes, that's were the problems start, because ideally we would want to store any such images in the other data repo, and not here. So for the sake of getting this reviewed and into 1.3, I would suggest to use just random np arrays, or three times one of the test images we already have floating around in the docs. Then over the weekend/next week we can still find some real images to put in there.

[bsipocz]

Also for the tests, we could use either or both pytest-mpl and pytest-arraydiff to test the content, but I would definitely leave it for another follow-up PR if the maintainers are happy with it.

[parejkoj]

I've pushed a short bit of documentation, and linked it from the visualization index.rst file.

Definitely, if the astropy folks think it's ok, l'd rather do the detailed content tests in a subsequent PR.

[bsipocz]

Any rst docs example should be a full script within the file, so you need to add import numpy as np here, too.

[parejkoj]

Done. Can I see what my rendered documentation looks like?

[bsipocz]

python setup.py sphinx_build, then it's got generated in docs/_build/html

[parejkoj]

Looks like I can't do matplotlib from within the docs.

[bsipocz]

You can, but you need to use the .. plot directive. Use the example from e.g. docs/visualization/wcsaxes/index.rst

For the plotting, though, you need to give a full standalone script every time. Adding :include-source: will repeat the code in the docs, too, so you don't need to have it in there twice

[parejkoj]

CircleCI gives the following error for import matplotlib.image:

UNEXPECTED EXCEPTION: ImportError('No module named matplotlib.image',)

Is the version of matplotlib that circleCI uses out of date?

[bsipocz]

I think optional dependencies are not installed on circleci, so you will need to skip all doctest where matplotlib is imported/used. In general use .. doctest-skip:: before the code block.

However, here you don't need the >>> before the lines, remove them from inside the plot directive, just as it's in the other files (e.g. docs/visualization/normalization.rst, as I realize the wcsaxes one is not yet in this branch)

[parejkoj]

Skipping the doctests got it to work. Once we get a better set of images (and a place to put them), we can insert them into the docs.

[parejkoj]

I don't know if travis is waiting because of the OSX problems, or something else, though the tests all pass on my machine. CircleCI passed, so there's now some minimal working documentation.

[parejkoj]

More comments are welcome at this point. I don't know if we'll actually be able to get it in before the freeze or not, given that I don't know the exact procedure there.

[parejkoj]

Just to note: there's one travis failure due to urllib.error.URLError: <urlopen error [Errno -3] Temporary failure in name resolution> in several places during the sphinx run.

[parejkoj]

Looks like all of the travis failures were due to timeouts of one sort or another.

[cdeil]

@parejkoj @bsipocz - Thanks for all the work on this!

IMO it would be better to not rush this into 1.3, but spend the time to polish it.

My main concern is that the end-user API hasn't seen a lot of discussion and most likely isn't final.
So merging it now would mean we're stuck with it forever, or immediately break it in 6 months in the next release.

To help discussion / review, it would be helpful to have a version of the HTML docs online.

I wanted to do this, but for me the Sphinx build is currently broken (see here), not sure yet if it's just for this PR or in master also. @bsipocz - have you seen this error? can you reproduce it with this branch?

So currently this PR adds the following functions / classes to the public API of astropy.visualization:

• makeRGB
• Mapping, LinearMapping, AsinhMapping, AsinhZScaleMapping

My suggestion would be to rename makeRGB to make_lupton_rgb or just lupton_rgb to make it more explicit, and easier to find.

And also in Astropy we follow PEP8 naming to use snake_case for functions and variables. That's something that should be cleaned up in all of the code.

For the "Mapping" classes I'm still a little fuzzy how they relate to the existing "Transform" or "Stretch" classes?
http://astropy.readthedocs.io/en/latest/visualization/index.html#class-inheritance-diagram
Do we need a new separate concept of "Mapping" or should those classes be rewritten as "Transform" or "Stretch" classes, so that the astropy.visualization API stays nice and coherent?

Concerning docs examples, maybe you could contact the authors of the paper and ask if they can still find one of the FITS files for the examples used in the paper?
Then those could be copied on data.astropy.org and used here, IMO that would be nicest if one of the Figures from the paper (e.g. the first) is reproduced.

Concerning license, you put it as BSD-3 here. But it's a copy of older code, right? What was the license there? The easiest would probably to get the author or the original code to leave a comment here that they agree to re-licensing this code as BSD-3 as done here, then it's fine.

There's several small things that should be cleaned up. Just one example: makeRGB has an argument fileName which isn't mentioned in the docstring. IMO it should be removed, because all it does is trigger to execute one line:

matplotlib.image.imsave(fileName, rgb)

which isn't a good way to do it.
If the user wants to write a file, they should just execute it as a separate line after makeRGB.
If they want to use
http://matplotlib.org/api/image_api.html#matplotlib.image.imsave
they can call it directly, and have access to the extra options like "format" or "dpi".

[cdeil]

I restarted the docs build, but there's again timeouts due to SESAME name lookups:
https://travis-ci.org/astropy/astropy/jobs/180892093#L2349

@bsipocz - we had the same problem in astroplan, and @bmorris3 mocked out the remote queries. Too many remote queries make continuous integration fail all the time. How about doing the same for Astropy core? (or just change the examples to not do name lookups? If someone agrees, please make a PR or split it out into a separate issue.

---

Concerning https://gist.github.com/cdeil/5cf16ba5146ed5308c122c82bf9d2020, I think that's actually two issues:

1. there's formatting issues in the RST in this pull request
2. there's a bug in the docs build, that makes it crash without a good error message in that case.

Again, @bsipocz or anyone that has time -- can you please have a look, and if there really is an issue with the docs build erroring out on invalid RST input without good error message, please split it out into a separate issue for Astropy or Astropy helpers.

[bsipocz]

@cdeil - How did you restart the build? It still says to be finished 'about 11 hours ago'. I think only people with commit rights can restart them (even though travis says the restart was successful, actually nothing happens when I try it).

Also there was a problem (probably still there is) with the astropy server at STScI, so many of the time outs comes from there. Changing the examples is something @adrn has to decide about, but I would be surprised if they are causing the issue and the gallery examples are not tested and only run during the docs building.

I'll try to look into the details of your gist tracelog later, but it's strange travis doesn't fails. I guess it may be a problem with python2 where we don't test. Have you tried to build the docs with python3?

[parejkoj]

My suggestion would be to rename makeRGB to make_lupton_rgb or just lupton_rgb to make it more explicit, and easier to find.

Done in latest commit.

And also in Astropy we follow PEP8 naming to use snake_case for functions and variables. That's something that should be cleaned up in all of the code.

I've de-cameled the interface.

Do we need a new separate concept of "Mapping" or should those classes be rewritten as "Transform" or "Stretch" classes, so that the astropy.visualization API stays nice and coherent?

I've looked at the Stretch classes, and they are similar, but have less functionality. I suppose we could lift these into Stretch, but these have some things that are quite specific to this algorithm (e.g. assuming there are 3 images). We could certainly rename *Mapping to something else, if you'd prefer.

Concerning license, you put it as BSD-3 here. But it's a copy of older code, right? What was the license there? The easiest would probably to get the author or the original code to leave a comment here that they agree to re-licensing this code as BSD-3 as done here, then it's fine.

The code this is converted from lived here, and all the LSST code is GPL3. I'd gotten permission from Paul Price and Robert Lutpon (the two primary contributors) to re-license, plus the LSST project managers support this re-licensing, so there's no worries there. If you'd rather something other than BSD3, we can certainly do that.

makeRGB has an argument fileName which isn't mentioned in the docstring. IMO it should be removed.

I've documented it. I disagree about removing it though: I suspect a significant use case of this code will be batch-producing a number of images, and it'll save the user from just re-doing that same matplotlib line. I think it's quite helpful.

Concerning docs examples, maybe you could contact the authors of the paper and ask if they can still find one of the FITS files for the examples used in the paper?

I've asked Mike Evans about whether the SDSS image use policy also applies to images generated from the data. I can easily reproduce e.g. Figure 1 of Lupton et al. from the SDSS data.

Finally, I don't understand why the sphinx build is failing. It fails for me locally when I run python setup.py build_docs, but the error message isn't helpful for me to debug it.

[bsipocz]

Having think about this a little more, I think we should only leave make_lupton_rgb here, so all the rest became quasi private, and thus easier to remove/rework later without a having to support the deprecated API. I think this would address some of the review comments.

[parejkoj]

That sounds good to me. The classes aren't that hard to use, but they're also quite specific use-cases that most people probably won't want.

[cdeil]

Left some more inline comments.

[cdeil]

These *** lines have to have the exact same length as the heading, or Sphinx will complain.
My guess is that if you fix this and potential other RST formatting issues, than python setup.py build_docs will work (because it does in master).

[parejkoj]

Got it, thanks. RST in general is not as picky so long as the header are at least as long as the text.

That said, the docs build is still failing on my local machine, with the same error.

[cdeil]

This will not generate a link in the HTML docs currently.
The only symbol Sphinx knows about is ~astropy.visualization.make_lupton_rgb, so you should link to that.

[parejkoj]

Done, thank you.

[cdeil]

The astropy.visualization.lupton_rgb module isn't part of the official API.
Only the astropy.visualization.make_lupton_rgb function is, so please change to:

from astropy.visualization import make_lupton_rgb

[parejkoj]

Done, thanks.

[cdeil]

Can you please link to the high-level docs section about this function from the docstring?

Here's an example of what I mean: a class Galactocentric which has a separate page in docs and the class docstring and separate page are cross-linked:
http://astropy.readthedocs.io/en/latest/coordinates/galactocentric.html

[parejkoj]

Thanks for giving me an example. Done.

[cdeil]

Remove np.asarray call on this line. Has already been called before.

[parejkoj]

Good catch. I've removed both of those asarrays, as the array sanitizing happens in the Mapping constructor.

[cdeil]

A saturated_border_width option that just raises NotImplementedError isn't useful.
I would suggest to just add the option in a future PR, if / when it's implemented.
(or just implement it here if you want to add it)

[parejkoj]

The code for this currently lives in some LSST C++ that would be a significant effort to convert, particularly because it assumes the inputs are afw.MaskedImages, and thus have the saturated pixels pre-masked.

Should I file a PR for this now?

[cdeil]

Is there any information about the data type and range of values the input image arrays should have and the output image array has?
Is it float or can it also be integers?
Does it have to be range 0 .. 1 or can input / output be any range?

If it's possible to quickly summarise how to use this function correctly in the docstring, so that users don't have to open up the paper or look at the code, I think that would be helpful.

[parejkoj]

Good point. I've updated the docs.

[cdeil]

Suggest to remove the example here.
A better example will be in the high-level docs page.
If you think another example in the "code" file is useful e.g. for interactive help lookup (I don't), then I'd suggest to add it to the function docstring.

[parejkoj]

Nope, that's fine. I put the example there because I didn't know that there would be docs elsewhere.

[cdeil]

I think this note referring to LSST objects "MaskedImages, Images" can just be removed. Here we use Numpy arrays and that's documented above.

[parejkoj]

Good catch! I thought I'd removed all the references to LSST Images. I found one more MaskedImage reference that I removed.

[parejkoj]

Docs still don't build, and I don't understand why.

I'm trying to reproduce Figure 1 of the paper using the SDSS FITS files, but I need to put the gri images in the same frame. I can't see quite how to do that with astropy.wcs: is there a "remap to sky" method?

[cdeil]

is there a "remap to sky" method?

There's https://reproject.readthedocs.io

[eteq]

Hooray! It all passed - merging! Thanks for all your work on this @parejkoj

[astrobot]

@eteq - thanks for merging this! However, I noticed the following issues with this pull request:

• The milestone has not been set
• Changelog entry not present (or pull request number missing) and neither the Affects-dev nor the no-changelog-entry-needed label are set

Would it be possible to fix these? Thanks!

This is an experimental bot being written by @astrofrog - let me know if the message above is incorrect!

[eteq]

Oh @astrobot, you are a lifesaver. Changelog added in c3832fd . But FYI, @parejkoj, in future contributions try to remember to do that ;)

Speaking of, @astrofrog, it might make sense to drop the "experimental" warning (just replace it with a "contact @astrofrog if there are problems")?

[pllim]

I would replace experimental with all seeing and powerful.

[eteq]

(I wish there was a +💯 option as a reaction, as I would give that for #5535 (comment))

[parejkoj]

@eteq Sorry, I totally forgot about the changelog! Thanks for taking care of it. Nobody noticed during review, though... ;-)

I see that the "contributing" docs say to put in the changelog note before the PR is opened, so I'll try to do that in the future.

[eteq]

@parejkoj - Yep that's true, but it's not quite so simple: the catch is that to do that you need to know the PR # - I do a trick where I prepare the commit and then quick check what the last issue # is and do that+1... but that's admittedly awkward. The hope is to develop a better changelog system for this in the near future, though!

[saimn]

Really great to have this in 1.3, thanks @parejkoj !

[larrybradley]

Thanks, @parejkoj!

[bsipocz]

This will need a what's new, too, won't it?

[astrofrog]

Um, I can't remember how I set up @astrobot. It's sentient now.

[cdeil]

@parejkoj - Thank you!

[cdeil]

I tried this out now, using the filename option to write a JPEG.
I think it's flipped in y-direction wrt the input FITS files if I open those in DS9!?

This should be changed, no?

[parejkoj]

Interesting. I think the solution would be to put origin='lower' or origin='upper' in the imsave() call. Could you try that, since you have an example ready to hand?

[cdeil]

@parejkoj - I tried to put origin='lower', but it doesn't work for me.

From a quick look at the MPL function, it seems that there's a bug where the origin option is not respected when I save to PNG!?

For now, I'm using np.flipup on the inputs as a workaround.

I don't have time to make a PR against Astropy or MPL this week.

@parejkoj - Can you check for the image example in the docs and compare to DS9, and make a PR if there is indeed an issue?

def imsave(fname, arr, vmin=None, vmax=None, cmap=None, format=None,
           origin=None, dpi=100):
    """
    Save an array as in image file.

    The output formats available depend on the backend being used.

    Arguments:
      *fname*:
        A string containing a path to a filename, or a Python file-like object.
        If *format* is *None* and *fname* is a string, the output
        format is deduced from the extension of the filename.
      *arr*:
        An MxN (luminance), MxNx3 (RGB) or MxNx4 (RGBA) array.
    Keyword arguments:
      *vmin*/*vmax*: [ None | scalar ]
        *vmin* and *vmax* set the color scaling for the image by fixing the
        values that map to the colormap color limits. If either *vmin*
        or *vmax* is None, that limit is determined from the *arr*
        min/max value.
      *cmap*:
        cmap is a colors.Colormap instance, e.g., cm.jet.
        If None, default to the rc image.cmap value.
      *format*:
        One of the file extensions supported by the active
        backend.  Most backends support png, pdf, ps, eps and svg.
      *origin*
        [ 'upper' | 'lower' ] Indicates where the [0,0] index of
        the array is in the upper left or lower left corner of
        the axes. Defaults to the rc image.origin value.
      *dpi*
        The DPI to store in the metadata of the file.  This does not affect the
        resolution of the output image.
    """
    from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
    from matplotlib.figure import Figure

    # Fast path for saving to PNG
    if (format == 'png' or format is None or
            isinstance(fname, six.string_types) and
            fname.lower().endswith('.png')):
        image = AxesImage(None, cmap=cmap, origin=origin)
        image.set_data(arr)
        image.set_clim(vmin, vmax)
        image.write_png(fname)
    else:
        fig = Figure(dpi=dpi, frameon=False)
        FigureCanvas(fig)
        fig.figimage(arr, cmap=cmap, vmin=vmin, vmax=vmax, origin=origin,
                     resize=True)
        fig.savefig(fname, dpi=dpi, format=format, transparent=True)

[parejkoj]

One question is whether we should assume that ds9 is "correct" here: there are a few possible conventions.

[cdeil]

One question is whether we should assume that ds9 is "correct" here: there are a few possible conventions.

Well, my recommendation was to not re-expose filename on lupton_rgb and call imsave, just to save one line. It's just an unnecessary complication, and an inflexible API if the imsave arguments like origin aren't re-exposed.

I agree there's no right and wrong, just different conventions.
But I would still claim this will "bite" users like it did me today if the orientation stays as-is, and the image is saved flipped up / down wrt. what DS9 shows.

So +1 to remove filename.
And if that doesn't find consensus, then to change to hard-code to origin="lower", and to file a MPL issue why that option doesn't work for PNG output (at least for me the output doesn't change).

[parejkoj]

Setting origin='lower' gives me an output jpegs that looks like ds9's default display.

[astrofrog]

I agree with @cdeil - the output PNG/JPEG files should be flipped vertically by default. FITS files are stored with the bottom pixels first and other formats have the top pixels first, so this is not just about DS9.