GalSim-developers / GalSim

The modular galaxy image simulation toolkit. Documentation:
http://galsim-developers.github.io/GalSim/
Other
223 stars 105 forks source link

`whitenNoise` not working properly with COSMOSCatalog real galaxies #1071

Open xuod opened 4 years ago

xuod commented 4 years ago

Hi, I'm trying to whiten the noise in images from the COSMOS catalog and it seems not to scale properly with the pixel scale and (maybe) the PSF fwhm I'm using. Below is a snippet reproducing what I'm finding.

@rmjarvis thinks that might be due to the noise object not scaling properly, so I tried to dilate the noise object by the ratio of pixel scales, which seems to improve things (not sure how to quantify that).

import galsim

nx = 256
pixel_scale = 0.2

cat = galsim.COSMOSCatalog()
gal = cat.makeGalaxy(index=245, gal_type='real', noise_pad_size=nx*pixel_scale)
PSF = galsim.Kolmogorov(fwhm=0.15)
gal_psf = galsim.Convolve([gal, PSF])

im1 = gal_psf.drawImage(scale=pixel_scale, nx=nx, ny=nx)
im2 = im1.copy()
im1.whitenNoise(gal_psf.noise)

fig, axes = plt.subplots(1, 3, figsize=(12,3))
im = axes[0].imshow(im2.array)
plt.colorbar(im, ax=axes[0])
im = axes[1].imshow(im1.array)
plt.colorbar(im, ax=axes[1])
im = axes[2].imshow(im1.array-im2.array)
plt.colorbar(im, ax=axes[2])

im1 = gal_psf.drawImage(scale=pixel_scale, nx=nx, ny=nx)
im2 = im1.copy()
im1.whitenNoise(gal_psf.noise.dilate(0.03/pixel_scale))

fig, axes = plt.subplots(1, 3, figsize=(12,3))
im = axes[0].imshow(im2.array)
plt.colorbar(im, ax=axes[0])
im = axes[1].imshow(im1.array)
plt.colorbar(im, ax=axes[1])
im = axes[2].imshow(im1.array-im2.array)
plt.colorbar(im, ax=axes[2])

image image

rmjarvis commented 4 years ago

I'll note that the unit tests for whitening are all using pretty low-level interface of drawing blank images with a particular correlated noise, and then whitening it. We don't seem to have any that use the high-level interface of getting a galaxy from COSMOSCatalog (or even directly from a RealGalaxyCatalog), convolving it with a PSF and then checking that it can properly whiten using the noise attribute of the resulting galaxy. We do use that interface in demo11, but there is nothing to check the correctness of the result. So I bet there is something wrong in that process that we never noticed due to lack of sufficient testing.

rmjarvis commented 4 years ago

I think we had visual checks that this used to work in version 1.0 (which was used for the Great3 simulations). My hypothesis is that some of the API changes we made to drawImage, like handing the WCS automatically and letting drawImage convolve by a pixel for you, messed up this functionality. Probably we need to port over some of those steps into the whiten function.

rmandelb commented 4 years ago

If nobody tackles this before ~2 weeks from now, I could play around with it then.

rmjarvis commented 4 years ago

I started looking into this, but the more I do, the more I think I don't actually understand how any of this CorrelatedNoise stuff works. :( I started with just an InterpolatedImage whose image is just noise, so I could track what the noise is when I draw it. I think it all works ok when the pixel scale of the image I'm drawing onto matches the pixel scale of the underlying InterpolatedImage, and if I draw with method='no_pixel'. The second part makes sense. We'll at least need to make it easier to handle the different draw methods that went into making the image when someone wants to whiten it. Especially the default method='auto' that convolves by a pixel for you. We'll need to do that to the noise object too. Probably by adding a method parameter to the whiten function. However, even without convolution with the pixel (or deconvolution or any other transformations), I can't figure out the right thing to do when the image has a different pixel scale than the InterpolatedImage. The resulting variance seems to scale as pixel_scale**3, which I don't understand yet. And whiten doesn't do the right thing in terms of figuring out what the existing noise is in the image, so it's not just me that doesn't understand -- it seems that the code doesn't get it right either.
So I think there is probably more we need to do with respect to keeping track of the wcs of the noise object and tracking what that means for images with different pixel scales.

rmjarvis commented 4 years ago

I pushed the test that I wrote, which is currently failing, in case you feel like taking a look. Maybe I'm just doing something stupid here. (I don't expect you too until next week ofc.)

rmjarvis commented 3 years ago

@rmandelb I know you're probably still super busy, but in case you feel like a change of pace over the holidays, this issue is still extant. :)