Incorrect image size for low surface brightness RealGalaxies

[msimet]

In doing some tests on #389, I've discovered a problem with some low-surface-brightness RealGalaxy images. Basically, the images drawn by RealGalaxy.original_image are much too small compared to the size of the image in the catalog. @rmandelb and @barnabytprowe have suggested this is due to a problem with the calculation of stepK in generating the InterpolatedImage, but in any case it could use some investigation (even if our only solution is "raise a warning when this happens").

For a test case, I will show the galaxy with index 1333 in the latest version of real_galaxy_catalog_23.5_sub.fits, which is the galaxy in HDU 333 of real_galaxy_images_23.5_n2.fits. Here's what the original image looks like, with an extent of 389x389 pixels:

And when I load it into GalSim and draw the original_image, I get this 18x18 image:

[rmandelb]

I think it's worth adding that the incidence of this is of order 0.1% (unless I missed something?).

[msimet]

Well, the incidence of this among things that FindAdaptiveMom() can measure at all is ~0.1% or less, yes. I don't know how large a fraction of the unmeasurable objects it is.

[rmjarvis]

Do we know if any of these pass our cuts on S/N and such for great3? This looks like a pretty marginal detection, so I hope we're not trying to use galaxies this faint.

[rmandelb]

They fall into the I<23.5 sample, but we of course have other cuts. There are a few options for what is going on, including (a) an error in the COSMOS catalog, (b) the object is more extensive than in that image, but there's an error in my routines to mask out nearby objects and am accidentally masking out part of a clumpy galaxy (basically shredding something complex), and I'm sure there are other options I'm missing.

Melanie, can you give me a few more ident values? I can quickly look at the non-masked postage stamps to check for option (b), but it might be better to check a few since we could always have a few effecfts going on here.

I agree we shouldn't be using things that are this marginal, but i think the best way to get rid of them is to figure out why we have them in the first place.

[msimet]

Melanie, can you give me a few more ident values?

1683, 2154, 2478, 2683, 2768. Also, a lot of these (esp 2154) seem to get MUCH noisier when noise_padded, which isn't true of the well-measured galaxies.

[rmandelb]

Just checking, are these ident or index in the catalog?

On Aug 14, 2013, at 6:02 PM, msimet notifications@github.com wrote:

Melanie, can you give me a few more ident values?

1683, 2154, 2478, 2683, 2768. Also, a lot of these (esp 2154) seem to get MUCH noisier when noise_padded, which isn't true of the well-measured galaxies.

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

---

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

[msimet]

Index in the catalog, sorry...

[rmandelb]

Okay, I should really be doing other things right now, but this was kind of intriguing so I did this instead.

The six indices you gave me correspond to the following ident values: (294931, 331253, 134424, 167472, 141727, 80509) I looked into four of them before deciding that it is complicated (no clear/neat explanation for all).

To investigate whether my masking routines (which identify nearby objects and mask them out) are the culprit, I made plots of (a) the original postage stamps Alexie sent me and (b) the masked ones that I use as the basis for these catalogs. In 2 of the four cases, 141727 and 167472, I found no clear explanation for the problem, in the sense that my routine was doing some masking but the resulting object - while kind of extended - should not be problematic for a stepK calculation. For example, here's the image for 141727 (left = original, right = masked, scales are linear with zscale set to the same scale parameters): You can see that above and slightly to the left of the central object there is something very bright that gets masked. And it's questionable as to whether some of the other things that got masked are part of the diffuse central object or truly are separate things. But still, the central object remains decently well-defined so I'm not sure why it would be problematic. The other case I put in this category, 167472, looks even better (i.e. I showed you the worse one).

I just noticed that if I draw ident=141727 out of the real galaxy catalog, draw its original_image into an image, and write it to file, it doesn't look bad - I mean, clearly stepK is somewhat large because it's not leaving much margin around the object, but it's not completely screwed up:

So I'm a little surprised that this one is giving trouble.

