LSSTDESC / descwl-shear-sims

simple simulations for testing weak lensing shear measurement
BSD 3-Clause "New" or "Revised" License
10 stars 9 forks source link

Add optimized drawing #117

Closed esheldon closed 3 years ago

esheldon commented 4 years ago

drawing in small postage stamps and adding to the image is much faster but we need to ensure we can still recover the shear, in other words there are no simulation artifacts associated with this

esheldon commented 4 years ago

There is a nice way to do this with newer galsim

image = ImageF(nx,ny, wcs=complex_wcs)
shear = Shear(g1=g1,g2=g2)
for obj, world_pos in zip(obj_list, world_pos_list):
    image_pos = (image.wcs.toImage(world_pos) - image.true_center).shear(shear) + image.true_center
    stamp = obj.shear(shear).drawImage(wcs=image.wcs, center=image_pos)
    b = stamp.bounds & image.bounds
    image[b] += stamp[b]

we currently calculated the offsets rather than getting world pos and then converting to offsets. Not sure which I would prefer in the long run, but the later is certainly easier if we have some full sky catalog.

Note I'm not sure if this will always work as well as the "Add all objects then draw" version. The latter adds all the objects together so it has a good idea of how bright each part of the image is. e.g. the value in a pixel is the sum over all objects and can be a lot higher than any individual object. Without that information the drawing might not set appropriate thresholds and, for example, make the stamp too small which could cause artifacts that will fool the shear code

esheldon commented 4 years ago

Also note that the shearing of the image position is not actually implemented yet :)

rmjarvis commented 3 years ago

FYI, shearing positions is implemented now in GalSim v2.3. Syntax is:

sheared_pos = pos.shear(shear)

cf. http://galsim-developers.github.io/GalSim/_build/html/pos.html

esheldon commented 3 years ago

How do I achieve this if I have a set of offsets in arcseconds rather than a set of world positions? Do I need to do the calculations by hand or is there a way to do this using galsim primitives?

beckermr commented 3 years ago

You want pos to be the distance to the center of the uv plane in arcseconds.

esheldon commented 3 years ago

my shifts are in uv.

To use mikes' formula I need to transform uv to pixels. Just wondering if there was built in way to do it or do I need to type out the transformation

esheldon commented 3 years ago

I think it would be this in uv

sheared_uv_shift = uv_shift.shear(shear)
pixel_shift = some_transformation(sheared_uv_shift)
image_pos = image_center + pixel_shift
beckermr commented 3 years ago

ahhhhh - this likely depends on how the coadd wcs is represented.

esheldon commented 3 years ago

This is for single epoch; but in any case the wcs is TanWCS

beckermr commented 3 years ago

