Object centering

[dkirkby]

I have noticed that the default centering of a galaxy is slightly different for images with an even or odd number of pixels along each side. Specifically, galaxies are centered on the central pixel for odd sizes, but centered on the pixel whose bottom-left corner is at the image center for even sizes.

I am not sure if this is a bug or a feature, but it at least needs to be documented. If it is not too hard (or too late), it might instead be worth changing the default for even sizes.

[rmjarvis]

I agree. I've also found that convention a bit confusing at times. The convention comes from Gary's original SBProfile code, and is basically because we use integer bounding boxes, and when we recenter it, we do so by declaring a particular pixel to be (0,0). But this is of course off-center for even-sized images.

It would not be too hard to change this. Especially after my recent updates to how the loop over pixels happens (using fillXValue rather than the old fillXGrid and fillXImage). The definition of what (x,y) each pixel is happens in one place now, so it would be easy to put in a 0.5 offset when appropriate.

The most painful part of the change would be updating all the unit tests that compare to a reference image, since pretty much all of these will now be wrong.

Any objections to this change?

[reikonakajima]

No objections from me.

[rmandelb]

I agree that this convention is confusing. The question in my mind is whether we should change the default, or make a keyword to let users choose?

On Mar 29, 2013, at 11:16 AM, Mike Jarvis notifications@github.com wrote:

I agree. I've also found that convention a bit confusing at times. The convention comes from Gary's original SBProfile code, and is basically because we use integer bounding boxes, and when we recenter it, we do so by declaring a particular pixel to be (0,0). But this is of course off-center for even-sized images.

It would not be too hard to change this. Especially after my recent updates to how the loop over pixels happens (using fillXValue rather than the old fillXGrid and fillXImage). The definition of what (x,y) each pixel is happens in one place now, so it would be easy to put in a 0.5 offset when appropriate.

The most painful part of the change would be updating all the unit tests that compare to a reference image, since pretty much all of these will now be wrong.

Any objections to this change?

— Reply to this email directly or view it on GitHub.

---

Rachel Mandelbaum rmandelb@andrew.cmu.edu http://www.andrew.cmu.edu/~rmandelb/

[rmjarvis]

We already have a way for users to choose. applyShift. So I don't think we need another way to specify that.

And technically, that's a workaround for the current situation if people want the profile to be really centered on an even image. But it isn't obvious a priori whether that should be (+0.5,+0.5) or (-0.5,-0.5) times the pixel_scale. Even I had to look at David's description above to figure out that it is (-0.5,-0.5), since I didn't remember which way it was off.

However, if the default were to use the real center, and they wanted it centered on one of the pixel centers, it would be clear (I think) how to do that by shifting (0.5,0.5) pixels in whatever direction they want.

[dkirkby]

The trouble with applying a shift is that the required amount is different for an even or odd image size (and depends on the pixel scale), e.g.

    shift :
        type: XY
        x: { type: Eval, str: '-0.5*scale if 0==size%2 else 0' }
        y: { type: Eval, str: '-0.5*scale if 0==size%2 else 0' }

[rmandelb]

I meant a keyword that doesn't make them have to investigate what shift is necessary. But if people want to change the default, I'm okay with that. Now is better than later...

On Mar 29, 2013, at 11:47 AM, Mike Jarvis notifications@github.com wrote:

We already have a way for users to choose. applyShift. So I don't think we need another way to specify that.

And technically, that's a workaround for the current situation if people want the profile to be really centered on an even image. But it isn't obvious a priori whether that should be (+0.5,+0.5) or (-0.5,-0.5) times the pixel_scale. Even I had to look at David's description above to figure out that it is (-0.5,-0.5), since I didn't remember which way it was off.

However, if the default were to use the real center, and they wanted it centered on one of the pixel centers, it would be clear (I think) how to do that by shifting (0.5,0.5) pixels in whatever direction they want.

— Reply to this email directly or view it on GitHub.

---

Rachel Mandelbaum rmandelb@andrew.cmu.edu http://www.andrew.cmu.edu/~rmandelb/

[rmjarvis]

The trouble with applying a shift is that the required amount is different for an even or odd image size (and depends on the pixel scale)