Anyway, moving on to cases where I might be culpable... :) 294931: looks like a possible superposition that has something near the edge, and it gets aggressively masked or something. The result is really low surface brightness. So probably I screwed up here?

And I don't think I can take responsibility for 331253, which is pretty screwed up. My code was expecting to get postage stamps that were centroided on decently bright (for COSMOS) objects, but you can see that in this case, the thing at the center is quite faint and low surface brightness rather near an extremely bright object. My masking code does a decent but not great job at masking the extremely bright nearby object, but what's left is the low surface brightness faint thing. Honestly, I don't think it belongs in our bright catalog - I think its photometry in the COSMOS catalog must have been screwed up by the nearby bright thing.

Bottom line: there is not one clear cause for these four cases - there are maybe two surprising ones (not sure they should be hard for GalSim to handle), and two that are messed up but for quite different reasons. Once there was no obvious trend I decided to stop. I think we should try to think of a way to just flag these based on the images themselves. Any thoughts?

[barnabytprowe]

Well Rachel you thoroughly beat me to this (although I will check to see if the case I found on https://github.com/rmandelb/great3-private/pull/50 is a similar beast, I think it's highly likely).

It looks to me like there are ~three things at play here:

1. Masking issues related to nearby objects (leading to super low SNR final images)
2. Spuriously high COSMOS catalog flux values, probably also related to the proximity of nearby objects and difficulties associated with deblending (leading to super low SNR final images).
3. Early termination of the stepK algorithm. It's my gut feeling that this related to the extremely centrally (if at all) peaked + faint-winged profiles of these objects. I would imagine that the algorithm we use simply increases its test radius until it finds no significant excess flux, then simply stops. These sorts of objects would trigger than easily.

The algorithm we use in 3. to get stepK is iterative and depends on its starting conditions, which might be why you saw the behaviour being better in the drawn original_image you tried Rachel. Getting this sort of thing to behave correctly in all cases often depends on the initial guesses, or on tweaks to the termination criteria, and this can be a black hole of work with no guarantee (ever) of total success.

However, it seems to me that in order to result in an 18x18 final image the iterative algorithm must be deciding that our object is only ~2-3 pixels in size, which makes it pretty small compared to the PSF. @rmjarvis : for the future, do you think it might be possible to encode a minimum size into the search algorithm for calculate_stepK, since we know that our objects cannot be smaller than the PSF?

All that said, clearly fixing 3. is not going to be achieved in the next 2-3 days, so we need a way to exclude these objects. The .stepK() method of GSObjects is one obvious way to do that - if that is stupidly large, corresponding to 2 pi / (stupidly small), we could ignore the object and reselect, This would, unfortunately, break the correspondence between our FITS catalogues and the actual galaxies used in GREAT3, but in future we could imagine iteratively making the FITS cats, building RealGalaxies, and remaking with stepK failures removed, to fix this minor annoyance.

Phew, all of that just on item 3!

For 1+2, I don't have much better to suggest than...

• We could try and remove objects with close or close+large neighbours a bit more stringently (I guess this is already done to some extent) to try to reduce the occurrence of these issues.
• We could do the absolute simplest aperture SNR estimate on the drawn image: put a circle of radius X pixels down over the centre, calculate the contained flux, and use our knowledge of the variance per pixel to get the SNR. If this comes up with a value less than ~5-10 (say), we ditch it. (If we wanted to calibrate these SNRs against our other definition it would be straightforward, but to remove these outliers I doubt it would be necessary.)

This is all stuff for the future, and I'm happy to be involved in it. I think Rachel you are right that there is probably not a huge amount we can get from the COSMOS catalogues themselves to indicate where this is an issue a priori, as I think it is profile-dependent and neighbour dependent, so it maybe we need a couple of passes of RealGalaxyCatalog -> RealGalaxy generation to weed out these trouble-makers.

[rmandelb]

Well Rachel you thoroughly beat me to this

It's not totally clear to me that this is the same problem, though that is one possible explanation. I think getting IDs for the problem cases over there is going to be very helpful.

Regarding the exclude if stepK >= (2pi/stupidly small) test, I think this is the best way forward for now. A starting point would be to do this test for the entire catalog and create a flag (to add to the catalog) for those that are problematic. We and others could then use this as a way to filter the sample that we already have. In future versions we could do this iteratively as you suggest to exclude them entirely.

[rmjarvis]

Early termination of the stepK algorithm. It's my gut feeling that this related to the extremely centrally (if at all) peaked + faint-winged profiles of these objects. I would imagine that the algorithm we use simply increases its test radius until it finds no significant excess flux, then simply stops.

I think that's pretty much correct. The calculation starts at the center and expands out square rings of pixels until the total flux reaches some percentage (e.g. 99.5%) of the total flux. The problem in these cases seems to be that the noise in the image means that the total can hit that threshold spuriously early, since there are enough negative noise fluctuations farther out that compensate. If there isn't much real flux farther out, then randomness means that sometimes the flux can be contained inside a radius that is really too small.

for the future, do you think it might be possible to encode a minimum size into the search algorithm for calculate_stepK, since we know that our objects cannot be smaller than the PSF?

Yes. That would be a good improvement. We can allow an argument to set a maximum value of stepk or minimum value of maxk for the two calculate functions when such values are known a priori. This would force the calculation to at least go out to a given radius before checking if it has most of the flux yet.

The algorithm was really designed with OpticalPSF in mind, where there is no noise (so no negative values!), where it was hard to determine a priori how big an image is appropriate. The noisy RealGalaxy images might require a completely different algorithm.

[rmandelb]

So to put these ideas together: we can change the stepK and maxK calculations to allow the optional input of a maximum value of stepk / minimum value of maxk. For RealGalaxy, we can first use the usual algorithm (no optional limits imposed) to calculate stepK, maxK for the HST PSF. Then we use that result as an input for the minimum value of stepk for the RealGalaxy image itself.

[rmjarvis]

This is done on #430. Please see if it helps. (I don't seem to have the 'real_galaxy_catalog_23.5_sub.fits' file in my DropBox folder, so I can't verify myself.)

[rmandelb]

That catalog is just a subset of the one you do have, real_galaxy_catalog_23.5_sub.fits. The "sub" catalog is the first 28k items from the full catalog, in the same order, so the indices and ident values are the same.

I'll also go take a look at #430, but I just wanted to clarify about those datasets.

[msimet]

This is done on #430. Please see if it helps.

It doesn't appear to have made a difference, unfortunately...the images are the same (smaller) size that they were before.

[rmandelb]

That is very strange, because it's definitely the case that for at least one of them, the stepK of the original_image is larger than that of the original_PSF (I confirmed this explicitly using the InterpolatedImages). So this code should enforce the former to be <= the latter.

[rmjarvis]

OK. I'll check it out, now that I know I can just ignore the 'sub' part of the filename.

[rmjarvis]

Actually, it looks like I don't have the images file for this catalog: 'real_galaxy_images_23.5_n2.fits'. Could you share that? Or point me to where I should get it? (I guess I've probably missed some of the banter about these new files...)

