134 """Run the Naive trailed source measurement algorithm.
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.
145 lsst.meas.base.SingleFramePlugin.measure
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)
161 Ixx, Iyy, Ixy = measRecord.getShape().getParameterVector()
166 a2 = 0.5 * (xpy + sqrt(xmy2 + 4.0*xy2))
167 b2 = 0.5 * (xpy - sqrt(xmy2 + 4.0*xy2))
171 if measRecord.get(
"base_SdssShape_flag_unweighted"):
172 self.
log.debug(
"Unweighted")
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)
184 theta = 0.5 * np.arctan2(2.0 * Ixy, xmy)
188 dydtheta = radius*np.cos(theta)
189 dxdtheta = radius*np.sin(theta)
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):
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)):
215 cutout = getMeasurementCutout(measRecord, exposure)
218 params = np.array([xc, yc, 1.0, length, theta])
219 model = VeresModel(cutout)
220 flux, gradFlux = model.computeFluxWithGradient(params)
223 if not np.isfinite(flux):
224 if np.isfinite(measRecord.getApInstFlux()):
225 flux = measRecord.getApInstFlux()
232 IxxErr2, IyyErr2, IxyErr2 = np.diag(measRecord.getShapeErr())
236 xcErr2, ycErr2 = np.diag(measRecord.getCentroidErr())
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)
249 dThetadIxx = -Ixy / (xmy2 + 4.0*xy2)
250 dThetadIxy = xmy / (xmy2 + 4.0*xy2)
251 thetaErr = sqrt(dThetadIxx*dThetadIxx*(IxxErr2 + IyyErr2) + dThetadIxy*dThetadIxy*IxyErr2)
254 dFdxc, dFdyc, _, dFdL, dFdTheta = gradFlux
255 fluxErr = sqrt(dFdL*dFdL*lengthErr*lengthErr + dFdTheta*dFdTheta*thetaErr*thetaErr
256 + dFdxc*dFdxc*xcErr2 + dFdyc*dFdyc*ycErr2)
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)
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)
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)