Right. It's a bit trickier in the config than in the python. I wasn't really proposing this as a good solution to the current situation. I think it makes more sense to change the default as you suggested.

Probably most people who want the profile to be centered on a pixel center will just use odd-sized images, rather than use the kind of tricky Eval that you wrote here. But if a user absolutely wants an even-sized image and wants the profile centered in the center of a particular pixel, I don't feel too bad about making them jump through that hoop.

[rmandelb]

Whether or not this is a temporary or permanent solution, I would recommend that (at least for now) this example should be given on the config wiki.

On Mar 29, 2013, at 12:45 PM, Mike Jarvis notifications@github.com wrote:

The trouble with applying a shift is that the required amount is different for an even or odd image size (and depends on the pixel scale)

Right. It's a bit trickier in the config than in the python. I wasn't really proposing this as a good solution to the current situation. I think it makes more sense to change the default as you suggested.

Probably most people who want the profile to be centered on a pixel center will just use odd-sized images, rather than use the kind of tricky Eval that you wrote here. But if a user absolutely wants an even-sized image and wants the profile centered in the center of a particular pixel, I don't feel too bad about making them jump through that hoop.

— Reply to this email directly or view it on GitHub.

---

Rachel Mandelbaum rmandelb@andrew.cmu.edu http://www.andrew.cmu.edu/~rmandelb/

[rmandelb]

@rmjarvis - were you volunteering to make this change, or suggesting that someone else do it? @barnabytprowe , any objections to this change in the default centering?

[barnabytprowe]

Ooops thanks for the ping, this fell off my "GalSim response todo list"! I have no objections to changing the default behaviour.

I will say thought that the current behaviour is in my experience very much the most convenient when working with Discrete Fourier transforms (where a certain pixel has to be kx,ky=(0,0)), so I think this also explains why SBProfile works in this way. If we change this at too deep a level I'm vaguely concerned that there is a risk that there are other parts of the C++ layer that may be implicitly assuming the old convention for its indubitable convenience... at best these will no longer match, at worst these we may then inadvertently break.

I'd advocate not making a root and branch overhaul of the entire code to make this change, or at least not right now while we have something of a deadline approaching. Perhaps this is default behaviour that could be implemented at the Python draw() level using a shift, I don't think that would be any more costly computationally and I think it is less fraught with unforseen potential consequences.

[rmjarvis]

I will say thought that the current behaviour is in my experience very much the most convenient when working with Discrete Fourier transforms (where a certain pixel has to be kx,ky=(0,0))

All the Fourier images would still keep this feature. The change is only meant to apply to the real-space images.

[dkirkby]

Another related issue I just noticed is that the the center keyword for a Scattered image seems to have a built in offset of (-1,-1) relative to the bottom left corner, so that you need to add one to x,y offsets to get the behavior I expected. The snippet below demonstrates the problem: a 7-pixel stamp whose flux max is in pixel [3,3](counting from zero) is centered at [4,1], but the max appears at [3,0] instead of [4,1]. The (-1,-1) offset seems to be the same for all combinations of even/odd image and stamp size.

gal :
    type : Gaussian
    half_light_radius : 2

image :
    type: Scattered
    size : 8
    stamp_size : 7
    pixel_scale : 1
    center : { type : XY, x : 4, y : 1 }
    nobjects : 1

[rmjarvis]

I've made good progress on this issue yesterday and today. I took Barney's suggestion to implement the shift in the python layer, which turned out to be pretty easy. I got all the unit tests working with the new centering behavior, and I discovered a bug in the code where shifts would also add a slight spurious shear (BAD!!!), which I fixed.

Another related issue I just noticed is that the the center keyword for a Scattered image seems to have a built in offset of (-1,-1) relative to the bottom left corner, so that you need to add one to x,y offsets to get the behavior I expected.

I'm now addressing this point raised by David. And I realized that this comes down to the standard convention question about 0-based vs 1-based counting. The standard fits numbering scheme is that the lower left corner of the image is (1,1), not (0,0). This is the convention we generally use with our Image class, and it is also the convention I followed when implementing the center positions in the Scattered image type. It is also what ds9 uses when it displays the image.

