Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

DM-41994: Use galsim primitives vs higher level interface. #48

Merged
merged 5 commits into from
Dec 7, 2023
Merged

DM-41994: Use galsim primitives vs higher level interface. #48

Changes from 1 commit
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
ddf09c7
Use galsim primitives vs higher level interface.
erykoff Dec 1, 2023
c2874cf
Update to use official galsim._Image API.
erykoff Dec 5, 2023
ec405c5
Change to direct bounds constructor.
erykoff Dec 5, 2023
8843d03
Add test for 100% masked image.
erykoff Dec 5, 2023
fe7a5ef
Add more optimized interfaces.
erykoff Dec 6, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
112 changes: 66 additions & 46 deletions python/lsst/meas/extensions/shapeHSM/_hsm_moments.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,23 @@
]


def _quickGalsimImageView(array, bounds):
match array.dtype:
case np.float64:
gcls = galsim._galsim.ImageViewD
case np.float32:
gcls = galsim._galsim.ImageViewF
case np.int32:
gcls = galsim._galsim.ImageViewI

return gcls(
array.__array_interface__["data"][0],
array.strides[1]//array.itemsize,
array.strides[0]//array.itemsize,
bounds._b,
)
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this actually faster than using galsim._Image(array, bounds, None)? That's the supported API that skips most of the sanity checks. This interface isn't guaranteed to be stable across versions.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks! I didn't know about this, and will take a look. It does look like it should do what we need, and will confirm.



class HsmMomentsConfig(measBase.SingleFramePluginConfig):
"""Base configuration for HSM adaptive moments measurement."""

Expand Down Expand Up @@ -91,12 +108,11 @@ def _calculate(
self,
record: afwTable.SourceRecord,
*,
image: galsim.Image,
weight: galsim.Image | None = None,
badpix: galsim.Image | None = None,
image: galsim._galsim.ImageViewD | galsim._galsim.ImageViewF,
weight_image: galsim._galsim.ImageViewI,
centroid: Point2D,
sigma: float = 5.0,
precision: float = 1.0e-6,
centroid: Point2D | None = None,
erykoff marked this conversation as resolved.
Show resolved Hide resolved
) -> None:
"""
Calculate adaptive moments using GalSim's HSM and modify the record in
Expand All @@ -106,22 +122,16 @@ def _calculate(
----------
record : `~lsst.afw.table.SourceRecord`
Record to store measurements.
image : `~galsim.Image`
image : `~galsim._galsim.ImageViewF` or `~galsim._galsim.ImageViewD`
Image on which to perform measurements.
weight : `~galsim.Image`, optional
The weight image for the galaxy being measured. Can be an int or a
float array. No weighting is done if None. Default is None.
badpix : `~galsim.Image`, optional
Image representing bad pixels, where zero indicates good pixels and
any nonzero value denotes a bad pixel. No bad pixel masking is done
if None. Default is None.
weight_image : `~galsim._galsim.ImageViewI`
The combined badpix/weight image for input to galsim HSM code.
centroid : `~lsst.geom.Point2D`
Centroid guess for HSM adaptive moments.
sigma : `float`, optional
Estimate of object's Gaussian sigma in pixels. Default is 5.0.
precision : `float`, optional
Precision for HSM adaptive moments. Default is 1.0e-6.
centroid : `~lsst.geom.Point2D`, optional
Centroid guess for HSM adaptive moments, defaulting to the image's
true center if None. Default is None.

Raises
------
Expand All @@ -130,21 +140,25 @@ def _calculate(
"""
# Convert centroid to GalSim's PositionD type.
guessCentroid = galsim.PositionD(centroid.x, centroid.y)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Make this _PositionD


