133 def measure(self, measRecord, exposure):
134 """Run the Naive trailed source measurement algorithm.
135
136 Parameters
137 ----------
138 measRecord : `lsst.afw.table.SourceRecord`
139 Record describing the object being measured.
140 exposure : `lsst.afw.image.Exposure`
141 Pixel data to be measured.
142
143 See also
144 --------
145 lsst.meas.base.SingleFramePlugin.measure
146 """
147
148
149
150 xc = measRecord.get("base_SdssShape_x")
151 yc = measRecord.get("base_SdssShape_y")
152 if not np.isfinite(xc) or not np.isfinite(yc):
153 xc, yc = self.centroidExtractor(measRecord, self.flagHandler)
154 self.flagHandler.setValue(measRecord, self.SAFE_CENTROID.number, True)
155 self.flagHandler.setValue(measRecord, self.FAILURE.number, True)
156 return
157
158 ra, dec = self.computeRaDec(exposure, xc, yc)
159
160
161 Ixx, Iyy, Ixy = measRecord.getShape().getParameterVector()
162 xmy = Ixx - Iyy
163 xpy = Ixx + Iyy
164 xmy2 = xmy*xmy
165 xy2 = Ixy*Ixy
166 a2 = 0.5 * (xpy + sqrt(xmy2 + 4.0*xy2))
167 b2 = 0.5 * (xpy - sqrt(xmy2 + 4.0*xy2))
168
169
170
171 if measRecord.get("base_SdssShape_flag_unweighted"):
172 self.log.debug("Unweighted")
173 length, gradLength = self.computeLength(a2, b2)
174 else:
175 self.log.debug("Weighted")
176 length, gradLength, results = self.findLength(a2, b2)
177 if not results.converged:
178 self.log.info("Results not converged: %s", results.flag)
179 self.flagHandler.setValue(measRecord, self.NO_CONVERGE.number, True)
180 self.flagHandler.setValue(measRecord, self.FAILURE.number, True)
181 return
182
183
184 theta = 0.5 * np.arctan2(2.0 * Ixy, xmy)
185
186
187 radius = length/2.0
188 dydtheta = radius*np.cos(theta)
189 dxdtheta = radius*np.sin(theta)
190 x0 = xc - dydtheta
191 y0 = yc - dxdtheta
192 x1 = xc + dydtheta
193 y1 = yc + dxdtheta
194
195
196 if not (exposure.getBBox().beginX <= x0 <= exposure.getBBox().endX
197 and exposure.getBBox().beginX <= x1 <= exposure.getBBox().endX
198 and exposure.getBBox().beginY <= y0 <= exposure.getBBox().endY
199 and exposure.getBBox().beginY <= y1 <= exposure.getBBox().endY):
200
201 self.flagHandler.setValue(measRecord, self.EDGE.number, True)
202
203 else:
204
205
206
207 if exposure.mask[Point2I(int(x0), int(y0))] and exposure.mask[Point2I(int(x1), int(y1))]:
208 if ((exposure.mask[Point2I(int(x0), int(y0))] & exposure.mask.getPlaneBitMask('EDGE') != 0)
209 or (exposure.mask[Point2I(int(x1), int(y1))]
210 & exposure.mask.getPlaneBitMask('EDGE') != 0)):
211
212 self.flagHandler.setValue(measRecord, self.EDGE.number, True)
213
214
215 cutout = getMeasurementCutout(measRecord, exposure)
216
217
218 params = np.array([xc, yc, 1.0, length, theta])
219 model = VeresModel(cutout)
220 flux, gradFlux = model.computeFluxWithGradient(params)
221
222
223 if not np.isfinite(flux):
224 if np.isfinite(measRecord.getApInstFlux()):
225 flux = measRecord.getApInstFlux()
226 else:
227 self.flagHandler.setValue(measRecord, self.NO_FLUX.number, True)
228 self.flagHandler.setValue(measRecord, self.FAILURE.number, True)
229 return
230
231
232 IxxErr2, IyyErr2, IxyErr2 = np.diag(measRecord.getShapeErr())
233
234
235
236 xcErr2, ycErr2 = np.diag(measRecord.getCentroidErr())
237
238
239 desc = sqrt(xmy2 + 4.0*xy2)
240 da2dIxx = 0.5*(1.0 + (xmy/desc))
241 da2dIyy = 0.5*(1.0 - (xmy/desc))
242 da2dIxy = 2.0*Ixy / desc
243 a2Err2 = IxxErr2*da2dIxx*da2dIxx + IyyErr2*da2dIyy*da2dIyy + IxyErr2*da2dIxy*da2dIxy
244 b2Err2 = IxxErr2*da2dIyy*da2dIyy + IyyErr2*da2dIxx*da2dIxx + IxyErr2*da2dIxy*da2dIxy
245 dLda2, dLdb2 = gradLength
246 lengthErr = np.sqrt(dLda2*dLda2*a2Err2 + dLdb2*dLdb2*b2Err2)
247
248
249 dThetadIxx = -Ixy / (xmy2 + 4.0*xy2)
250 dThetadIxy = xmy / (xmy2 + 4.0*xy2)
251 thetaErr = sqrt(dThetadIxx*dThetadIxx*(IxxErr2 + IyyErr2) + dThetadIxy*dThetadIxy*IxyErr2)
252
253
254 dFdxc, dFdyc, _, dFdL, dFdTheta = gradFlux
255 fluxErr = sqrt(dFdL*dFdL*lengthErr*lengthErr + dFdTheta*dFdTheta*thetaErr*thetaErr
256 + dFdxc*dFdxc*xcErr2 + dFdyc*dFdyc*ycErr2)
257
258
259 dxdradius = np.cos(theta)
260 dydradius = np.sin(theta)
261 radiusErr2 = lengthErr*lengthErr/4.0
262 xErr2 = sqrt(xcErr2 + radiusErr2*dxdradius*dxdradius + thetaErr*thetaErr*dxdtheta*dxdtheta)
263 yErr2 = sqrt(ycErr2 + radiusErr2*dydradius*dydradius + thetaErr*thetaErr*dydtheta*dydtheta)
264 x0Err = sqrt(xErr2)
265 y0Err = sqrt(yErr2)
266
267
268 measRecord.set(self.keyRa, ra)
269 measRecord.set(self.keyDec, dec)
270 measRecord.set(self.keyX0, x0)
271 measRecord.set(self.keyY0, y0)
272 measRecord.set(self.keyX1, x1)
273 measRecord.set(self.keyY1, y1)
274 measRecord.set(self.keyFlux, flux)
275 measRecord.set(self.keyLength, length)
276 measRecord.set(self.keyAngle, theta)
277 measRecord.set(self.keyX0Err, x0Err)
278 measRecord.set(self.keyY0Err, y0Err)
279 measRecord.set(self.keyX1Err, x0Err)
280 measRecord.set(self.keyY1Err, y0Err)
281 measRecord.set(self.keyFluxErr, fluxErr)
282 measRecord.set(self.keyLengthErr, lengthErr)
283 measRecord.set(self.keyAngleErr, thetaErr)
284