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