LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
Baseline.cc
Go to the documentation of this file.
1 #include <list>
2 
4 #include "lsst/pex/logging.h"
5 #include "lsst/pex/exceptions.h"
6 #include "lsst/afw/geom/Box.h"
7 
8 #include <boost/math/special_functions/round.hpp>
9 
10 using boost::math::iround;
11 
12 namespace image = lsst::afw::image;
13 namespace det = lsst::afw::detection;
14 namespace deblend = lsst::meas::deblender;
15 namespace geom = lsst::afw::geom;
16 namespace pexLog = lsst::pex::logging;
17 
18 
19 static bool span_ptr_compare(PTR(det::Span) sp1, PTR(det::Span) sp2) {
20  return (*sp1 < *sp2);
21 }
22 
36 template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
37 void
39 medianFilter(ImageT const& img,
40  ImageT & out,
41  int halfsize) {
42  int S = halfsize*2 + 1;
43  int SS = S*S;
44  typedef typename ImageT::xy_locator xy_loc;
45  xy_loc pix = img.xy_at(halfsize,halfsize);
46  std::vector<typename xy_loc::cached_location_t> locs;
47  for (int i=0; i<S; ++i) {
48  for (int j=0; j<S; ++j) {
49  locs.push_back(pix.cache_location(j-halfsize, i-halfsize));
50  }
51  }
52  int W = img.getWidth();
53  int H = img.getHeight();
54  ImagePixelT vals[S*S];
55  for (int y=halfsize; y<H-halfsize; ++y) {
56  xy_loc inpix = img.xy_at(halfsize, y), end = img.xy_at(W-halfsize, y);
57  for (typename ImageT::x_iterator optr = out.row_begin(y) + halfsize;
58  inpix != end; ++inpix.x(), ++optr) {
59  for (int i=0; i<SS; ++i)
60  vals[i] = inpix[locs[i]];
61  std::nth_element(vals, vals+SS/2, vals+SS);
62  *optr = vals[SS/2];
63  }
64  }
65 
66  // grumble grumble margins
67  for (int y=0; y<2*halfsize; ++y) {
68  int iy = y;
69  if (y >= halfsize)
70  iy = H - 1 - (y-halfsize);
71  typename ImageT::x_iterator optr = out.row_begin(iy);
72  typename ImageT::x_iterator iptr = img.row_begin(iy), end=img.row_end(iy);
73  for (; iptr != end; ++iptr,++optr)
74  *optr = *iptr;
75  }
76  for (int y=halfsize; y<H-halfsize; ++y) {
77  typename ImageT::x_iterator optr = out.row_begin(y);
78  typename ImageT::x_iterator iptr = img.row_begin(y), end=img.row_begin(y)+halfsize;
79  for (; iptr != end; ++iptr,++optr)
80  *optr = *iptr;
81  iptr = img.row_begin(y) + ((W-1) - halfsize);
82  end = img.row_begin(y) + (W-1);
83  optr = out.row_begin(y) + ((W-1) - halfsize);
84  for (; iptr != end; ++iptr,++optr)
85  *optr = *iptr;
86  }
87 
88 }
89 
114 template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
115 void
118  ImageT & img,
119  det::PeakRecord const& peak) {
120 
121  int cx = peak.getIx();
122  int cy = peak.getIy();
123  int ix0 = img.getX0();
124  int iy0 = img.getY0();
125  int iW = img.getWidth();
126  int iH = img.getHeight();
127 
128  ImagePtrT shadowingImg = ImagePtrT(new ImageT(img, true));
129 
130  int DW = std::max(cx - img.getX0(), img.getX0() + img.getWidth() - cx);
131  int DH = std::max(cy - img.getY0(), img.getY0() + img.getHeight() - cy);
132 
133  const int S = 5;
134 
135  // Work out from the peak in chunks of "S" pixels.
136  int s;
137  for (s = 0; s < std::max(DW,DH); s += S) {
138  int p;
139  for (p=0; p<S; p++) {
140  // visit pixels with L_inf distance = s + p from the
141  // center (ie, the s+p'th square ring of pixels)
142  // L is the half-length of the ring (box).
143  int L = s+p;
144  int x = L, y = -L;
145  int dx = 0, dy = 0; // initialized here to satisfy the
146  // compiler; initialized for real
147  // below (first time through loop)
148  /*
149  int i;
150  int leg;
151  for (i=0; i<(8*L); i++, x += dx, y += dy) {
152  if (i % (2*L) == 0) {
153  leg = (i/(2*L));
154  dx = ( leg % 2) * (-1 + 2*(leg/2));
155  dy = ((leg+1) % 2) * ( 1 - 2*(leg/2));
156  }
157  */
158 
159  /*
160  We visit pixels in a box of "radius" L, in this order:
161 
162  L=1:
163 
164  4 3 2
165  5 1
166  6 7 0
167 
168  L=2:
169 
170  8 7 6 5 4
171  9 3
172  10 2
173  11 1
174  12 13 14 15 0
175 
176  Note that the number of pixel visited is 8*L, and that we
177  change "dx" or "dy" each "2*L" steps.
178  */
179  for (int i=0; i<(8*L); i++, x += dx, y += dy) {
180  // time to change directions? (Note that this runs
181  // the first time through the loop, initializing dx,dy
182  // appropriately.)
183  if (i % (2*L) == 0) {
184  int leg = (i/(2*L));
185  // dx = [ 0, -1, 0, 1 ][leg]
186  dx = ( leg % 2) * (-1 + 2*(leg/2));
187  // dy = [ 1, 0, -1, 0 ][leg]
188  dy = ((leg+1) % 2) * ( 1 - 2*(leg/2));
189  }
190  //printf(" i=%i, leg=%i, dx=%i, dy=%i, x=%i, y=%i\n", i, leg, dx, dy, x, y);
191  int px = cx + x - ix0;
192  int py = cy + y - iy0;
193  // If the shadowing pixel is out of bounds, nothing to do.
194  if (px < 0 || px >= iW || py < 0 || py >= iH)
195  continue;
196  // The pixel casting the shadow
197  ImagePixelT pix = (*shadowingImg)(px,py);
198 
199  // Cast this pixel's shadow S pixels long in a cone.
200  // We compute the range of slopes (or inverse-slopes)
201  // shadowed by the pixel, [ds0,ds1]
202  double ds0, ds1;
203  // Range of slopes shadowed
204  const double A = 0.3;
205  int shx, shy;
206  int psx, psy;
207  // Are we traversing a vertical edge of the box?
208  if (dx == 0) {
209  // (if so, then "x" is +- L, so no div-by-zero)
210  ds0 = (double(y) / double(x)) - A;
211  ds1 = ds0 + 2.0 * A;
212  // cast the shadow on column x + sign(x)*shx
213  for (shx=1; shx<=S; shx++) {
214  int xsign = (x>0?1:-1);
215  // the column being shadowed
216  psx = cx + x + (xsign*shx) - ix0;
217  if (psx < 0 || psx >= iW)
218  continue;
219  // shadow covers a range of y values based on slope
220  for (shy = iround(shx * ds0);
221  shy <= iround(shx * ds1); shy++) {
222  psy = cy + y + xsign*shy - iy0;
223  if (psy < 0 || psy >= iH)
224  continue;
225  img(psx, psy) = std::min(img(psx, psy), pix);
226  }
227  }
228 
229  } else {
230  // We're traversing a horizontal edge of the box; y = +-L
231  ds0 = (double(x) / double(y)) - A;
232  ds1 = ds0 + 2.0 * A;
233  // Cast shadow on row y + sign(y)*shy
234  for (shy=1; shy<=S; shy++) {
235  int ysign = (y>0?1:-1);
236  psy = cy + y + (ysign*shy) - iy0;
237  if (psy < 0 || psy >= iH)
238  continue;
239  // shadow covers a range of x vals based on slope
240  for (shx = iround(shy * ds0);
241  shx <= iround(shy * ds1); shx++) {
242  psx = cx + x + ysign*shx - ix0;
243  if (psx < 0 || psx >= iW)
244  continue;
245  img(psx, psy) = std::min(img(psx, psy), pix);
246  }
247  }
248  }
249  }
250  }
251  //*shadowingImg <<= *img;
252  shadowingImg->operator<<=(img);
253  }
254 }
255 
256 static double _get_contrib_r_to_footprint(int x, int y,
257  PTR(det::Footprint) tfoot) {
258  typedef det::Footprint::SpanList SpanList;
259  double minr2 = 1e12;
260  SpanList const& tspans = tfoot->getSpans();
261  for (SpanList::const_iterator ts = tspans.begin(); ts < tspans.end(); ++ts) {
262  geom::Span* sp = ts->get();
263  int mindx;
264  // span is to right of pixel?
265  int dx = sp->getX0() - x;
266  if (dx >= 0) {
267  mindx = dx;
268  } else {
269  // span is to left of pixel?
270  dx = x - sp->getX1();
271  if (dx >= 0) {
272  mindx = dx;
273  } else {
274  // span contains pixel (in x direction)
275  mindx = 0;
276  }
277  }
278  int dy = sp->getY() - y;
279  minr2 = std::min(minr2, (double)(mindx*mindx + dy*dy));
280  }
281  //printf("stray flux at (%i,%i): dist to t %i is %g\n", x, y, (int)i, sqrt(minr2));
282  return 1. / (1. + minr2);
283 }
284 
285 
286 template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
287 void
290  ImagePtrT tsum,
291  MaskedImageT const& img,
292  int strayFluxOptions,
293  std::vector<PTR(det::Footprint)> tfoots,
294  std::vector<bool> const& ispsf,
295  std::vector<int> const& pkx,
296  std::vector<int> const& pky,
297  double clipStrayFluxFraction,
298  std::vector<boost::shared_ptr<typename det::HeavyFootprint<ImagePixelT,MaskPixelT,VariancePixelT> > > & strays
299  ) {
300 
301  typedef typename det::Footprint::SpanList SpanList;
302  typedef typename det::HeavyFootprint<ImagePixelT, MaskPixelT, VariancePixelT> HeavyFootprint;
303  typedef typename boost::shared_ptr< HeavyFootprint > HeavyFootprintPtrT;
304 
305  // when doing stray flux: the footprints and pixels, which we'll
306  // combine into the return 'strays' HeavyFootprint at the end.
307  std::vector<PTR(det::Footprint) > strayfoot;
308  std::vector<std::vector<ImagePixelT> > straypix;
309  std::vector<std::vector<MaskPixelT> > straymask;
310  std::vector<std::vector<VariancePixelT> > strayvar;
311 
312  int ix0 = img.getX0();
313  int iy0 = img.getY0();
314  geom::Box2I sumbb = tsum->getBBox();
315  int sumx0 = sumbb.getMinX();
316  int sumy0 = sumbb.getMinY();
317 
318  for (size_t i=0; i<tfoots.size(); ++i) {
319  strayfoot.push_back(PTR(det::Footprint)());
320  straypix.push_back(std::vector<ImagePixelT>());
321  straymask.push_back(std::vector<MaskPixelT>());
322  strayvar.push_back(std::vector<VariancePixelT>());
323  }
324 
325  bool always = (strayFluxOptions & STRAYFLUX_TO_POINT_SOURCES_ALWAYS);
326 
327  typedef boost::uint16_t itype;
328  PTR(image::Image<itype>) nearest;
329 
330  if (strayFluxOptions & STRAYFLUX_NEAREST_FOOTPRINT) {
331  // Compute the map of which footprint is closest to each
332  // pixel in the bbox.
333  typedef boost::uint16_t dtype;
334  PTR(image::Image<dtype>) dist(new image::Image<dtype>(sumbb));
335  nearest = PTR(image::Image<itype>)(new image::Image<itype>(sumbb));
336 
337  std::vector<PTR(det::Footprint)> templist;
338  std::vector<PTR(det::Footprint)>* footlist = &tfoots;
339 
340  if (!always && ispsf.size()) {
341  // create a temp list that has empty footprints in place
342  // of all the point sources.
343  PTR(det::Footprint) empty(new det::Footprint(foot.getPeaks().getSchema()));
344  for (size_t i=0; i<tfoots.size(); ++i) {
345  if (ispsf[i]) {
346  templist.push_back(empty);
347  } else {
348  templist.push_back(tfoots[i]);
349  }
350  }
351  footlist = &templist;
352  }
353  nearestFootprint(*footlist, nearest, dist);
354  }
355 
356  // Go through the (parent) Footprint looking for stray flux:
357  // pixels that are not claimed by any template, and positive.
358  const SpanList& spans = foot.getSpans();
359  for (SpanList::const_iterator s = spans.begin(); s < spans.end(); s++) {
360  int y = (*s)->getY();
361  int x0 = (*s)->getX0();
362  int x1 = (*s)->getX1();
363  typename ImageT::x_iterator tsum_it =
364  tsum->row_begin(y - sumy0) + (x0 - sumx0);
365  typename MaskedImageT::x_iterator in_it =
366  img.row_begin(y - iy0) + (x0 - ix0);
367  double contrib[tfoots.size()];
368 
369  for (int x = x0; x <= x1; ++x, ++tsum_it, ++in_it) {
370  // Skip pixels that are covered by at least one
371  // template (*tsum_it > 0) or the input is not
372  // positive (*in_it <= 0).
373  if ((*tsum_it > 0) || (*in_it).image() <= 0) {
374  continue;
375  }
376 
377  if (strayFluxOptions & STRAYFLUX_R_TO_FOOTPRINT) {
378  // we'll compute these just-in-time
379  for (size_t i=0; i<tfoots.size(); ++i) {
380  contrib[i] = -1.0;
381  }
382  } else if (strayFluxOptions & STRAYFLUX_NEAREST_FOOTPRINT) {
383  for (size_t i=0; i<tfoots.size(); ++i) {
384  contrib[i] = 0.0;
385  }
386  int i = nearest->get0(x, y);
387  contrib[i] = 1.0;
388  } else {
389  // R_TO_PEAK
390  for (size_t i=0; i<tfoots.size(); ++i) {
391  // Split the stray flux by 1/(1+r^2) to peaks
392  int dx, dy;
393  dx = pkx[i] - x;
394  dy = pky[i] - y;
395  contrib[i] = 1. / (1. + dx*dx + dy*dy);
396  }
397  }
398 
399  // Round 1: skip point sources unless STRAYFLUX_TO_POINT_SOURCES_ALWAYS
400  // are we going to assign stray flux to ptsrcs?
401  bool ptsrcs = always;
402  double csum = 0.;
403  for (size_t i=0; i<tfoots.size(); ++i) {
404  // if we're skipping point sources and this is a point source...
405  if ((!ptsrcs) && ispsf.size() && ispsf[i]) {
406  continue;
407  }
408  if (contrib[i] == -1.0) {
409  contrib[i] = _get_contrib_r_to_footprint(x, y, tfoots[i]);
410  }
411  csum += contrib[i];
412  }
413  if ((csum == 0.) &&
414  (strayFluxOptions &
415  STRAYFLUX_TO_POINT_SOURCES_WHEN_NECESSARY)) {
416  // No extended sources -- assign to pt sources
417  //log.debugf("necessary to assign stray flux to point sources");
418  ptsrcs = true;
419  for (size_t i=0; i<tfoots.size(); ++i) {
420  if (contrib[i] == -1.0) {
421  contrib[i] = _get_contrib_r_to_footprint(x, y, tfoots[i]);
422  }
423  csum += contrib[i];
424  }
425  }
426 
427  // Drop small contributions...
428  double strayclip = (clipStrayFluxFraction * csum);
429  csum = 0.;
430  for (size_t i=0; i<tfoots.size(); ++i) {
431  // skip ptsrcs?
432  if ((!ptsrcs) && ispsf.size() && ispsf[i]) {
433  contrib[i] = 0.;
434  continue;
435  }
436  // skip small contributions
437  if (contrib[i] < strayclip) {
438  contrib[i] = 0.;
439  continue;
440  }
441  csum += contrib[i];
442  }
443 
444  for (size_t i=0; i<tfoots.size(); ++i) {
445  if (contrib[i] == 0.) {
446  continue;
447  }
448  // the stray flux to give to template i
449  double p = (contrib[i] / csum) * (*in_it).image();
450  if (!strayfoot[i]) {
451  strayfoot[i] = PTR(det::Footprint)(new det::Footprint(foot.getPeaks().getSchema()));
452  }
453  strayfoot[i]->addSpanInSeries(y, x, x); // XXX this may be rather heavy in memory use
454  straypix[i].push_back(p);
455  straymask[i].push_back((*in_it).mask());
456  strayvar[i].push_back((*in_it).variance());
457  }
458  }
459  }
460 
461  // Store the stray flux in HeavyFootprints
462  for (size_t i=0; i<tfoots.size(); ++i) {
463  if (!strayfoot[i]) {
464  strays.push_back(HeavyFootprintPtrT());
465  } else {
469  HeavyFootprintPtrT heavy(new HeavyFootprint(*strayfoot[i]));
470  ndarray::Array<ImagePixelT,1,1> himg = heavy->getImageArray();
471  typename std::vector<ImagePixelT>::const_iterator spix;
472  typename std::vector<MaskPixelT>::const_iterator smask;
473  typename std::vector<VariancePixelT>::const_iterator svar;
477 
478  assert((size_t)strayfoot[i]->getNpix() == straypix[i].size());
479 
480  for (spix = straypix[i].begin(),
481  smask = straymask[i].begin(),
482  svar = strayvar [i].begin(),
483  hpix = himg.begin(),
484  mpix = heavy->getMaskArray().begin(),
485  vpix = heavy->getVarianceArray().begin();
486  spix != straypix[i].end();
487  ++spix, ++smask, ++svar, ++hpix, ++mpix, ++vpix) {
488  *hpix = *spix;
489  *mpix = *smask;
490  *vpix = *svar;
491  }
492  strays.push_back(heavy);
493  }
494  }
495 }
496 
497 template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
498 void
500 _sum_templates(std::vector<ImagePtrT> timgs,
501  ImagePtrT tsum) {
502  geom::Box2I sumbb = tsum->getBBox();
503  int sumx0 = sumbb.getMinX();
504  int sumy0 = sumbb.getMinY();
505 
506  // Compute tsum = the sum of templates
507  for (size_t i=0; i<timgs.size(); ++i) {
508  ImagePtrT timg = timgs[i];
509  geom::Box2I tbb = timg->getBBox();
510  int tx0 = tbb.getMinX();
511  int ty0 = tbb.getMinY();
512  // To handle "ramped" templates that can extend outside the
513  // parent, clip the bbox. Note that we saved tx0,ty0 BEFORE
514  // doing this!
515  tbb.clip(sumbb);
516  int copyx0 = tbb.getMinX();
517  // Here we iterate over the template bbox -- we could instead
518  // iterate over the "tfoot"s.
519  for (int y=tbb.getMinY(); y<=tbb.getMaxY(); ++y) {
520  typename ImageT::x_iterator in_it = timg->row_begin(y - ty0) + (copyx0 - tx0);
521  typename ImageT::x_iterator inend = in_it + tbb.getWidth();
522  typename ImageT::x_iterator tsum_it =
523  tsum->row_begin(y - sumy0) + (copyx0 - sumx0);
524  for (; in_it != inend; ++in_it, ++tsum_it) {
525  *tsum_it += std::max((ImagePixelT)0., static_cast<ImagePixelT>(*in_it));
526  }
527  }
528  }
529 
530 }
531 
580 template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
581 std::vector<typename PTR(image::MaskedImage<ImagePixelT, MaskPixelT, VariancePixelT>)>
583 apportionFlux(MaskedImageT const& img,
584  det::Footprint const& foot,
585  std::vector<ImagePtrT> timgs,
586  std::vector<PTR(det::Footprint)> tfoots,
587  ImagePtrT tsum,
588  std::vector<bool> const& ispsf,
589  std::vector<int> const& pkx,
590  std::vector<int> const& pky,
591  std::vector<boost::shared_ptr<typename det::HeavyFootprint<ImagePixelT,MaskPixelT,VariancePixelT> > > & strays,
592  int strayFluxOptions,
593  double clipStrayFluxFraction
594  ) {
595  typedef typename det::Footprint::SpanList SpanList;
596 
597  if (timgs.size() != tfoots.size()) {
598  throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
599  (boost::format("Template images must be the same length as template footprints (%d vs %d)")
600  % timgs.size() % tfoots.size()).str());
601  }
602 
603  for (size_t i=0; i<timgs.size(); ++i) {
604  if (!timgs[i]->getBBox().contains(tfoots[i]->getBBox())) {
605  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
606  "Template image MUST contain template footprint");
607  }
608  if (!img.getBBox().contains(foot.getBBox())) {
609  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
610  "Image bbox MUST contain parent footprint");
611  }
612  // template bounding-boxes *can* extend outside the parent
613  // footprint if we are ramping templates with significant flux
614  // at the edges. We handle this below.
615  }
616 
617  // the apportioned flux return value
618  std::vector<MaskedImagePtrT> portions;
619 
621  "lsst.meas.deblender.apportionFlux");
622  bool findStrayFlux = (strayFluxOptions & ASSIGN_STRAYFLUX);
623 
624  int ix0 = img.getX0();
625  int iy0 = img.getY0();
626  geom::Box2I fbb = foot.getBBox();
627 
628  if (!tsum) {
629  tsum = ImagePtrT(new ImageT(fbb.getDimensions()));
630  tsum->setXY0(fbb.getMinX(), fbb.getMinY());
631  }
632 
633  if (!tsum->getBBox().contains(foot.getBBox())) {
634  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
635  "Template sum image MUST contain parent footprint");
636  }
637 
638  geom::Box2I sumbb = tsum->getBBox();
639  int sumx0 = sumbb.getMinX();
640  int sumy0 = sumbb.getMinY();
641 
642  _sum_templates(timgs, tsum);
643 
644  // Compute flux portions
645  for (size_t i=0; i<timgs.size(); ++i) {
646  ImagePtrT timg = timgs[i];
647  // Initialize return value:
648  MaskedImagePtrT port(new MaskedImageT(timg->getDimensions()));
649  port->setXY0(timg->getXY0());
650  portions.push_back(port);
651 
652  // Split flux = image * template / tsum
653  geom::Box2I tbb = timg->getBBox();
654  int tx0 = tbb.getMinX();
655  int ty0 = tbb.getMinY();
656  // As above
657  tbb.clip(sumbb);
658  int copyx0 = tbb.getMinX();
659  for (int y=tbb.getMinY(); y<=tbb.getMaxY(); ++y) {
660  typename MaskedImageT::x_iterator in_it =
661  img.row_begin(y - iy0) + (copyx0 - ix0);
662  typename ImageT::x_iterator tptr =
663  timg->row_begin(y - ty0) + (copyx0 - tx0);
664  typename ImageT::x_iterator tend = tptr + tbb.getWidth();
665  typename ImageT::x_iterator tsum_it =
666  tsum->row_begin(y - sumy0) + (copyx0 - sumx0);
667  typename MaskedImageT::x_iterator out_it =
668  port->row_begin(y - ty0) + (copyx0 - tx0);
669  for (; tptr != tend; ++tptr, ++in_it, ++out_it, ++tsum_it) {
670  if (*tsum_it == 0) {
671  continue;
672  }
673  double frac = std::max((ImagePixelT)0., static_cast<ImagePixelT>(*tptr)) / (*tsum_it);
674  //if (frac == 0) {
675  // treat mask planes differently?
676  // }
677  out_it.mask() = (*in_it).mask();
678  out_it.variance() = (*in_it).variance()*frac;
679  out_it.image() = (*in_it).image() * frac;
680  }
681  }
682  }
683 
684  if (findStrayFlux) {
685  if ((ispsf.size() > 0) && (ispsf.size() != timgs.size())) {
686  throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
687  (boost::format("'ispsf' must be the same length as templates (%d vs %d)")
688  % ispsf.size() % timgs.size()).str());
689  }
690  if ((pkx.size() != timgs.size()) || (pky.size() != timgs.size())) {
691  throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
692  (boost::format("'pkx' and 'pky' must be the same length as templates (%d,%d vs %d)")
693  % pkx.size() % pky.size() % timgs.size()).str());
694  }
695  _find_stray_flux(foot, tsum, img, strayFluxOptions, tfoots,
696  ispsf, pkx, pky, clipStrayFluxFraction, strays);
697  }
698  return portions;
699 }
700 
701 
711 public:
713 
715 
716  RelativeSpanIterator(SpanList::const_iterator const& real,
717  SpanList const& arr,
718  int cx, int cy, bool forward=true)
719  : _real(real), _cx(cx), _cy(cy), _forward(forward)
720  {
721  if (_forward) {
722  _end = arr.end();
723  } else {
724  _end = arr.begin();
725  }
726  }
727 
728  bool operator==(const SpanList::const_iterator & other) {
729  return _real == other;
730  }
731  bool operator!=(const SpanList::const_iterator & other) {
732  return _real != other;
733  }
734  bool operator<=(const SpanList::const_iterator & other) {
735  return _real <= other;
736  }
737  bool operator<(const SpanList::const_iterator & other) {
738  return _real < other;
739  }
740  bool operator>=(const SpanList::const_iterator & other) {
741  return _real >= other;
742  }
743  bool operator>(const SpanList::const_iterator & other) {
744  return _real > other;
745  }
746 
748  return (_real == other._real) &&
749  (_cx == other._cx) && (_cy == other._cy) &&
750  (_forward == other._forward);
751  }
753  return !(*this == other);
754  }
755 
757  if (_forward) {
758  _real++;
759  } else {
760  _real--;
761  }
762  return *this;
763  }
765  RelativeSpanIterator result = *this;
766  ++(*this);
767  return result;
768  }
769 
770  bool notDone() {
771  if (_forward) {
772  return _real < _end;
773  } else {
774  return _real >= _end;
775  }
776  }
777 
778  int dxlo() {
779  if (_forward) {
780  return (*_real)->getX0() - _cx;
781  } else {
782  return _cx - (*_real)->getX1();
783  }
784  }
785  int dxhi() {
786  if (_forward) {
787  return (*_real)->getX1() - _cx;
788  } else {
789  return _cx - (*_real)->getX0();
790  }
791  }
792  int dy() {
793  return std::abs((*_real)->getY() - _cy);
794  }
795  int x0() {
796  return (*_real)->getX0();
797  }
798  int x1() {
799  return (*_real)->getX1();
800  }
801  int y() {
802  return (*_real)->getY();
803  }
804 
805 private:
806  SpanList::const_iterator _real;
807  SpanList::const_iterator _end;
808  int _cx;
809  int _cy;
810  bool _forward;
811 };
812 
813 
814 /*
815  // Check symmetrizeFootprint by computing truth naively.
816  // compute correct answer dumbly
817  det::Footprint truefoot;
818  geom::Box2I bbox = foot.getBBox();
819  for (int y=bbox.getMinY(); y<=bbox.getMaxY(); y++) {
820  for (int x=bbox.getMinX(); x<=bbox.getMaxX(); x++) {
821  int dy = y - cy;
822  int dx = x - cx;
823  int x2 = cx - dx;
824  int y2 = cy - dy;
825  if (foot.contains(geom::Point2I(x, y)) &&
826  foot.contains(geom::Point2I(x2, y2))) {
827  truefoot.addSpan(y, x, x);
828  }
829  }
830  }
831  truefoot.normalize();
832 
833  SpanList sp1 = truefoot.getSpans();
834  SpanList sp2 = sfoot->getSpans();
835  for (SpanList::const_iterator spit1 = sp1.begin(),
836  spit2 = sp2.begin();
837  spit1 != sp1.end() && spit2 != sp2.end();
838  spit1++, spit2++) {
839  //printf("\n");
840  printf(" true y %i, x [%i, %i]\n", (*spit1)->getY(), (*spit1)->getX0(), (*spit1)->getX1());
841  printf(" sfoot y %i, x [%i, %i]\n", (*spit2)->getY(), (*spit2)->getX0(), (*spit2)->getX1());
842  if (((*spit1)->getY() != (*spit2)->getY()) ||
843  ((*spit1)->getX0() != (*spit2)->getX0()) ||
844  ((*spit1)->getX1() != (*spit2)->getX1())) {
845  printf("*******************************************\n");
846  }
847  }
848  assert(sp1.size() == sp2.size());
849  for (SpanList::const_iterator spit1 = sp1.begin(),
850  spit2 = sp2.begin();
851  spit1 != sp1.end() && spit2 != sp2.end();
852  spit1++, spit2++) {
853  assert((*spit1)->getY() == (*spit2)->getY());
854  assert((*spit1)->getX0() == (*spit2)->getX0());
855  assert((*spit1)->getX1() == (*spit2)->getX1());
856  }
857  */
858 
864 template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
866 deblend::BaselineUtils<ImagePixelT,MaskPixelT,VariancePixelT>::
867 symmetrizeFootprint(
868  det::Footprint const& foot,
869  int cx, int cy) {
870 
871  typedef typename det::Footprint::SpanList SpanList;
872 
873  PTR(det::Footprint) sfoot(new det::Footprint(foot.getPeaks().getSchema()));
874  const SpanList spans = foot.getSpans();
875  assert(foot.isNormalized());
876 
878  "lsst.meas.deblender.symmetrizeFootprint");
879 
880  // Find the Span containing the peak.
881  PTR(det::Span) target(new det::Span(cy, cx, cx));
882  SpanList::const_iterator peakspan =
883  std::upper_bound(spans.begin(), spans.end(), target, span_ptr_compare);
884  // upper_bound returns the last position where "target" could be inserted;
885  // ie, the first Span larger than "target". The Span containing "target"
886  // should be peakspan-1 or (if the peak is on the first pixel in the span,
887  // peakspan.
888  PTR(det::Span) sp;
889  if (peakspan == spans.begin()) {
890  sp = *peakspan;
891  if (!sp->contains(cx, cy)) {
892  log.warnf(
893  "Failed to find span containing (%i,%i): before the beginning of this footprint", cx, cy);
894  return PTR(det::Footprint)();
895  }
896  } else {
897  --peakspan;
898  sp = *peakspan;
899 
900  if (!sp->contains(cx, cy)) {
901  ++peakspan;
902  sp = *peakspan;
903  if (!sp->contains(cx, cy)) {
904  geom::Box2I fbb = foot.getBBox();
905  log.warnf("Failed to find span containing (%i,%i): nearest is %i, [%i,%i]. "
906  "Footprint bbox is [%i,%i],[%i,%i]",
907  cx, cy, sp->getY(), sp->getX0(), sp->getX1(),
908  fbb.getMinX(), fbb.getMaxX(), fbb.getMinY(), fbb.getMaxY());
909  return PTR(det::Footprint)();
910  }
911  }
912  }
913  log.debugf("Span containing (%i,%i): (x=[%i,%i], y=%i)",
914  cx, cy, sp->getX0(), sp->getX1(), sp->getY());
915 
916  // The symmetric templates are essentially an AND of the footprint
917  // pixels and its 180-degree-rotated self, rotated around the
918  // peak (cx,cy).
919  //
920  // We iterate forward and backward simultaneously, starting from
921  // the span containing the peak and moving out, row by row.
922  //
923  // In the loop below, we search for the next pair of Spans that
924  // overlap (in "dx" from the center), output the overlapping
925  // portion of the Spans, and advance either the "fwd" or "back"
926  // iterator. When we fail to find an overlapping pair of Spans,
927  // we move on to the next row.
928  //
929  // [The following paragraph is somewhat obsoleted by the
930  // RelativeSpanIterator class, which performs some of the renaming
931  // and the dx,dy coords.]
932  //
933  // '''In reading the code, "forward", "advancing", etc, are all
934  // from the perspective of the "fwd" iterator (the one going
935  // forward through the Span list, from low to high Y and then low
936  // to high X). It will help to imagine making a copy of the
937  // footprint and rotating it around the center pixel by 180
938  // degrees, so that "fwd" and "back" are both iterating the same
939  // direction; we're then just finding the AND of those two
940  // iterators, except we have to work in dx,dy coordinates rather
941  // than original x,y coords, and the accessors for "back" are
942  // opposite.'''
943 
944  RelativeSpanIterator fwd (peakspan, spans, cx, cy, true);
945  RelativeSpanIterator back(peakspan, spans, cx, cy, false);
946 
947  int dy = 0;
948  while (fwd.notDone() && back.notDone()) {
949  // forward and backward "y"; just symmetric around cy
950  int fy = cy + dy;
951  int by = cy - dy;
952  // delta-x of the beginnings of the spans, for "fwd" and "back"
953  int fdxlo = fwd.dxlo();
954  int bdxlo = back.dxlo();
955 
956  // First find:
957  // fend -- first span in the next row, or end(); ie,
958  // the end of this row in the forward direction
959  // bend -- the end of this row in the backward direction
960  RelativeSpanIterator fend, bend;
961  for (fend = fwd; fend.notDone(); ++fend) {
962  if (fend.dy() != dy)
963  break;
964  }
965  for (bend = back; bend.notDone(); ++bend) {
966  if (bend.dy() != dy)
967  break;
968  }
969 
970  log.debugf("dy=%i, fy=%i, fx=[%i, %i], by=%i, fx=[%i, %i], fdx=%i, bdx=%i",
971  dy, fy, fwd.x0(), fwd.x1(), by, back.x0(), back.x1(),
972  fdxlo, bdxlo);
973 
974  // Find possibly-overlapping span
975  if (bdxlo > fdxlo) {
976  log.debugf("Advancing forward.");
977  // While the "forward" span is entirely to the "left" of the "backward" span,
978  // (in dx coords), ie, |---fwd---X X---back---|
979  // and we are comparing the edges marked X
980  while ((fwd != fend) && (fwd.dxhi() < bdxlo)) {
981  fwd++;
982  if (fwd == fend) {
983  log.debugf("Reached fend");
984  } else {
985  log.debugf("Advanced to forward span %i, [%i, %i]",
986  fy, fwd.x0(), fwd.x1());
987  }
988  }
989  } else if (fdxlo > bdxlo) {
990  log.debugf("Advancing backward.");
991  // While the "backward" span is entirely to the "left" of the "foreward" span,
992  // (in dx coords), ie, |---back---X X---fwd---|
993  // and we are comparing the edges marked X
994  while ((back != bend) && (back.dxhi() < fdxlo)) {
995  back++;
996  if (back == bend) {
997  log.debugf("Reached bend");
998  } else {
999  log.debugf("Advanced to backward span %i, [%i, %i]",
1000  by, back.x0(), back.x1());
1001  }
1002  }
1003  }
1004 
1005  if ((back == bend) || (fwd == fend)) {
1006  // We reached the end of the row without finding spans that could
1007  // overlap. Move onto the next dy.
1008  if (back == bend) {
1009  log.debugf("Reached bend");
1010  }
1011  if (fwd == fend) {
1012  log.debugf("Reached fend");
1013  }
1014  back = bend;
1015  fwd = fend;
1016  dy++;
1017  continue;
1018  }
1019 
1020  // Spans may overlap -- find the overlapping part.
1021  int dxlo = std::max(fwd.dxlo(), back.dxlo());
1022  int dxhi = std::min(fwd.dxhi(), back.dxhi());
1023  if (dxlo <= dxhi) {
1024  log.debugf("Adding span fwd %i, [%i, %i], back %i, [%i, %i]",
1025  fy, cx+dxlo, cx+dxhi, by, cx-dxhi, cx-dxlo);
1026  sfoot->addSpan(fy, cx + dxlo, cx + dxhi);
1027  sfoot->addSpan(by, cx - dxhi, cx - dxlo);
1028  }
1029 
1030  // Advance the one whose "hi" edge is smallest
1031  if (fwd.dxhi() < back.dxhi()) {
1032  fwd++;
1033  if (fwd == fend) {
1034  log.debugf("Stepped to fend");
1035  } else {
1036  log.debugf("Stepped forward to span %i, [%i, %i]",
1037  fwd.y(), fwd.x0(), fwd.x1());
1038  }
1039  } else {
1040  back++;
1041  if (back == bend) {
1042  log.debugf("Stepped to bend");
1043  } else {
1044  log.debugf("Stepped backward to span %i, [%i, %i]",
1045  back.y(), back.x0(), back.x1());
1046  }
1047  }
1048 
1049  if ((back == bend) || (fwd == fend)) {
1050  // Reached the end of the row. On to the next dy!
1051  if (back == bend) {
1052  log.debugf("Reached bend");
1053  }
1054  if (fwd == fend) {
1055  log.debugf("Reached fend");
1056  }
1057  back = bend;
1058  fwd = fend;
1059  dy++;
1060  continue;
1061  }
1062 
1063  }
1064  sfoot->normalize();
1065  return sfoot;
1066 }
1067 
1080 template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
1081 std::pair<typename PTR(lsst::afw::image::Image<ImagePixelT>),
1085  MaskedImageT const& img,
1086  det::Footprint const& foot,
1087  det::PeakRecord const& peak,
1088  double sigma1,
1089  bool minZero,
1090  bool patchEdge,
1091  bool* patchedEdges) {
1092 
1093  typedef typename MaskedImageT::const_xy_locator xy_loc;
1094  typedef typename det::Footprint::SpanList SpanList;
1095 
1096  *patchedEdges = false;
1097 
1098  int cx = peak.getIx();
1099  int cy = peak.getIy();
1100 
1102  "lsst.meas.deblender.symmetricFootprint");
1103 
1104  if (!img.getBBox(image::PARENT).contains(foot.getBBox())) {
1105  throw LSST_EXCEPT(lsst::pex::exceptions::LengthError, "Image too small for footprint");
1106  }
1107 
1108  FootprintPtrT sfoot = symmetrizeFootprint(foot, cx, cy);
1109  if (!sfoot) {
1110  return std::pair<ImagePtrT, FootprintPtrT>(ImagePtrT(), sfoot);
1111  }
1112 
1113  if (!img.getBBox(image::PARENT).contains(sfoot->getBBox())) {
1114  throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
1115  "Image too small for symmetrized footprint");
1116  }
1117  const SpanList spans = sfoot->getSpans();
1118 
1119  // does this footprint touch an EDGE?
1120  bool touchesEdge = false;
1121  if (patchEdge) {
1122  log.debugf("Checking footprint for EDGE bits");
1123  MaskPtrT mask = img.getMask();
1124  bool edge = false;
1125  MaskPixelT edgebit = mask->getPlaneBitMask("EDGE");
1126  for (SpanList::const_iterator fwd=spans.begin();
1127  fwd != spans.end(); ++fwd) {
1128  int x0 = (*fwd)->getX0();
1129  int x1 = (*fwd)->getX1();
1130  typename MaskT::x_iterator xiter =
1131  mask->x_at(x0 - mask->getX0(), (*fwd)->getY() - mask->getY0());
1132  for (int x=x0; x<=x1; ++x, ++xiter) {
1133  if ((*xiter) & edgebit) {
1134  edge = true;
1135  break;
1136  }
1137  }
1138  if (edge)
1139  break;
1140  }
1141  if (edge) {
1142  log.debugf("Footprint includes an EDGE pixel.");
1143  touchesEdge = true;
1144  }
1145  }
1146 
1147  // The result image:
1148  ImagePtrT targetimg(new ImageT(sfoot->getBBox()));
1149 
1150  SpanList::const_iterator fwd = spans.begin();
1151  SpanList::const_iterator back = spans.end()-1;
1152 
1153  ImagePtrT theimg = img.getImage();
1154 
1155  for (; fwd <= back; fwd++, back--) {
1156  int fy = (*fwd)->getY();
1157  int by = (*back)->getY();
1158 
1159  for (int fx=(*fwd)->getX0(), bx=(*back)->getX1();
1160  fx <= (*fwd)->getX1();
1161  fx++, bx--) {
1162  // FIXME -- CURRENTLY WE IGNORE THE MASK PLANE! options
1163  // include ORing the mask bits, or being clever about
1164  // ignoring some masked pixels, or copying the mask bits
1165  // of the min pixel
1166 
1167  // We have already checked the bounding box, so this should always be satisfied
1168  assert(theimg->getBBox(image::PARENT).contains(geom::Point2I(fx, fy)));
1169  assert(theimg->getBBox(image::PARENT).contains(geom::Point2I(bx, by)));
1170 
1171  // FIXME -- we could do this with image iterators instead.
1172  // But first profile to show that it's necessary and an
1173  // improvement.
1174  ImagePixelT pixf = theimg->get0(fx, fy);
1175  ImagePixelT pixb = theimg->get0(bx, by);
1176  ImagePixelT pix = std::min(pixf, pixb);
1177  if (minZero) {
1178  pix = std::max(pix, static_cast<ImagePixelT>(0));
1179  }
1180  targetimg->set0(fx, fy, pix);
1181  targetimg->set0(bx, by, pix);
1182 
1183  }
1184  }
1185 
1186 
1187  if (touchesEdge) {
1188  // Find spans whose mirrors fall outside the image bounds,
1189  // grow the footprint to include those spans, and plug in
1190  // their pixel values.
1191  geom::Box2I bb = sfoot->getBBox();
1192 
1193  // Actually, it's not necessarily the IMAGE bounds that count
1194  //-- the footprint may not go right to the image edge.
1195  //geom::Box2I imbb = img.getBBox();
1196  geom::Box2I imbb = foot.getBBox();
1197 
1198  log.debugf("Footprint touches EDGE: start bbox [%i,%i],[%i,%i]",
1199  bb.getMinX(), bb.getMaxX(), bb.getMinY(), bb.getMaxY());
1200  // original footprint spans
1201  const SpanList ospans = foot.getSpans();
1202  for (fwd = ospans.begin(); fwd != ospans.end(); ++fwd) {
1203  int y = (*fwd)->getY();
1204  int x = (*fwd)->getX0();
1205  // mirrored coords
1206  int ym = cy + (cy - y);
1207  int xm = cx + (cx - x);
1208  if (!imbb.contains(geom::Point2I(xm, ym))) {
1209  bb.include(geom::Point2I(x, y));
1210  }
1211  x = (*fwd)->getX1();
1212  xm = cx + (cx - x);
1213  if (!imbb.contains(geom::Point2I(xm, ym))) {
1214  bb.include(geom::Point2I(x, y));
1215  }
1216  }
1217  log.debugf("Footprint touches EDGE: grown bbox [%i,%i],[%i,%i]",
1218  bb.getMinX(), bb.getMaxX(), bb.getMinY(), bb.getMaxY());
1219 
1220  // New template image
1221  ImagePtrT targetimg2(new ImageT(bb));
1222  copyWithinFootprint(*sfoot, targetimg, targetimg2);
1223 
1224  log.debugf("Symmetric footprint spans:");
1225  const SpanList sspans = sfoot->getSpans();
1226  for (fwd = sspans.begin(); fwd != sspans.end(); ++fwd) {
1227  log.debugf(" %s", (*fwd)->toString().c_str());
1228  }
1229 
1230  // copy original 'img' pixels for the portion of spans whose
1231  // mirrors are out of bounds.
1232  for (fwd = ospans.begin(); fwd != ospans.end(); ++fwd) {
1233  int y = (*fwd)->getY();
1234  int x0 = (*fwd)->getX0();
1235  int x1 = (*fwd)->getX1();
1236  // mirrored coords
1237  int ym = cy + (cy - y);
1238  int xm0 = cx + (cx - x0);
1239  int xm1 = cx + (cx - x1);
1240  bool in0 = imbb.contains(geom::Point2I(xm0, ym));
1241  bool in1 = imbb.contains(geom::Point2I(xm1, ym));
1242  if (in0 && in1) {
1243  // both endpoints of the symmetric span are in bounds; nothing to do
1244  continue;
1245  }
1246  // clip to the part of the span where the mirror is out of bounds
1247  if (in0) {
1248  // the mirror of x0 is in-bounds; move x0 to be the first pixel
1249  // whose mirror would be out-of-bounds
1250  x0 = cx + (cx - (imbb.getMinX() - 1));
1251  }
1252  if (in1) {
1253  x1 = cx + (cx - (imbb.getMaxX() + 1));
1254  }
1255  log.debugf("Span y=%i, x=[%i,%i] has mirror (%i,[%i,%i]) out-of-bounds; clipped to %i,[%i,%i]",
1256  y, (*fwd)->getX0(), (*fwd)->getX1(), ym, xm1, xm0, y, x0, x1);
1257  typename MaskedImageT::x_iterator initer =
1258  img.x_at(x0 - img.getX0(), y - img.getY0());
1259  typename ImageT::x_iterator outiter =
1260  targetimg2->x_at(x0 - targetimg2->getX0(), y - targetimg2->getY0());
1261  for (int x=x0; x<=x1; ++x, ++outiter, ++initer) {
1262  *outiter = initer.image();
1263  }
1264  sfoot->addSpan(y, x0, x1);
1265  }
1266  sfoot->normalize();
1267  targetimg = targetimg2;
1268  }
1269 
1270  *patchedEdges = touchesEdge;
1271  return std::pair<ImagePtrT, FootprintPtrT>(targetimg, sfoot);
1272 }
1273 
1278 template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
1279 bool
1282  PTR(det::Footprint) sfoot,
1283  ImagePixelT thresh) {
1284  typedef typename det::Footprint::SpanList SpanList;
1285 
1287  "lsst.meas.deblender.hasSignificantFluxAtEdge");
1288 
1289  // Find edge template pixels with significant flux -- perhaps
1290  // because their symmetric pixels were outside the footprint?
1291  // (clipped by an image edge, etc)
1292  PTR(det::Footprint) edge = sfoot->findEdgePixels();
1293 
1294  const SpanList spans = edge->getSpans();
1295  for (SpanList::const_iterator sp = spans.begin(); sp != spans.end(); ++sp) {
1296  int const y = (*sp)->getY();
1297  int const x0 = (*sp)->getX0();
1298  int const x1 = (*sp)->getX1();
1299  int x;
1300  typename ImageT::const_x_iterator xiter;
1301  for (xiter = img->x_at(x0 - img->getX0(), y - img->getY0()), x=x0; x<=x1; ++x, ++xiter) {
1302  if (*xiter >= thresh) {
1303  return true;
1304  }
1305  }
1306  }
1307  return false;
1308 }
1309 
1314 template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
1315 boost::shared_ptr<det::Footprint>
1318  PTR(det::Footprint) sfoot,
1319  ImagePixelT thresh) {
1321  "lsst.meas.deblender.getSignificantEdgePixels");
1322 
1323  sfoot->normalize(); // Algorithms require spans sorted by y
1324 
1325  PTR(det::Footprint) significant(new det::Footprint(sfoot->getPeaks().getSchema()));
1326 
1327  int const x0 = img->getX0(), y0 = img->getY0();
1328  CONST_PTR(det::Footprint) const edges = sfoot->findEdgePixels();
1329  for (det::Footprint::SpanList::const_iterator ss = edges->getSpans().begin();
1330  ss != edges->getSpans().end(); ++ss) {
1331  geom::Span const& span = **ss;
1332  int const y = span.getY();
1333  int x = span.getX0();
1334  typename ImageT::const_x_iterator iter = img->x_at(x - x0, y - y0);
1335  bool onSpan = false; // Are we in a span of interest
1336  int xSpan; // Starting x of span
1337  for (; x <= span.getX1(); ++x, ++iter) {
1338  if (*iter >= thresh) {
1339  onSpan = true;
1340  xSpan = x;
1341  } else if (onSpan) {
1342  onSpan = false;
1343  significant->addSpanInSeries(y, xSpan, x - 1);
1344  }
1345  }
1346  if (onSpan) {
1347  significant->addSpanInSeries(y, xSpan, span.getX1());
1348  }
1349  }
1350 
1351  return significant;
1352 }
1353 
1354 
1355 // Instantiate
1356 template class deblend::BaselineUtils<float>;
int y
int iter
void nearestFootprint(std::vector< boost::shared_ptr< Footprint >> const &foots, lsst::afw::image::Image< boost::uint16_t >::Ptr argmin, lsst::afw::image::Image< boost::uint16_t >::Ptr dist)
ImagePtr getImage(bool const noThrow=false) const
Return a (Ptr to) the MaskedImage&#39;s image.
Definition: MaskedImage.h:869
void clip(Box2I const &other)
Shrink this to ensure that other.contains(*this).
bool operator>=(const SpanList::const_iterator &other)
Definition: Baseline.cc:740
#define PTR(...)
Definition: base.h:41
x_iterator row_begin(int y) const
Definition: Image.h:319
x_iterator row_end(int y) const
Return an x_iterator to the end of the y&#39;th row.
Definition: Image.h:324
RelativeSpanIterator(SpanList::const_iterator const &real, SpanList const &arr, int cx, int cy, bool forward=true)
Definition: Baseline.cc:716
A range of pixels within one row of an Image.
Definition: Span.h:54
bool operator<(const SpanList::const_iterator &other)
Definition: Baseline.cc:737
#define CONST_PTR(...)
Definition: base.h:47
bool operator==(RelativeSpanIterator &other)
Definition: Baseline.cc:747
lsst::afw::detection::Footprint Footprint
Definition: Source.h:61
Ref< ImagePixelT >::type image()
Return (a reference to) the image part of the Pixel pointed at by the iterator.
Definition: MaskedImage.h:147
PeakCatalog & getPeaks()
Definition: Footprint.h:123
bool operator==(const SpanList::const_iterator &other)
Definition: Baseline.cc:728
int getMinY() const
Definition: Box.h:125
boost::shared_ptr< lsst::afw::detection::Footprint > FootprintPtrT
Definition: Baseline.h:33
int getMaxX() const
Definition: Box.h:128
int getIy() const
Convenience accessors for the keys in the minimal schema.
Definition: Peak.h:210
boost::shared_ptr< Footprint > findEdgePixels() const
Find edge pixels on the footprint.
a place to record messages and descriptions of the state of processing.
Definition: Log.h:154
static Log & getDefaultLog()
int const x0
Definition: saturated.cc:45
An integer coordinate rectangle.
Definition: Box.h:53
A set of pixels in an Image, including those pixels&#39; actual values.
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
Ref< MaskPixelT >::type mask()
Return (a reference to) the mask part of the Pixel pointed at by the iterator.
Definition: MaskedImage.h:152
Iterator begin() const
Return an iterator to the first pixel in the Span.
Definition: Span.h:72
An iterator to the MaskedImage.
Definition: MaskedImage.h:109
x_iterator x_at(int x, int y) const
Return an x_iterator at the point (x, y)
Definition: MaskedImage.h:1006
def log
Definition: log.py:85
static void medianFilter(ImageT const &img, ImageT &outimg, int halfsize)
Definition: Baseline.cc:39
int getWidth() const
Return the number of columns in the image.
Definition: Image.h:237
std::vector< Span::Ptr > SpanList
The Footprint&#39;s Span list.
Definition: Footprint.h:71
_const_view_t::x_iterator const_x_iterator
A const iterator for traversing the pixels in a row.
Definition: Image.h:158
Extent2I const getDimensions() const
Definition: Box.h:153
_view_t::xy_locator xy_locator
An xy_locator.
Definition: Image.h:139
geom::Box2I getBBox() const
Return the Footprint&#39;s bounding box.
Definition: Footprint.h:206
Schema getSchema() const
Return the schema associated with the catalog&#39;s table.
Definition: Catalog.h:114
boost::shared_ptr< lsst::afw::image::Mask< MaskPixelT > > MaskPtrT
Definition: Baseline.h:30
geom::Box2I getBBox(ImageOrigin const origin=PARENT) const
Definition: MaskedImage.h:905
void debugf(const char *fmt,...)
Definition: Log.h:527
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:77
SpanList::const_iterator _real
Definition: Baseline.cc:806
det::Footprint::SpanList SpanList
Definition: Baseline.cc:712
static bool hasSignificantFluxAtEdge(ImagePtrT, boost::shared_ptr< lsst::afw::detection::Footprint >, ImagePixelT threshold)
Definition: Baseline.cc:1281
bool contains(Point2I const &point) const
Return true if the box contains the point.
def deblend
Deblend a single footprint in a maskedImage.
Definition: baseline.py:233
bool operator<=(const SpanList::const_iterator &other)
Definition: Baseline.cc:734
int getMaxY() const
Definition: Box.h:129
int getMinX() const
Definition: Box.h:124
A set of pixels in an Image.
Definition: Footprint.h:62
void include(Point2I const &point)
Expand this to ensure that this-&gt;contains(point).
int getY0() const
Definition: Image.h:255
int x
static void makeMonotonic(ImageT &img, lsst::afw::detection::PeakRecord const &pk)
Definition: Baseline.cc:117
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
Iterator end() const
Return an iterator to one past the last pixel in the Span.
Definition: Span.h:75
static void _find_stray_flux(lsst::afw::detection::Footprint const &foot, ImagePtrT tsum, MaskedImageT const &img, int strayFluxOptions, std::vector< boost::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< boost::shared_ptr< typename lsst::afw::detection::HeavyFootprint< ImagePixelT, MaskPixelT, VariancePixelT > > > &strays)
Definition: Baseline.cc:289
MaskPixelT mask() const
Return the mask part of a Pixel.
Definition: Pixel.h:201
static boost::shared_ptr< lsst::afw::detection::Footprint > getSignificantEdgePixels(ImagePtrT, boost::shared_ptr< lsst::afw::detection::Footprint >, ImagePixelT threshold)
Definition: Baseline.cc:1317
SpanList::const_iterator _end
Definition: Baseline.cc:807
MaskPtr getMask(bool const noThrow=false) const
Return a (Ptr to) the MaskedImage&#39;s mask.
Definition: MaskedImage.h:879
afw::table::Key< double > sigma1
static void _sum_templates(std::vector< ImagePtrT > timgs, ImagePtrT tsum)
Definition: Baseline.cc:500
_view_t::x_iterator x_iterator
An iterator for traversing the pixels in a row.
Definition: Image.h:151
bool operator>(const SpanList::const_iterator &other)
Definition: Baseline.cc:743
int getY() const
Return the y-value.
Definition: Span.h:81
bool operator!=(RelativeSpanIterator &other)
Definition: Baseline.cc:752
xy_locator xy_at(int x, int y) const
Definition: Image.h:351
x_iterator row_begin(int y) const
Return an x_iterator to the start of the image.
Definition: MaskedImage.cc:742
ExpressionTraits< Derived >::Iterator Iterator
Nested expression or element iterator.
int getHeight() const
Return the number of rows in the image.
Definition: Image.h:239
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)
Definition: Baseline.cc:1084
bool operator!=(const SpanList::const_iterator &other)
Definition: Baseline.cc:731
RelativeSpanIterator operator++()
Definition: Baseline.cc:756
int const y0
Definition: saturated.cc:45
RelativeSpanIterator operator++(int dummy)
Definition: Baseline.cc:764
Record class that represents a peak in a Footprint.
Definition: Peak.h:40
boost::shared_ptr< lsst::afw::image::MaskedImage< ImagePixelT > MaskPixelT, VariancePixelT > MaskedImagePtrT
Definition: Baseline.h:26
int getX1() const
Return the ending x-value.
Definition: Span.h:79
T real(const T &x)
Definition: Integrate.h:193
int getWidth() const
Definition: Box.h:154
boost::shared_ptr< lsst::afw::detection::HeavyFootprint< ImagePixelT > MaskPixelT, VariancePixelT > HeavyFootprintPtrT
Definition: Baseline.h:36
int getIx() const
Convenience accessors for the keys in the minimal schema.
Definition: Peak.h:209
void copyWithinFootprint(Footprint const &foot, boost::shared_ptr< ImageOrMaskedImageT > const input, boost::shared_ptr< ImageOrMaskedImageT > output)
Include files required for standard LSST Exception handling.
boost::shared_ptr< lsst::afw::image::Image< ImagePixelT > > ImagePtrT
Definition: Baseline.h:28
int getX0() const
Definition: Image.h:247
A const locator for the MaskedImage.
Definition: MaskedImage.h:115
int getX0() const
Return the starting x-value.
Definition: Span.h:77
Iterator begin() const
Return an Iterator to the beginning of the array.
Definition: ArrayBase.h:99