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
Public Types | Static Public Member Functions | Static Public Attributes | List of all members
lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT > Class Template Reference

#include <Baseline.h>

Public Types

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

Static Public Member Functions

static boost::shared_ptr
< lsst::afw::detection::Footprint
symmetrizeFootprint (lsst::afw::detection::Footprint const &foot, int cx, int cy)
 
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)
 
static void medianFilter (ImageT const &img, ImageT &outimg, int halfsize)
 
static void makeMonotonic (ImageT &img, lsst::afw::detection::PeakRecord const &pk)
 
static std::vector< typename
boost::shared_ptr
< lsst::afw::image::MaskedImage
< ImagePixelT > MaskPixelT,
VariancePixelT >)> 
apportionFlux (MaskedImageT const &img, lsst::afw::detection::Footprint const &foot, std::vector< typename boost::shared_ptr< lsst::afw::image::Image< ImagePixelT > >> templates, std::vector< boost::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< boost::shared_ptr< typename lsst::afw::detection::HeavyFootprint< ImagePixelT, MaskPixelT, VariancePixelT > > > &strays, int strayFluxOptions, double clipStrayFluxFraction)
 
static bool hasSignificantFluxAtEdge (ImagePtrT, boost::shared_ptr< lsst::afw::detection::Footprint >, ImagePixelT threshold)
 
static boost::shared_ptr
< lsst::afw::detection::Footprint
getSignificantEdgePixels (ImagePtrT, boost::shared_ptr< lsst::afw::detection::Footprint >, ImagePixelT threshold)
 
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< 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)
 

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 Baseline.h.

Member Typedef Documentation

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

Definition at line 33 of file Baseline.h.

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 Baseline.h.

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

Definition at line 36 of file Baseline.h.

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 Baseline.h.

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

Definition at line 28 of file Baseline.h.

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 Baseline.h.

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

Definition at line 26 of file Baseline.h.

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 Baseline.h.

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

Definition at line 30 of file Baseline.h.

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 Baseline.h.

Member Function Documentation

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< 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 
)
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 289 of file Baseline.cc.

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 &
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 }
int y
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)
#define PTR(...)
Definition: base.h:41
for(FootprintList::const_iterator ptr=feet->begin(), end=feet->end();ptr!=end;++ptr)
Definition: saturated.cc:82
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:578
static const int STRAYFLUX_TO_POINT_SOURCES_WHEN_NECESSARY
Definition: Baseline.h:63
lsst::afw::detection::Footprint Footprint
Definition: Source.h:61
int getMinY() const
Definition: Box.h:125
int const x0
Definition: saturated.cc:45
PixelConstReference get0(int x, int y) const
Definition: Image.h:223
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
static const int STRAYFLUX_NEAREST_FOOTPRINT
Definition: Baseline.h:68
std::vector< Span::Ptr > SpanList
The Footprint&#39;s Span list.
Definition: Footprint.h:71
if(width!=gim.getWidth()||height!=gim.getHeight()||x0!=gim.getX0()||y0!=gim.getY0())
Definition: saturated.cc:47
int getMinX() const
Definition: Box.h:124
A set of pixels in an Image.
Definition: Footprint.h:62
int x
_view_t::x_iterator x_iterator
An iterator for traversing the pixels in a row.
Definition: Image.h:151
ExpressionTraits< Derived >::Iterator Iterator
Nested expression or element iterator.
static const int STRAYFLUX_TO_POINT_SOURCES_ALWAYS
Definition: Baseline.h:64
A class to represent a 2-dimensional array of pixels.
Definition: PSF.h:43
static const int STRAYFLUX_R_TO_FOOTPRINT
Definition: Baseline.h:67
boost::shared_ptr< lsst::afw::detection::HeavyFootprint< ImagePixelT > MaskPixelT, VariancePixelT > HeavyFootprintPtrT
Definition: Baseline.h:36
Iterator begin() const
Return an Iterator to the beginning of the array.
Definition: ArrayBase.h:99
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 500 of file Baseline.cc.

