242 ):
243 r"""Fit a PSF + smooth background model (linear) to a small region
244 around a peak.
245
246 See fitPsfs for a more thorough description, including all
247 parameters not described below.
248
249 Parameters
250 ----------
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.
275
276 Results
277 -------
278 ispsf: `bool`
279 Whether or not the peak matches a PSF model.
280 """
281 import lsstDebug
282
283
286
287
288
289 R0 = int(np.ceil(psffwhm*1.))
290
291 R1 = int(np.ceil(psffwhm*1.5))
292 cx, cy = pkF.getX(), pkF.getY()
293 psfimg = psf.computeImage(cx, cy)
294
295 R2 = R1 +
min(psfimg.getWidth(), psfimg.getHeight())/2.
296
297 pbb = psfimg.getBBox()
298 pbb.clip(fbb)
299 px0, py0 = psfimg.getX0(), psfimg.getY0()
300
301
302
304 pkres.setOutOfBounds()
305 return
306
307
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))
313 stampbb.clip(fbb)
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()
319 return
320
321
322 if min(stampbb.getWidth(), stampbb.getHeight()) <=
max(tinyFootprintSize, 2):
323
324
325 log.trace('Skipping this peak: tiny footprint / close to edge')
326 pkres.setTinyFootprint()
327 return
328
329
330 otherpeaks = []
331 for pk2, pkF2 in zip(peaks, peaksF):
332 if pk2 == pk:
333 continue
334 if pkF.distanceSquared(pkF2) > R2**2:
335 continue
336 opsfimg = psf.computeImage(pkF2.getX(), pkF2.getY())
337 if not opsfimg.getBBox().overlaps(stampbb):
338 continue
339 otherpeaks.append(opsfimg)
340 log.trace('%i other peaks within range', len(otherpeaks))
341
342
343
344
345
346 NT1 = 4 + len(otherpeaks)
347
348 NT2 = NT1 + 2
349
350 NP = (1 + yhi - ylo)*(1 + xhi - xlo)
351
352 I_psf = 0
353 I_sky = 1
354 I_sky_ramp_x = 2
355 I_sky_ramp_y = 3
356
357 I_opsf = 4
358 I_dx = NT1 + 0
359 I_dy = NT1 + 1
360
361
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]
369
370
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()
375
376
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)
382 NP = valid.sum()
383
384 if NP == 0:
385 log.warning('Skipping peak at (%.1f, %.1f): no unmasked pixels nearby', cx, cy)
386 pkres.setNoValidPixels()
387 return
388
389
390 XX, YY = np.meshgrid(xx, yy)
391 ipixes = np.vstack((XX[valid] - xlo, YY[valid] - ylo)).T
392
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)
398
399 del inpsfx
400 del inpsfy
401
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)
405 Xlo = xloclamp - xlo
406 xhiclamp =
min(xhi, xmax)
407 Xhi = Xlo + (xhiclamp - xloclamp)
408 assert xloclamp >= 0
409 assert Xlo >= 0
410 return (xloclamp, xhiclamp+1, Xlo, Xhi+1)
411
412 A = np.zeros((NP, NT2))
413
414 A[:, I_sky] = 1.
415
416 A[:, I_sky_ramp_x] = ipixes[:, 0] + (xlo-cx)
417 A[:, I_sky_ramp_y] = ipixes[:, 1] + (ylo-cy)
418
419
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]
430
431
432
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]
439
440 (sx1, sx2, sx3, sx4) = oldsx
441
442
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]
448
449
450 for j, opsf in enumerate(otherpeaks):
451 obb = opsf.getBBox()
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]
461
462 b = img_sub[valid]
463
464
465
466 rw = np.ones_like(RR)
467 ii = (RR > R0**2)
468 rr = np.sqrt(RR[ii])
469 rw[ii] = np.maximum(0, 1. - ((rr - R0)/(R1 - R0)))
470 w = np.sqrt(rw[valid]/var_sub[valid])
471
472 sumr = np.sum(rw[valid])
473 log.debug('sumr = %g', sumr)
474
475 del ii
476
477 Aw = A*w[:, np.newaxis]
478 bw = b*w
479
480 if debugPlots:
481 import pylab as plt
482 plt.clf()
483 N = NT2 + 2
484 R, C = 2, (N+1)/2
485 for i in range(NT2):
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')
498 plt.savefig('A.png')
499
500
501
502
503
504
505
506 try:
507
508
509 X1, r1, rank1, s1 = np.linalg.lstsq(Aw[:, :NT1], bw, rcond=-1)
510
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()
515 return
516
517 log.debug('r1 r2 %s %s', r1, r2)
518
519
520
521 if len(r1) > 0:
522 chisq1 = r1[0]
523 else:
524 chisq1 = 1e30
525 if len(r2) > 0:
526 chisq2 = r2[0]
527 else:
528 chisq2 = 1e30
529 dof1 = sumr - len(X1)
530 dof2 = sumr - len(X2)
531 log.debug('dof1, dof2 %g %g', dof1, dof2)
532
533
534 if dof1 <= 0 or dof2 <= 0:
535 log.trace('Skipping this peak: bad DOF %g, %g', dof1, dof2)
536 pkres.setBadPsfDof()
537 return
538
539 q1 = chisq1/dof1
540 q2 = chisq2/dof2
541 log.trace('PSF fits: chisq/dof = %g, %g', q1, q2)
542 ispsf1 = (q1 < psfChisqCut1)
543 ispsf2 = (q2 < psfChisqCut2)
544
545 pkres.psfFit1 = (chisq1, dof1)
546 pkres.psfFit2 = (chisq2, dof2)
547
548
549 if ispsf2:
550 fdx, fdy = X2[I_dx], X2[I_dy]
551 f0 = X2[I_psf]
552
553 dx = fdx/f0
554 dy = fdy/f0
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))
557 if not ispsf2:
558 pkres.psfFitBigDecenter = True
559
560
561
562 if ispsf2:
563 psfimg2 = psf.computeImage(cx + dx, cy + dy)
564
565 pbb2 = psfimg2.getBBox()
566 pbb2.clip(fbb)
567
568
569
570 if not pbb2.contains(
geom.Point2I(int(cx + dx), int(cy + dy))):
571 ispsf2 = False
572 else:
573
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()
579
580
581 Ab = A[:, :NT1]
582
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]
591
592 Aw = Ab*w[:, np.newaxis]
593
594 Xb, rb, rankb, sb = np.linalg.lstsq(Aw, bw, rcond=-1)
595 if len(rb) > 0:
596 chisqb = rb[0]
597 else:
598 chisqb = 1e30
599 dofb = sumr - len(Xb)
600 qb = chisqb/dofb
601 ispsf2 = (qb < psfChisqCut2b)
602 q2 = qb
603 X2 = Xb
604 log.trace('shifted PSF: new chisq/dof = %g; good? %s', qb, ispsf2)
605 pkres.psfFit3 = (chisqb, dofb)
606
607
608 if (((ispsf1 and ispsf2) and (q2 < q1))
609 or (ispsf2 and not ispsf1)):
610 Xpsf = X2
611 chisq = chisq2
612 dof = dof2
613 log.debug('dof %g', dof)
614 log.trace('Keeping shifted-PSF model')
615 cx += dx
616 cy += dy
617 pkres.psfFitWithDecenter = True
618 else:
619
620 Xpsf = X1
621 chisq = chisq1
622 dof = dof1
623 log.debug('dof %g', dof)
624 log.trace('Keeping unshifted PSF model')
625
626 ispsf = (ispsf1 or ispsf2)
627
628
629 if debugPsf:
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))
644 for ii in range(NP):
645 x, y = ipixes[ii, :]
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))
650 modelfp.normalize()
651
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)
660 ww[valid] = w
661 pkres.psfFitDebugWeight = ww
662 pkres.psfFitDebugRampWeight = rw
663
664
665 pkres.psfFitR0 = R0
666 pkres.psfFitR1 = R1
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)
674
675 if ispsf:
676 pkres.setDeblendedAsPsf()
677
678
679
680 log.trace('Deblending as PSF; setting template to PSF model')
681
682
683 psfimg = psf.computeImage(cx, cy)
684
685 psfimg *= Xpsf[I_psf]
686 psfimg = psfimg.convertF()
687
688
690 psfbb = psfimg.getBBox()
691 fpcopy.clipTo(psfbb)
692 bb = fpcopy.getBBox()
693
694
695 psfmod = afwImage.ImageF(bb)
696 fpcopy.spans.copyImage(psfimg, psfmod)
697
698 clipFootprintToNonzeroImpl(fpcopy, psfmod)
699 pkres.setTemplate(psfmod, fpcopy)
700
701
702 pkres.setPsfTemplate(psfmod, fpcopy)
703
704 return ispsf
705
706
An integer coordinate rectangle.