LSST Applications g0b6bd0c080+a72a5dd7e6,g1182afd7b4+2a019aa3bb,g17e5ecfddb+2b8207f7de,g1d67935e3f+06cf436103,g38293774b4+ac198e9f13,g396055baef+6a2097e274,g3b44f30a73+6611e0205b,g480783c3b1+98f8679e14,g48ccf36440+89c08d0516,g4b93dc025c+98f8679e14,g5c4744a4d9+a302e8c7f0,g613e996a0d+e1c447f2e0,g6c8d09e9e7+25247a063c,g7271f0639c+98f8679e14,g7a9cd813b8+124095ede6,g9d27549199+a302e8c7f0,ga1cf026fa3+ac198e9f13,ga32aa97882+7403ac30ac,ga786bb30fb+7a139211af,gaa63f70f4e+9994eb9896,gabf319e997+ade567573c,gba47b54d5d+94dc90c3ea,gbec6a3398f+06cf436103,gc6308e37c7+07dd123edb,gc655b1545f+ade567573c,gcc9029db3c+ab229f5caf,gd01420fc67+06cf436103,gd877ba84e5+06cf436103,gdb4cecd868+6f279b5b48,ge2d134c3d5+cc4dbb2e3f,ge448b5faa6+86d1ceac1d,gecc7e12556+98f8679e14,gf3ee170dca+25247a063c,gf4ac96e456+ade567573c,gf9f5ea5b4d+ac198e9f13,gff490e6085+8c2580be5c,w.2022.27
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"
10
11using std::lround;
12
13namespace image = lsst::afw::image;
14namespace det = lsst::afw::detection;
16namespace afwGeom = lsst::afw::geom;
17namespace geom = lsst::geom;
18
19template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
21
22template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
23const int
25
26template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
28
29template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
31
32template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
34
35template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
37
38static bool span_compare(afwGeom::Span const & sp1,
39 afwGeom::Span const & sp2) {
40 return (sp1 < sp2);
41}
42
43namespace {
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) {
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
145template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
146void
148medianFilter(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
223template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
224void
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
364static 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
391template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
392void
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
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) {
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
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]) {
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
607template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
608void
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
690template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
693apportionFlux(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
819public:
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
837 return _real == other;
838 }
840 return _real != other;
841 }
843 return _real <= other;
844 }
845 bool operator<(const SpanSet::const_iterator & other) {
846 return _real < other;
847 }
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
913private:
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
972template<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;
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
1186template<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
1383template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
1384bool
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
1416template<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();
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
1457template 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
lsst::geom::Box2I getBBox() const
Return the Footprint's bounding box.
Definition: Footprint.h:208
std::shared_ptr< geom::SpanSet > getSpans() const
Return a shared pointer to the SpanSet.
Definition: Footprint.h:115
PeakCatalog & getPeaks()
Return the Peaks contained in this Footprint.
Definition: Footprint.h:129
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< ImagePixelT >::type image()
Return (a reference to) the image part of the Pixel pointed at by the iterator.
Definition: MaskedImage.h:139
Ref< MaskPixelT >::type mask()
Return (a reference to) the mask part of the Pixel pointed at by the iterator.
Definition: MaskedImage.h:142
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 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::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 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)