LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
Public Types | Static Public Member Functions | Static Public Attributes | List of all members
lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT > Class Template Reference

#include <BaselineUtils.h>

Public Types

typedef lsst::afw::image::MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > MaskedImageT
 
typedef std::shared_ptr< lsst::afw::image::MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > > MaskedImagePtrT
 
typedef lsst::afw::image::Image< ImagePixelT > ImageT
 
typedef std::shared_ptr< lsst::afw::image::Image< ImagePixelT > > ImagePtrT
 
typedef lsst::afw::image::Mask< MaskPixelT > MaskT
 
typedef std::shared_ptr< lsst::afw::image::Mask< MaskPixelT > > MaskPtrT
 
typedef lsst::afw::detection::Footprint FootprintT
 
typedef std::shared_ptr< lsst::afw::detection::FootprintFootprintPtrT
 
typedef lsst::afw::detection::HeavyFootprint< ImagePixelT, MaskPixelT, VariancePixelT > HeavyFootprintT
 
typedef std::shared_ptr< lsst::afw::detection::HeavyFootprint< ImagePixelT, MaskPixelT, VariancePixelT > > HeavyFootprintPtrT
 

Static Public Member Functions

static std::shared_ptr< lsst::afw::detection::FootprintsymmetrizeFootprint (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 twofold rotational symmetry) – the AND of the two symmetric halves. More...
 
