diff --git a/include/lsst/meas/algorithms/SpanSetMoments.h b/include/lsst/meas/algorithms/SpanSetMoments.h new file mode 100644 index 000000000..195dc06e2 --- /dev/null +++ b/include/lsst/meas/algorithms/SpanSetMoments.h @@ -0,0 +1,132 @@ +// -*- lsst-c++ -*- +/* + * LSST Data Management System + * Copyright 2008-2013 LSST Corporation. + * + * This product includes software developed by the + * LSST Project (http://www.lsst.org/). + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the LSST License Statement and + * the GNU General Public License along with this program. If not, + * see . + */ + +#ifndef LSST_MEAS_ALGORITHMS_SpanSetMoments_h_INCLUDED +#define LSST_MEAS_ALGORITHMS_SpanSetMoments_h_INCLUDED + +#include "ndarray.h" +#include "lsst/geom/Point.h" +#include "lsst/afw/image/MaskedImage.h" +#include "lsst/afw/geom/ellipses/Quadrupole.h" +#include "lsst/afw/geom/SpanSet.h" +#include "lsst/shapelet/ShapeletFunction.h" + +namespace lsst { +namespace meas { +namespace algorithms { + +/** + * A struct that computes the unweighted moments of the pixels in an + * `lsst.afw.geom.SpanSet`. + * + * This class provides low-level pixel-processing code for + * `ComputeRoughPsfShapeletsTask`. + */ +struct SpanSetMoments { + + /// Total flux within the SpanSet. + double flux; + + /// Total variance within the SpanSet. + double variance; + + /// Center derived from the unweighted first moment of the image. + geom::Point2D center; + + /// Shape derived from the unweighted second moment of the image. + afw::geom::ellipses::Quadrupole shape; + + /// The pixels actually used to compute the moments. + std::shared_ptr spans; + + /// Flag set if there were too many bad pixels to compute the moments. + bool too_many_bad_pixels = false; + + /// Flag set if the center did not lie within the SpanSet. + bool center_out_of_bounds = false; + + /// Flag set if there was a bad pixel too close to the center. + bool bad_pixel_in_center = false; + + /// Flag set if the second moments did not resolve to an ellipse. + bool singular_second_moments = false; + + /// Test whether any failure flag is set. + bool any_flags_set() const { + return too_many_bad_pixels || center_out_of_bounds || bad_pixel_in_center || singular_second_moments; + } + + /// Return a flattened array of the x coordinates in `spans`. + ndarray::Array get_x_array() const; + + /// Return a flattened array of the y coordinates in `spans`. + ndarray::Array get_y_array() const; + + /** + * Compute the unweighted moments of an image within a SpanSet. + * + * @param[in] spans Pixel region to use. + * @param[in] masked_image Image to measure. + * @param[in] bad_bitmask Mask of bad pixels to remove from `spans` + * before computing the moments. + * @param[in] bad_pixel_mask_fraction Maximum fraction of the pixels + * in `spans` that can be bad + * before giving up. + * @param[in] bad_pixel_exlusion_radius Radius around the estimated + * center where the present of a + * bad pixels will cause the + * algorithm to give up. + */ + static std::shared_ptr compute(afw::geom::SpanSet const& spans, + afw::image::MaskedImage const& masked_image, + afw::image::MaskPixel bad_bitmask, + double bad_pixel_max_fraction, + double bad_pixel_exclusion_radius); + + /** + * Fit a common shapelet expansion to multiple sources whose moments have + * already been computed. + * + * The expected use case is a sample of stars with similar but not + * identical moments, which might plausibly have an approximately common + * shapelet representation while benefitting from allowing the ellipse + * used for the expansion to vary (i.e. to partially capture spatial + * variation in the PSF). + * + * @param[in] masked_image Image ot measure. + * @param[in] moments Moments already measured on stars. + * @param[in] order Order of the shapelet expansion. + * @param[in] scale Factor to scale the moments ellipses by to + * yield the ellipse for the shapelet fit. + */ + static shapelet::ShapeletFunction fit_shapelets( + afw::image::MaskedImage const& masked_image, + std::vector> const& moments, + int order, double scale); +}; + +} // namespace algorithms +} // namespace meas +} // namespace lsst + +#endif // !LSST_MEAS_ALGORITHMS_SpanSetMoments_h_INCLUDED diff --git a/python/lsst/meas/algorithms/SConscript b/python/lsst/meas/algorithms/SConscript index ffa853e99..4f285f874 100644 --- a/python/lsst/meas/algorithms/SConscript +++ b/python/lsst/meas/algorithms/SConscript @@ -15,5 +15,6 @@ scripts.BasicSConscript.python(['_algorithmsLib'], [ 'pcaPsf.cc', 'psfCandidate/psfCandidate.cc', 'singleGaussianPsf.cc', + 'spanSetMoments.cc', 'spatialModelPsf.cc', 'warpedPsf.cc']) diff --git a/python/lsst/meas/algorithms/_algorithmsLib.cc b/python/lsst/meas/algorithms/_algorithmsLib.cc index fe70c2035..9b0f940c0 100644 --- a/python/lsst/meas/algorithms/_algorithmsLib.cc +++ b/python/lsst/meas/algorithms/_algorithmsLib.cc @@ -45,6 +45,7 @@ void wrapCoaddPsf(WrapperCollection &wrappers); void wrapCoaddBoundedField(WrapperCollection &wrappers); void wrapCr(WrapperCollection &wrappers); void wrapCloughTocher2DInterpolatorUtils(WrapperCollection &wrappers); +void wrapSpanSetMoments(WrapperCollection & wrappers); PYBIND11_MODULE(_algorithmsLib, mod) { WrapperCollection wrappers(mod, "lsst.meas.algorithms"); @@ -67,9 +68,10 @@ WrapperCollection wrappers(mod, "lsst.meas.algorithms"); wrapCoaddBoundedField(wrappers); wrapCr(wrappers); wrapCloughTocher2DInterpolatorUtils(wrappers); + wrapSpanSetMoments(wrappers); wrappers.finish(); } } // algorithms } // meas -} // lsst \ No newline at end of file +} // lsst diff --git a/python/lsst/meas/algorithms/computePsfChiSq.py b/python/lsst/meas/algorithms/computePsfChiSq.py new file mode 100644 index 000000000..5d8d8d38c --- /dev/null +++ b/python/lsst/meas/algorithms/computePsfChiSq.py @@ -0,0 +1,90 @@ +# This file is part of meas_algorithms. +# +# Developed for the LSST Data Management System. +# This product includes software developed by the LSST Project +# (https://www.lsst.org). +# See the COPYRIGHT file at the top-level directory of this distribution +# for details of code ownership. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +__all__ = ["ComputePsfChiSqTask", "ComputePsfChiSqConfig"] + +import numpy as np + +from lsst.afw.detection import Psf +from lsst.afw.image import ImageF, MaskedImageF +from lsst.geom import Point2D +from lsst.pex.config import Config, Field, ListField +from lsst.pipe.base import Struct, Task + + +class ComputePsfChiSqConfig(Config): + padding = Field( + "Number of pixels to grow the PSF model bbox by on all sides.", + dtype=int, + default=10, + ) + bad_mask_planes = ListField( + "Mask planes to identify pixels to drop from the calculation.", + dtype=str, + default=["SAT", "SUSPECT"], + ) + + +class ComputePsfChiSqTask(Task): + ConfigClass = ComputePsfChiSqConfig + _DefaultName = "computePsfChiSq" + config: ComputePsfChiSqConfig + + def run( + self, *, masked_image: MaskedImageF, psf: Psf, x: np.ndarray, y: np.ndarray + ) -> Struct: + bitmask = masked_image.mask.getPlaneBitMask(self.config.bad_mask_planes) + result = Struct( + chisq=np.zeros(x.shape, dtype=np.float32), + snr=np.zeros(x.shape, dtype=np.float32), + n_missing_pixels=np.zeros(x.shape, dtype=np.uint16), + ) + for i, (sx, sy) in enumerate(zip(x, y)): + inner_model = psf.computeImage(Point2D(sx, sy)) + inner_bbox = inner_model.getBBox() + outer_bbox = inner_bbox.dilatedBy(self.config.padding) + n_clipped_pixels = 0 + if not masked_image.getBBox().contains(outer_bbox): + original_bbox_area = outer_bbox.getArea() + outer_bbox.clip(masked_image.getBBox()) + n_clipped_pixels = original_bbox_area - outer_bbox.getArea() + outer_model = ImageF(outer_bbox) + outer_model.array[:, :] = 0 + outer_model[inner_bbox] = inner_model.convertF() + outer_data = masked_image.image[outer_bbox] + outer_variance = masked_image.variance[outer_bbox] + mask = np.logical_not(masked_image.mask[outer_bbox].array & bitmask) + flat_model = outer_model.array[mask] + flat_data = outer_data.array[mask] + flat_variance = outer_variance.array[mask] + sigma = flat_variance**0.5 + weighted_model = flat_model / sigma + weighted_data = flat_data / sigma + inv_var_alpha = np.dot(weighted_model, weighted_model) + alpha = np.dot(weighted_model, weighted_data) / inv_var_alpha + weighted_residuals = weighted_model * alpha + weighted_residuals -= weighted_data + result.chisq[i] = sum(weighted_residuals**2) / np.sum(mask) + result.snr[i] = alpha * inv_var_alpha**0.5 + result.n_missing_pixels[i] = ( + outer_bbox.getArea() - np.sum(mask) + n_clipped_pixels + ) + return result diff --git a/python/lsst/meas/algorithms/computeRoughPsfShapelets.py b/python/lsst/meas/algorithms/computeRoughPsfShapelets.py new file mode 100644 index 000000000..23a8946cc --- /dev/null +++ b/python/lsst/meas/algorithms/computeRoughPsfShapelets.py @@ -0,0 +1,785 @@ +# This file is part of meas_algorithms. +# +# Developed for the LSST Data Management System. +# This product includes software developed by the LSST Project +# (https://www.lsst.org). +# See the COPYRIGHT file at the top-level directory of this distribution +# for details of code ownership. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +from __future__ import annotations + +__all__ = ["ComputeRoughPsfShapeletsTask", "ComputeRoughPsfShapeletsConfig"] + +import itertools +from typing import TYPE_CHECKING, Any, Literal + +import numpy as np +import scipy.signal +import scipy.stats +from sklearn.covariance import MinCovDet +from sklearn.neighbors import KernelDensity + +from lsst.afw.geom import ellipses +from lsst.afw.image import ImageD, ImageF, MaskedImageF +from lsst.afw.table import Point2DKey, QuadrupoleKey, Schema, SourceCatalog +from lsst.geom import Box2D +from lsst.pex.config import Config, Field, ListField +from lsst.pipe.base import AlgorithmError, Struct, Task +from lsst.shapelet import LAGUERRE, ShapeletFunction, computeOffset + +from ._algorithmsLib import SpanSetMoments + +if TYPE_CHECKING: + import matplotlib.axes + import matplotlib.figure + import matplotlib.image + import matplotlib.patches + + +class NoStarsForShapeletsError(AlgorithmError): + """Exception raised when ComputeRawPsfMomentsTask fails to find any usable + stars. + """ + + @property + def metadata(self) -> dict[str, Any]: + return {} + + +class ComputeRoughPsfShapeletsConfig(Config): + bad_mask_planes = ListField( + "Mask planes to identify pixels to drop from the calculation.", + dtype=str, + default=["SAT", "SUSPECT", "INTRP"], + ) + bad_pixel_max_fraction = Field( + "Maximum fraction of a footprint's pixels that can be bad (according " + "to bad_mask_planes) before that footprint is fully dropped.", + dtype=float, + default=0.25, + ) + bad_pixel_exclusion_radius = Field( + "If a bad pixel (according to bad_mask_planes) falls within this " + "radius of the unweighted centroid of a footprint, drop that footprint " + "entirely.", + dtype=float, + default=2.0, + ) + max_footprint_area = Field( + "Footprints larger than this threshold are dropped before computing moments.", + dtype=int, + default=10000, + ) + min_snr = Field( + "Mininum flux S/N for inclusion in the star sample.", + dtype=float, + default=50.0, + ) + max_radius_factor = Field( + "Maximium multiple of the mode of the radius distribution for inclusion in the star sample.", + dtype=float, + default=2.0, + ) + max_shape_distance = Field( + "Maximum Mahalanobis distance (distance from the center of the shape " + "distribution in elliptical sigma units) to select a star from the candidate " + "sample, comparing the shape of that source vs. a robust estimate of the " + "distribution of the shapes of all sources.", + dtype=float, + default=3.0, + ) + min_n_stars = Field( + "Minimum number of stars to select. The S/N, radius, and shape distance thresholds are " + "relaxed as needed to meet this target.", + dtype=int, + default=10, + ) + max_n_stars = Field( + "Maximum number of stars to select. High shape-distance sources " + "are dropped as needed to meet this target.", + dtype=int, + default=20, + ) + logarithmic_shapes = Field( + "If True, transform the (xx, yy, xy) moments to conformal shear " + "and log trace radius when selecting stars in order to map the ellipse " + "parameter to a space with (-inf, inf) bounds on all quantities (but " + "a less-linear relationship to the pixel data).", + dtype=bool, + default=False, + ) + radius_mode_min = Field( + "Minimum radius at which to start searching for the first mode. " + "This just needs be large enough to avoid unphysically small garbage (e.g. CRs).", + dtype=float, + default=1.0, + ) + radius_kde_bandwidth = Field( + "Bandwidth of the Gaussian kernel density estimator used to find the " + "first mode of the radius distribution.", + dtype=float, + default=1.0, + ) + shapelet_order = Field( + "Order of the shapelet expansion fit to the stars.", + dtype=int, + default=4, + ) + shapelet_scale_factor = Field( + "Scale fractor to apply to the moments ellipse when computing the ellipse for the shapelet basis.", + dtype=float, + default=1.0, + ) + + def validate(self) -> None: + if self.min_n_stars > self.max_n_stars: + raise ValueError( + f"min_n_stars={self.min_n_stars} is greater than max_n_stars={self.max_n_stars}." + ) + if self.shapelet_order < 0: + raise ValueError(f"shapelet order {self.shapelet_order} must be nonnegative.") + if self.shapelet_scale_factor <= 0.0: + raise ValueError(f"shapelet scale factor {self.shapelet_scale_factor} must be positive.") + + +class ComputeRoughPsfShapeletsTask(Task): + """A task that computes a rough shapelet expansion of the PSF from a set + of high S/N detections. + + Notes + ----- + This task is expected to be run early in single-epoch processing - just + after background subtraction and an initial high S/N detection phase, and + before any deblending or measurement - in order to identify out-of-focus + or otherwise bad PSFs. + + Given a background-subtracted `lsst.afw.image.MaskedImage`, an + `lsst.afw.table.SourceCatalog` with footprints attached, and a random + number generator seed, the `run` method will: + + - Compute the *unweighted* 0th-2nd moments of every non-child source over + the footprint (except certain configurable masked pixels). This is + delegated to the `compute_raw_moments` method (which uses the C++ + `SpanSetMoments` class for the pixel-level processing). Unweighted + moments are used to avoid "latching onto" a small piece of PSF + substructure, but can be much noiser than the Gaussian-weighed moments + we usually use. + + - Select a "candidate" sample of sources with successfully measured + moments that satisfy a S/N cut and a radius cut (determined from the + first mode of the radius distribution, via kernel density estimation), + and then use a robust covariance estimator (`scikit_learn.MinCovDet`) + to select presumed isolated stars that are close to the center of that + distribution, in 3-parameter shape space. This is delegated to the + `select_stars` method. + + - Fit a single shapelet expansion to the selected stars. This is + mostly delegated to the `SpanSetMoments.fit_shapelets` method. + + The radial shapelet terms at 0th, 2nd, and 4th order are expected to form + a space in which donut-shaped PSFs are well-separated from those with + monotonic profiles. Other terms *may* be useful in identifying other + kinds of undesirable PSF structure. + """ + + ConfigClass = ComputeRoughPsfShapeletsConfig + _DefaultName = "computeRoughPsfShapelets" + config: ComputeRoughPsfShapeletsConfig + + def __init__( + self, + config: ComputeRoughPsfShapeletsConfig | None = None, + *, + schema: Schema, + **kwargs: Any, + ): + super().__init__(config=config, **kwargs) + self.schema = schema + self._flux_key = schema.addField("RawPsfMoments_flux", type=float, doc="Unweighted zeroth moment.") + self._flux_err_key = schema.addField( + "RoughPsfShapelets_fluxErr", + type=float, + doc="Uncertainty on the unweighted zeroth moment.", + ) + self._center_key = Point2DKey.addFields( + schema, "RoughPsfShapelets", "Center from unweighted first moments.", "pixels" + ) + self._shape_key = QuadrupoleKey.addFields( + schema, "RoughPsfShapelets", "Shape from unweighted second moments." + ) + self._flag_key = schema.addField( + "RoughPsfShapelets_flag", + type="Flag", + doc="Flag set if the raw PSF moments were not computed.", + ) + self._candidate_key = schema.addField( + "RoughPsfShapelets_candidate", + type="Flag", + doc="Flag set if this source passed the radius_fraction cut (see configuration).", + ) + self._used_key = schema.addField( + "RoughPsfShapelets_used", + type="Flag", + doc=( + "Flag set if this source passed the radius_fraction and shape_distance cuts " + "(see configuration) and was used to fit the shapelet expansion." + ), + ) + + def run(self, *, masked_image: MaskedImageF, catalog: SourceCatalog, seed: int) -> Struct: + """Compute raw moments, select stars, and fit a shapelet expansion to + them. + + Parameters + ---------- + masked_image + Masked image to measure on. Must be background-subtracted. + catalog + Catalog of detections to extract footprints from and fill output + columns of. Its schema must be a superset of ``self.schema``. + seed + A random-number generator seed, used for the robust covariance + estimator. + + Returns + ------- + `lsst.pipe.base.Struct` + A struct of results containing: + + - ``shapelet`` (`lsst.shapelet.ShapeletFunction`): A + Gauss-Laguerre (polar shaplet) expansion of the PSF. + - ``radial`` (`list` [`float`]): the purely radial coefficients + of the shapelet expansion. + - all attributes returned by the `select_stars` method. + """ + moments = self.compute_raw_moments(masked_image=masked_image, catalog=catalog) + result = self.select_stars(catalog, seed=seed) + star_moments = [moments[star_id] for star_id in result.star_ids] + result.shapelet = SpanSetMoments.fit_shapelets( + masked_image, + star_moments, + self.config.shapelet_order, + self.config.shapelet_scale_factor, + ) + result.shapelet.getEllipse().setCore(result.mean_shape) + result.shapelet.changeBasisType(LAGUERRE) + result.radial = result.shapelet.getCoefficients()[ + [computeOffset(i) for i in range(0, self.config.shapelet_order + 1, 2)] + ] + self.log.info("Rough PSF shapelet radial terms: %s.", result.radial) + return result + + def compute_raw_moments( + self, *, masked_image: MaskedImageF, catalog: SourceCatalog + ) -> dict[int, SpanSetMoments]: + """Compute the unweighted moments of the footprints in a catalog. + + Parameters + ---------- + masked_image + Masked image to measure on. Must be background-subtracted. + catalog + Catalog of detections to extract footprints from and fill output + columns of. Its schema must be a superset of ``self.schema``. + + Returns + ------- + `dict` [`int`, `SpanSetMoments`] + Objects used to construct and hold the unweighted moments and the + pixel region used to computed them, keyed by source ID. + """ + bitmask = masked_image.mask.getPlaneBitMask(self.config.bad_mask_planes) + all_moments: dict[int, SpanSetMoments] = {} + for record in catalog: + if record.getParent() != 0: + record.set(self._flag_key, True) + self.log.debug("Skipping child source %s", record.getId()) + continue + footprint_spans = record.getFootprint().getSpans() + if footprint_spans.getArea() > self.config.max_footprint_area: + record.set(self._flag_key, True) + self.log.debug( + "Skipping source %s with footprint area %d > %d.", + record.getId(), + footprint_spans.getArea(), + self.config.max_footprint_area, + ) + continue + moments = SpanSetMoments.compute( + record.getFootprint().getSpans(), + masked_image=masked_image, + bad_bitmask=bitmask, + bad_pixel_max_fraction=self.config.bad_pixel_max_fraction, + bad_pixel_exclusion_radius=self.config.bad_pixel_exclusion_radius, + ) + record.set(self._flux_key, moments.flux) + record.set(self._flux_err_key, moments.variance**0.5) + record.set(self._center_key, moments.center) + record.set(self._shape_key, moments.shape) + record.set(self._flag_key, moments.any_flags_set) + if not moments.any_flags_set: + all_moments[record.getId()] = moments + if all_moments: + self.log.verbose( + "Successfully measured raw moments for %d of %d sources.", + len(all_moments), + len(catalog), + ) + else: + raise NoStarsForShapeletsError("No raw moments could be measured.") + return all_moments + + def select_stars(self, catalog: SourceCatalog, seed: int) -> Struct: + """Select probable stars from the distribution of second moments. + + Parameters + ---------- + catalog + Catalog of detections to extract footprints from and fill output + columns of. Its schema must be a superset of ``self.schema``. + seed + A random-number generator seed, used for the robust covariance + estimator. + + Returns + ------- + `lsst.pipe.base.Struct` + A struct of results containing: + + - ``star_ids`` (`numpy.ndarray`): the source IDs that are expected + to be stars. + - ``mean_shape`` (`lsst.afw.geom.ellipses.BaseCore`): the mean of + the shape distribution. + - ``shape_covariance`` (`numpy.ndarray`): the covariance of the + distribution of shapes; a 3x3 matrix. This uses the same + parameterization of the shapes as ``mean_shape``. + - ``radius_cut`` (`float`): the indended radius cut (i.e. the + mode of the radius distribution multipled by the + ``radius_factor`` configuration option). + - ``radius_kde`` (`sklearn.neighbors.KernelDensity`): kernel + density estimator on the radius distribution, used to determine + the radius cut. + """ + # Cut on flags and SNR first. + indices = np.arange(len(catalog), dtype=int)[np.logical_not(catalog[self._flag_key])] + indices = indices[ + self._threshold_with_bounds( + catalog[self._flux_key][indices] / catalog[self._flux_err_key][indices], + threshold=self.config.min_snr, + min_count=self.config.min_n_stars, + max_count=len(catalog), + name="S/N", + kind=">", + ) + ] + # Cut on radius next. + radii = np.zeros(indices.size, dtype=np.float64) + for n, index in enumerate(indices): + record = catalog[index] + shape = record.get(self._shape_key) + radii[n] = shape.getTraceRadius() + radius_mode, radius_kde = self._find_first_radius_mode(radii) + radius_cut = self.config.max_radius_factor * radius_mode + indices = indices[ + self._threshold_with_bounds( + radii, + threshold=radius_cut, + min_count=self.config.min_n_stars, + max_count=len(catalog), + name="radius", + kind="<", + ) + ] + shape_data = np.zeros((len(indices), 3), dtype=np.float64) + for n, index in enumerate(indices): + record = catalog[index] + record.set(self._candidate_key, True) + shape = record.get(self._shape_key) + if self.config.logarithmic_shapes: + shape = ellipses.SeparableConformalShearLogTraceRadius(shape) + shape_data[n, :] = shape.getParameterVector() + shape_dist = MinCovDet(random_state=seed).fit(shape_data) + m_distances = shape_dist.mahalanobis(shape_data) + indices = indices[ + self._threshold_with_bounds( + m_distances, + threshold=self.config.max_shape_distance, + min_count=self.config.min_n_stars, + max_count=self.config.max_n_stars, + name="Mahalanobis distance", + kind="<", + ) + ] + for index in indices: + catalog[index].set(self._used_key, True) + star_ids = catalog["id"][indices] + if self.config.logarithmic_shapes: + mean_shape = ellipses.SeparableConformalShearLogTraceRadius(shape_dist.location_) + else: + mean_shape = ellipses.Quadrupole(shape_dist.location_) + return Struct( + star_ids=star_ids, + mean_shape=mean_shape, + shape_covariance=shape_dist.covariance_, + radius_cut=radius_cut, + radius_kde=radius_kde, + ) + + def plot_selection( + self, figure: matplotlib.figure.Figure, *, catalog: SourceCatalog, results: Struct + ) -> None: + """Create plots of the shape distribution space used to select stars. + + Parameters + ---------- + figure + Matplotlib figure to plot to. + catalog + Catalog of sources with columns populated by the `run` method (at + least through the `select_stars` step). + results + Result struct returned by `run` or `select_stars`. + """ + from matplotlib.lines import Line2D + + shape_data = np.zeros((len(catalog), 3), dtype=np.float64) + radii = np.zeros(len(catalog), dtype=np.float64) + for n, record in enumerate(catalog): + if record[self._flag_key]: + continue + shape = record[self._shape_key] + if self.config.logarithmic_shapes: + shape = ellipses.SeparableConformalShearLogTraceRadius(shape) + shape_data[n, :] = shape.getParameterVector() + radii[n] = shape.getTraceRadius() + used_mask = catalog[self._used_key] + candidate_mask = np.logical_and(catalog[self._candidate_key], np.logical_not(used_mask)) + measured_mask = np.logical_and( + np.logical_not(catalog[self._flag_key]), np.logical_not(catalog[self._candidate_key]) + ) + # Set up the axes. + axes = figure.subplot_mosaic( + [ + ["radius", "radius", "radius"], + ["hist0", ".", "."], + ["scatter01", "hist1", "."], + ["scatter02", "scatter12", "hist2"], + ], + gridspec_kw=dict(bottom=0.2, top=1.0), + ) + axes["scatter01"].sharex(axes["hist0"]) + axes["scatter02"].sharex(axes["hist0"]) + axes["scatter12"].sharex(axes["hist1"]) + axes["scatter12"].sharey(axes["scatter02"]) + for tk in itertools.chain( + axes["hist0"].get_xticklabels(), + axes["scatter01"].get_xticklabels(), + axes["hist1"].get_xticklabels(), + axes["scatter12"].get_yticklabels(), + ): + tk.set_visible(False) + # Move hist y axes to the outside. + for ax in [axes["hist0"], axes["hist1"], axes["hist2"]]: + ax.yaxis.set_label_position("right") + ax.yaxis.tick_right() + # Add labels to outer axes. + names = ["Ixx", "Iyy", "Ixy"] if not self.config.logarithmic_shapes else ["𝜂1", "𝜂2", "ln(r)"] + axes["scatter02"].set_xlabel(names[0]) + axes["scatter12"].set_xlabel(names[1]) + axes["hist2"].set_xlabel(names[2]) + axes["scatter01"].set_ylabel(names[1]) + axes["scatter02"].set_ylabel(names[2]) + # Make the plots. + mu = results.mean_shape.getParameterVector() + sigma = np.diagonal(results.shape_covariance) ** 0.5 + lower_bounds = [max(mu[i] - 3 * sigma[i], min(shape_data[:, i])) for i in range(3)] + upper_bounds = [min(mu[i] + 3 * sigma[i], max(shape_data[:, i])) for i in range(3)] + grids = [np.linspace(lower_bounds[i], upper_bounds[i], 50) for i in range(3)] + axes["radius"].axvline(results.radius_cut, color="k") + for color, mask, alpha in [ + ("grey", measured_mask, 0.5), + ("blue", candidate_mask, 0.75), + ("green", used_mask, 1.0), + ]: + axes["radius"].hist( + radii[mask], + color=color, + alpha=alpha, + bins=16, + histtype="step", + range=(radii.min(), 15.0), + density=True, + ) + for i in range(3): + axes[f"hist{i}"].hist( + shape_data[mask, i], + bins=16, + range=(lower_bounds[i], upper_bounds[i]), + density=True, + color=color, + histtype="step", + alpha=alpha, + ) + for j in range(3): + if (ax := axes.get(f"scatter{i}{j}")) is not None: + ax.scatter( + shape_data[mask, i], + shape_data[mask, j], + c=color, + s=4, + alpha=alpha, + edgecolors=None, + ) + for i in range(3): + axes[f"hist{i}"].plot( + grids[i], scipy.stats.norm.pdf(grids[i], loc=mu[i], scale=sigma[i]), "k", alpha=0.5 + ) + axes[f"hist{i}"].set_xlim(lower_bounds[i], upper_bounds[i]) + for j in range(3): + if (ax := axes.get(f"scatter{i}{j}")) is not None: + sigma_ellipse = ellipses.Quadrupole( + results.shape_covariance[i, i], + results.shape_covariance[j, j], + results.shape_covariance[i, j], + ) + for factor in [1, 2, 3]: + self._draw_ellipse( + ax, + sigma_ellipse, + x=mu[i], + y=mu[j], + scale=factor, + fill=False, + edgecolor="k", + alpha=0.5, + ) + ax.set_xlim(lower_bounds[i], upper_bounds[i]) + ax.set_ylim(lower_bounds[j], upper_bounds[j]) + figure.legend( + [ + Line2D([], [], color="green", alpha=1.0), + Line2D([], [], color="blue", alpha=0.75), + Line2D([], [], color="gray", alpha=0.5), + ], + [ + f"RoughPsfShapelets_used ({used_mask.sum()})", + f"RoughPsfShapelets_candidate & ~RoughPsfShapelets_used ({candidate_mask.sum()})", + f"~RoughPsfShapelets_flag & ~RoughPsfShapelets_candidate ({measured_mask.sum()})", + ], + loc="lower center", + ) + + def plot_shapelets( + self, + figure: matplotlib.figure.Figure, + *, + image: ImageF, + catalog: SourceCatalog, + results: Struct, + n_stars: int = 3, + stamp_size: float = 2.0, + ) -> None: + """Create data/model/residual plots of stars and the shapelet model. + + Parameters + ---------- + figure + Matplotlib figure to plot to. + image + The image the stars were measured on. + catalog + Catalog of sources with columns populated by the `run` method . + results + Result struct returned by `run`. + n_stars + Number of stars to include. + stamp_size + Stamp size in inches. + """ + from matplotlib.colors import Normalize + + width = stamp_size * 3 + 1.5 + figure.set_size_inches(w=width, h=stamp_size * n_stars) + axes = figure.subplot_mosaic( + [ + ["image_cbar", f"d{star_id}", f"m{star_id}", f"r{star_id}", "res_cbar"] + for star_id in results.star_ids[:n_stars] + ], + gridspec_kw=dict( + wspace=0.01, hspace=0.01, left=0.5 / width, right=1.0 - 0.5 / width, bottom=0.01, top=0.99 + ), + width_ratios=[0.25, stamp_size, stamp_size, stamp_size, 0.25], + ) + for name, ax in axes.items(): + if not name.endswith("cbar"): + ax.axis("off") + norm: Normalize | None = None + res_norm: Normalize | None = None + for star_id in results.star_ids[:n_stars]: + record = catalog.find(star_id) + star_bbox = record.getFootprint().getBBox() + star_model = ImageD(star_bbox) + star_center = record[self._center_key] + star_ellipse = ellipses.Ellipse(ellipses.Axes(record[self._shape_key]), star_center) + star_shapelet = ShapeletFunction(results.shapelet) + star_shapelet.setEllipse(star_ellipse) + star_shapelet.evaluate().addToImage(star_model) + if norm is None: + norm = Normalize(vmin=star_model.array.min(), vmax=star_model.array.max()) + star_image = image[star_bbox].clone() + star_image /= record[self._flux_key] + self._draw_image(axes[f"d{star_id}"], star_image, norm=norm, cmap="YlGnBu") + self._draw_ellipse(axes[f"d{star_id}"], star_ellipse, fill=False, edgecolor="blue", alpha=0.5) + image_plot = self._draw_image(axes[f"m{star_id}"], star_model, norm=norm, cmap="YlGnBu") + self._draw_ellipse(axes[f"m{star_id}"], star_ellipse, fill=False, edgecolor="blue", alpha=0.5) + star_image -= star_model.convertF() + amax = np.abs(star_image.array).max() + if res_norm is None: + res_norm = Normalize(vmin=-amax, vmax=amax) + res_plot = self._draw_image(axes[f"r{star_id}"], star_image, norm=res_norm, cmap="RdBu") + self._draw_ellipse(axes[f"r{star_id}"], star_ellipse, fill=False, edgecolor="blue", alpha=0.5) + figure.colorbar(image_plot, cax=axes["image_cbar"], location="left") + figure.colorbar(res_plot, cax=axes["res_cbar"], location="right") + + def _threshold_with_bounds( + self, + values: np.ndarray, + threshold: float, + min_count: int, + max_count: int, + name: str, + kind: Literal["<", ">"], + ) -> np.ndarray: + """Return the indices of an array that satisfy an inequality + and/or lower and upper bounds on the number of indices returned. + + Parameters + ---------- + values + Array of values to threshold on. + threshold + Threshold value that selected elements must be above or below. + min_count + The minimum number of indices returned. When thresholding would + yield fewer than this number, the threshold is ignored. Note that + the number of indices may still be less than this if the size of + ``values`` is less than this. + max_count + The maximum number of indices returned. + name + Name of the quantity being thresholded, for log messages. + kind + Whether the threshold is a upper bound (``<``) or lower bound + (``>``). This also sets how values are ranked when they are added + or dropped to satisfy the count constraints. + + Returns + ------- + indices + Indices into ``values``. + """ + sorter = values.argsort() + n = np.searchsorted(values[sorter], threshold) + if kind == ">": + sorter = sorter[::-1] + n = len(sorter) - n + if n < min_count: + self.log.verbose( + "Applying a %s %s %f cut yields only %d sources; keeping the top %d (%s %s %f) instead.", + name, + kind, + threshold, + n, + min_count, + name, + kind, + values[sorter[min_count - 1]], + ) + n = min_count + elif n > max_count: + self.log.verbose( + "%d sources have %s %s %f; keeping only the top %d (%s %s %f) instead.", + n, + name, + kind, + threshold, + max_count, + name, + kind, + values[sorter[max_count - 1]], + ) + n = max_count + else: + self.log.verbose("Keeping %d sources with %s %s %f.", n, name, kind, threshold) + return sorter[:n] + + def _find_first_radius_mode(self, radii: np.ndarray) -> tuple[float, KernelDensity]: + """Find the first peak in a 1-d distribution of radii.""" + kde = KernelDensity(bandwidth=self.config.radius_kde_bandwidth).fit(radii.reshape(-1, 1)) + sorted_radii = radii.copy() + sorted_radii.sort() + sorted_radii = sorted_radii[sorted_radii.searchsorted(self.config.radius_mode_min):] + scores = kde.score_samples(radii.reshape(-1, 1)) + peaks, _ = scipy.signal.find_peaks(scores) + if not peaks.size: + raise NoStarsForShapeletsError("Radius distribute has no mode.") + return sorted_radii[peaks.min()], kde + + @staticmethod + def _draw_ellipse( + axes: matplotlib.axes.Axes, + ellipse: ellipses.BaseCore | ellipses.Ellipse, + *, + x: float | None = None, + y: float | None = None, + scale: float = 1.0, + **kwargs: Any, + ) -> matplotlib.patches.Ellipse: + from matplotlib.patches import Ellipse as EllipsePatch + + if isinstance(ellipse, ellipses.Ellipse): + if x is None: + x = ellipse.getCenter().getX() + if y is None: + y = ellipse.getCenter().getY() + ellipse = ellipse.getCore() + else: + if x is None: + x = 0.0 + if y is None: + y = 0.0 + ellipse = ellipses.Axes(ellipse) + patch = EllipsePatch( + (x, y), + ellipse.getA() * 2 * scale, # factor of 2 for radius->diameter + ellipse.getB() * 2 * scale, + angle=ellipse.getTheta() * 180.0 / np.pi, + **kwargs, + ) + axes.add_patch(patch) + return patch + + @staticmethod + def _draw_image(axes: matplotlib.axes.Axes, image: ImageF, **kwargs: Any) -> matplotlib.image.AxesImage: + fp_bbox = Box2D(image.getBBox()) + return axes.imshow( + image.array, + interpolation="nearest", + origin="lower", + aspect="equal", + extent=(fp_bbox.x.min, fp_bbox.x.max, fp_bbox.y.min, fp_bbox.y.max), + **kwargs, + ) diff --git a/python/lsst/meas/algorithms/spanSetMoments.cc b/python/lsst/meas/algorithms/spanSetMoments.cc new file mode 100644 index 000000000..fc72571cb --- /dev/null +++ b/python/lsst/meas/algorithms/spanSetMoments.cc @@ -0,0 +1,67 @@ +/* + * LSST Data Management System + * + * This product includes software developed by the + * LSST Project (http://www.lsst.org/). + * See the COPYRIGHT file + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the LSST License Statement and + * the GNU General Public License along with this program. If not, + * see . + */ +#include "pybind11/pybind11.h" +#include "pybind11/stl.h" +#include "lsst/cpputils/python.h" + +#include "lsst/meas/algorithms/SpanSetMoments.h" + +namespace py = pybind11; +using namespace pybind11::literals; + +namespace lsst { +namespace meas { +namespace algorithms { + +void wrapSpanSetMoments(lsst::cpputils::python::WrapperCollection& wrappers) { + using PySpanSetMoments = py::classh; + auto clsSpanSetMoments = + wrappers.wrapType(PySpanSetMoments(wrappers.module, "SpanSetMoments"), [](auto& mod, auto& cls) { + cls.def_readonly("flux", &SpanSetMoments::flux); + cls.def_readonly("variance", &SpanSetMoments::variance); + // def_readonly doesn't ever seem to do the right thing w.r.t. + // memory management, even when explicit return value policies. + // Copying in a getter seems to be the only safe way to avoid + // issues. + cls.def_property_readonly( + "center", [](SpanSetMoments const& self) -> geom::Point2D { return self.center; }); + cls.def_property_readonly("shape", + [](SpanSetMoments const& self) -> afw::geom::ellipses::Quadrupole { + return self.shape; + }); + cls.def_readonly("spans", &SpanSetMoments::spans); + cls.def_readonly("too_many_bad_pixels", &SpanSetMoments::too_many_bad_pixels); + cls.def_readonly("center_out_of_bounds", &SpanSetMoments::center_out_of_bounds); + cls.def_readonly("bad_pixel_in_center", &SpanSetMoments::bad_pixel_in_center); + cls.def_readonly("singular_second_moments", &SpanSetMoments::singular_second_moments); + cls.def_property_readonly("any_flags_set", &SpanSetMoments::any_flags_set); + cls.def("get_x_array", &SpanSetMoments::get_x_array); + cls.def("get_y_array", &SpanSetMoments::get_y_array); + cls.def_static("compute", &SpanSetMoments::compute, "spans"_a, "masked_image"_a, + "bad_bitmask"_a, "bad_pixel_max_fraction"_a, "bad_pixel_exclusion_radius"_a); + cls.def_static("fit_shapelets", &SpanSetMoments::fit_shapelets, "masked_image"_a, "moments"_a, + "order"_a, "scale"_a); + }); +} +} // namespace algorithms +} // namespace meas +} // namespace lsst diff --git a/src/SpanSetMoments.cc b/src/SpanSetMoments.cc new file mode 100644 index 000000000..8b8a17bff --- /dev/null +++ b/src/SpanSetMoments.cc @@ -0,0 +1,172 @@ +// -*- LSST-C++ -*- + +/* + * LSST Data Management System + * Copyright 2008, 2009, 2010 LSST Corporation. + * + * This product includes software developed by the + * LSST Project (http://www.lsst.org/). + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the LSST License Statement and + * the GNU General Public License along with this program. If not, + * see . + */ + +#include +#include "Eigen/Dense" +#include "lsst/shapelet/constants.h" +#include "ndarray/eigen.h" +#include "lsst/meas/algorithms/SpanSetMoments.h" +#include +#include "lsst/shapelet/MatrixBuilder.h" + +namespace lsst { +namespace meas { +namespace algorithms { + +ndarray::Array SpanSetMoments::get_x_array() const { + ndarray::Array result = ndarray::allocate(spans->getArea()); + spans->applyFunctor([](geom::Point2I const& point, float& out) { out = point.getX(); }, + ndarray::ndFlat(result)); + return result; +} + +ndarray::Array SpanSetMoments::get_y_array() const { + ndarray::Array result = ndarray::allocate(spans->getArea()); + spans->applyFunctor([](geom::Point2I const& point, float& out) { out = point.getY(); }, + ndarray::ndFlat(result)); + return result; +} + +std::shared_ptr SpanSetMoments::compute(afw::geom::SpanSet const& spans, + afw::image::MaskedImage const& masked_image, + afw::image::MaskPixel bad_bitmask, + double bad_pixel_max_fraction, + double bad_pixel_exclusion_radius) { + double const nan = std::numeric_limits::quiet_NaN(); + std::shared_ptr result( + new SpanSetMoments{nan, nan, geom::Point2D(), afw::geom::ellipses::Quadrupole()}); + auto area = spans.getArea(); + auto bbox = spans.getBBox(); + auto bad_spans = afw::geom::SpanSet::fromMask((*masked_image.getMask())[bbox], bad_bitmask); + result->spans = spans.intersectNot(*bad_spans); + if ((area - result->spans->getArea()) > (bad_pixel_max_fraction * area)) { + result->too_many_bad_pixels = true; + return result; + } + // Compute the variance of the flux. + result->variance = 0.0; + result->spans->applyFunctor( + [&result](geom::Point2I const& point, float data) { result->variance += data; }, + *masked_image.getVariance()); + // Compute the flux (zeroth moment) and center (first moments). + result->flux = 0.0; + double qx = 0.0; + double qy = 0.0; + result->spans->applyFunctor( + [&result, &qx, &qy](geom::Point2I const& point, float data) { + result->flux += data; + qx += point.getX() * data; + qy += point.getY() * data; + }, + *masked_image.getImage()); + result->center.setX(qx / result->flux); + result->center.setY(qy / result->flux); + // Check to see if the center is plausible, and that it is not too close + // to a bad pixel. + if (!spans.contains(geom::Point2I(result->center))) { + result->center_out_of_bounds = true; + return result; + } + double min_bad_pixel_center_distance_squared = std::numeric_limits::infinity(); + for (auto const& span : *bad_spans) { + for (auto point : span) { + min_bad_pixel_center_distance_squared = + std::min(geom::Point2D(point).distanceSquared(result->center), + min_bad_pixel_center_distance_squared); + } + } + if (min_bad_pixel_center_distance_squared < bad_pixel_exclusion_radius * bad_pixel_exclusion_radius) { + result->bad_pixel_in_center = true; + return result; + } + // Compute the second moments. + double qxx = 0.0; + double qyy = 0.0; + double qxy = 0.0; + result->spans->applyFunctor( + [&result, &qxx, &qyy, &qxy](geom::Point2I const& point, float data) { + double dx = point.getX() - result->center.getX(); + double dy = point.getY() - result->center.getY(); + qxx += dx * dx * data; + qyy += dy * dy * data; + qxy += dx * dy * data; + }, + *masked_image.getImage()); + result->shape.setIxx(qxx / result->flux); + result->shape.setIyy(qyy / result->flux); + result->shape.setIxy(qxy / result->flux); + if (result->shape.getDeterminant() <= 0.0 || !std::isfinite(result->shape.getTraceRadius())) { + result->singular_second_moments = true; + return result; + } + return result; +} + +shapelet::ShapeletFunction SpanSetMoments::fit_shapelets( + afw::image::MaskedImage const& masked_image, + std::vector> const& moments, int order, double scale) { + std::size_t n_pixels = std::accumulate(moments.begin(), moments.end(), 0, + [](std::size_t n, std::shared_ptr const& moments) { + return n + moments->spans->getArea(); + }); + std::size_t n_coeff = shapelet::computeSize(order); + ndarray::Array data_array = ndarray::allocate(n_pixels); + data_array.deep() = 0.0; + // We allocate the design matrix via its transpose, because that yields + // the memory order expected by the shapelet library and the better + // memory order for the least-squares solver. + ndarray::Array design_matrix = + ndarray::Array(ndarray::allocate(n_coeff, n_pixels)).transpose(); + design_matrix.deep() = 0.0; + std::size_t start = 0; + for (auto const& source_moments : moments) { + shapelet::MatrixBuilder builder(source_moments->get_x_array(), source_moments->get_y_array(), order); + std::size_t stop = start + builder.getDataSize(); + auto source_data_array = data_array[ndarray::view(start, stop)].shallow(); + ndarray::Array source_weight_array = ndarray::allocate(builder.getDataSize()); + source_moments->spans->applyFunctor( + [](geom::Point2I const& point, float& data, float& weight, float img, float var) { + weight = 1.0 / std::sqrt(var); + data = img * weight; + }, + ndFlat(source_data_array), ndFlat(source_weight_array), *masked_image.getImage(), + *masked_image.getVariance()); + auto source_design_matrix = design_matrix[ndarray::view(start, stop)()].shallow(); + afw::geom::ellipses::Ellipse ellipse(source_moments->shape, source_moments->center); + ellipse.scale(scale); + builder(source_design_matrix, ellipse); + asEigenArray(source_design_matrix) *= source_moments->flux; + asEigenArray(source_design_matrix).colwise() *= asEigenArray(source_weight_array); + start = stop; + } + shapelet::ShapeletFunction result(order, shapelet::HERMITE); + asEigenMatrix(result.getCoefficients()) = asEigenMatrix(design_matrix) + .bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV) + .solve(asEigenMatrix(data_array)).cast(); + return result; +} + +} // namespace algorithms +} // namespace meas +} // namespace lsst diff --git a/ups/meas_algorithms.cfg b/ups/meas_algorithms.cfg index 2e2c1802d..ebc9fc583 100644 --- a/ups/meas_algorithms.cfg +++ b/ups/meas_algorithms.cfg @@ -4,7 +4,7 @@ import lsst.sconsUtils dependencies = { "required": ["astshim", "cpputils", "geom", "afw", "boost_math", "pex_config", "meas_base", - "minuit2"], + "minuit2", "shapelet"], "buildRequired": ["boost_test", "pybind11"], } diff --git a/ups/meas_algorithms.table b/ups/meas_algorithms.table index 5f574125e..57cd69619 100644 --- a/ups/meas_algorithms.table +++ b/ups/meas_algorithms.table @@ -12,6 +12,7 @@ setupRequired(utils) setupRequired(cpputils) setupRequired(kht) setupRequired(obs_base) +setupRequired(shapelet) envPrepend(LD_LIBRARY_PATH, ${PRODUCT_DIR}/lib) envPrepend(DYLD_LIBRARY_PATH, ${PRODUCT_DIR}/lib)