So David, when I used your config specification above and displayed the image in ds9, I found the maximum pixel located at (4,1), which is the 4th column, bottom row. However, I see why you might think this should be considered (3,0) instead.

So I guess I can claim that it is at least not an inadvertent bug, but rather an intentional feature. So the question (to everyone: @rmandelb @barnabytprowe @reikonakajima @msimet) is, is it a good feature, or should we change it?

Either way, I'm sure David's confusion is indicating that we at least need better documentation about the convention we decide on.

[rmandelb]

Given that it's the standard fits numbering convention, I think that we should use it, but document it better.

One option to consider if we want to be really flexible is to add some kind of keyword like zero_index that means all pixel values (e.g., for the center specification) get interpreted as being zero-indexed?

On another note related to indexing: I noticed that if I make a 10x10 galsim.ImageD, and try to use the setValue method to set some x/y value outside of the range from 1-10, it doesn't complain. That doesn't entirely make sense to me... what do you think?

[rmjarvis]

One option to consider if we want to be really flexible is to add some kind of keyword like zero_index that means all pixel values (e.g., for the center specification) get interpreted as being zero-indexed?

I like this idea of a keyword to switch from 1-based to 0-based positions. That wouldn't be too hard to implement I think. David, does this sound reasonable to you?

Some suggestions:

zero_index = True or False
zero_convention = True or False
index_convention = 0 or 1
index_convention = 'C' or 'Fortran'
index_convention = 'python' or 'FITS'

We could even allow all the index_convention possibilities as synonyms (0 == 'C' == 'python'). I think this would be my preference, but I'm open to further suggestions or statements of preference.

[dkirkby]

Hi Mike - I agree that there is no bug in the implementation of the Scattered image center parameter. My confusion arose because I expected that the same numerical values would have the same effect when used as a galaxy shift or a scattered image center. But these are not generally the same since a shift is a difference between pixel coords and so independent of the indexing convention, so I was implicitly assuming the zero-offset convention.

I would be a bit hesitant to introduce an index_convention parameter since it makes it a bit harder to read scripts (because the meaning of one line can now depend on preceding lines that may be far away in a big script), and I wouldn't use it myself. There's a reason why compilers generally don't provide a flag or #pragma to allow you switch between the two array conventions :-)

Another related confusion I experienced with gal/shift and scattered/center is that they use different units (arcsecs vs pixels) when I would have expected them to be the same. Is there a reason for this? Did you consider treating image separations similar to how angle_value is currently treated, where units (Pixels, ArcSecs) must always be provided explicitly. At first, I found this a bit pedantic for angles but it is now growing on me and certainly makes scripts easier to read.

[rmjarvis]

I would be a bit hesitant to introduce an index_convention parameter since it makes it a bit harder to read scripts (because the meaning of one line can now depend on preceding lines that may be far away in a big script), and I wouldn't use it myself.

I only added it to the config stuff, not to the python layer. I agree that it would be a bad idea to change the convention in the python code, but it seems like it could be useful to be able to declare in the config file what value you want the image origin to have. In the python code, the same behavior is achieved with im.setOrigin(0,0).

Also, I think adding the index_convention to the documentation of config.image will at least serve the purpose of documenting what the default convention is.

Another related confusion I experienced with gal/shift and scattered/center is that they use different units (arcsecs vs pixels) when I would have expected them to be the same.

The logic here is that the galaxy itself has nothing to do with pixels. Everything about it always uses world coordinates (normally arcsec). This includes the size, the shift in applyShift, etc.

However, the image is intrinsically a pixel creature, so positions on it are defined in pixels. So image.shift(dx,dy) takes offsets in pixels, as does image.setCenter(x0,y0). In the config stuff, this includes the image size, the stamp sizes, and the center item for Scattered images.

However, I can see why it can be confusing for the center values, since they are kind of related to both the image and the galaxy. And depending on the use case, one or the other might be more natural. e.g. in your case, it seems that arcsec would be the more natural unit, since you are trying to get two galaxies separated by a particular distance, which makes sense to be in arcsec.

So maybe this would be the better thing to be able to switch between instead of index_convention. We could add a config item for Scattered that would allow a switch between using pixels (1-based convention) to using world coordinates (perhaps relative to the image center). This would probably be more convenient for you. And it would still serve the documentation purpose I mentioned above. What do you think? Any preference about a name?