static std::pair< ImagePtrT, FootprintPtrTbuildSymmetricTemplate (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 MaskedImage and Footprint describing a footprint where: output pixels (cx + dx, cy + dy) and (cx - dx, cy - dy) = min(input pixels (cx + dx, cy + dy) and (cx - dx, cy - dy)) More...
 
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. More...
 
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 values smaller than those close to the peak; make the profile monotonic-decreasing. More...
 
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,tfoots. More...
 
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. More...
 
static std::shared_ptr< lsst::afw::detection::FootprintgetSignificantEdgePixels (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, above threshold thresh. More...
 
static void _sum_templates (std::vector< ImagePtrT > timgs, ImagePtrT tsum)
 
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 Public Attributes

static const int ASSIGN_STRAYFLUX = 0x1
 
static const int STRAYFLUX_TO_POINT_SOURCES_WHEN_NECESSARY = 0x2
 
static const int STRAYFLUX_TO_POINT_SOURCES_ALWAYS = 0x4
 
static const int STRAYFLUX_R_TO_FOOTPRINT = 0x8
 
static const int STRAYFLUX_NEAREST_FOOTPRINT = 0x10
 
static const int STRAYFLUX_TRIM = 0x20
 

Detailed Description

template<typename ImagePixelT, typename MaskPixelT = lsst::afw::image::MaskPixel, typename VariancePixelT = lsst::afw::image::VariancePixel>
class lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >

Definition at line 22 of file BaselineUtils.h.

Member Typedef Documentation

◆ FootprintPtrT

template<typename ImagePixelT , typename MaskPixelT = lsst::afw::image::MaskPixel, typename VariancePixelT = lsst::afw::image::VariancePixel>
typedef std::shared_ptr<lsst::afw::detection::Footprint> lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::FootprintPtrT

Definition at line 33 of file BaselineUtils.h.

◆ FootprintT

template<typename ImagePixelT , typename MaskPixelT = lsst::afw::image::MaskPixel, typename VariancePixelT = lsst::afw::image::VariancePixel>
typedef lsst::afw::detection::Footprint lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::FootprintT

Definition at line 32 of file BaselineUtils.h.

◆ HeavyFootprintPtrT

template<typename ImagePixelT , typename MaskPixelT = lsst::afw::image::MaskPixel, typename VariancePixelT = lsst::afw::image::VariancePixel>
typedef std::shared_ptr<lsst::afw::detection::HeavyFootprint<ImagePixelT, MaskPixelT, VariancePixelT> > lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::HeavyFootprintPtrT

Definition at line 36 of file BaselineUtils.h.

◆ HeavyFootprintT

template<typename ImagePixelT , typename MaskPixelT = lsst::afw::image::MaskPixel, typename VariancePixelT = lsst::afw::image::VariancePixel>
typedef lsst::afw::detection::HeavyFootprint<ImagePixelT, MaskPixelT, VariancePixelT> lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::HeavyFootprintT

Definition at line 34 of file BaselineUtils.h.

◆ ImagePtrT

template<typename ImagePixelT , typename MaskPixelT = lsst::afw::image::MaskPixel, typename VariancePixelT = lsst::afw::image::VariancePixel>
typedef std::shared_ptr<lsst::afw::image::Image<ImagePixelT> > lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::ImagePtrT

Definition at line 28 of file BaselineUtils.h.

◆ ImageT

template<typename ImagePixelT , typename MaskPixelT = lsst::afw::image::MaskPixel, typename VariancePixelT = lsst::afw::image::VariancePixel>
typedef lsst::afw::image::Image<ImagePixelT> lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::ImageT

Definition at line 27 of file BaselineUtils.h.

◆ MaskedImagePtrT

template<typename ImagePixelT , typename MaskPixelT = lsst::afw::image::MaskPixel, typename VariancePixelT = lsst::afw::image::VariancePixel>
typedef std::shared_ptr<lsst::afw::image::MaskedImage<ImagePixelT, MaskPixelT, VariancePixelT> > lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::MaskedImagePtrT

Definition at line 26 of file BaselineUtils.h.

◆ MaskedImageT

template<typename ImagePixelT , typename MaskPixelT = lsst::afw::image::MaskPixel, typename VariancePixelT = lsst::afw::image::VariancePixel>
typedef lsst::afw::image::MaskedImage<ImagePixelT, MaskPixelT, VariancePixelT> lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::MaskedImageT

Definition at line 25 of file BaselineUtils.h.

◆ MaskPtrT

template<typename ImagePixelT , typename MaskPixelT = lsst::afw::image::MaskPixel, typename VariancePixelT = lsst::afw::image::VariancePixel>
typedef std::shared_ptr<lsst::afw::image::Mask<MaskPixelT> > lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::MaskPtrT

Definition at line 30 of file BaselineUtils.h.

◆ MaskT

template<typename ImagePixelT , typename MaskPixelT = lsst::afw::image::MaskPixel, typename VariancePixelT = lsst::afw::image::VariancePixel>
typedef lsst::afw::image::Mask<MaskPixelT> lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::MaskT

Definition at line 29 of file BaselineUtils.h.

Member Function Documentation

◆ _find_stray_flux()

template<typename ImagePixelT , typename MaskPixelT = lsst::afw::image::MaskPixel, typename VariancePixelT = lsst::afw::image::VariancePixel>
void lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::_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

Hmm, this is a little bit dangerous: we're assuming that the HeavyFootprint stores its pixels in the same order that we iterate over them above (ie, lexicographic).

Definition at line 393 of file BaselineUtils.cc.

404  {
405 
406  typedef typename det::HeavyFootprint<ImagePixelT, MaskPixelT, VariancePixelT> HeavyFootprint;
408 
409  // when doing stray flux: the footprints and pixels, which we'll
410  // combine into the return 'strays' HeavyFootprint at the end.
412  std::vector<std::vector<afwGeom::Span>> straySpans(tfoots.size());
416 
417  int ix0 = img.getX0();
418  int iy0 = img.getY0();
419  geom::Box2I sumbb = tsum->getBBox();
420  int sumx0 = sumbb.getMinX();
421  int sumy0 = sumbb.getMinY();
422 
423  for (size_t i=0; i<tfoots.size(); ++i) {
426  straymask.push_back(std::vector<MaskPixelT>());
428  }
429 
430  bool always = (strayFluxOptions & STRAYFLUX_TO_POINT_SOURCES_ALWAYS);
431 
432  typedef std::uint16_t itype;
434 
435  if (strayFluxOptions & STRAYFLUX_NEAREST_FOOTPRINT) {
436  // Compute the map of which footprint is closest to each
437  // pixel in the bbox.
438  typedef std::uint16_t dtype;
441 
443  std::vector<std::shared_ptr<det::Footprint>>* footlist = &tfoots;
444 
445  if (!always && ispsf.size()) {
446  // create a temp list that has empty footprints in place
447  // of all the point sources.
448  auto empty = std::make_shared<det::Footprint>();
449  empty->setPeakSchema(foot.getPeaks().getSchema());
450  for (size_t i=0; i<tfoots.size(); ++i) {
451  if (ispsf[i]) {
452  templist.push_back(empty);
453  } else {
454  templist.push_back(tfoots[i]);
455  }
456  }
457  footlist = &templist;
458  }
459  nearestFootprint(*footlist, nearest, dist);
460  }
461 
462  // Go through the (parent) Footprint looking for stray flux:
463  // pixels that are not claimed by any template, and positive.
464  for (afwGeom::Span const & s : *foot.getSpans()) {
465  int y = s.getY();
466  int x0 = s.getX0();
467  int x1 = s.getX1();
468  typename ImageT::x_iterator tsum_it =
469  tsum->row_begin(y - sumy0) + (x0 - sumx0);
470  typename MaskedImageT::x_iterator in_it =
471  img.row_begin(y - iy0) + (x0 - ix0);
472  double contrib[tfoots.size()];
473 
474  for (int x = x0; x <= x1; ++x, ++tsum_it, ++in_it) {
475  // Skip pixels that are covered by at least one
476  // template (*tsum_it > 0) or the input is not
477  // positive (*in_it <= 0).
478  if ((*tsum_it > 0) || (*in_it).image() <= 0) {
479  continue;
480  }
481 
482  if (strayFluxOptions & STRAYFLUX_R_TO_FOOTPRINT) {
483  // we'll compute these just-in-time
484  for (size_t i=0; i<tfoots.size(); ++i) {
485  contrib[i] = -1.0;
486  }
487  } else if (strayFluxOptions & STRAYFLUX_NEAREST_FOOTPRINT) {
488  for (size_t i=0; i<tfoots.size(); ++i) {
489  contrib[i] = 0.0;
490  }
491  int i = (*nearest)[geom::Point2I(x, y)];
492  contrib[i] = 1.0;
493  } else {
494  // R_TO_PEAK
495  for (size_t i=0; i<tfoots.size(); ++i) {
496  // Split the stray flux by 1/(1+r^2) to peaks
497  int dx, dy;
498  dx = pkx[i] - x;
499  dy = pky[i] - y;
500  contrib[i] = 1. / (1. + dx*dx + dy*dy);
501  }
502  }
503 
504  // Round 1: skip point sources unless STRAYFLUX_TO_POINT_SOURCES_ALWAYS
505  // are we going to assign stray flux to ptsrcs?
506  bool ptsrcs = always;
507  double csum = 0.;
508  for (size_t i=0; i<tfoots.size(); ++i) {
509  // if we're skipping point sources and this is a point source...
510  if ((!ptsrcs) && ispsf.size() && ispsf[i]) {
511  continue;
512  }
513  if (contrib[i] == -1.0) {
514  contrib[i] = _get_contrib_r_to_footprint(x, y, tfoots[i]);
515  }
516  csum += contrib[i];
517  }
518  if ((csum == 0.) &&
519  (strayFluxOptions &
521  // No extended sources -- assign to pt sources
522  //LOGL_DEBUG(_log, "necessary to assign stray flux to point sources");
523  ptsrcs = true;
524  for (size_t i=0; i<tfoots.size(); ++i) {
525  if (contrib[i] == -1.0) {
526  contrib[i] = _get_contrib_r_to_footprint(x, y, tfoots[i]);
527  }
528  csum += contrib[i];
529  }
530  }
531 
532  // Drop small contributions...
533  double strayclip = (clipStrayFluxFraction * csum);
534  csum = 0.;
535  for (size_t i=0; i<tfoots.size(); ++i) {
536  // skip ptsrcs?
537  if ((!ptsrcs) && ispsf.size() && ispsf[i]) {
538  contrib[i] = 0.;
539  continue;
540  }
541  // skip small contributions
542  if (contrib[i] < strayclip) {
543  contrib[i] = 0.;
544  continue;
545  }
546  csum += contrib[i];
547  }
548 
549  for (size_t i=0; i<tfoots.size(); ++i) {
550  if (contrib[i] == 0.) {
551  continue;
552  }
553  // the stray flux to give to template i
554  double p = (contrib[i] / csum) * (*in_it).image();
555 
556  if (!strayfoot[i]) {
557  strayfoot[i] = std::make_shared<det::Footprint>();
558  strayfoot[i]->setPeakSchema(foot.getPeaks().getSchema());
559  }
560  straySpans[i].push_back(afwGeom::Span(y, x, x));
561  straypix[i].push_back(p);
562  straymask[i].push_back((*in_it).mask());
563  strayvar[i].push_back((*in_it).variance());
564  }
565  }
566  }
567 
568  // Store the stray flux in HeavyFootprints
569  for (size_t i=0; i<tfoots.size(); ++i) {
570  if (strayfoot[i]) {
571  strayfoot[i]->setSpans(std::make_shared<afwGeom::SpanSet>(straySpans[i]));
572  }
573  if (!strayfoot[i]) {
574  strays.push_back(HeavyFootprintPtrT());
575  } else {
579  HeavyFootprintPtrT heavy(new HeavyFootprint(*strayfoot[i]));
580  ndarray::Array<ImagePixelT,1,1> himg = heavy->getImageArray();
587 
588  assert((size_t)strayfoot[i]->getArea() == straypix[i].size());
589 
590  for (spix = straypix[i].begin(),
591  smask = straymask[i].begin(),
592  svar = strayvar [i].begin(),
593  hpix = himg.begin(),
594  mpix = heavy->getMaskArray().begin(),
595  vpix = heavy->getVarianceArray().begin();
596  spix != straypix[i].end();
597  ++spix, ++smask, ++svar, ++hpix, ++mpix, ++vpix) {
598  *hpix = *spix;
599  *mpix = *smask;
600  *vpix = *svar;
601  }
602  strays.push_back(heavy);
603  }
604  }
605 }
double x
int y
Definition: SpanSet.cc:48
T begin(T... args)
A set of pixels in an Image, including those pixels' actual values.
A range of pixels within one row of an Image.
Definition: Span.h:47
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
MaskedImageIterator< typename Image::x_iterator, typename Mask::x_iterator, typename Variance::x_iterator > x_iterator
An iterator to a row of a MaskedImage.
Definition: MaskedImage.h:540
An integer coordinate rectangle.
Definition: Box.h:55
int getMinY() const noexcept
Definition: Box.h:158
int getMinX() const noexcept
Definition: Box.h:157
static const int STRAYFLUX_TO_POINT_SOURCES_WHEN_NECESSARY
Definition: BaselineUtils.h:63
static const int STRAYFLUX_NEAREST_FOOTPRINT
Definition: BaselineUtils.h:68
std::shared_ptr< lsst::afw::detection::HeavyFootprint< ImagePixelT, MaskPixelT, VariancePixelT > > HeavyFootprintPtrT
Definition: BaselineUtils.h:36
static const int STRAYFLUX_TO_POINT_SOURCES_ALWAYS
Definition: BaselineUtils.h:64
Point< int, 2 > Point2I
Definition: Point.h:321
FastFinder::Iterator Iterator
Definition: FastFinder.cc:178
T push_back(T... args)
T size(T... args)

◆ _sum_templates()

template<typename ImagePixelT , typename MaskPixelT , typename VariancePixelT >
void lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::_sum_templates ( std::vector< ImagePtrT timgs,
ImagePtrT  tsum 
)
static

Definition at line 609 of file BaselineUtils.cc.

611  {
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 }
void clip(Box2I const &other) noexcept
Shrink this to ensure that other.contains(*this).
Definition: Box.cc:189
int getWidth() const noexcept
Definition: Box.h:187
int getMaxY() const noexcept
Definition: Box.h:162
std::shared_ptr< lsst::afw::image::Image< ImagePixelT > > ImagePtrT
Definition: BaselineUtils.h:28
T max(T... args)

◆ apportionFlux()

template<typename ImagePixelT , typename MaskPixelT = lsst::afw::image::MaskPixel, typename VariancePixelT = lsst::afw::image::VariancePixel>
std::vector< typename std::shared_ptr< image::MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > > > lsst::meas::deblender::BaselineUtils< 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 
)
static

Splits flux in a given image img, within a given footprint foot, among a number of templates timgs,tfoots.

This is where actual "deblending" takes place.

timgs* and tfoots MUST be the same length.

Flux is assigned to templates according to their relative heights at each pixel.

If strayFluxOptions includes ASSIGN_STRAYFLUX, then "stray flux" – flux in the parent footprint that is not covered by any of the template footprints – is assigned to templates based on their 1/(1+r^2) distance.

If strayFluxOptions includes STRAYFLUX_R_TO_FOOTPRINT, the stray flux is distributed to the footprints based on 1/(1+r^2) of the minimum distance from the stray flux to footprint.

If strayFluxOptions includes "STRAYFLUX_NEAREST_FOOTPRINT*, the stray flux is assigned to the footprint with lowest L-1 (Manhattan) distance to the stray flux.

Otherwise, stray flux is assigned based on (1/(1+r^2) from the peaks.

If strayFluxOptions includes STRAYFLUX_TO_POINT_SOURCES_ALWAYS, then point sources are always included in the 1/(1+r^2) splitting. Otherwise, if STRAYFLUX_TO_POINT_SOURCES_WHEN_NECESSARY, point sources are included only if there are no extended sources nearby.

If any stray-flux portion is less than clipStrayFluxFraction, it is clipped to zero.

When doing stray flux, the "strays" arg is used as an extra return value, the stray flux assigned to each template.

When doing stray flux, the ispsf, pkx, and pky arrays are required. They give the peak x,y coords plus whether the peak is believed (by the deblender) to be a point source. pkx and pky MUST be the same length as timgs. If ispsf has nonzero length, it MUST be the same length as timgs.

If tsum is given, is it set to the sum of max(0, template).

The return value is a vector of MaskedImages containing the flux assigned to each template.

Definition at line 692 of file BaselineUtils.cc.

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 }
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
#define LOG_GET(logger)
Returns a Log object associated with logger.
Definition: Log.h:75
#define LOG_LOGGER
Definition: Log.h:714
Extent2I const getDimensions() const noexcept
Definition: Box.h:186
std::shared_ptr< lsst::afw::image::MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > > MaskedImagePtrT
Definition: BaselineUtils.h:26
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)
lsst::afw::image::MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > MaskedImageT
Definition: BaselineUtils.h:25
lsst::afw::image::Image< ImagePixelT > ImageT
Definition: BaselineUtils.h:27
static void _sum_templates(std::vector< ImagePtrT > timgs, ImagePtrT tsum)
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
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174

◆ buildSymmetricTemplate()

template<typename ImagePixelT , typename MaskPixelT , typename VariancePixelT >
std::pair< typename std::shared_ptr< lsst::afw::image::Image< ImagePixelT > >, typename std::shared_ptr< lsst::afw::detection::Footprint > > lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::buildSymmetricTemplate ( MaskedImageT const &  img,
lsst::afw::detection::Footprint const &  foot,
lsst::afw::detection::PeakRecord const &  pk,
double  sigma1,
bool  minZero,
bool  patchEdge,
bool *  patchedEdges 
)
static

Given an img, footprint foot, and peak, creates a symmetric template around the peak; produce a MaskedImage and Footprint describing a footprint where: output pixels (cx + dx, cy + dy) and (cx - dx, cy - dy) = min(input pixels (cx + dx, cy + dy) and (cx - dx, cy - dy))

If patchEdge is true and if the footprint touches pixels with the EDGE bit set, then for spans whose symmetric mirror are outside the image, the symmetric footprint is grown to include them and their pixel values are stored.

Definition at line 1189 of file BaselineUtils.cc.

1197  {
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 }
afw::table::Key< afw::table::Array< MaskPixelT > > mask
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
Definition: Log.h:515
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
const_MaskedImageLocator< typename Image::xy_locator, typename Mask::xy_locator, typename Variance::xy_locator > const_xy_locator
A const_locator for a MaskedImage.
Definition: MaskedImage.h:555
void include(Point2I const &point)
Expand this to ensure that this->contains(point).
Definition: Box.cc:152
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
std::shared_ptr< lsst::afw::detection::Footprint > FootprintPtrT
Definition: BaselineUtils.h:33
std::shared_ptr< lsst::afw::image::Mask< MaskPixelT > > MaskPtrT
Definition: BaselineUtils.h:30
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...
T min(T... args)
T move(T... args)

◆ getSignificantEdgePixels()

template<typename ImagePixelT , typename MaskPixelT , typename VariancePixelT >
std::shared_ptr< det::Footprint > lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::getSignificantEdgePixels ( ImagePtrT  img,
std::shared_ptr< lsst::afw::detection::Footprint sfoot,
ImagePixelT  threshold 
)
static

Returns a list of pixels that are on the edge of the given Footprint sfoot* in image img, above threshold thresh.

Definition at line 1418 of file BaselineUtils.cc.

1421  {
1422  LOG_LOGGER _log = LOG_GET("meas.deblender.getSignificantEdgePixels");
1423 
1424 
1425  auto significant = std::make_shared<det::Footprint>();
1426  significant->setPeakSchema(sfoot->getPeaks().getSchema());
1427 
1428  int const x0 = img->getX0(), y0 = img->getY0();
1429  std::shared_ptr<afwGeom::SpanSet> edgeSpans = sfoot->getSpans()->findEdgePixels();
1430  std::vector<afwGeom::Span> tmpSpans;
1431  for (afwGeom::SpanSet::const_iterator ss = edgeSpans->begin(); ss != edgeSpans->end(); ++ss) {
1432  afwGeom::Span const& span = *ss;
1433  int const y = span.getY();
1434  int x = span.getX0();
1435  typename ImageT::const_x_iterator iter = img->x_at(x - x0, y - y0);
1436  bool onSpan = false; // Are we in a span of interest
1437  int xSpan; // Starting x of span
1438  for (; x <= span.getX1(); ++x, ++iter) {
1439  if (*iter >= thresh) {
1440  onSpan = true;
1441  xSpan = x;
1442  } else if (onSpan) {
1443  onSpan = false;
1444  tmpSpans.push_back(afwGeom::Span(y, xSpan, x - 1));
1445  }
1446  }
1447  if (onSpan) {
1448  tmpSpans.push_back(afwGeom::Span(y, xSpan, span.getX1()));
1449  }
1450  }
1451  significant->setSpans(std::make_shared<afwGeom::SpanSet>(std::move(tmpSpans)));
1452  return significant;
1453 }
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
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
typename _const_view_t::x_iterator const_x_iterator
A const iterator for traversing the pixels in a row.
Definition: ImageBase.h:141
Schema getSchema() const
Return the schema associated with the catalog's table.
Definition: Catalog.h:118

◆ hasSignificantFluxAtEdge()

template<typename ImagePixelT , typename MaskPixelT , typename VariancePixelT >
bool lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::hasSignificantFluxAtEdge ( ImagePtrT  img,
std::shared_ptr< lsst::afw::detection::Footprint sfoot,
ImagePixelT  threshold 
)
static

Returns true if the given Footprint sfoot in image img has flux above value thresh at its edge.

Definition at line 1385 of file BaselineUtils.cc.

1388  {
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 }

◆ makeMonotonic()

template<typename ImagePixelT , typename MaskPixelT , typename VariancePixelT >
void lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::makeMonotonic ( ImageT img,
lsst::afw::detection::PeakRecord const &  pk 
)
static

Given an image mimg and Peak location peak, overwrite mimg so that pixels further from the peak have values smaller than those close to the peak; make the profile monotonic-decreasing.

The exact algorithm is a little more complicated than that. The basic idea is of "casting a shadow" from a pixel to pixels farther from the peak in the same direction. Done naively, this results in very narrow "shadows" and ragged profiles. A tweak is to make the shadows "fatter" – make a pixel shadow a wedge of pixels – but if one does this naively, the wedge gets wider and wider too quickly. The algorithm works out from the peak in square "rings" of pixels, so if a pixel shadows a wedge 30 degrees wide, in the next ring of pixels the shadowed pixel at largest angle from the shadowing pixel will shade a yet-larger wedge, expanding the shadowing angle. To reduce this effect, we work in chunks of 5 pixels in radius, only copying the intermediate pixels to the "shadowing" image at the end of each chunk.

Currently the mask and variance planes of the input are totally ignored.

For illustration, run tests/monotonic.py and look at im*.png

Definition at line 225 of file BaselineUtils.cc.

228  {
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 }
T lround(T... args)

◆ medianFilter()

template<typename ImagePixelT , typename MaskPixelT , typename VariancePixelT >
void lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::medianFilter ( ImageT const &  img,
ImageT out,
int  halfsize 
)
static

Run a spatial median filter over the given input img, writing the results to out.

halfsize is half the box size of the filter; ie, a halfsize of 50 means that each output ixel will be the median of the pixels in a 101 x 101-pixel box in the input image.

Margins are handled horribly: the median is computed only for pixels more than halfsize away from the edges; pixels near the edges are simply copied from the input img to out.

Mask and variance planes are, likewise, simply copied from img to out*.

Definition at line 147 of file BaselineUtils.cc.

150  {
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 }
int end
typename _view_t::xy_locator xy_locator
An xy_locator.
Definition: ImageBase.h:121
T nth_element(T... args)

◆ symmetrizeFootprint()

template<typename ImagePixelT , typename MaskPixelT , typename VariancePixelT >
std::shared_ptr< lsst::afw::detection::Footprint > lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::symmetrizeFootprint ( lsst::afw::detection::Footprint const &  foot,
int  cx,
int  cy 
)
static

Given a Footprint foot and peak cx,cy, returns a Footprint that is symmetric around the peak (with twofold rotational symmetry) – the AND of the two symmetric halves.

Definition at line 974 of file BaselineUtils.cc.

977  {
978 
979  auto sfoot = std::make_shared<det::Footprint>();
980  sfoot->setPeakSchema(foot.getPeaks().getSchema());
981  afwGeom::SpanSet const & spans = *foot.getSpans();
982 
983  LOG_LOGGER _log = LOG_GET("lsst.meas.deblender.symmetrizeFootprint");
984 
985  // Find the Span containing the peak.
986  afwGeom::Span target(cy, cx, cx);
988  std::upper_bound(spans.begin(), spans.end(), target, span_compare);
989  // upper_bound returns the last position where "target" could be inserted;
990  // ie, the first Span larger than "target". The Span containing "target"
991  // should be peakspan-1 or (if the peak is on the first pixel in the span,
992  // peakspan.
993  afwGeom::Span sp;
994  if (peakspan == spans.begin()) {
995  sp = *peakspan;
996  if (!sp.contains(cx, cy)) {
997  LOGL_WARN(_log,
998  "Failed to find span containing (%i,%i): before the beginning of this footprint", cx, cy);
1000  }
1001  } else {
1002  --peakspan;
1003  sp = *peakspan;
1004 
1005  if (!sp.contains(cx, cy)) {
1006  ++peakspan;
1007  sp = *peakspan;
1008  if (!sp.contains(cx, cy)) {
1009  geom::Box2I fbb = foot.getBBox();
1010  LOGL_WARN(_log, "Failed to find span containing (%i,%i): nearest is %i, [%i,%i]. "
1011  "Footprint bbox is [%i,%i],[%i,%i]",
1012  cx, cy, sp.getY(), sp.getX0(), sp.getX1(),
1013  fbb.getMinX(), fbb.getMaxX(), fbb.getMinY(), fbb.getMaxY());
1015  }
1016  }
1017  }
1018  LOGL_DEBUG(_log, "Span containing (%i,%i): (x=[%i,%i], y=%i)",
1019  cx, cy, sp.getX0(), sp.getX1(), sp.getY());
1020 
1021  // The symmetric templates are essentially an AND of the footprint
1022  // pixels and its 180-degree-rotated self, rotated around the
1023  // peak (cx,cy).
1024  //
1025  // We iterate forward and backward simultaneously, starting from
1026  // the span containing the peak and moving out, row by row.
1027  //
1028  // In the loop below, we search for the next pair of Spans that
1029  // overlap (in "dx" from the center), output the overlapping
1030  // portion of the Spans, and advance either the "fwd" or "back"
1031  // iterator. When we fail to find an overlapping pair of Spans,
1032  // we move on to the next row.
1033  //
1034  // [The following paragraph is somewhat obsoleted by the
1035  // RelativeSpanIterator class, which performs some of the renaming
1036  // and the dx,dy coords.]
1037  //
1038  // '''In reading the code, "forward", "advancing", etc, are all
1039  // from the perspective of the "fwd" iterator (the one going
1040  // forward through the Span list, from low to high Y and then low
1041  // to high X). It will help to imagine making a copy of the
1042  // footprint and rotating it around the center pixel by 180
1043  // degrees, so that "fwd" and "back" are both iterating the same
1044  // direction; we're then just finding the AND of those two
1045  // iterators, except we have to work in dx,dy coordinates rather
1046  // than original x,y coords, and the accessors for "back" are
1047  // opposite.'''
1048 
1049  RelativeSpanIterator fwd (peakspan, spans, cx, cy, true);
1050  RelativeSpanIterator back(peakspan, spans, cx, cy, false);
1051 
1052  int dy = 0;
1053  std::vector<afwGeom::Span> tmpSpans;
1054  while (fwd.notDone() && back.notDone()) {
1055  // forward and backward "y"; just symmetric around cy
1056  int fy = cy + dy;
1057  int by = cy - dy;
1058  // delta-x of the beginnings of the spans, for "fwd" and "back"
1059  int fdxlo = fwd.dxlo();
1060  int bdxlo = back.dxlo();
1061 
1062  // First find:
1063  // fend -- first span in the next row, or end(); ie,
1064  // the end of this row in the forward direction
1065  // bend -- the end of this row in the backward direction
1066  RelativeSpanIterator fend, bend;
1067  for (fend = fwd; fend.notDone(); ++fend) {
1068  if (fend.dy() != dy)
1069  break;
1070  }
1071  for (bend = back; bend.notDone(); ++bend) {
1072  if (bend.dy() != dy)
1073  break;
1074  }
1075 
1076  LOGL_DEBUG(_log, "dy=%i, fy=%i, fx=[%i, %i], by=%i, fx=[%i, %i], fdx=%i, bdx=%i",
1077  dy, fy, fwd.x0(), fwd.x1(), by, back.x0(), back.x1(),
1078  fdxlo, bdxlo);
1079 
1080  // Find possibly-overlapping span
1081  if (bdxlo > fdxlo) {
1082  LOGL_DEBUG(_log, "Advancing forward.");
1083  // While the "forward" span is entirely to the "left" of the "backward" span,
1084  // (in dx coords), ie, |---fwd---X X---back---|
1085  // and we are comparing the edges marked X
1086  while ((fwd != fend) && (fwd.dxhi() < bdxlo)) {
1087  fwd++;
1088  if (fwd == fend) {
1089  LOGL_DEBUG(_log, "Reached fend");
1090  } else {
1091  LOGL_DEBUG(_log, "Advanced to forward span %i, [%i, %i]",
1092  fy, fwd.x0(), fwd.x1());
1093  }
1094  }
1095  } else if (fdxlo > bdxlo) {
1096  LOGL_DEBUG(_log, "Advancing backward.");
1097  // While the "backward" span is entirely to the "left" of the "foreward" span,
1098  // (in dx coords), ie, |---back---X X---fwd---|
1099  // and we are comparing the edges marked X
1100  while ((back != bend) && (back.dxhi() < fdxlo)) {
1101  back++;
1102  if (back == bend) {
1103  LOGL_DEBUG(_log, "Reached bend");
1104  } else {
1105  LOGL_DEBUG(_log, "Advanced to backward span %i, [%i, %i]",
1106  by, back.x0(), back.x1());
1107  }
1108  }
1109  }
1110 
1111  if ((back == bend) || (fwd == fend)) {
1112  // We reached the end of the row without finding spans that could
1113  // overlap. Move onto the next dy.
1114  if (back == bend) {
1115  LOGL_DEBUG(_log, "Reached bend");
1116  }
1117  if (fwd == fend) {
1118  LOGL_DEBUG(_log, "Reached fend");
1119  }
1120  back = bend;
1121  fwd = fend;
1122  dy++;
1123  continue;
1124  }
1125 
1126  // Spans may overlap -- find the overlapping part.
1127  int dxlo = std::max(fwd.dxlo(), back.dxlo());
1128  int dxhi = std::min(fwd.dxhi(), back.dxhi());
1129  if (dxlo <= dxhi) {
1130  LOGL_DEBUG(_log, "Adding span fwd %i, [%i, %i], back %i, [%i, %i]",
1131  fy, cx+dxlo, cx+dxhi, by, cx-dxhi, cx-dxlo);
1132  tmpSpans.push_back(afwGeom::Span(fy, cx + dxlo, cx + dxhi));
1133  tmpSpans.push_back(afwGeom::Span(by, cx - dxhi, cx - dxlo));
1134  }
1135 
1136  // Advance the one whose "hi" edge is smallest
1137  if (fwd.dxhi() < back.dxhi()) {
1138  fwd++;
1139  if (fwd == fend) {
1140  LOGL_DEBUG(_log, "Stepped to fend");
1141  } else {
1142  LOGL_DEBUG(_log, "Stepped forward to span %i, [%i, %i]",
1143  fwd.y(), fwd.x0(), fwd.x1());
1144  }
1145  } else {
1146  back++;
1147  if (back == bend) {
1148  LOGL_DEBUG(_log, "Stepped to bend");
1149  } else {
1150  LOGL_DEBUG(_log, "Stepped backward to span %i, [%i, %i]",
1151  back.y(), back.x0(), back.x1());
1152  }
1153  }
1154 
1155  if ((back == bend) || (fwd == fend)) {
1156  // Reached the end of the row. On to the next dy!
1157  if (back == bend) {
1158  LOGL_DEBUG(_log, "Reached bend");
1159  }
1160  if (fwd == fend) {
1161  LOGL_DEBUG(_log, "Reached fend");
1162  }
1163  back = bend;
1164  fwd = fend;
1165  dy++;
1166  continue;
1167  }
1168 
1169  }
1170  sfoot->setSpans(std::make_shared<afwGeom::SpanSet>(std::move(tmpSpans)));
1171  return sfoot;
1172 }
Key< Flag > const & target
#define LOGL_WARN(logger, message...)
Log a warn-level message using a varargs/printf style interface.
Definition: Log.h:547
This is a convenience class used in symmetrizeFootprint, wrapping the idea of iterating through a Spa...
bool contains(int x) const noexcept
Definition: Span.h:95
T upper_bound(T... args)

Member Data Documentation

◆ ASSIGN_STRAYFLUX

template<typename ImagePixelT , typename MaskPixelT , typename VariancePixelT >
const int lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::ASSIGN_STRAYFLUX = 0x1
static

Definition at line 62 of file BaselineUtils.h.

◆ STRAYFLUX_NEAREST_FOOTPRINT

template<typename ImagePixelT , typename MaskPixelT , typename VariancePixelT >
const int lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::STRAYFLUX_NEAREST_FOOTPRINT = 0x10
static

Definition at line 68 of file BaselineUtils.h.

◆ STRAYFLUX_R_TO_FOOTPRINT

template<typename ImagePixelT , typename MaskPixelT , typename VariancePixelT >
const int lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::STRAYFLUX_R_TO_FOOTPRINT = 0x8
static

Definition at line 67 of file BaselineUtils.h.

◆ STRAYFLUX_TO_POINT_SOURCES_ALWAYS

template<typename ImagePixelT , typename MaskPixelT , typename VariancePixelT >
const int lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::STRAYFLUX_TO_POINT_SOURCES_ALWAYS = 0x4
static

Definition at line 64 of file BaselineUtils.h.

◆ STRAYFLUX_TO_POINT_SOURCES_WHEN_NECESSARY

template<typename ImagePixelT , typename MaskPixelT , typename VariancePixelT >
const int lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::STRAYFLUX_TO_POINT_SOURCES_WHEN_NECESSARY = 0x2
static

Definition at line 63 of file BaselineUtils.h.

◆ STRAYFLUX_TRIM

template<typename ImagePixelT , typename MaskPixelT , typename VariancePixelT >
const int lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::STRAYFLUX_TRIM = 0x20
static

Definition at line 69 of file BaselineUtils.h.


The documentation for this class was generated from the following files: