Indexing afw.image objects

Data Management
, ,
#1

A recent discussion with Astropy developers on improving interoperability between LSST’s image classes and Astropy’s NDData classes led to a set of questions about how our classes behave (or should behave) that @parejkoj captured in the #dm-astropy Slack channel.

I’ve tried to answer those below. Opinions about how they ought to behave are my own, though I believe in many cases they’re shared by others (or at least there’s broad agreement that the current interface is not ideal). I certainly hope anyone who disagrees will speak up on this thread, and I hope this discussion might spawn some RFCs to fix some of the problems; I think any meaningful fixes will be backwards incompatible at a level that will have to involve a frightening number of changes throughout the stack, but I don’t like the alternative of living with these problems either.

How do users access pixels?

Many of our algorithms only access pixels in C++, using standard C++ iterator patterns or a small extension to them.

Those that do access pixels in Python typically do it through NumPy views: calling getArray() (or, after DM-10541, accessing the .array attribute) returns a writeable, (height, width) array view into an image. These views ignore xy0, which I think is desirable; our Image class is conceptually just a 2-d NumPy array (guaranteed to have contiguous values within a row) combined with the xy0 offset, so accessing the array view naturally discards that offset.

Images also support direct pixel access via image.get(x, y) and image.set(x, y, value) (which we could and perhaps should expose as image[y, x]), but using those is discouraged for the same performance reasons that discourage Python programmers from writing for loops over array elements in general. Unfortunately, our get and set methods currently do not pay attention to xy0. They absolutely should.

Our image classes also support in-place binary operators (+=, -=, etc.) with other images, scalar values, and instances of our Function class. We’ve deliberately not implemented the versions of these operators that return a new object to force users to think before blindly allocating memory for a temporary. I’m ambivalent as to whether that was a good idea; it has no doubt prevented some bad code, but I imagine it’s also led to some confusion and frustration.

What happens if you try to add two MaskedImages with different xy0

Because we only support in-place addition, the symmetry between the lhs (input/output) and rhs (input only) is broken. That sidesteps some issues that would otherwise appear.

What does happen: the two images are added, and the xy0 of the image being updated is not changed.

What should happen: If the bounding box (including xy0) of the rhs image is wholly contained by the bounding box of the lhs image, add the images in the overlap region, modifying only that part of lhs. If this is not true, raise an exception.

If we were to support regular (non in-place) binary operators, I think it’d probably be best if those required the bounding boxes (including xy0) to be identical, and raised an exception if they were not.

Can I slice my_image[5:10, 17:42] in an…

…Exposure?

Yes, and the result will be another Exposure that is a view into the original. The ExposureInfo is transferred without modification (more on that below).

What does happen: the slice indices ignore my_image.xy0, they are interpreted as x, y, and the returned view has x0 = 5, y0 = 17.

What should happen: the slice indices account for my_image.xy0, they are interpreted as y, x, and the returned view has x0 = (my_image.x0 + 17), y0 = (my_image.y0 + 5). This is what happens if you pass a bbox argument to the Exposure constructor, and IMO that’s the right way to get a subimage at present; slice indexing is too broken to be safe.

…ExposureInfo?

This is not possible at present, and I’m not sure it needs to be.

The only reason I think it would need to be sliceable is if it contains components that both depend on the coordinate system of the image (e.g. PSF models and WCS transforms) and have a bounding box that needs to modified to match the resulting subimage’s bounding box (as opposed to just indicating a possibly-larger region of validity for the component). Depending on how we want to interpret the bounding boxes on those components, then, I think you could argue that there’s no need to ever slice an ExposureInfo. And it’s desirable not to, because it avoids a lot of copying and/or proxy objects that just impose a new bounding box on their referent.

And that’s one big advantage of having xy0 in Exposure: if you do slice it (with a bbox argument to the constructor, not through the broken slice syntax; see above), you don’t have to copy or modify the components, because the subimage and the original Exposure have the same coordinate system. The big advantage, of course, is that this means subimages also remain valid for objects (like source catalogs, model-fitting parameters, or Footprints) that are tied to the original image’s coordinate system but aren’t a part of the Exposure object at all.

…MaskedImage?

Yes, with the same qualifications about the actual slice syntax being broken as Exposure. The result is another MaskedImage.

…MaskedImage.data?

I’m not sure what this means; MaskedImage doesn’t have a .data, and its arrays aren’t allocated as a single block, so there’s nothing I could even imagine filling that role.

If you mean e.g. MaskedImage.image.array, that will do regular NumPy slicing on the array view, which of course returns a regular NumPy array and deliberately does not take xy0 into account.

What is the result of slicing those objects?

(included in answers above)

How does a user get world coordinates from pixel coordinates? How much linkage needs to be present for that to work, i.e. which components of Exposure are necessary to calculate a WCS?

Because we’ve defined our pixel coordinates such that the lower-left corner of the image is xy0, not (0, 0), all you need to go from pixel coordinates to sky coordinates is exposure.getWcs().pixelToSky(position). All you really need is the Wcs object, and you can use the same Wcs object to perform coordinate transformations when working with a subimage of the original Exposure the Wcs was associated with (see above discussion on slicing ExposureInfo). It does not have to be attached to the Exposure at all to work.

If you do care about array coordinates (i.e. indexing and slicing Image.array), you’d need to utilize both the Wcs object and the xy0 offset, of course. But the whole philosophy of this approach is that you don’t want to care about array coordinates: you put the offset between pixel and array coordinates in your image objects, and that saves you from putting (and constantly updating) an offset in all of your other objects.

1 Like
#2

Why should images be indexed in the order y, x? That seems like a recipe for bugs.

#3

Very well said.

It is clear that if AstroPy is to adopt any of this then they will want to support binary operators that return new image-like objects. I agree that for such operations it makes sense to insist that xy0 and dimensions both match. From that perspective I would argue that we should do the same for our in-place operations. it is a bit clumsier than your suggestion of only setting pixels in the matching subregion, but avoids unwanted surprises (only setting a subregion when one was expecting to set all pixels).

I wonder if it makes sense to try to add a bounding-box like object to AstroPy? We find it very convenient, e.g. for specifying a subregion or comparing two subregions (a single equality operation).

As to bbox in ExposureInfo…a timely discussion, since adding a bbox to ExposureInfo is an outstanding request. I had hoped that ExposureInfo would contain the bbox of the subregion to which it is attached. That would allow such operations as searching a reference catalog for sources that overlap a subimage (such as a coadd patch) without unpersisting the pixels; instead one could unpersist just the ExposureInfo. You make a good argument for keeping the parent bbox, but either we need the ExposureInfo to contain both, or we need an efficient way to unpersist the bbox of a subimage. I suspect the latter would suffice, leaving us with only the parent bbox in the ExposureInfo.

#4

There’s no good choice for this convention; people say both “x, y” and “row, column”, which are inconsistent. The reason I think “y, x” is the least bad answer in this case is because it’s very much the NumPy convention.

For me, adding a subimage is such a useful operation that I’d like to make it as simple and readable as possible, even at the expense of an asymmetry with the non-in-place operators.

I certainly think it’s useful. Before recommending ours to Astropy I’d want to think about making it immutable and/or adopting the sphgeom naming conventions for mutators.

I’m not certain what’s best here either. I actually do think that if we do add a bbox to ExposureInfo, it should probably be the exact bbox of the Exposure and hence have to be copied and updated when subimages are created - what I really want to avoid is updating bboxes in all of the things ExposureInfo holds.

But I do think we need to separate the need for unpersisting a bbox independently of loading an Exposure from the need to add bbox to ExposureInfo. The former is obviously necessary, but I’m still not sure what we really get from adding the bbox to ExposureInfo. There is nothing that says we need to add something to ExposureInfo to load it independently of Exposure; on the contrary, if we ever support reading ExposureInfo itself without reading the full Exposure, I imagine we’d do it by loading all of the things in an ExposureInfo individually and attaching them to a newly-constructed ExposureInfo.

#5

Thinking about this a bit more, it may not be trivial to include a full parent bbox in ExposureInfo. The problem is that we save xy0 in our images, but we do not save the dimensions of the parent image. As long as we build the ExposureInfo for the parent image to start with then we’ll be OK, but I doubt we always work that way – certainly it’s quite easy to make an Exposure with a WCS and a specified xy0.

#6

I think the more important issue is that we don’t have a use case for putting a parent bbox in ExposureInfo. I could imagine having the exact bbox there being useful, just because it’s one more property of the image that you could pass around by passing around an ExposureInfo. The question is whether that bundling-for-convenience is worth the pain of having to slice an ExposureInfo whenever you slice its Exposure.

The only parent bboxes I envision us having in ExposureInfo are stored indirectly - they’re the regions of validity of something like a PhotoCalib, which would only be coincidentally the same as the bbox of the parent image because that might be the region over which the PhotoCalib is considered valid. But that bbox could also cover a whole focal plane if we have a single PhotoCalib whose region of validity extends that far, and then we wouldn’t have any image with the same bbox.

#7

Thanks for the detailed answers – sorry about the delay in responding. One high-level take-away from this seems to be that some of the details of Exposure and MaskedImage are still being worked out (at least with regards to slicing bounding boxes).

In astropy there are (broadly speaking) two classes in the nddata part of astropy. One is NDData, which is really just intended as a simple container to provide a very basic, uniform interface to gridded data (the interface itself is defined in NDDataBase). NDData should be fairly easy to integrate with LSST.

The closest equivalent to Exposure in nddata is NDDataArray, which adds slicing, some arithmetic methods, and a StdDevUncertainty which is propagated correctly by the arithmetic methods. This object is intended to hold image data (i.e. a two-dimensional array of pixel values).

There are some important differences in implementation between NDDataArray and Exposure; perhaps the largest is that NDDataArray doesn’t keep track of an origin like xy0, so when an image is sliced, the WCS is updated to reflect the change in origin. I also don’t believe the arithmetic operators are overridden (neither the binary + nor the in-place +=); instead you add im1 and im2 with im1.add(im2). The motivation for that decision was the same as yours, as I understand it, in not arithmetic operations on Exposures: it is not clear how to add metadata. Sticking with a method allows the user to suggest strategies for merging or replacing the metadata.

One thing I like about the LSST implementation is that you’ve bundled the data, mask and uncertainty together – that is not the case in NDDataArray, and it makes the class much more complicated.

I’d like to move the discussion over to the astropy proposal for enhancement (APE) related to this: https://github.com/astropy/astropy-APEs/pull/14

My next step there will be to put together a table summarizing the differences between LSST images and the current nddata classes, and will try to update the APE on Monday.

It would be very nice to have the MaskedImage object available in astropy – we briefly discussed this option at pyastro17 and my impression is that it would be hard to do in the short run. If it happens down the road, though, we could get rid of a few hundred lines of code, at least, in astropy.

#8

@jbosch – one follow-up question: is array a property on MaskedImage or on another object? I took a quick look at the JIRA issue you linked to and it mentions Image and Mask , so I’m guessing it is something like MaskedImage.image.array and MaskedImage.mask.array.

#9

Yes, that’s correct: there’s only a .array on Image and Mask. Now that you mention it, I believe we have a .getArrays() on MaskedImage that returns a tuple of the three arrays, but that’s so rarely used that I actually forgot about it when adding the properties, and I imagine no one would notice if we got rid of it.

#10

I think you already know this, but just to be sure: MaskedImage is a C++ object with a pybind11 interface. I think we’d be happy to try to move it into astropy, but would this be the first piece of C++ in the guts of astropy?

#11

I think you already know this, but just to be sure: MaskedImage is a C++ object with a pybind11 interface. I think we’d be happy to try to move it into astropy, but would this be the first piece of C++ in the guts of astropy?

It might be, though I wouldn’t see that as an impediment necessarily. It might make more sense for MakedImage to be in a separate package on which astropy depends, but which is smaller than afw. My memory from the discussions a few weeks ago is that the problem with adding afw as a dependency was its size (don’t remember who raised that concern).

I’m working on a large-ish PR to the original NDData/LSST APE right now, hopefully ready to share in ~1-2 hours.

#12

That could be a very nice approach. It wouldn’t be a very large package (MaskedImage doesn’t have that many dependencies)

Size, complexity, and it containing a bunch of things that are either already in astropy, or that you probably don’t want (at least, you don’t want this implementation of them).

#13

Please take a look at the NDData/LSST APE and comment when you have a chance. It has been updated based on the discussions at pyastro17 and the responses to the questions raised there.

Looking forward to your comments!