LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
BaselineUtils.cc
Go to the documentation of this file.
1 #include <list>
2 #include <cmath>
3 #include <cstdint>
4 
5 #include "lsst/log/Log.h"
6 #include "lsst/geom.h"
8 #include "lsst/pex/exceptions.h"
9 #include "lsst/afw/geom/Span.h"
10 
11 using std::lround;
12 
13 namespace image = lsst::afw::image;
14 namespace det = lsst::afw::detection;
15 namespace deblend = lsst::meas::deblender;
16 namespace afwGeom = lsst::afw::geom;
17 namespace geom = lsst::geom;
18 
19 template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
21 
22 template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
23 const int
25 
26 template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
28 
29 template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
31 
32 template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
34 
35 template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
37 
38 static bool span_compare(afwGeom::Span const & sp1,
39  afwGeom::Span const & sp2) {
40  return (sp1 < sp2);
41 }
42 
43 namespace {
44  void nearestFootprint(std::vector<std::shared_ptr<det::Footprint>> const& foots,
47  {
48  /*
49  * insert the footprints into an image, set all the pixels to the
50  * Manhattan distance from the nearest set pixel.
51  */
52  typedef std::uint16_t dtype;
53  typedef std::uint16_t itype;
54 
55  const itype nil = 0xffff;
56 
57  *argmin = 0;
58  *dist = 0;
59 
60  {
61  int i = 0;
62  for (std::shared_ptr<det::Footprint> const & foot : foots) {
63  foot->getSpans()->setImage(*argmin, static_cast<itype>(i));
64  foot->getSpans()->setImage(*dist, static_cast<dtype>(1));
65  ++i;
66  }
67  }
68 
69  int const height = dist->getHeight();
70  int const width = dist->getWidth();
71 
72  // traverse from bottom left to top right
73  for (int y = 0; y != height; ++y) {
74  image::Image<dtype>::xy_locator dim = dist->xy_at(0, y);
75  image::Image<itype>::xy_locator aim = argmin->xy_at(0, y);
76  for (int x = 0; x != width; ++x, ++dim.x(), ++aim.x()) {
77  if (dim(0, 0) == 1) {
78  // first pass and pixel was on, it gets a zero
79  dim(0, 0) = 0;
80  // its argmin is already set
81  } else {
82  // pixel was off. It is at most the sum of lengths of
83  // the array away from a pixel that is on
84  dim(0, 0) = width + height;
85  aim(0, 0) = nil;
86  // or one more than the pixel to the north
87  if (y > 0) {
88  dtype ndist = dim(0,-1) + 1;
89  if (ndist < dim(0,0)) {
90  dim(0,0) = ndist;
91  aim(0,0) = aim(0,-1);
92  }
93  }
94  // or one more than the pixel to the west
95  if (x > 0) {
96  dtype ndist = dim(-1,0) + 1;
97  if (ndist < dim(0,0)) {
98  dim(0,0) = ndist;
99  aim(0,0) = aim(-1,0);
100  }
101  }
102  }
103  }
104  }
105  // traverse from top right to bottom left
106  for (int y = height - 1; y >= 0; --y) {
107  image::Image<dtype>::xy_locator dim = dist->xy_at(width-1, y);
108  image::Image<itype>::xy_locator aim = argmin->xy_at(width-1, y);
109  for (int x = width - 1; x >= 0; --x, --dim.x(), --aim.x()) {
110  // either what we had on the first pass or one more than
111  // the pixel to the south
112  if (y + 1 < height) {
113  dtype ndist = dim(0,1) + 1;
114  if (ndist < dim(0,0)) {
115  dim(0,0) = ndist;
116  aim(0,0) = aim(0,1);
117  }
118  }
119  // or one more than the pixel to the east
120  if (x + 1 < width) {
121  dtype ndist = dim(1,0) + 1;
122  if (ndist < dim(0,0)) {
123  dim(0,0) = ndist;
124  aim(0,0) = aim(1,0);
125  }
126  }
127  }
128  }
129  }
130 } // end anonymous namespace
131 
145 template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
146 void
148 medianFilter(ImageT const& img,
149  ImageT & out,
150  int halfsize) {
151  int S = halfsize*2 + 1;
152  int SS = S*S;
153  typedef typename ImageT::xy_locator xy_loc;
154  xy_loc pix = img.xy_at(halfsize,halfsize);
156  for (int i=0; i<S; ++i) {
157  for (int j=0; j<S; ++j) {
158  locs.push_back(pix.cache_location(j-halfsize, i-halfsize));
159  }
160  }
161  int W = img.getWidth();
162  int H = img.getHeight();
163  ImagePixelT vals[S*S];
164  for (int y=halfsize; y<H-halfsize; ++y) {
165  xy_loc inpix = img.xy_at(halfsize, y), end = img.xy_at(W-halfsize, y);
166  for (typename ImageT::x_iterator optr = out.row_begin(y) + halfsize;
167  inpix != end; ++inpix.x(), ++optr) {
168  for (int i=0; i<SS; ++i)
169  vals[i] = inpix[locs[i]];
170  std::nth_element(vals, vals+SS/2, vals+SS);
171  *optr = vals[SS/2];
172  }
173  }
174 
175  // grumble grumble margins
176  for (int y=0; y<2*halfsize; ++y) {
177  int iy = y;
178  if (y >= halfsize)
179  iy = H - 1 - (y-halfsize);
180  typename ImageT::x_iterator optr = out.row_begin(iy);
181  typename ImageT::x_iterator iptr = img.row_begin(iy), end=img.row_end(iy);
182  for (; iptr != end; ++iptr,++optr)
183  *optr = *iptr;
184  }
185  for (int y=halfsize; y<H-halfsize; ++y) {
186  typename ImageT::x_iterator optr = out.row_begin(y);
187  typename ImageT::x_iterator iptr = img.row_begin(y), end=img.row_begin(y)+halfsize;
188  for (; iptr != end; ++iptr,++optr)
189  *optr = *iptr;
190  iptr = img.row_begin(y) + ((W-1) - halfsize);
191  end = img.row_begin(y) + (W-1);
192  optr = out.row_begin(y) + ((W-1) - halfsize);
193  for (; iptr != end; ++iptr,++optr)
194  *optr = *iptr;
195  }
196 
197 }
198 
223 template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
224 void
227  ImageT & img,
228  det::PeakRecord const& peak) {
229 
230  int cx = peak.getIx();
231  int cy = peak.getIy();
232  int ix0 = img.getX0();
233  int iy0 = img.getY0();
234  int iW = img.getWidth();
235  int iH = img.getHeight();
236 
237  ImagePtrT shadowingImg = ImagePtrT(new ImageT(img, true));
238 
239  int DW = std::max(cx - img.getX0(), img.getX0() + img.getWidth() - cx);
240  int DH = std::max(cy - img.getY0(), img.getY0() + img.getHeight() - cy);
241 
242  const int S = 5;
243 
244  // Work out from the peak in chunks of "S" pixels.
245  int s;
246  for (s = 0; s < std::max(DW,DH); s += S) {
247  int p;
248  for (p=0; p<S; p++) {
249  // visit pixels with L_inf distance = s + p from the
250  // center (ie, the s+p'th square ring of pixels)
251  // L is the half-length of the ring (box).
252  int L = s+p;
253  int x = L, y = -L;
254  int dx = 0, dy = 0; // initialized here to satisfy the
255  // compiler; initialized for real
256  // below (first time through loop)
257  /*
258  int i;
259  int leg;
260  for (i=0; i<(8*L); i++, x += dx, y += dy) {
261  if (i % (2*L) == 0) {
262  leg = (i/(2*L));
263  dx = ( leg % 2) * (-1 + 2*(leg/2));
264  dy = ((leg+1) % 2) * ( 1 - 2*(leg/2));
265  }
266  */
267 
268  /*
269  We visit pixels in a box of "radius" L, in this order:
270 
271  L=1:
272 
273  4 3 2
274  5 1
275  6 7 0
276 
277  L=2:
278 
279  8 7 6 5 4
280  9 3
281  10 2
282  11 1
283  12 13 14 15 0
284 
285  Note that the number of pixel visited is 8*L, and that we
286  change "dx" or "dy" each "2*L" steps.
287  */
288  for (int i=0; i<(8*L); i++, x += dx, y += dy) {
289  // time to change directions? (Note that this runs
290  // the first time through the loop, initializing dx,dy
291  // appropriately.)
292  if (i % (2*L) == 0) {
293  int leg = (i/(2*L));
294  // dx = [ 0, -1, 0, 1 ][leg]
295  dx = ( leg % 2) * (-1 + 2*(leg/2));
296  // dy = [ 1, 0, -1, 0 ][leg]
297  dy = ((leg+1) % 2) * ( 1 - 2*(leg/2));
298  }
299  //printf(" i=%i, leg=%i, dx=%i, dy=%i, x=%i, y=%i\n", i, leg, dx, dy, x, y);
300  int px = cx + x - ix0;
301  int py = cy + y - iy0;
302  // If the shadowing pixel is out of bounds, nothing to do.
303  if (px < 0 || px >= iW || py < 0 || py >= iH)
304  continue;
305  // The pixel casting the shadow
306  ImagePixelT pix = (*shadowingImg)(px,py);
307 
308  // Cast this pixel's shadow S pixels long in a cone.
309  // We compute the range of slopes (or inverse-slopes)
310  // shadowed by the pixel, [ds0,ds1]
311  double ds0, ds1;
312  // Range of slopes shadowed
313  const double A = 0.3;
314  int shx, shy;
315  int psx, psy;
316  // Are we traversing a vertical edge of the box?
317  if (dx == 0) {
318  // (if so, then "x" is +- L, so no div-by-zero)
319  ds0 = (double(y) / double(x)) - A;
320  ds1 = ds0 + 2.0 * A;
321  // cast the shadow on column x + sign(x)*shx
322  for (shx=1; shx<=S; shx++) {
323  int xsign = (x>0?1:-1);
324  // the column being shadowed
325  psx = cx + x + (xsign*shx) - ix0;
326  if (psx < 0 || psx >= iW)
327  continue;
328  // shadow covers a range of y values based on slope
329  for (shy = lround(shx * ds0);
330  shy <= lround(shx * ds1); shy++) {
331  psy = cy + y + xsign*shy - iy0;
332  if (psy < 0 || psy >= iH)
333  continue;
334  img(psx, psy) = std::min(img(psx, psy), pix);
335  }
336  }
337 
338  } else {
339  // We're traversing a horizontal edge of the box; y = +-L
340  ds0 = (double(x) / double(y)) - A;
341  ds1 = ds0 + 2.0 * A;
342  // Cast shadow on row y + sign(y)*shy
343  for (shy=1; shy<=S; shy++) {
344  int ysign = (y>0?1:-1);
345  psy = cy + y + (ysign*shy) - iy0;
346  if (psy < 0 || psy >= iH)
347  continue;
348  // shadow covers a range of x vals based on slope
349  for (shx = lround(shy * ds0);
350  shx <= lround(shy * ds1); shx++) {
351  psx = cx + x + ysign*shx - ix0;
352  if (psx < 0 || psx >= iW)
353  continue;
354  img(psx, psy) = std::min(img(psx, psy), pix);
355  }
356  }
357  }
358  }
359  }
360  shadowingImg->assign(img);
361  }
362 }
363 
364 static double _get_contrib_r_to_footprint(int x, int y,
366  double minr2 = 1e12;
367  for (afwGeom::Span const & sp : *(tfoot->getSpans())) {
368  int mindx;
369  // span is to right of pixel?
370  int dx = sp.getX0() - x;
371  if (dx >= 0) {
372  mindx = dx;
373  } else {
374  // span is to left of pixel?
375  dx = x - sp.getX1();
376  if (dx >= 0) {
377  mindx = dx;
378  } else {
379  // span contains pixel (in x direction)
380  mindx = 0;
381  }
382  }
383  int dy = sp.getY() - y;
384  minr2 = std::min(minr2, (double)(mindx*mindx + dy*dy));
385  }
386  //printf("stray flux at (%i,%i): dist to t %i is %g\n", x, y, (int)i, sqrt(minr2));
387  return 1. / (1. + minr2);
388 }
389 
390 
391 template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
392 void
394 _find_stray_flux(det::Footprint const& foot,
395  ImagePtrT tsum,
396  MaskedImageT const& img,
397  int strayFluxOptions,
399  std::vector<bool> const& ispsf,
400  std::vector<int> const& pkx,
401  std::vector<int> const& pky,
402  double clipStrayFluxFraction,
404  ) {
405 
406  typedef typename det::HeavyFootprint<ImagePixelT, MaskPixelT, VariancePixelT> HeavyFootprint;
408 
409  // when doing stray flux: the footprints and pixels, which we'll
410  // combine into the return 'strays' HeavyFootprint at the end.
412  std::vector<std::vector<afwGeom::Span>> straySpans(tfoots.size());
416 
417  int ix0 = img.getX0();
418  int iy0 = img.getY0();
419  geom::Box2I sumbb = tsum->getBBox();
420  int sumx0 = sumbb.getMinX();
421  int sumy0 = sumbb.getMinY();
422 
423  for (size_t i=0; i<tfoots.size(); ++i) {
426  straymask.push_back(std::vector<MaskPixelT>());
428  }
429 
430  bool always = (strayFluxOptions & STRAYFLUX_TO_POINT_SOURCES_ALWAYS);
431 
432  typedef std::uint16_t itype;
434 
435  if (strayFluxOptions & STRAYFLUX_NEAREST_FOOTPRINT) {
436  // Compute the map of which footprint is closest to each
437  // pixel in the bbox.
438  typedef std::uint16_t dtype;
441 
443  std::vector<std::shared_ptr<det::Footprint>>* footlist = &tfoots;
444 
445  if (!always && ispsf.size()) {
446  // create a temp list that has empty footprints in place
447  // of all the point sources.
448  auto empty = std::make_shared<det::Footprint>();
449  empty->setPeakSchema(foot.getPeaks().getSchema());
450  for (size_t i=0; i<tfoots.size(); ++i) {
451  if (ispsf[i]) {
452  templist.push_back(empty);
453  } else {
454  templist.push_back(tfoots[i]);
455  }
456  }
457  footlist = &templist;
458  }
459  nearestFootprint(*footlist, nearest, dist);
460  }
461 
462  // Go through the (parent) Footprint looking for stray flux:
463  // pixels that are not claimed by any template, and positive.
464  for (afwGeom::Span const & s : *foot.getSpans()) {
465  int y = s.getY();
466  int x0 = s.getX0();
467  int x1 = s.getX1();
468  typename ImageT::x_iterator tsum_it =
469  tsum->row_begin(y - sumy0) + (x0 - sumx0);
470  typename MaskedImageT::x_iterator in_it =
471  img.row_begin(y - iy0) + (x0 - ix0);
472  double contrib[tfoots.size()];
473 
474  for (int x = x0; x <= x1; ++x, ++tsum_it, ++in_it) {
475  // Skip pixels that are covered by at least one
476  // template (*tsum_it > 0) or the input is not
477  // positive (*in_it <= 0).
478  if ((*tsum_it > 0) || (*in_it).image() <= 0) {
479  continue;
480  }
481 
482  if (strayFluxOptions & STRAYFLUX_R_TO_FOOTPRINT) {
483  // we'll compute these just-in-time
484  for (size_t i=0; i<tfoots.size(); ++i) {
485  contrib[i] = -1.0;
486  }
487  } else if (strayFluxOptions & STRAYFLUX_NEAREST_FOOTPRINT) {
488  for (size_t i=0; i<tfoots.size(); ++i) {
489  contrib[i] = 0.0;
490  }
491  int i = (*nearest)[geom::Point2I(x, y)];
492  contrib[i] = 1.0;
493  } else {
494  // R_TO_PEAK
495  for (size_t i=0; i<tfoots.size(); ++i) {
496  // Split the stray flux by 1/(1+r^2) to peaks
497  int dx, dy;
498  dx = pkx[i] - x;
499  dy = pky[i] - y;
500  contrib[i] = 1. / (1. + dx*dx + dy*dy);
501  }
502  }
503 
504  // Round 1: skip point sources unless STRAYFLUX_TO_POINT_SOURCES_ALWAYS
505  // are we going to assign stray flux to ptsrcs?
506  bool ptsrcs = always;
507  double csum = 0.;
508  for (size_t i=0; i<tfoots.size(); ++i) {
509  // if we're skipping point sources and this is a point source...
510  if ((!ptsrcs) && ispsf.size() && ispsf[i]) {
511  continue;
512  }
513  if (contrib[i] == -1.0) {
514  contrib[i] = _get_contrib_r_to_footprint(x, y, tfoots[i]);
515  }
516  csum += contrib[i];
517  }
518  if ((csum == 0.) &&
519  (strayFluxOptions &
520  STRAYFLUX_TO_POINT_SOURCES_WHEN_NECESSARY)) {
521  // No extended sources -- assign to pt sources
522  //LOGL_DEBUG(_log, "necessary to assign stray flux to point sources");
523  ptsrcs = true;
524  for (size_t i=0; i<tfoots.size(); ++i) {
525  if (contrib[i] == -1.0) {
526  contrib[i] = _get_contrib_r_to_footprint(x, y, tfoots[i]);
527  }
528  csum += contrib[i];
529  }
530  }
531 
532  // Drop small contributions...
533  double strayclip = (clipStrayFluxFraction * csum);
534  csum = 0.;
535  for (size_t i=0; i<tfoots.size(); ++i) {
536  // skip ptsrcs?
537  if ((!ptsrcs) && ispsf.size() && ispsf[i]) {
538  contrib[i] = 0.;
539  continue;
540  }
541  // skip small contributions
542  if (contrib[i] < strayclip) {
543  contrib[i] = 0.;
544  continue;
545  }
546  csum += contrib[i];
547  }
548 
549  for (size_t i=0; i<tfoots.size(); ++i) {
550  if (contrib[i] == 0.) {
551  continue;
552  }
553  // the stray flux to give to template i
554  double p = (contrib[i] / csum) * (*in_it).image();
555 
556  if (!strayfoot[i]) {
557  strayfoot[i] = std::make_shared<det::Footprint>();
558  strayfoot[i]->setPeakSchema(foot.getPeaks().getSchema());
559  }
560  straySpans[i].push_back(afwGeom::Span(y, x, x));
561  straypix[i].push_back(p);
562  straymask[i].push_back((*in_it).mask());
563  strayvar[i].push_back((*in_it).variance());
564  }
565  }
566  }
567 
568  // Store the stray flux in HeavyFootprints
569  for (size_t i=0; i<tfoots.size(); ++i) {
570  if (strayfoot[i]) {
571  strayfoot[i]->setSpans(std::make_shared<afwGeom::SpanSet>(straySpans[i]));
572  }
573  if (!strayfoot[i]) {
574  strays.push_back(HeavyFootprintPtrT());
575  } else {
579  HeavyFootprintPtrT heavy(new HeavyFootprint(*strayfoot[i]));
580  ndarray::Array<ImagePixelT,1,1> himg = heavy->getImageArray();
587 
588  assert((size_t)strayfoot[i]->getArea() == straypix[i].size());
589 
590  for (spix = straypix[i].begin(),
591  smask = straymask[i].begin(),
592  svar = strayvar [i].begin(),
593  hpix = himg.begin(),
594  mpix = heavy->getMaskArray().begin(),
595  vpix = heavy->getVarianceArray().begin();
596  spix != straypix[i].end();
597  ++spix, ++smask, ++svar, ++hpix, ++mpix, ++vpix) {
598  *hpix = *spix;
599  *mpix = *smask;
600  *vpix = *svar;
601  }
602  strays.push_back(heavy);
603  }
604  }
605 }
606 
607 template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
608 void
611  ImagePtrT tsum) {
612  geom::Box2I sumbb = tsum->getBBox();
613  int sumx0 = sumbb.getMinX();
614  int sumy0 = sumbb.getMinY();
615 
616  // Compute tsum = the sum of templates
617  for (size_t i=0; i<timgs.size(); ++i) {
618  ImagePtrT timg = timgs[i];
619  geom::Box2I tbb = timg->getBBox();
620  int tx0 = tbb.getMinX();
621  int ty0 = tbb.getMinY();
622  // To handle "ramped" templates that can extend outside the
623  // parent, clip the bbox. Note that we saved tx0,ty0 BEFORE
624  // doing this!
625  tbb.clip(sumbb);
626  int copyx0 = tbb.getMinX();
627  // Here we iterate over the template bbox -- we could instead
628  // iterate over the "tfoot"s.
629  for (int y=tbb.getMinY(); y<=tbb.getMaxY(); ++y) {
630  typename ImageT::x_iterator in_it = timg->row_begin(y - ty0) + (copyx0 - tx0);
631  typename ImageT::x_iterator inend = in_it + tbb.getWidth();
632  typename ImageT::x_iterator tsum_it =
633  tsum->row_begin(y - sumy0) + (copyx0 - sumx0);
634  for (; in_it != inend; ++in_it, ++tsum_it) {
635  *tsum_it += std::max((ImagePixelT)0., static_cast<ImagePixelT>(*in_it));
636  }
637  }
638  }
639 
640 }
641 
690 template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
693 apportionFlux(MaskedImageT const& img,
694  det::Footprint const& foot,
697  ImagePtrT tsum,
698  std::vector<bool> const& ispsf,
699  std::vector<int> const& pkx,
700  std::vector<int> const& pky,
702  int strayFluxOptions,
703  double clipStrayFluxFraction
704  ) {
705 
706  if (timgs.size() != tfoots.size()) {
708  (boost::format("Template images must be the same length as template footprints (%d vs %d)")
709  % timgs.size() % tfoots.size()).str());
710  }
711 
712  for (size_t i=0; i<timgs.size(); ++i) {
713  if (!timgs[i]->getBBox().contains(tfoots[i]->getBBox())) {
715  "Template image MUST contain template footprint");
716  }
717  if (!img.getBBox().contains(foot.getBBox())) {
719  "Image bbox MUST contain parent footprint");
720  }
721  // template bounding-boxes *can* extend outside the parent
722  // footprint if we are ramping templates with significant flux
723  // at the edges. We handle this below.
724  }
725 
726  // the apportioned flux return value
728 
729  LOG_LOGGER _log = LOG_GET("lsst.meas.deblender.apportionFlux");
730  bool findStrayFlux = (strayFluxOptions & ASSIGN_STRAYFLUX);
731 
732  int ix0 = img.getX0();
733  int iy0 = img.getY0();
734  geom::Box2I fbb = foot.getBBox();
735 
736  if (!tsum) {
737  tsum = ImagePtrT(new ImageT(fbb.getDimensions()));
738  tsum->setXY0(fbb.getMinX(), fbb.getMinY());
739  }
740 
741  if (!tsum->getBBox().contains(foot.getBBox())) {
743  "Template sum image MUST contain parent footprint");
744  }
745 
746  geom::Box2I sumbb = tsum->getBBox();
747  int sumx0 = sumbb.getMinX();
748  int sumy0 = sumbb.getMinY();
749 
750  _sum_templates(timgs, tsum);
751 
752  // Compute flux portions
753  for (size_t i=0; i<timgs.size(); ++i) {
754  ImagePtrT timg = timgs[i];
755  // Initialize return value:
756  MaskedImagePtrT port(new MaskedImageT(timg->getDimensions()));
757  port->setXY0(timg->getXY0());
758  portions.push_back(port);
759 
760  // Split flux = image * template / tsum
761  geom::Box2I tbb = timg->getBBox();
762  int tx0 = tbb.getMinX();
763  int ty0 = tbb.getMinY();
764  // As above
765  tbb.clip(sumbb);
766  int copyx0 = tbb.getMinX();
767  for (int y=tbb.getMinY(); y<=tbb.getMaxY(); ++y) {
768  typename MaskedImageT::x_iterator in_it =
769  img.row_begin(y - iy0) + (copyx0 - ix0);
770  typename ImageT::x_iterator tptr =
771  timg->row_begin(y - ty0) + (copyx0 - tx0);
772  typename ImageT::x_iterator tend = tptr + tbb.getWidth();
773  typename ImageT::x_iterator tsum_it =
774  tsum->row_begin(y - sumy0) + (copyx0 - sumx0);
775  typename MaskedImageT::x_iterator out_it =
776  port->row_begin(y - ty0) + (copyx0 - tx0);
777  for (; tptr != tend; ++tptr, ++in_it, ++out_it, ++tsum_it) {
778  if (*tsum_it == 0) {
779  continue;
780  }
781  double frac = std::max((ImagePixelT)0., static_cast<ImagePixelT>(*tptr)) / (*tsum_it);
782  //if (frac == 0) {
783  // treat mask planes differently?
784  // }
785  out_it.mask() = (*in_it).mask();
786  out_it.variance() = (*in_it).variance();
787  out_it.image() = (*in_it).image() * frac;
788  }
789  }
790  }
791 
792  if (findStrayFlux) {
793  if ((ispsf.size() > 0) && (ispsf.size() != timgs.size())) {
795  (boost::format("'ispsf' must be the same length as templates (%d vs %d)")
796  % ispsf.size() % timgs.size()).str());
797  }
798  if ((pkx.size() != timgs.size()) || (pky.size() != timgs.size())) {
800  (boost::format("'pkx' and 'pky' must be the same length as templates (%d,%d vs %d)")
801  % pkx.size() % pky.size() % timgs.size()).str());
802  }
803  _find_stray_flux(foot, tsum, img, strayFluxOptions, tfoots,
804  ispsf, pkx, pky, clipStrayFluxFraction, strays);
805  }
806  return portions;
807 }
808 
809 
819 public:
821 
823 
825  SpanSet const& arr,
826  int cx, int cy, bool forward=true)
827  : _real(real), _cx(cx), _cy(cy), _forward(forward)
828  {
829  if (_forward) {
830  _end = arr.end();
831  } else {
832  _end = arr.begin();
833  }
834  }
835 
836  bool operator==(const SpanSet::const_iterator & other) {
837  return _real == other;
838  }
839  bool operator!=(const SpanSet::const_iterator & other) {
840  return _real != other;
841  }
842  bool operator<=(const SpanSet::const_iterator & other) {
843  return _real <= other;
844  }
845  bool operator<(const SpanSet::const_iterator & other) {
846  return _real < other;
847  }
848  bool operator>=(const SpanSet::const_iterator & other) {
849  return _real >= other;
850  }
851  bool operator>(const SpanSet::const_iterator & other) {
852  return _real > other;
853  }
854 
856  return (_real == other._real) &&
857  (_cx == other._cx) && (_cy == other._cy) &&
858  (_forward == other._forward);
859  }
861  return !(*this == other);
862  }
863 
865  if (_forward) {
866  _real++;
867  } else {
868  _real--;
869  }
870  return *this;
871  }
874  ++(*this);
875  return result;
876  }
877 
878  bool notDone() {
879  if (_forward) {
880  return _real < _end;
881  } else {
882  return _real >= _end;
883  }
884  }
885 
886  int dxlo() {
887  if (_forward) {
888  return _real->getX0() - _cx;
889  } else {
890  return _cx - _real->getX1();
891  }
892  }
893  int dxhi() {
894  if (_forward) {
895  return _real->getX1() - _cx;
896  } else {
897  return _cx - _real->getX0();
898  }
899  }
900  int dy() {
901  return std::abs(_real->getY() - _cy);
902  }
903  int x0() {
904  return _real->getX0();
905  }
906  int x1() {
907  return _real->getX1();
908  }
909  int y() {
910  return _real->getY();
911  }
912 
913 private:
916  int _cx;
917  int _cy;
918  bool _forward;
919 };
920 
921 
922 /*
923  // Check symmetrizeFootprint by computing truth naively.
924  // compute correct answer dumbly
925  det::Footprint truefoot;
926  geom::Box2I bbox = foot.getBBox();
927  for (int y=bbox.getMinY(); y<=bbox.getMaxY(); y++) {
928  for (int x=bbox.getMinX(); x<=bbox.getMaxX(); x++) {
929  int dy = y - cy;
930  int dx = x - cx;
931  int x2 = cx - dx;
932  int y2 = cy - dy;
933  if (foot.contains(geom::Point2I(x, y)) &&
934  foot.contains(geom::Point2I(x2, y2))) {
935  truefoot.addSpan(y, x, x);
936  }
937  }
938  }
939  truefoot.normalize();
940 
941  SpanList sp1 = truefoot.getSpans();
942  SpanList sp2 = sfoot->getSpans();
943  for (SpanList::const_iterator spit1 = sp1.begin(),
944  spit2 = sp2.begin();
945  spit1 != sp1.end() && spit2 != sp2.end();
946  spit1++, spit2++) {
947  //printf("\n");
948  printf(" true y %i, x [%i, %i]\n", (*spit1)->getY(), (*spit1)->getX0(), (*spit1)->getX1());
949  printf(" sfoot y %i, x [%i, %i]\n", (*spit2)->getY(), (*spit2)->getX0(), (*spit2)->getX1());
950  if (((*spit1)->getY() != (*spit2)->getY()) ||
951  ((*spit1)->getX0() != (*spit2)->getX0()) ||
952  ((*spit1)->getX1() != (*spit2)->getX1())) {
953  printf("*******************************************\n");
954  }
955  }
956  assert(sp1.size() == sp2.size());
957  for (SpanList::const_iterator spit1 = sp1.begin(),
958  spit2 = sp2.begin();
959  spit1 != sp1.end() && spit2 != sp2.end();
960  spit1++, spit2++) {
961  assert((*spit1)->getY() == (*spit2)->getY());
962  assert((*spit1)->getX0() == (*spit2)->getX0());
963  assert((*spit1)->getX1() == (*spit2)->getX1());
964  }
965  */
966 
972 template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
976  det::Footprint const& foot,
977  int cx, int cy) {
978 
979  auto sfoot = std::make_shared<det::Footprint>();
980  sfoot->setPeakSchema(foot.getPeaks().getSchema());
981  afwGeom::SpanSet const & spans = *foot.getSpans();
982 
983  LOG_LOGGER _log = LOG_GET("lsst.meas.deblender.symmetrizeFootprint");
984 
985  // Find the Span containing the peak.
986  afwGeom::Span target(cy, cx, cx);
988  std::upper_bound(spans.begin(), spans.end(), target, span_compare);
989  // upper_bound returns the last position where "target" could be inserted;
990  // ie, the first Span larger than "target". The Span containing "target"
991  // should be peakspan-1 or (if the peak is on the first pixel in the span,
992  // peakspan.
993  afwGeom::Span sp;
994  if (peakspan == spans.begin()) {
995  sp = *peakspan;
996  if (!sp.contains(cx, cy)) {
997  LOGL_WARN(_log,
998  "Failed to find span containing (%i,%i): before the beginning of this footprint", cx, cy);
1000  }
1001  } else {
1002  --peakspan;
1003  sp = *peakspan;
1004 
1005  if (!sp.contains(cx, cy)) {
1006  ++peakspan;
1007  sp = *peakspan;
1008  if (!sp.contains(cx, cy)) {
1009  geom::Box2I fbb = foot.getBBox();
1010  LOGL_WARN(_log, "Failed to find span containing (%i,%i): nearest is %i, [%i,%i]. "
1011  "Footprint bbox is [%i,%i],[%i,%i]",
1012  cx, cy, sp.getY(), sp.getX0(), sp.getX1(),
1013  fbb.getMinX(), fbb.getMaxX(), fbb.getMinY(), fbb.getMaxY());
1015  }
1016  }
1017  }
1018  LOGL_DEBUG(_log, "Span containing (%i,%i): (x=[%i,%i], y=%i)",
1019  cx, cy, sp.getX0(), sp.getX1(), sp.getY());
1020 
1021  // The symmetric templates are essentially an AND of the footprint
1022  // pixels and its 180-degree-rotated self, rotated around the
1023  // peak (cx,cy).
1024  //
1025  // We iterate forward and backward simultaneously, starting from
1026  // the span containing the peak and moving out, row by row.
1027  //
1028  // In the loop below, we search for the next pair of Spans that
1029  // overlap (in "dx" from the center), output the overlapping
1030  // portion of the Spans, and advance either the "fwd" or "back"
1031  // iterator. When we fail to find an overlapping pair of Spans,
1032  // we move on to the next row.
1033  //
1034  // [The following paragraph is somewhat obsoleted by the
1035  // RelativeSpanIterator class, which performs some of the renaming
1036  // and the dx,dy coords.]
1037  //
1038  // '''In reading the code, "forward", "advancing", etc, are all
1039  // from the perspective of the "fwd" iterator (the one going
1040  // forward through the Span list, from low to high Y and then low
1041  // to high X). It will help to imagine making a copy of the
1042  // footprint and rotating it around the center pixel by 180
1043  // degrees, so that "fwd" and "back" are both iterating the same
1044  // direction; we're then just finding the AND of those two
1045  // iterators, except we have to work in dx,dy coordinates rather
1046  // than original x,y coords, and the accessors for "back" are
1047  // opposite.'''
1048 
1049  RelativeSpanIterator fwd (peakspan, spans, cx, cy, true);
1050  RelativeSpanIterator back(peakspan, spans, cx, cy, false);
1051 
1052  int dy = 0;
1053  std::vector<afwGeom::Span> tmpSpans;
1054  while (fwd.notDone() && back.notDone()) {
1055  // forward and backward "y"; just symmetric around cy
1056  int fy = cy + dy;
1057  int by = cy - dy;
1058  // delta-x of the beginnings of the spans, for "fwd" and "back"
1059  int fdxlo = fwd.dxlo();
1060  int bdxlo = back.dxlo();
1061 
1062  // First find:
1063  // fend -- first span in the next row, or end(); ie,
1064  // the end of this row in the forward direction
1065  // bend -- the end of this row in the backward direction
1066  RelativeSpanIterator fend, bend;
1067  for (fend = fwd; fend.notDone(); ++fend) {
1068  if (fend.dy() != dy)
1069  break;
1070  }
1071  for (bend = back; bend.notDone(); ++bend) {
1072  if (bend.dy() != dy)
1073  break;
1074  }
1075 
1076  LOGL_DEBUG(_log, "dy=%i, fy=%i, fx=[%i, %i], by=%i, fx=[%i, %i], fdx=%i, bdx=%i",
1077  dy, fy, fwd.x0(), fwd.x1(), by, back.x0(), back.x1(),
1078  fdxlo, bdxlo);
1079 
1080  // Find possibly-overlapping span
1081  if (bdxlo > fdxlo) {
1082  LOGL_DEBUG(_log, "Advancing forward.");
1083  // While the "forward" span is entirely to the "left" of the "backward" span,
1084  // (in dx coords), ie, |---fwd---X X---back---|
1085  // and we are comparing the edges marked X
1086  while ((fwd != fend) && (fwd.dxhi() < bdxlo)) {
1087  fwd++;
1088  if (fwd == fend) {
1089  LOGL_DEBUG(_log, "Reached fend");
1090  } else {
1091  LOGL_DEBUG(_log, "Advanced to forward span %i, [%i, %i]",
1092  fy, fwd.x0(), fwd.x1());
1093  }
1094  }
1095  } else if (fdxlo > bdxlo) {
1096  LOGL_DEBUG(_log, "Advancing backward.");
1097  // While the "backward" span is entirely to the "left" of the "foreward" span,
1098  // (in dx coords), ie, |---back---X X---fwd---|
1099  // and we are comparing the edges marked X
1100  while ((back != bend) && (back.dxhi() < fdxlo)) {
1101  back++;
1102  if (back == bend) {
1103  LOGL_DEBUG(_log, "Reached bend");
1104  } else {
1105  LOGL_DEBUG(_log, "Advanced to backward span %i, [%i, %i]",
1106  by, back.x0(), back.x1());
1107  }
1108  }
1109  }
1110 
1111  if ((back == bend) || (fwd == fend)) {
1112  // We reached the end of the row without finding spans that could
1113  // overlap. Move onto the next dy.
1114  if (back == bend) {
1115  LOGL_DEBUG(_log, "Reached bend");
1116  }
1117  if (fwd == fend) {
1118  LOGL_DEBUG(_log, "Reached fend");
1119  }
1120  back = bend;
1121  fwd = fend;
1122  dy++;
1123  continue;
1124  }
1125 
1126  // Spans may overlap -- find the overlapping part.
1127  int dxlo = std::max(fwd.dxlo(), back.dxlo());
1128  int dxhi = std::min(fwd.dxhi(), back.dxhi());
1129  if (dxlo <= dxhi) {
1130  LOGL_DEBUG(_log, "Adding span fwd %i, [%i, %i], back %i, [%i, %i]",
1131  fy, cx+dxlo, cx+dxhi, by, cx-dxhi, cx-dxlo);
1132  tmpSpans.push_back(afwGeom::Span(fy, cx + dxlo, cx + dxhi));
1133  tmpSpans.push_back(afwGeom::Span(by, cx - dxhi, cx - dxlo));
1134  }
1135 
1136  // Advance the one whose "hi" edge is smallest
1137  if (fwd.dxhi() < back.dxhi()) {
1138  fwd++;
1139  if (fwd == fend) {
1140  LOGL_DEBUG(_log, "Stepped to fend");
1141  } else {
1142  LOGL_DEBUG(_log, "Stepped forward to span %i, [%i, %i]",
1143  fwd.y(), fwd.x0(), fwd.x1());
1144  }
1145  } else {
1146  back++;
1147  if (back == bend) {
1148  LOGL_DEBUG(_log, "Stepped to bend");
1149  } else {
1150  LOGL_DEBUG(_log, "Stepped backward to span %i, [%i, %i]",
1151  back.y(), back.x0(), back.x1());
1152  }
1153  }
1154 
1155  if ((back == bend) || (fwd == fend)) {
1156  // Reached the end of the row. On to the next dy!
1157  if (back == bend) {
1158  LOGL_DEBUG(_log, "Reached bend");
1159  }
1160  if (fwd == fend) {
1161  LOGL_DEBUG(_log, "Reached fend");
1162  }
1163  back = bend;
1164  fwd = fend;
1165  dy++;
1166  continue;
1167  }
1168 
1169  }
1170  sfoot->setSpans(std::make_shared<afwGeom::SpanSet>(std::move(tmpSpans)));
1171  return sfoot;
1172 }
1173 
1186 template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
1191  MaskedImageT const& img,
1192  det::Footprint const& foot,
1193  det::PeakRecord const& peak,
1194  double sigma1,
1195  bool minZero,
1196  bool patchEdge,
1197  bool* patchedEdges) {
1198 
1199  typedef typename MaskedImageT::const_xy_locator xy_loc;
1200 
1201  *patchedEdges = false;
1202 
1203  int cx = peak.getIx();
1204  int cy = peak.getIy();
1205 
1206  LOG_LOGGER _log = LOG_GET("lsst.meas.deblender.symmetricFootprint");
1207 
1208  if (!img.getBBox(image::PARENT).contains(foot.getBBox())) {
1209  throw LSST_EXCEPT(lsst::pex::exceptions::LengthError, "Image too small for footprint");
1210  }
1211 
1212  FootprintPtrT sfoot = symmetrizeFootprint(foot, cx, cy);
1213 
1214  if (!sfoot) {
1216  }
1217 
1218  if (!img.getBBox(image::PARENT).contains(sfoot->getBBox())) {
1220  "Image too small for symmetrized footprint");
1221  }
1222  afwGeom::SpanSet const & spans = *sfoot->getSpans();
1223 
1224  // does this footprint touch an EDGE?
1225  bool touchesEdge = false;
1226  if (patchEdge) {
1227  LOGL_DEBUG(_log, "Checking footprint for EDGE bits");
1228  MaskPtrT mask = img.getMask();
1229  bool edge = false;
1230  MaskPixelT edgebit = mask->getPlaneBitMask("EDGE");
1231  for (afwGeom::SpanSet::const_iterator fwd=spans.begin();
1232  fwd != spans.end(); ++fwd) {
1233  int x0 = fwd->getX0();
1234  int x1 = fwd->getX1();
1235  typename MaskT::x_iterator xiter =
1236  mask->x_at(x0 - mask->getX0(), fwd->getY() - mask->getY0());
1237  for (int x=x0; x<=x1; ++x, ++xiter) {
1238  if ((*xiter) & edgebit) {
1239  edge = true;
1240  break;
1241  }
1242  }
1243  if (edge)
1244  break;
1245  }
1246  if (edge) {
1247  LOGL_DEBUG(_log, "Footprint includes an EDGE pixel.");
1248  touchesEdge = true;
1249  }
1250  }
1251 
1252  // The result image:
1253  ImagePtrT targetimg(new ImageT(sfoot->getBBox()));
1254 
1256  afwGeom::SpanSet::const_iterator back = spans.end()-1;
1257 
1258  ImagePtrT theimg = img.getImage();
1259 
1260  for (; fwd <= back; fwd++, back--) {
1261  int fy = fwd->getY();
1262  int by = back->getY();
1263 
1264  for (int fx=fwd->getX0(), bx=back->getX1();
1265  fx <= fwd->getX1();
1266  fx++, bx--) {
1267  // FIXME -- CURRENTLY WE IGNORE THE MASK PLANE! options
1268  // include ORing the mask bits, or being clever about
1269  // ignoring some masked pixels, or copying the mask bits
1270  // of the min pixel
1271 
1272  // We have already checked the bounding box, so this should always be satisfied
1273  assert(theimg->getBBox(image::PARENT).contains(geom::Point2I(fx, fy)));
1274  assert(theimg->getBBox(image::PARENT).contains(geom::Point2I(bx, by)));
1275 
1276  // FIXME -- we could do this with image iterators instead.
1277  // But first profile to show that it's necessary and an
1278  // improvement.
1279  ImagePixelT pixf = (*theimg)[geom::Point2I(fx, fy)];
1280  ImagePixelT pixb = (*theimg)[geom::Point2I(bx, by)];
1281  ImagePixelT pix = std::min(pixf, pixb);
1282  if (minZero) {
1283  pix = std::max(pix, static_cast<ImagePixelT>(0));
1284  }
1285  (*targetimg)[geom::Point2I(fx, fy)] = pix;
1286  (*targetimg)[geom::Point2I(bx, by)] = pix;
1287 
1288  }
1289  }
1290 
1291  if (touchesEdge) {
1292  // Find spans whose mirrors fall outside the image bounds,
1293  // grow the footprint to include those spans, and plug in
1294  // their pixel values.
1295  geom::Box2I bb = sfoot->getBBox();
1296 
1297  // Actually, it's not necessarily the IMAGE bounds that count
1298  //-- the footprint may not go right to the image edge.
1299  //geom::Box2I imbb = img.getBBox();
1300  geom::Box2I imbb = foot.getBBox();
1301 
1302  LOGL_DEBUG(_log, "Footprint touches EDGE: start bbox [%i,%i],[%i,%i]",
1303  bb.getMinX(), bb.getMaxX(), bb.getMinY(), bb.getMaxY());
1304  // original footprint spans
1305  const afwGeom::SpanSet & ospans = *foot.getSpans();
1306  for (fwd = ospans.begin(); fwd != ospans.end(); ++fwd) {
1307  int y = fwd->getY();
1308  int x = fwd->getX0();
1309  // mirrored coords
1310  int ym = cy + (cy - y);
1311  int xm = cx + (cx - x);
1312  if (!imbb.contains(geom::Point2I(xm, ym))) {
1313  bb.include(geom::Point2I(x, y));
1314  }
1315  x = fwd->getX1();
1316  xm = cx + (cx - x);
1317  if (!imbb.contains(geom::Point2I(xm, ym))) {
1318  bb.include(geom::Point2I(x, y));
1319  }
1320  }
1321  LOGL_DEBUG(_log, "Footprint touches EDGE: grown bbox [%i,%i],[%i,%i]",
1322  bb.getMinX(), bb.getMaxX(), bb.getMinY(), bb.getMaxY());
1323 
1324  // New template image
1325  ImagePtrT targetimg2(new ImageT(bb));
1326  sfoot->getSpans()->copyImage(*targetimg, *targetimg2);
1327 
1328  LOGL_DEBUG(_log, "Symmetric footprint spans:");
1329  const afwGeom::SpanSet & sspans = *sfoot->getSpans();
1330  for (fwd = sspans.begin(); fwd != sspans.end(); ++fwd) {
1331  LOGL_DEBUG(_log, " %s", fwd->toString().c_str());
1332  }
1333 
1334  // copy original 'img' pixels for the portion of spans whose
1335  // mirrors are out of bounds.
1336  std::vector<afwGeom::Span> newSpans(sfoot->getSpans()->begin(), sfoot->getSpans()->end());
1337  for (fwd = ospans.begin(); fwd != ospans.end(); ++fwd) {
1338  int y = fwd->getY();
1339  int x0 = fwd->getX0();
1340  int x1 = fwd->getX1();
1341  // mirrored coords
1342  int ym = cy + (cy - y);
1343  int xm0 = cx + (cx - x0);
1344  int xm1 = cx + (cx - x1);
1345  bool in0 = imbb.contains(geom::Point2I(xm0, ym));
1346  bool in1 = imbb.contains(geom::Point2I(xm1, ym));
1347  if (in0 && in1) {
1348  // both endpoints of the symmetric span are in bounds; nothing to do
1349  continue;
1350  }
1351  // clip to the part of the span where the mirror is out of bounds
1352  if (in0) {
1353  // the mirror of x0 is in-bounds; move x0 to be the first pixel
1354  // whose mirror would be out-of-bounds
1355  x0 = cx + (cx - (imbb.getMinX() - 1));
1356  }
1357  if (in1) {
1358  x1 = cx + (cx - (imbb.getMaxX() + 1));
1359  }
1360  LOGL_DEBUG(_log, "Span y=%i, x=[%i,%i] has mirror (%i,[%i,%i]) out-of-bounds; clipped to %i,[%i,%i]",
1361  y, fwd->getX0(), fwd->getX1(), ym, xm1, xm0, y, x0, x1);
1362  typename MaskedImageT::x_iterator initer =
1363  img.x_at(x0 - img.getX0(), y - img.getY0());
1364  typename ImageT::x_iterator outiter =
1365  targetimg2->x_at(x0 - targetimg2->getX0(), y - targetimg2->getY0());
1366  for (int x=x0; x<=x1; ++x, ++outiter, ++initer) {
1367  *outiter = initer.image();
1368  }
1369  newSpans.push_back(afwGeom::Span(y, x0, x1));
1370  }
1371  sfoot->setSpans(std::make_shared<afwGeom::SpanSet>(std::move(newSpans)));
1372  targetimg = targetimg2;
1373  }
1374 
1375  *patchedEdges = touchesEdge;
1376  return std::pair<ImagePtrT, FootprintPtrT>(targetimg, sfoot);
1377 }
1378 
1383 template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
1384 bool
1388  ImagePixelT thresh) {
1389 
1390  LOG_LOGGER _log = LOG_GET("lsst.meas.deblender.hasSignificantFluxAtEdge");
1391 
1392  // Find edge template pixels with significant flux -- perhaps
1393  // because their symmetric pixels were outside the footprint?
1394  // (clipped by an image edge, etc)
1395  std::shared_ptr<afwGeom::SpanSet> spans = sfoot->getSpans()->findEdgePixels();
1396 
1397  for (afwGeom::SpanSet::const_iterator sp = spans->begin(); sp != spans->end(); ++sp) {
1398  int const y = sp->getY();
1399  int const x0 = sp->getX0();
1400  int const x1 = sp->getX1();
1401  int x;
1402  typename ImageT::const_x_iterator xiter;
1403  for (xiter = img->x_at(x0 - img->getX0(), y - img->getY0()), x=x0; x<=x1; ++x, ++xiter) {
1404  if (*xiter >= thresh) {
1405  return true;
1406  }
1407  }
1408  }
1409  return false;
1410 }
1411 
1416 template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
1421  ImagePixelT thresh) {
1422  LOG_LOGGER _log = LOG_GET("meas.deblender.getSignificantEdgePixels");
1423 
1424 
1425  auto significant = std::make_shared<det::Footprint>();
1426  significant->setPeakSchema(sfoot->getPeaks().getSchema());
1427 
1428  int const x0 = img->getX0(), y0 = img->getY0();
1429  std::shared_ptr<afwGeom::SpanSet> edgeSpans = sfoot->getSpans()->findEdgePixels();
1430  std::vector<afwGeom::Span> tmpSpans;
1431  for (afwGeom::SpanSet::const_iterator ss = edgeSpans->begin(); ss != edgeSpans->end(); ++ss) {
1432  afwGeom::Span const& span = *ss;
1433  int const y = span.getY();
1434  int x = span.getX0();
1435  typename ImageT::const_x_iterator iter = img->x_at(x - x0, y - y0);
1436  bool onSpan = false; // Are we in a span of interest
1437  int xSpan; // Starting x of span
1438  for (; x <= span.getX1(); ++x, ++iter) {
1439  if (*iter >= thresh) {
1440  onSpan = true;
1441  xSpan = x;
1442  } else if (onSpan) {
1443  onSpan = false;
1444  tmpSpans.push_back(afwGeom::Span(y, xSpan, x - 1));
1445  }
1446  }
1447  if (onSpan) {
1448  tmpSpans.push_back(afwGeom::Span(y, xSpan, span.getX1()));
1449  }
1450  }
1451  significant->setSpans(std::make_shared<afwGeom::SpanSet>(std::move(tmpSpans)));
1452  return significant;
1453 }
1454 
1455 
1456 // Instantiate
1457 template class deblend::BaselineUtils<float>;
py::object result
Definition: _schema.cc:429
Key< Flag > const & target
int end
double x
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
table::Key< double > sigma1
afw::table::Key< afw::table::Array< ImagePixelT > > image
afw::table::Key< afw::table::Array< MaskPixelT > > mask
LSST DM logging module built on log4cxx.
#define LOGL_WARN(logger, message...)
Log a warn-level message using a varargs/printf style interface.
Definition: Log.h:547
#define LOG_GET(logger)
Returns a Log object associated with logger.
Definition: Log.h:75
#define LOG_LOGGER
Definition: Log.h:714
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
Definition: Log.h:515
int y
Definition: SpanSet.cc:48
This is a convenience class used in symmetrizeFootprint, wrapping the idea of iterating through a Spa...
bool operator<(const SpanSet::const_iterator &other)
RelativeSpanIterator(SpanSet::const_iterator const &real, SpanSet const &arr, int cx, int cy, bool forward=true)
bool operator!=(const SpanSet::const_iterator &other)
RelativeSpanIterator operator++(int dummy)
RelativeSpanIterator operator++()
bool operator!=(RelativeSpanIterator &other)
bool operator>(const SpanSet::const_iterator &other)
bool operator==(RelativeSpanIterator &other)
bool operator<=(const SpanSet::const_iterator &other)
bool operator>=(const SpanSet::const_iterator &other)
bool operator==(const SpanSet::const_iterator &other)
Class to describe the properties of a detected object from an image.
Definition: Footprint.h:63
std::shared_ptr< geom::SpanSet > getSpans() const
Return a shared pointer to the SpanSet.
Definition: Footprint.h:115
lsst::geom::Box2I getBBox() const
Return the Footprint's bounding box.
Definition: Footprint.h:208
PeakCatalog & getPeaks()
Return the Peaks contained in this Footprint.
Definition: Footprint.h:129
void setSpans(std::shared_ptr< geom::SpanSet > otherSpanSet)
Sets the shared pointer to the SpanSet in the Footprint.
Definition: Footprint.cc:45
A set of pixels in an Image, including those pixels' actual values.
Record class that represents a peak in a Footprint.
Definition: Peak.h:42
A range of pixels within one row of an Image.
Definition: Span.h:47
bool contains(int x) const noexcept
Definition: Span.h:95
int getY() const noexcept
Return the y-value.
Definition: Span.h:81
int getX0() const noexcept
Return the starting x-value.
Definition: Span.h:76
int getX1() const noexcept
Return the ending x-value.
Definition: Span.h:78
A compact representation of a collection of pixels.
Definition: SpanSet.h:78
const_iterator end() const
Definition: SpanSet.h:89
const_iterator begin() const
Definition: SpanSet.h:88
std::vector< Span >::const_iterator const_iterator
Definition: SpanSet.h:80
int getX0() const
Return the image's column-origin.
Definition: ImageBase.h:306
xy_locator xy_at(int x, int y) const
Return an xy_locator at the point (x, y) in the image.
Definition: ImageBase.h:425
typename _view_t::xy_locator xy_locator
An xy_locator.
Definition: ImageBase.h:121
typename _const_view_t::x_iterator const_x_iterator
A const iterator for traversing the pixels in a row.
Definition: ImageBase.h:141
int getWidth() const
Return the number of columns in the image.
Definition: ImageBase.h:294
int getY0() const
Return the image's row-origin.
Definition: ImageBase.h:314
int getHeight() const
Return the number of rows in the image.
Definition: ImageBase.h:296
x_iterator row_begin(int y) const
Return an x_iterator to the start of the y'th row.
Definition: ImageBase.h:401
x_iterator row_end(int y) const
Return an x_iterator to the end of the y'th row.
Definition: ImageBase.h:404
typename _view_t::x_iterator x_iterator
An iterator for traversing the pixels in a row.
Definition: ImageBase.h:133
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:51
Ref< MaskPixelT >::type mask()
Return (a reference to) the mask part of the Pixel pointed at by the iterator.
Definition: MaskedImage.h:142
Ref< ImagePixelT >::type image()
Return (a reference to) the image part of the Pixel pointed at by the iterator.
Definition: MaskedImage.h:139
Ref< VariancePixelT >::type variance()
Return (a reference to) the variance part of the Pixel pointed at by the iterator.
Definition: MaskedImage.h:145
A const locator for the MaskedImage.
Definition: MaskedImage.h:512
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:73
int getX0() const
Return the image's column-origin.
Definition: MaskedImage.h:1066
lsst::geom::Box2I getBBox(ImageOrigin const origin=PARENT) const
Definition: MaskedImage.h:1058
int getY0() const
Return the image's row-origin.
Definition: MaskedImage.h:1074
x_iterator row_begin(int y) const
Return an x_iterator to the start of the image.
Definition: MaskedImage.cc:621
x_iterator x_at(int x, int y) const
Return an x_iterator at the point (x, y)
Definition: MaskedImage.h:1180
MaskPtr getMask() const
Return a (shared_ptr to) the MaskedImage's mask.
Definition: MaskedImage.h:1030
ImagePtr getImage() const
Return a (shared_ptr to) the MaskedImage's image.
Definition: MaskedImage.h:1018
MaskPixelT mask() const
Return the mask part of a Pixel.
Definition: Pixel.h:219
ImagePixelT image() const
Return the image part of a Pixel.
Definition: Pixel.h:217
VariancePixelT variance() const
Return the variance part of a Pixel.
Definition: Pixel.h:221
Schema getSchema() const
Return the schema associated with the catalog's table.
Definition: Catalog.h:118
An integer coordinate rectangle.
Definition: Box.h:55
void clip(Box2I const &other) noexcept
Shrink this to ensure that other.contains(*this).
Definition: Box.cc:189
int getMinY() const noexcept
Definition: Box.h:158
void include(Point2I const &point)
Expand this to ensure that this->contains(point).
Definition: Box.cc:152
int getMinX() const noexcept
Definition: Box.h:157
int getWidth() const noexcept
Definition: Box.h:187
int getMaxX() const noexcept
Definition: Box.h:161
bool contains(Point2I const &point) const noexcept
Return true if the box contains the point.
Definition: Box.cc:114
int getMaxY() const noexcept
Definition: Box.h:162
Extent2I const getDimensions() const noexcept
Definition: Box.h:186
static void medianFilter(ImageT const &img, ImageT &outimg, int halfsize)
Run a spatial median filter over the given input img, writing the results to out.
static std::vector< typename std::shared_ptr< lsst::afw::image::MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > > > apportionFlux(MaskedImageT const &img, lsst::afw::detection::Footprint const &foot, std::vector< typename std::shared_ptr< lsst::afw::image::Image< ImagePixelT >>> templates, std::vector< std::shared_ptr< lsst::afw::detection::Footprint > > templ_footprints, ImagePtrT templ_sum, std::vector< bool > const &ispsf, std::vector< int > const &pkx, std::vector< int > const &pky, std::vector< std::shared_ptr< typename lsst::afw::detection::HeavyFootprint< ImagePixelT, MaskPixelT, VariancePixelT > > > &strays, int strayFluxOptions, double clipStrayFluxFraction)
Splits flux in a given image img, within a given footprint foot, among a number of templates timgs,...
static void _find_stray_flux(lsst::afw::detection::Footprint const &foot, ImagePtrT tsum, MaskedImageT const &img, int strayFluxOptions, std::vector< std::shared_ptr< lsst::afw::detection::Footprint > > tfoots, std::vector< bool > const &ispsf, std::vector< int > const &pkx, std::vector< int > const &pky, double clipStrayFluxFraction, std::vector< std::shared_ptr< typename lsst::afw::detection::HeavyFootprint< ImagePixelT, MaskPixelT, VariancePixelT > > > &strays)
static std::shared_ptr< lsst::afw::detection::Footprint > getSignificantEdgePixels(ImagePtrT, std::shared_ptr< lsst::afw::detection::Footprint >, ImagePixelT threshold)
Returns a list of pixels that are on the edge of the given Footprint sfoot* in image img,...
static std::pair< ImagePtrT, FootprintPtrT > buildSymmetricTemplate(MaskedImageT const &img, lsst::afw::detection::Footprint const &foot, lsst::afw::detection::PeakRecord const &pk, double sigma1, bool minZero, bool patchEdges, bool *patchedEdges)
Given an img, footprint foot, and peak, creates a symmetric template around the peak; produce a Maske...
static void makeMonotonic(ImageT &img, lsst::afw::detection::PeakRecord const &pk)
Given an image mimg and Peak location peak, overwrite mimg so that pixels further from the peak have ...
static std::shared_ptr< lsst::afw::detection::Footprint > symmetrizeFootprint(lsst::afw::detection::Footprint const &foot, int cx, int cy)
Given a Footprint foot and peak cx,cy, returns a Footprint that is symmetric around the peak (with tw...
static void _sum_templates(std::vector< ImagePtrT > timgs, ImagePtrT tsum)
static bool hasSignificantFluxAtEdge(ImagePtrT, std::shared_ptr< lsst::afw::detection::Footprint >, ImagePixelT threshold)
Returns true if the given Footprint sfoot in image img has flux above value thresh at its edge.
Reports attempts to exceed implementation-defined length limits for some classes.
Definition: Runtime.h:76
Reports errors that are due to events beyond the control of the program.
Definition: Runtime.h:104
T end(T... args)
T max(T... args)
T min(T... args)
T move(T... args)
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
T real(const T &x)
Definition: Integrate.h:165
Point< int, 2 > Point2I
Definition: Point.h:321
FastFinder::Iterator Iterator
Definition: FastFinder.cc:178
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
Angle abs(Angle const &a)
Definition: Angle.h:106
T nth_element(T... args)
T push_back(T... args)
T lround(T... args)
T size(T... args)
T upper_bound(T... args)