239def _fitPsf(fp, fmask, pk, pkF, pkres, fbb, peaks, peaksF, log, psf, psffwhm,
240 img, varimg, psfChisqCut1, psfChisqCut2, psfChisqCut2b,
243 r"""Fit a PSF + smooth background model (linear) to a small region
246 See fitPsfs for a more thorough description, including all
247 parameters not described below.
251 fp: `afw.detection.Footprint`
252 Footprint containing the Peaks to model.
253 fmask: `afw.image.Mask`
254 The Mask plane for pixels in the Footprint
255 pk: `afw.detection.PeakRecord`
256 The peak within the Footprint that we are going to fit with PSF model
257 pkF: `afw.geom.Point2D`
258 Floating point coordinates of the peak.
259 pkres: `meas.deblender.DeblendedPeak`
260 Peak results object that will hold the results.
261 fbb: `afw.geom.Box2I`
262 Bounding box of ``fp``
263 peaks: `afw.detection.PeakCatalog`
264 Catalog of peaks contained in the parent footprint.
265 peaksF: list of `afw.geom.Point2D`
266 List of floating point coordinates of all of the peaks.
267 psf: list of `afw.detection.Psf`\ s
268 Psf of the ``maskedImage`` for each band.
269 psffwhm: list pf `float`\ s
270 FWHM of the ``maskedImage``\ 's ``psf`` in each band.
271 img: `afw.image.ImageF`
272 The image that contains the footprint.
273 varimg: `afw.image.ImageF`
274 The variance of the image that contains the footprint.
279 Whether or not the peak matches a PSF model.
289 R0 = int(np.ceil(psffwhm*1.))
291 R1 = int(np.ceil(psffwhm*1.5))
292 cx, cy = pkF.getX(), pkF.getY()
293 psfimg = psf.computeImage(cx, cy)
295 R2 = R1 +
min(psfimg.getWidth(), psfimg.getHeight())/2.
297 pbb = psfimg.getBBox()
299 px0, py0 = psfimg.getX0(), psfimg.getY0()
304 pkres.setOutOfBounds()
308 xlo = int(np.floor(cx - R1))
309 ylo = int(np.floor(cy - R1))
310 xhi = int(np.ceil(cx + R1))
311 yhi = int(np.ceil(cy + R1))
314 xlo, xhi = stampbb.getMinX(), stampbb.getMaxX()
315 ylo, yhi = stampbb.getMinY(), stampbb.getMaxY()
316 if xlo > xhi
or ylo > yhi:
317 log.trace(
'Skipping this peak: out of bounds')
318 pkres.setOutOfBounds()
322 if min(stampbb.getWidth(), stampbb.getHeight()) <=
max(tinyFootprintSize, 2):
325 log.trace(
'Skipping this peak: tiny footprint / close to edge')
326 pkres.setTinyFootprint()
331 for pk2, pkF2
in zip(peaks, peaksF):
334 if pkF.distanceSquared(pkF2) > R2**2:
336 opsfimg = psf.computeImage(pkF2.getX(), pkF2.getY())
337 if not opsfimg.getBBox().overlaps(stampbb):
339 otherpeaks.append(opsfimg)
340 log.trace(
'%i other peaks within range', len(otherpeaks))
346 NT1 = 4 + len(otherpeaks)
350 NP = (1 + yhi - ylo)*(1 + xhi - xlo)
362 ix0, iy0 = img.getX0(), img.getY0()
363 fx0, fy0 = fbb.getMinX(), fbb.getMinY()
364 fslice = (slice(ylo-fy0, yhi-fy0+1), slice(xlo-fx0, xhi-fx0+1))
365 islice = (slice(ylo-iy0, yhi-iy0+1), slice(xlo-ix0, xhi-ix0+1))
366 fmask_sub = fmask .getArray()[fslice]
367 var_sub = varimg.getArray()[islice]
368 img_sub = img.getArray()[islice]
371 psfarr = psfimg.getArray()[pbb.getMinY()-py0: 1+pbb.getMaxY()-py0,
372 pbb.getMinX()-px0: 1+pbb.getMaxX()-px0]
373 px0, px1 = pbb.getMinX(), pbb.getMaxX()
374 py0, py1 = pbb.getMinY(), pbb.getMaxY()
377 valid = (fmask_sub > 0)
378 xx, yy = np.arange(xlo, xhi+1), np.arange(ylo, yhi+1)
379 RR = ((xx - cx)**2)[np.newaxis, :] + ((yy - cy)**2)[:, np.newaxis]
380 valid *= (RR <= R1**2)
381 valid *= (var_sub > 0)
385 log.warning(
'Skipping peak at (%.1f, %.1f): no unmasked pixels nearby', cx, cy)
386 pkres.setNoValidPixels()
390 XX, YY = np.meshgrid(xx, yy)
391 ipixes = np.vstack((XX[valid] - xlo, YY[valid] - ylo)).T
393 inpsfx = (xx >= px0)*(xx <= px1)
394 inpsfy = (yy >= py0)*(yy <= py1)
395 inpsf = np.outer(inpsfy, inpsfx)
396 indx = np.outer(inpsfy, (xx > px0)*(xx < px1))
397 indy = np.outer((yy > py0)*(yy < py1), inpsfx)
402 def _overlap(xlo, xhi, xmin, xmax):
403 assert (xlo <= xmax)
and (xhi >= xmin)
and (xlo <= xhi)
and (xmin <= xmax)
404 xloclamp =
max(xlo, xmin)
406 xhiclamp =
min(xhi, xmax)
407 Xhi = Xlo + (xhiclamp - xloclamp)
410 return (xloclamp, xhiclamp+1, Xlo, Xhi+1)
412 A = np.zeros((NP, NT2))
416 A[:, I_sky_ramp_x] = ipixes[:, 0] + (xlo-cx)
417 A[:, I_sky_ramp_y] = ipixes[:, 1] + (ylo-cy)
420 px0, px1 = pbb.getMinX(), pbb.getMaxX()
421 py0, py1 = pbb.getMinY(), pbb.getMaxY()
422 sx1, sx2, sx3, sx4 = _overlap(xlo, xhi, px0, px1)
423 sy1, sy2, sy3, sy4 = _overlap(ylo, yhi, py0, py1)
424 dpx0, dpy0 = px0 - xlo, py0 - ylo
425 psf_y_slice = slice(sy3 - dpy0, sy4 - dpy0)
426 psf_x_slice = slice(sx3 - dpx0, sx4 - dpx0)
427 psfsub = psfarr[psf_y_slice, psf_x_slice]
428 vsub = valid[sy1-ylo: sy2-ylo, sx1-xlo: sx2-xlo]
429 A[inpsf[valid], I_psf] = psfsub[vsub]
433 oldsx = (sx1, sx2, sx3, sx4)
434 sx1, sx2, sx3, sx4 = _overlap(xlo, xhi, px0+1, px1-1)
435 psfsub = (psfarr[psf_y_slice, sx3 - dpx0 + 1: sx4 - dpx0 + 1]
436 - psfarr[psf_y_slice, sx3 - dpx0 - 1: sx4 - dpx0 - 1])/2.
437 vsub = valid[sy1-ylo: sy2-ylo, sx1-xlo: sx2-xlo]
438 A[indx[valid], I_dx] = psfsub[vsub]
440 (sx1, sx2, sx3, sx4) = oldsx
443 sy1, sy2, sy3, sy4 = _overlap(ylo, yhi, py0+1, py1-1)
444 psfsub = (psfarr[sy3 - dpy0 + 1: sy4 - dpy0 + 1, psf_x_slice]
445 - psfarr[sy3 - dpy0 - 1: sy4 - dpy0 - 1, psf_x_slice])/2.
446 vsub = valid[sy1-ylo: sy2-ylo, sx1-xlo: sx2-xlo]
447 A[indy[valid], I_dy] = psfsub[vsub]
450 for j, opsf
in enumerate(otherpeaks):
452 ino = np.outer((yy >= obb.getMinY())*(yy <= obb.getMaxY()),
453 (xx >= obb.getMinX())*(xx <= obb.getMaxX()))
454 dpx0, dpy0 = obb.getMinX() - xlo, obb.getMinY() - ylo
455 sx1, sx2, sx3, sx4 = _overlap(xlo, xhi, obb.getMinX(), obb.getMaxX())
456 sy1, sy2, sy3, sy4 = _overlap(ylo, yhi, obb.getMinY(), obb.getMaxY())
457 opsfarr = opsf.getArray()
458 psfsub = opsfarr[sy3 - dpy0: sy4 - dpy0, sx3 - dpx0: sx4 - dpx0]
459 vsub = valid[sy1-ylo: sy2-ylo, sx1-xlo: sx2-xlo]
460 A[ino[valid], I_opsf + j] = psfsub[vsub]
466 rw = np.ones_like(RR)
469 rw[ii] = np.maximum(0, 1. - ((rr - R0)/(R1 - R0)))
470 w = np.sqrt(rw[valid]/var_sub[valid])
472 sumr = np.sum(rw[valid])
473 log.debug(
'sumr = %g', sumr)
477 Aw = A*w[:, np.newaxis]
486 im1 = np.zeros((1+yhi-ylo, 1+xhi-xlo))
487 im1[ipixes[:, 1], ipixes[:, 0]] = A[:, i]
488 plt.subplot(R, C, i+1)
489 plt.imshow(im1, interpolation=
'nearest', origin=
'lower')
490 plt.subplot(R, C, NT2+1)
491 im1 = np.zeros((1+yhi-ylo, 1+xhi-xlo))
492 im1[ipixes[:, 1], ipixes[:, 0]] = b
493 plt.imshow(im1, interpolation=
'nearest', origin=
'lower')
494 plt.subplot(R, C, NT2+2)
495 im1 = np.zeros((1+yhi-ylo, 1+xhi-xlo))
496 im1[ipixes[:, 1], ipixes[:, 0]] = w
497 plt.imshow(im1, interpolation=
'nearest', origin=
'lower')
509 X1, r1, rank1, s1 = np.linalg.lstsq(Aw[:, :NT1], bw, rcond=-1)
511 X2, r2, rank2, s2 = np.linalg.lstsq(Aw, bw, rcond=-1)
512 except np.linalg.LinAlgError
as e:
513 log.warning(
"Failed to fit PSF to child: %s", e)
514 pkres.setPsfFitFailed()
517 log.debug(
'r1 r2 %s %s', r1, r2)
529 dof1 = sumr - len(X1)
530 dof2 = sumr - len(X2)
531 log.debug(
'dof1, dof2 %g %g', dof1, dof2)
534 if dof1 <= 0
or dof2 <= 0:
535 log.trace(
'Skipping this peak: bad DOF %g, %g', dof1, dof2)
541 log.trace(
'PSF fits: chisq/dof = %g, %g', q1, q2)
542 ispsf1 = (q1 < psfChisqCut1)
543 ispsf2 = (q2 < psfChisqCut2)
545 pkres.psfFit1 = (chisq1, dof1)
546 pkres.psfFit2 = (chisq2, dof2)
550 fdx, fdy = X2[I_dx], X2[I_dy]
555 ispsf2 = ispsf2
and (abs(dx) < 1.
and abs(dy) < 1.)
556 log.trace(
'isPSF2 -- checking derivatives: dx,dy = %g, %g -> %s', dx, dy, str(ispsf2))
558 pkres.psfFitBigDecenter =
True
563 psfimg2 = psf.computeImage(cx + dx, cy + dy)
565 pbb2 = psfimg2.getBBox()
570 if not pbb2.contains(
geom.Point2I(int(cx + dx), int(cy + dy))):
574 px0, py0 = psfimg2.getX0(), psfimg2.getY0()
575 psfarr = psfimg2.getArray()[pbb2.getMinY()-py0:1+pbb2.getMaxY()-py0,
576 pbb2.getMinX()-px0:1+pbb2.getMaxX()-px0]
577 px0, py0 = pbb2.getMinX(), pbb2.getMinY()
578 px1, py1 = pbb2.getMaxX(), pbb2.getMaxY()
583 sx1, sx2, sx3, sx4 = _overlap(xlo, xhi, px0, px1)
584 sy1, sy2, sy3, sy4 = _overlap(ylo, yhi, py0, py1)
585 dpx0, dpy0 = px0 - xlo, py0 - ylo
586 psfsub = psfarr[sy3-dpy0:sy4-dpy0, sx3-dpx0:sx4-dpx0]
587 vsub = valid[sy1-ylo:sy2-ylo, sx1-xlo:sx2-xlo]
588 xx, yy = np.arange(xlo, xhi+1), np.arange(ylo, yhi+1)
589 inpsf = np.outer((yy >= py0)*(yy <= py1), (xx >= px0)*(xx <= px1))
590 Ab[inpsf[valid], I_psf] = psfsub[vsub]
592 Aw = Ab*w[:, np.newaxis]
594 Xb, rb, rankb, sb = np.linalg.lstsq(Aw, bw, rcond=-1)
599 dofb = sumr - len(Xb)
601 ispsf2 = (qb < psfChisqCut2b)
604 log.trace(
'shifted PSF: new chisq/dof = %g; good? %s', qb, ispsf2)
605 pkres.psfFit3 = (chisqb, dofb)
608 if (((ispsf1
and ispsf2)
and (q2 < q1))
609 or (ispsf2
and not ispsf1)):
613 log.debug(
'dof %g', dof)
614 log.trace(
'Keeping shifted-PSF model')
617 pkres.psfFitWithDecenter =
True
623 log.debug(
'dof %g', dof)
624 log.trace(
'Keeping unshifted PSF model')
626 ispsf = (ispsf1
or ispsf2)
630 SW, SH = 1+xhi-xlo, 1+yhi-ylo
631 psfmod = afwImage.ImageF(SW, SH)
632 psfmod.setXY0(xlo, ylo)
633 psfderivmodm = afwImage.MaskedImageF(SW, SH)
634 psfderivmod = psfderivmodm.getImage()
635 psfderivmod.setXY0(xlo, ylo)
636 model = afwImage.ImageF(SW, SH)
637 model.setXY0(xlo, ylo)
638 for i
in range(len(Xpsf)):
639 for (x, y), v
in zip(ipixes, A[:, i]*Xpsf[i]):
640 ix, iy = int(x), int(y)
641 model.set(ix, iy, model.get(ix, iy) + float(v))
642 if i
in [I_psf, I_dx, I_dy]:
643 psfderivmod.set(ix, iy, psfderivmod.get(ix, iy) + float(v))
646 psfmod.set(int(x), int(y), float(A[ii, I_psf]*Xpsf[I_psf]))
648 for (x, y)
in ipixes:
649 modelfp.addSpan(int(y+ylo), int(x+xlo), int(x+xlo))
652 pkres.psfFitDebugPsf0Img = psfimg
653 pkres.psfFitDebugPsfImg = psfmod
654 pkres.psfFitDebugPsfDerivImg = psfderivmod
655 pkres.psfFitDebugPsfModel = model
656 pkres.psfFitDebugStamp = img.Factory(img, stampbb,
True)
657 pkres.psfFitDebugValidPix = valid
658 pkres.psfFitDebugVar = varimg.Factory(varimg, stampbb,
True)
659 ww = np.zeros(valid.shape, np.float64)
661 pkres.psfFitDebugWeight = ww
662 pkres.psfFitDebugRampWeight = rw
667 pkres.psfFitStampExtent = (xlo, xhi, ylo, yhi)
668 pkres.psfFitCenter = (cx, cy)
669 log.debug(
'saving chisq,dof %g %g', chisq, dof)
670 pkres.psfFitBest = (chisq, dof)
671 pkres.psfFitParams = Xpsf
672 pkres.psfFitFlux = Xpsf[I_psf]
673 pkres.psfFitNOthers = len(otherpeaks)
676 pkres.setDeblendedAsPsf()
680 log.trace(
'Deblending as PSF; setting template to PSF model')
683 psfimg = psf.computeImage(cx, cy)
685 psfimg *= Xpsf[I_psf]
686 psfimg = psfimg.convertF()
690 psfbb = psfimg.getBBox()
692 bb = fpcopy.getBBox()
695 psfmod = afwImage.ImageF(bb)
696 fpcopy.spans.copyImage(psfimg, psfmod)
699 pkres.setTemplate(psfmod, fpcopy)
702 pkres.setPsfTemplate(psfmod, fpcopy)