Right. The last time I tried to solve this problem I introduced that bug. :(

esheldon commented 3 years ago

Do I need to get the jacobian and run uvToxy ?

esheldon commented 3 years ago

(I give that as an example of something that doesn't seem right)

beckermr commented 3 years ago

The celestial wcs has uv handling stuff internally iirc but idk what precisely to call

esheldon commented 3 years ago

I'm stumped then. Maybe @rmjarvis has a suggestion.

Another way forward would be to generate ra, dec positions rather than offsets so we can use mike's formulae

beckermr commented 3 years ago

I am 99% sure there are private methods that have what need. I think we should make a PR to galsim to expose them.

beckermr commented 3 years ago

There is an _uv method in the celestial wcs classes. This is the coordinates right before deprojecting to the sphere. I think this is what we need (or it's inverse)

esheldon commented 3 years ago

That is relative to crpix, but the center of shearing can be different

beckermr commented 3 years ago

The u,v plane is flat, so we shift centers there I think

esheldon commented 3 years ago

That's what I mean, we would need to copy that code out and have our own function for this which I don't want to do

beckermr commented 3 years ago

I think we should make a PR to galsim to expose this publicly. Then we can write code against that new API.

esheldon commented 3 years ago

Separate issue, apparently we need to specify the image size for the stamp

stamp = obj.shear(shear).drawImage(wcs=image.wcs, center=image_pos)
galsim.errors.GalSimIncompatibleValuesError: Cannot provide non-local wcs with automatically sized image Values {'wcs': galsim.GSFitsWCS(_data = ['TAN', array([100.0, 100.0]), array([[-7.305555555555556e-05, -0.0], [0.0, 7.305555555555556e-05]]), coord.CelestialCoord(coord.Angle(3.490658503988659, coord.radians), coord.Angle(0.0, coord.radians)), None, None, None]), 'image': None, 'bounds': galsim.BoundsI()}
esheldon commented 3 years ago

On Fri, Jul 30, 2021 at 2:43 PM Matthew R. Becker @.***> wrote:

I think we should make a PR to galsim to expose this publicly. Then we can write code against that new API.

Do you mean with an api to re-center?

-- Erin Scott Sheldon Brookhaven National Laboratory erin dot sheldon at gmail dot com

esheldon commented 3 years ago

If I do send dims for the stamp, I get a different error

galsim.errors.GalSimBoundsError: Attempt to access subImage not (fully) in image galsim.BoundsI() not in galsim.BoundsI(1,200,1,200)

I thought the b = stamp.bounds & image.bounds was supposed to deal with this

esheldon commented 3 years ago

I can avoid this error if I check that b.area() > 0

beckermr commented 3 years ago

I mean if we have coordinates in uv, we can compute the uv center like and shear about that center.

esheldon commented 3 years ago

I don't think _uv does that so we would have to add a keyword to use a different center

beckermr commented 3 years ago

we can call _uv(point) - _uv(center)

rmjarvis commented 3 years ago

You guys have been busy here. OK, a few things to respond to. First, your current method isn't really doing a non-uniform WCS. When you add up all the objects and draw them all together, it draws them with the local WCS at the center of the image. So if you're happy with that functionality, I'd avoid using CelestialWCS for any of this. It's an extra layer of complication, which is not likely to add any useful realism to the sim. So, with that slight adjustment to the spec sheet, I think I have a script that runs things both ways so you can see how it works. I'll send the full script on Slack, but here are the relevant snippets:

#
# Method 1.  Add up all objects into a single uber object and draw that.
#

# Position are all shifts relative to the image center in (u,v) coordinates.
uber_obj = galsim.Add([ obj.shift(pos) for obj,pos in zip(all_obj, all_pos) ])
uber_obj = uber_obj.shear(applied_shear)

image1 = galsim.ImageD(bounds=bounds, wcs=wcs)
uber_obj.drawImage(image1, center=image1.true_center)
#
# Method 2.  Draw each object separately and add to the final image
#

image2 = galsim.ImageD(bounds=bounds, wcs=wcs)

xy_center = image2.true_center
uv_center = wcs.toWorld(xy_center)

for obj,pos in zip(all_obj, all_pos):
    sheared_obj = obj.shear(applied_shear)
    sheared_pos = pos.shear(applied_shear)
    image_pos = wcs.toImage(uv_center + sheared_pos)
    stamp = sheared_obj.drawImage(center=image_pos, wcs=wcs, dtype=float)
    b = stamp.bounds & image2.bounds
    image2[b] += stamp[b]
rmjarvis commented 3 years ago

I get that the two images are equal to <1.e-3 of the peak. I'm a little surprised that it wasn't better than this. I think I expected more like 1.e-4, but it's the same level of accuracy when I use a trivial WCS and no applied shear, so I don't think it's an error in any of these steps. Maybe just an underlying precision level of the ability of GalSim to apply a large shift (for method 1) versus drawing something centered in the stamp (method 2).

beckermr commented 3 years ago

Maybe just an underlying precision level of the ability of GalSim to apply a large shift (for method 1) versus drawing something centered in the stamp (method 2).

When I implemented this in our old code base (https://github.com/beckermr/metadetect-sims/blob/master/mdetsims/sim_utils.py#L362), I noticed similar variations in the accuracy.

rmjarvis commented 3 years ago

I thought the b = stamp.bounds & image.bounds was supposed to deal with this

Re this error, the b.area() > 0 check is fine. The canonical way would be

if b.isDefined():
   image2[b] += stamp[b]
esheldon commented 3 years ago

from mike: "First, your current method isn't really doing a non-uniform WCS"

This is surprising. Does this invalidate the tests we have done since it is essentially the same as what metacal does?

esheldon commented 3 years ago

(Note a key part of this sim program was to test non uniform wcs)

rmjarvis commented 3 years ago

With the draw everything on separate stamps method, we could get that to do non-uniform WCS properly. But what you're doing now is definitely not testing that.

esheldon commented 3 years ago

Could this explain why we saw few parts in a thousand bias with Matt's method drawing in stamps?

beckermr commented 3 years ago

What aspects of the WCS we are testing is a bit complicated. For sure when we render SE images, we are using the jacobian approximation. However, the coadding process uses the full WCS. We then use a jacobian approximation at the center of the cell-based coadd. So we are somewhere in the middle.

However, I don't think any of this actually matters that much. In the data, we only care about WCS variations over an arcminute in the coordinate system we choose to make the cell-based coadd. We can easily get rid of those by using a different coordinate system with the center of the projection at the center of the patch instead of at the center of the tract.

beckermr commented 3 years ago

Could this explain why we saw few parts in a thousand bias with Matt's method drawing in stamps?

Yes this is a possibility. We never got to the bottom of it.

beckermr commented 3 years ago

For the record, my rendering technique was trying to

  1. shear the position of the object in the uv-plane about the center of the image
  2. render the sheared object with the correct jacobian using the full non-uniform WCS

This introduces slightly different shearing effects for the positions and the shapes. That could break metadetect maybe.