try:
# Attempt to compute HSM moments.
shape = galsim.hsm.FindAdaptiveMom(

# Use galsim c++/python interface directly.
shape = galsim.hsm.ShapeData()
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pass in arguments here: observed_shape=galsim._Shear(), psf_shape=galsim._Shear(), moments_centroid=galsim._PositionD(). These avoid initializing objects underneath without the overhead.

hsmparams = galsim.hsm.HSMParams.default

galsim._galsim.FindAdaptiveMomView(
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Even this isn't really guaranteed to be stable, but it should be very close to equivalent to image.FindAdaptiveMom(...).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So there was actually a significant amount of time spent in _convertMask which was some of the motivation to use the primitive. If there was an officially supported way of bypassing this I think we could use the official API.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh right. The direct Python API calls the _convertMask function which has a bunch of galsim.Image instantiations making it slow as well. I could open an issue on the GalSim repo if there's no immediate supported API.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's fair. What mask do you use? None? Or something real?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We supply a value for badpix that masks some pixels, so very much real. I think the overheads there can be reduced by using galsim._ImageI instead of galsim.ImageI in _convertMask.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

GalSim-developers/GalSim#1264

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Right now, this is the only sticking point I believe. Either we use the FindAdaptiveMom call or have a TODO comment here referring to a ticket that should replace it with FindAdaptiveMom as soon as we have a version of GalSim with the above issue implemented.

shape._data,
image,
weight=weight,
badpix=badpix,
guess_sig=sigma,
precision=precision,
guess_centroid=guessCentroid,
strict=True,
round_moments=self.config.roundMoments,
hsmparams=None,
weight_image,
float(sigma),
float(precision),
guessCentroid._p,
bool(self.config.roundMoments),
erykoff marked this conversation as resolved.
Show resolved Hide resolved
hsmparams._hsmp,
)
except galsim.hsm.GalSimHSMError as error:

except RuntimeError as error:
raise measBase.MeasurementError(str(error), self.GALSIM.number)

# Retrieve computed moments sigma and centroid.
Expand Down Expand Up @@ -278,19 +292,20 @@ def measure(self, record, exposure):
# Create a GalSim image using the extracted array.
# NOTE: GalSim's HSM uses the FITS convention of 1,1 for the
# lower-left corner.
image = galsim.Image(imageArray, bounds=bounds, copy=False)
_image = _quickGalsimImageView(imageArray, bounds)
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here, I would recommend image = galsim._Image(imageArray, bounds, None).


# Convert the mask of bad pixels to a format suitable for galsim.
# NOTE: galsim.Image will match whatever dtype the input array is
# (here int32).
badpix = galsim.Image(badpix, bounds=bounds, copy=False)
gd = (badpix == 0)
badpix[gd] = 1
badpix[~gd] = 0

_weight_image = _quickGalsimImageView(badpix, bounds)

# Call the internal method to calculate adaptive moments using GalSim.
self._calculate(
record,
image=image,
weight=None,
badpix=badpix,
image=_image,
weight_image=_weight_image,
sigma=2.5 * psfSigma,
precision=1.0e-6,
centroid=center,
Expand Down Expand Up @@ -430,7 +445,17 @@ def measure(self, record, exposure):
imageArray = psfImage.array

# Create a GalSim image using the PSF image array.
image = galsim.Image(imageArray, bounds=bounds, copy=False)
_image = _quickGalsimImageView(imageArray, bounds)
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also here. image = galsim._Image(imageArray, bounds, None).


if badpix is not None:
gd = (badpix == 0)
badpix[gd] = 1
badpix[~gd] = 0

_weight_image = _quickGalsimImageView(badpix, bounds)
else:
arr = np.ones(imageArray.shape, dtype=np.int32)
_weight_image = _quickGalsimImageView(arr, bounds)

# Decide on the centroid position based on configuration.
if self.config.useSourceCentroidOffset:
Expand All @@ -444,9 +469,8 @@ def measure(self, record, exposure):
# Call the internal method to calculate adaptive moments using GalSim.
self._calculate(
record,
image=image,
weight=None,
badpix=badpix,
image=_image,
weight_image=_weight_image,
sigma=psfSigma,
centroid=centroid,
)
Expand Down Expand Up @@ -499,7 +523,7 @@ def _adjustNoise(
exposure: afwImage.Exposure,
record: afwTable.SourceRecord,
bounds: galsim.bounds.BoundsI,
) -> galsim.Image:
) -> np.ndarray:
"""
Adjusts noise in the PSF image and updates the bad pixel mask based on
exposure data. This method modifies `psfImage` in place and returns a
Expand All @@ -522,9 +546,9 @@ def _adjustNoise(

Returns
-------
badpix : `~galSim.Image`
Image representing bad pixels, where zero indicates good pixels and
any nonzero value denotes a bad pixel.
badpix : `~np.ndarray`
Numpy image array (np.int32) representing bad pixels, where zero
indicates good pixels and any nonzero value denotes a bad pixel.

Raises
------
Expand Down Expand Up @@ -589,14 +613,10 @@ def _adjustNoise(
overlap = badpix.getBBox()
overlap.clip(exposure.getBBox())
badpix[overlap] = exposure.mask[overlap]
# Pull out the numpy view of the badpix mask image.
badpix = badpix.array

bitValue = exposure.mask.getPlaneBitMask(self.config.badMaskPlanes)
badpix &= bitValue

# Convert the mask of bad pixels to a format suitable for galsim.
# NOTE: galsim.Image will match whatever dtype the input array is
# (here int32).
badpix = galsim.Image(badpix, bounds=bounds, copy=False)

return badpix