Calculate centers

[esheldon]

This adds code to calculate centers from shot photons

These positions will be accurate in all cases, whereas the current available x, y are from converting ra, dec to x, y using the batoid wcs. That wcs does not include all physics (e.g. DCR) making the WCS inaccurate in some cases.

Currently this is "always on" and uses a fixed number of photons.

I need help to understand how to make this configurable within the galsim config system, as imsim interfaces with it.

[esheldon]

The tests such as test_lsst_silicon_builder_passes_correct_photon_ops_to_drawImage don't actually result in any flux in the image, so the code is currently failing. Was no flux drawn the intention?

[esheldon]

Anyone know how to make the tests draw a non zero flux?

[rmjarvis]

Rather than change what that test does, probably better to make this centroid calculation optional, which we need to do anyway. We definitely don't want it turned on by default, since it's much too slow for objects that have very few photons. (Which are the vast majority of objects in typical use cases.)

[esheldon]

The recent commits make it optional

[esheldon]

I can't make any progress on this PR until I can get the tests to actually produce photons. Currently the tests as set up produce zero flux drawn in the stamp but I don't understand why. Does anyone know how to modify the test to produce flux?

[rmjarvis]

Don't use mock. Have this call the real drawImage function, and then it will produce an actual output image. The above tests that use mock were testing something very specific about how parameters were being passed to drawImage, which mock lets you test. I think you want to actually call the drawImage function and check that the centroid values get computed correctly.

[esheldon]

Yeah, I didn't know what that context did, but I found I had to use it because otherwise I got this error:

imsim/stamp.py:778: in draw
    image = gal.drawImage(bp_for_drawImage,
../../miniforge3/envs/stack/lib/python3.11/site-packages/galsim/chromatic.py:2136: in drawImage
    image = ChromaticObject.drawImage(self, bandpass, image, integrator, **kwargs)
../../miniforge3/envs/stack/lib/python3.11/site-packages/galsim/chromatic.py:485: in drawImage
    image = prof0.drawImage(image=image, **kwargs)
../../miniforge3/envs/stack/lib/python3.11/site-packages/galsim/gsobject.py:1775: in drawImage
    added_photons, photons = prof.drawPhot(imview, gain, add_to_image,
../../miniforge3/envs/stack/lib/python3.11/site-packages/galsim/gsobject.py:2383: in drawPhot
    op.applyTo(photons, local_wcs, rng)
imsim/photon_ops.py:90: in applyTo
    v = self.photon_velocity(photon_array, local_wcs, rng)
imsim/photon_ops.py:189: in photon_velocity
    return self.rubin_diffraction.photon_velocity(
imsim/photon_ops.py:267: in photon_velocity
    v = photon_velocity(
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

photon_array = galsim.PhotonArray(100, x=array([-4.95703380547105, -1.7774137125821585, 5.4679574148088355, 1.3853752380310527, -11.0...28, 567.8319745320754, 610.7805936794955, 588.7526624330898, 631.0916929335601, 659.
8606846139876, 662.0577956513312]))
xy_to_v = <imsim.photon_ops.XyToV object at 0x7165d0d3b350>, get_n = <bound method Medium.getN of ConstMedium(1.0)>

    def photon_velocity(photon_array, xy_to_v: "XyToV", get_n) -> np.ndarray:
        """Computes the velocity of a photon array."""
    
>       assert photon_array.hasAllocatedPupil()
E       AssertionError

imsim/photon_ops.py:123: AssertionError

[rmjarvis]

You need to either use an atmosphere (as PSF or photon_op) or PupilAnnulusSampler. The RubinDiffractionOptics photon_op requires that something else has already populated the pupil plane positions. In normal imsim usage, it's the AtmosphericPSF. But you can also do it manually with PupilAnnulusSampler (or PupilImageSampler works too, but you probably don't want to use that one.)

[esheldon]

I also had to add a TimeSampler

[rmjarvis]

Yeah, that sounds right.

[esheldon]

I think this is ready for review.

[rmjarvis]

The code looks fine for what it is, but I don't love this API for turning on the calculation. I think I'd rather have two optional parameters in the stamp: calculate_centroids: bool and centroid_n_photons: int. The latter should have a documented default value (maybe 1000 is good, though you mention 100_000 in the doc above). Then to turn this on, most people would just set calculate_centroids: True. And if they want to adjust the precision, they could also set centroid_n_photons to something.

[esheldon]

I have my own issues with this PR, mainly that it is calculating a centroid for every object which is not necessary.

For images with a lot of galaxies this slows down the code significantly.

We only need a subset of the objects, say a random 100 or 200, with which to fit a WCS and validate it.