501  {
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 }
int y
void clip(Box2I const &other)
Shrink this to ensure that other.contains(*this).
int getMinY() const
Definition: Box.h:125
An integer coordinate rectangle.
Definition: Box.h:53
int getMaxY() const
Definition: Box.h:129
int getMinX() const
Definition: Box.h:124
_view_t::x_iterator x_iterator
An iterator for traversing the pixels in a row.
Definition: Image.h:151
int getWidth() const
Definition: Box.h:154
boost::shared_ptr< lsst::afw::image::Image< ImagePixelT > > ImagePtrT
Definition: Baseline.h:28
template<typename ImagePixelT , typename MaskPixelT = lsst::afw::image::MaskPixel, typename VariancePixelT = lsst::afw::image::VariancePixel>
std::vector< typename boost::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 boost::shared_ptr< lsst::afw::image::Image< ImagePixelT > >>  templates,
std::vector< boost::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< boost::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 583 of file Baseline.cc.

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 }
int y
void clip(Box2I const &other)
Shrink this to ensure that other.contains(*this).
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:578
lsst::afw::image::Image< ImagePixelT > ImageT
Definition: Baseline.h:27
int getMinY() const
Definition: Box.h:125
lsst::afw::image::MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > MaskedImageT
Definition: Baseline.h:25
a place to record messages and descriptions of the state of processing.
Definition: Log.h:154
static Log & getDefaultLog()
An integer coordinate rectangle.
Definition: Box.h:53
def log
Definition: log.py:85
std::vector< Span::Ptr > SpanList
The Footprint&#39;s Span list.
Definition: Footprint.h:71
Extent2I const getDimensions() const
Definition: Box.h:153
static const int ASSIGN_STRAYFLUX
Definition: Baseline.h:62
int getMaxY() const
Definition: Box.h:129
int getMinX() const
Definition: Box.h:124
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
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
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
boost::shared_ptr< lsst::afw::image::MaskedImage< ImagePixelT > MaskPixelT, VariancePixelT > MaskedImagePtrT
Definition: Baseline.h:26
int getWidth() const
Definition: Box.h:154
boost::shared_ptr< lsst::afw::image::Image< ImagePixelT > > ImagePtrT
Definition: Baseline.h:28
template<typename ImagePixelT , typename MaskPixelT , typename VariancePixelT >
std::pair< typename boost::shared_ptr< lsst::afw::image::Image< ImagePixelT > >, typename boost::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 1084 of file Baseline.cc.

