LSST Applications g2079a07aa2+86d27d4dc4,g2305ad1205+a659bff248,g2bbee38e9b+3c60f8fe34,g337abbeb29+3c60f8fe34,g33d1c0ed96+3c60f8fe34,g3502564af9+d77d6d1350,g3a166c0a6a+3c60f8fe34,g487adcacf7+25d9892218,g4be5004598+d77d6d1350,g50ff169b8f+96c6868917,g52b1c1532d+585e252eca,g591dd9f2cf+4d81263f9a,g5cd07815a0+980d2b1c3b,g607f77f49a+d77d6d1350,g858d7b2824+d77d6d1350,g88963caddf+83e433e629,g99cad8db69+a4d3c48eeb,g9ddcbc5298+9a081db1e4,ga1e77700b3+bcf1af89ad,ga57fefb910+9a39d7b2d7,gae0086650b+585e252eca,gb065fddaf9+4f9fd82a2c,gb0e22166c9+60f28cb32d,gb363559e06+d84b1d3d07,gb3b7280ab2+4563d032e1,gb4b16eec92+babe958938,gba4ed39666+c2a2e4ac27,gbb8dafda3b+ed6854b564,gc120e1dc64+b72d212f87,gc28159a63d+3c60f8fe34,gc3e9b769f7+921dbcd359,gcf0d15dbbd+9a39d7b2d7,gdaeeff99f8+f9a426f77a,gddc38dedce+585e252eca,ge79ae78c31+3c60f8fe34,w.2024.21
LSST Data Management Base Package
Loading...
Searching...
No Matches
_compensatedGaussian.py
Go to the documentation of this file.
1# This file is part of meas_base.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
21
22from __future__ import annotations
23
24__all__ = (
25 "SingleFrameCompensatedGaussianFluxConfig",
26 "SingleFrameCompensatedGaussianFluxPlugin",
27)
28
29import math
30import numpy as np
31import scipy.stats as sps
32
33from lsst.pex.config import Field, ListField
34from lsst.geom import Point2I
35
36from ..sfm import SingleFramePlugin, SingleFramePluginConfig
37from ..pluginRegistry import register
38
39from .._measBaseLib import _compensatedGaussianFiltInnerProduct, FlagHandler, FlagDefinitionList
40
41
43 kernel_widths = ListField(
44 doc="The widths (in pixels) of the kernels for which to measure compensated apertures.",
45 dtype=int,
46 minLength=1,
47 default=[3, 5]
48 )
49
50 t = Field(
51 doc="Scale parameter of outer Gaussian compared to inner Gaussian.",
52 dtype=float,
53 default=2.0,
54 )
55
56
57@register("base_CompensatedGaussianFlux")
59 ConfigClass = SingleFrameCompensatedGaussianFluxConfig
60
61 @classmethod
63 return cls.FLUX_ORDER
64
66 self,
67 config: SingleFrameCompensatedGaussianFluxConfig,
68 name: str,
69 schema,
70 metadata,
71 logName=None,
72 **kwds,
73 ):
74 super().__init__(config, name, schema, metadata, logName, **kwds)
75
76 flagDefs = FlagDefinitionList()
77
78 self.width_keys = {}
79 self._rads = {}
82 self._t = config.t
83 for width in config.kernel_widths:
84 base_key = f"{name}_{width}"
85
86 # flux
87 flux_str = f"{base_key}_instFlux"
88 flux_key = schema.addField(
89 flux_str,
90 type="D",
91 doc="Compensated Gaussian flux measurement.",
92 units="count",
93 )
94
95 # flux error
96 err_str = f"{base_key}_instFluxErr"
97 err_key = schema.addField(
98 err_str,
99 type="D",
100 doc="Compensated Gaussian flux error.",
101 units="count",
102 )
103
104 # mask bits
105 mask_str = f"{base_key}_mask_bits"
106 mask_key = schema.addField(mask_str, type=np.int32, doc="Mask bits set within aperture.")
107
108 # failure flags
109 failure_flag = flagDefs.add(f"{width}_flag", "Compensated Gaussian measurement failed")
110 oob_flag = flagDefs.add(f"{width}_flag_bounds", "Compensated Gaussian out-of-bounds")
111
112 self.width_keys[width] = (flux_key, err_key, mask_key, failure_flag, oob_flag)
113 self._rads[width] = math.ceil(sps.norm.ppf((0.995,), scale=width * config.t)[0])
114
115 self.flagHandler = FlagHandler.addFields(schema, name, flagDefs)
116 self._max_rad = max(self._rads)
117
118 def fail(self, measRecord, error=None):
119 """Record failure
120
121 See also
122 --------
123 lsst.meas.base.SingleFramePlugin.fail
124 """
125 if error is None:
126 self.flagHandler.handleFailure(measRecord)
127 else:
128 self.flagHandler.handleFailure(measRecord, error.cpp)
129
130 def measure(self, measRecord, exposure):
131 center = measRecord.getCentroid()
132 bbox = exposure.getBBox()
133
134 y = center.getY() - bbox.beginY
135 x = center.getX() - bbox.beginX
136
137 y_floor = math.floor(y)
138 x_floor = math.floor(x)
139
140 for width, (flux_key, err_key, mask_key, failure_flag, oob_flag) in self.width_keys.items():
141 rad = self._rads[width]
142
143 if Point2I(center) not in exposure.getBBox().erodedBy(rad):
144 self.flagHandler.setValue(measRecord, failure_flag.number, True)
145 self.flagHandler.setValue(measRecord, oob_flag.number, True)
146 continue
147
148 y_slice = slice(y_floor - rad, y_floor + rad + 1, 1)
149 x_slice = slice(x_floor - rad, x_floor + rad + 1, 1)
150 y_mean = y - y_floor + rad
151 x_mean = x - x_floor + rad
152
153 sub_im = exposure.image.array[y_slice, x_slice]
154 sub_var = exposure.variance.array[y_slice, x_slice]
155
156 if sub_im.size == 0 or sub_im.shape[0] != sub_im.shape[1] or (sub_im.shape[0] % 2) == 0:
157 self.flagHandler.setValue(measRecord, failure_flag.number, True)
158 self.flagHandler.setValue(measRecord, oob_flag.number, True)
159 continue
160
161 flux, var = _compensatedGaussianFiltInnerProduct(
162 sub_im,
163 sub_var,
164 x_mean,
165 y_mean,
166 width,
167 self._t,
168 )
169
170 measRecord.set(flux_key, flux)
171 measRecord.set(err_key, np.sqrt(var))
172 measRecord.set(mask_key, np.bitwise_or.reduce(exposure.mask.array[y_slice, x_slice], axis=None))
std::vector< SchemaItem< Flag > > * items
int max
__init__(self, SingleFrameCompensatedGaussianFluxConfig config, str name, schema, metadata, logName=None, **kwds)