[rmandelb]

Leader board website, "data" tab. http://great3.projects.phys.ucl.ac.uk/leaderboard/board if you're not an admin, there is a tarball with the "sub" catalog and all images. If you're an admin, you'll see the full catalog and all images. Actually, the n2 file alone is small enough that I can just stick it in the great3_fits_data dropbox for you - will do that in a moment.

On Aug 15, 2013, at 7:13 PM, Mike Jarvis notifications@github.com wrote:

Actually, it looks like I don't have the images file for this catalog: 'real_galaxy_images_23.5_n2.fits'. Could you share that? Or point me to where I should get it? (I guess I've probably missed some of the banter about these new files...)

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

---

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

[rmjarvis]

Thanks. I got it.

The trickiness with this image (1333, that is) is that the noise is sufficient that the total flux is much less than the central flux. The sum of all the pixel values is 0.89, while the total noise (sqrt(noise_variance * npix)) is 1.59, making the total S/N < 1. However, the central region has a total flux of around 10 with a noise in that region of only 0.5. So near the center, we have S/N ~ 20.

In fact, the accumulated flux steadily decreases from around 10 down to 0.89, so we could be looking at a background calibration problem here. The average background seems to be negative.

Anyway, when SBInterpolatedImage calculates the image flux, it gets 0.89, not 10. So the radius enclosing 99.5% of the flux is pretty much anything. All radii enclose more than the total flux. So when I set the minimum radius to be that of the PSF, it still stops immediately at that point, since that radius also encloses much more than the total flux.