1091  {
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 }
int y
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:578
lsst::afw::image::Image< ImagePixelT > ImageT
Definition: Baseline.h:27
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
static boost::shared_ptr< lsst::afw::detection::Footprint > symmetrizeFootprint(lsst::afw::detection::Footprint const &foot, int cx, int cy)
Definition: Baseline.cc:867
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
def log
Definition: log.py:85
std::vector< Span::Ptr > SpanList
The Footprint&#39;s Span list.
Definition: Footprint.h:71
boost::shared_ptr< lsst::afw::image::Mask< MaskPixelT > > MaskPtrT
Definition: Baseline.h:30
bool contains(Point2I const &point) const
Return true if the box contains the point.
int getMaxY() const
Definition: Box.h:129
int getMinX() const
Definition: Box.h:124
void include(Point2I const &point)
Expand this to ensure that this-&gt;contains(point).
int x
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
_view_t::x_iterator x_iterator
An iterator for traversing the pixels in a row.
Definition: Image.h:151
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:597
void copyWithinFootprint(Footprint const &foot, boost::shared_ptr< ImageOrMaskedImageT > const input, boost::shared_ptr< ImageOrMaskedImageT > output)
boost::shared_ptr< lsst::afw::image::Image< ImagePixelT > > ImagePtrT
Definition: Baseline.h:28
template<typename ImagePixelT , typename MaskPixelT = lsst::afw::image::MaskPixel, typename VariancePixelT = lsst::afw::image::VariancePixel>
boost::shared_ptr< det::Footprint > lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::getSignificantEdgePixels ( ImagePtrT  ,
boost::shared_ptr< lsst::afw::detection::Footprint ,
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 1317 of file Baseline.cc.

1319  {
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 }
int y
int iter
#define PTR(...)
Definition: base.h:41
for(FootprintList::const_iterator ptr=feet->begin(), end=feet->end();ptr!=end;++ptr)
Definition: saturated.cc:82
A range of pixels within one row of an Image.
Definition: Span.h:54
#define CONST_PTR(...)
Definition: base.h:47
lsst::afw::detection::Footprint Footprint
Definition: Source.h:61
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
def log
Definition: log.py:85
_const_view_t::x_iterator const_x_iterator
A const iterator for traversing the pixels in a row.
Definition: Image.h:158
A set of pixels in an Image.
Definition: Footprint.h:62
int x
int getY() const
Return the y-value.
Definition: Span.h:81
int const y0
Definition: saturated.cc:45
int getX1() const
Return the ending x-value.
Definition: Span.h:79
int getX0() const
Return the starting x-value.
Definition: Span.h:77
template<typename ImagePixelT , typename MaskPixelT = lsst::afw::image::MaskPixel, typename VariancePixelT = lsst::afw::image::VariancePixel>
bool lsst::meas::deblender::BaselineUtils< ImagePixelT, MaskPixelT, VariancePixelT >::hasSignificantFluxAtEdge ( ImagePtrT  ,
boost::shared_ptr< lsst::afw::detection::Footprint ,
ImagePixelT  threshold 
)
static

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

Definition at line 1281 of file Baseline.cc.

1283  {
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 }
int y
#define PTR(...)
Definition: base.h:41
for(FootprintList::const_iterator ptr=feet->begin(), end=feet->end();ptr!=end;++ptr)
Definition: saturated.cc:82
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
def log
Definition: log.py:85
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
A set of pixels in an Image.
Definition: Footprint.h:62
int x
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 117 of file Baseline.cc.

119  {
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 }
int y
lsst::afw::image::Image< ImagePixelT > ImageT
Definition: Baseline.h:27
int x
boost::shared_ptr< lsst::afw::image::Image< ImagePixelT > > ImagePtrT
Definition: Baseline.h:28
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 39 of file Baseline.cc.

41  {
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 }
int y
_view_t::xy_locator xy_locator
An xy_locator.
Definition: Image.h:139
_view_t::x_iterator x_iterator
An iterator for traversing the pixels in a row.
Definition: Image.h:151
template<typename ImagePixelT , typename MaskPixelT , typename VariancePixelT >
boost::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 867 of file Baseline.cc.

869  {
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 
877  pexLog::Log log(pexLog::Log::getDefaultLog(),
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 }
#define PTR(...)
Definition: base.h:41
lsst::afw::detection::Footprint Footprint
Definition: Source.h:61
int getMinY() const
Definition: Box.h:125
int getMaxX() const
Definition: Box.h:128
static boost::shared_ptr< lsst::afw::detection::Footprint > symmetrizeFootprint(lsst::afw::detection::Footprint const &foot, int cx, int cy)
Definition: Baseline.cc:867
An integer coordinate rectangle.
Definition: Box.h:53
def log
Definition: log.py:85
std::vector< Span::Ptr > SpanList
The Footprint&#39;s Span list.
Definition: Footprint.h:71
if(width!=gim.getWidth()||height!=gim.getHeight()||x0!=gim.getX0()||y0!=gim.getY0())
Definition: saturated.cc:47
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

Member Data Documentation

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

Definition at line 62 of file Baseline.h.

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

Definition at line 68 of file Baseline.h.

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

Definition at line 67 of file Baseline.h.

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

Definition at line 64 of file Baseline.h.

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

Definition at line 63 of file Baseline.h.

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

Definition at line 69 of file Baseline.h.


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