[dkirkby]

So maybe this would be the better thing to be able to switch between instead of index_convention. We could add a config item for Scattered that would allow a switch between using pixels (1-based convention) to using world coordinates (perhaps relative to the image center). This would probably be more convenient for you. And it would still serve the documentation purpose I mentioned above. What do you think? Any preference about a name?

This sounds reasonable, but then there are other parameters where the same logic could apply. For example, I am using stamp.shift(dx,dy) to add a galaxy stamp to a field in its correct catalog position, which is most naturally expressed in arcsecs, not pixels. I would either leave things as they are now, with clearly documented arcsec/pixel conventions for each parameter (even if this is necessarily a bit arbitrary), or else acknowledge that any parameter that describes a separation on the sky or image can validly be expressed in arscecs or pixels, and introduce a new separation_type similar to angle_type to support this, e.g.

shift: { type: XY, x: { type: ArcSec, distance: 1.2 }, y: { type: ArcSec, distance: 1.3 } }

center: { type: RTheta, r: { type: Pixel, distance: 5.5 }, theta: { type: Deg, theta: 45 } }

[rmjarvis]

I think it might help to think in terms of a more general WCS, which we haven't yet implemented in galsim, but which we plan to. In general arcsec and pixels aren't just a unit convention. They are really different things, which can have a complicated functional to convert between them.

So for a galaxy shift (gal.shift in config, gal.applyShift in python), the galaxy is translated in the sky, and the natural units are arcsec.

However, shifting an image is really a different thing. Shifting some number of pixels can in general be a different number of arcsec depending on which direction you go (i.e. if the WCS includes a distortion). We don't have this in config, but that's what stamp.shift(dx,dy) does in python.

So I think what we want is an alternate way to place the stamp on the larger image, where the position is specified in world coordinates rather than chip coordinates. e.g. sky_center, and maybe change the name of the current center to chip_center to be more explicit. Or perhaps stamp_chip_pos and stamp_sky_pos would be clearer that this is talking about the position of the (center of the) stamp on the larger image, rather than the center of the full image.

[rmjarvis]

So I think what we want is an alternate way to place the stamp on the larger image, where the position is specified in world coordinates rather than chip coordinates. e.g. sky_center, and maybe change the name of the current center to chip_center to be more explicit. Or perhaps stamp_chip_pos and stamp_sky_pos would be clearer that this is talking about the position of the (center of the) stamp on the larger image, rather than the center of the full image.

David,

Any comments on this idea to get rid of center and instead have stamp_sky_pos and stamp_chip_pos that would be in the two coordinate systems? I think this would aptly address the confusion you were having about the units of the center position for Scattered images.

[dkirkby]

Hi Mike - your proposal solves the problem at hand (although you might want to replace "chip" with "image" or "camera") but I still wonder if the problem is more general and may eventually need a more general solution once you support non-trivial WCS transforms.

Imagine that you have encapsulated all knowledge of the relevant transforms into an instance of some WCS class, then how do you update the set of methods that currently accept either arcsecs or pixels? I see two options: either you make a pair of methods (as you are proposing here) where you think it makes sense, or else you encapsulate the world/camera distinction into a generic separation_type class that all these methods accept. I think the second approach is architecturally cleaner and leads to more self-documenting code, but the first method is easier to phase in gradually and more appropriate if the set of functions where both world and camera coords conceivably make sense is small enough (or there are some methods where using either world or camera coords would be conceptually wrong).

[rmjarvis]

There are currently three items that take a position value in the config stuff: image.center, gal.shift, and input.nfw_halo.pos.

I wasn't planning on having the galaxy know about the WCS, so I'd like to have gal.shift only be in sky coordinates.

The NFW halo position is also given in sky coords (relative to the image center), since it seems like that's also something that properly lives in the sky coordinate system. But it wouldn't be conceptually hard to add an option to have it in image coordinates, given where in the code the value is set. (i.e. We will know the WCS at that point.)

So given that there are really only two items where we would need this distinction, I think it would be easier to just have two keywords in each case. But we can keep open the idea of combining them when I do eventually add the WCS stuff.