So I'm not sure what the best solution is for this image. I'm inclined to push back and say the background calibration should be fixed. But that might be a bigger project than we want to tackle before Great3. (Certainly not before next week!) So if someone has a suggested solution in the meanwhile, please let me know.

[rmandelb]

If you want a solution for the immediate future, just to keep this from causing holes in GREAT3 images, is there anything wrong with the solution I proposed to identify objects with stepK of galaxy image >~ c * stepK of PSF image (c is some fudge factor to account for noise, e.g., 1.2 or 1.5 or something like that) and just toss them out?

Did you check out whether htis helps with some of the other objects noted in this issue, the ones that are not as obviously weird as this one? (ident = 141727 and 167472)

Background calibration is a longer term project.

[msimet]

For the ones that are not as notably weird, the amount of clipping is also much less--they're 80% of the original size or more. But that was, apparently, enough that they got flagged as being above the 100sigma cut in at least 5 of the tests where they were measurable (as opposed to eg 1333, which was thrown out 20 times based on that cut--I note that this was not quite all of the tests that I ran!).

[rmandelb]

Good point, Melanie. Also I should clarify that my suggestion to check whether this helps with those objects was a suggestion to check whether the new code on #430 helps with the stepK calculation. It wasn't a suggestion to check whether those objects fail the "stepK of galaxy image >~ c * stepK of PSF image" test (I already checked, and they don't fail it - apparently only the ones that look visually awful fail it).

Can you point us to any that fail a lot of the time, like the one with index 1333?

[msimet]

Okay. Indices (not idents, sorry, didn't write out that info) with >=5 failures. (Note that the max number of failures is much larger than 20 now--I must have made that calculation with only some of the fits being done, but since I'm not sure which ones, this is now from all fits including sizes).

98: 1333 80: 4271 76: 4686 72: 3932 54: 2154 52: 5639 50: 1683, 3252 42: 3559 34: 3665 26: 5202 22: 3016, 3755 20: 3329, 5143 16: 3784, 5896 14: 2768, 3706, 5638 12: 2683 10: 112, 2478, 3298 8: 465, 1724, 3106, 3444 6: 85, 310, 1712, 2392, 4332, 4587, 5137

[rmjarvis]

At least some of these seem to have a low nominal S/N ratio. That is, the total flux in the image divided by the total noise, so not trying to shrink the aperture to something more relevant. Here are the snr values for the ones with at least 50 failures from Melanie's list above:

snr =  0.56079090832
snr =  52.1575623322
snr =  63.9439126901
snr =  13.4308546507
snr =  107.521022754
snr =  32.8602442018
snr =  6.90641380765
snr =  25.7757166609

By contrast, the first 20 indices (which I take to be a random sample of objects that don't have problems) have the following snr values:

snr =  167.686095995
snr =  132.144419567
snr =  94.028567234
snr =  101.536624313
snr =  96.8352622825
snr =  107.100143202
snr =  63.9062429332
snr =  124.505033336
snr =  92.8020207715
snr =  145.850335564
snr =  105.43480089
snr =  102.319414439
snr =  458.124993859
snr =  46.4018338045
snr =  268.456546928
snr =  250.562976888
snr =  72.1179733674
snr =  71.5630982699
snr =  50.5220084313
snr =  137.161943336

Only 1/20 less than 50 compared to 5/8. So I think S/N might be a good (one of several probably) cuts to impose to remove the problem images from the sample. I don't think stepK is necessarily a good one, partly because I like the new code that fixes the problem, so the resulting stepk's aren't diagnostic anymore.

Finally, I looked at the 5th item in detail (index 2154) that had a high snr, and I'm not sure why it failed so often for you. Here is a screen shot. On the left is what I get when I draw the RealGalaxy. On the right is the original image in the RealGalaxyCatalog. What nature of failures did it have for you?

[rmandelb]

At least some of these seem to have a low nominal S/N ratio.

Yes. One thing Barney and I were just discussing is that given the low surface brightness and our current S/N definition (intensity-weighted using the ideal profile from Claire's fits), these objects could be 1-sigma above the noise is a large number of pixels, and so they are barely passing a S/N>20 cut for the fits, while still looking like crap in the real image. In the long term, we may wish to filter using a different S/N definition and/or one that uses the real image.

I could see adaptive moment failures happening for the object you pictured because of the blobbiness getting the iterative process confused about the direction of the elliptical gaussian (along the elongated part, or guided by the blobs).

[rmjarvis]

I don't think we want to use objects that have an original S/N near 20. I thought we were planning to add noise to these images to bring them down to S/N = 20 at the lowest. If the original image has S/N ~ 20, then we can't add any more. We are stuck with the noise that is already frozen into the image.

I'd feel more comfortable having a cut of S/N > 40 or 50 for the originals so that most of the noise in the final images would be added by us.

[rmjarvis]

I could see adaptive moment failures happening for the object you pictured because of the blobbiness getting the iterative process confused about the direction of the elliptical gaussian (along the elongated part, or guided by the blobs).

I see. Probably worth leaving this in the sample then. More of a problem with the HSM code than with the image per se.

[rmandelb]

Probably worth leaving this in the sample then.

We can't. We need to be able to measure the shape if we want to do B-mode shape noise.

I don't think we want to use objects that have an original S/N near 20.

I understand this, and I agree. Claire's code fits this with something that is very very extended and low surface brightness. Supposedly if we do the intensity-weighted S/N calculation for that fit, we get a sum(I^2) values that would give us S/N~20 for our added noise level. For such objects that is not very useful, and you're right that if we use the object itself the S/N is too low. So we have to add a cut to our list of cuts to apply, based on the observed image itself. That will have to happen after the initial release.

[rmjarvis]

Probably worth leaving this in the sample then.

We can't. We need to be able to measure the shape if we want to do B-mode shape noise.

Ah. I see. I hadn't appreciated that aspect. OK, as long as these are rare, they probably don't bias the sample too badly. (Although it would be nice to be able to include such tricky objects for people to try to measure shear on!)

[rmandelb]

When I ran the script to estimate the shears, I think they were 1.5%. (We only estimate shears that are better-fit by 2 components; if they are better-fit by a single Sersic, we just use the Sersic galaxy shape directly.) Irregular galaxies are a much higher fraction of the sample than that, so even though the galaxies we're removing are likely irregular galaxies, we're not removing most of the irregulars in the sample with this cut.

[rmjarvis]

Should this issue be closed now that #430 is done? Or is there more to do about these images?

[rmandelb]

Ah, it's really convenient that you just posted here because I had been looking for the issue where you mentioned something about SNR in the original images, and i couldn't find it. But it's this comment on this issue.

So how exactly did you get these SNRs? I do agree we may need to cut on them, but first I have to be able to calculate them in some sensible way.

And to answer your question, I'm not sure there's anything else to do here, but I will leave it to @msimet since she opened the issue.

[rmjarvis]

I don't seem to have saved the script I used to calculate those values. (Probably it was an interactive session.) But I think I just did something like

snr = im.array.sum() / np.sqrt(im.array.variance())

This is very rough of course, since it counts the real flux as part of the noise. So you can certainly do better. I was just trying to get something approximate that might give some indication of what the problem was on these images. Not be super precise.

[msimet]

I haven't seen these failures for a while, so I think this is solved (at least at a level that's below other areas of concern!).

[rmandelb]

Okay, let's close this issue and Mike, I will bring the S/N discussion elsewhere.

[rmandelb]

Closing #454.