LSST Applications g0265f82a02+c6dfa2ddaf,g1162b98a3f+b2075782a9,g2079a07aa2+1b2e822518,g2bbee38e9b+c6dfa2ddaf,g337abbeb29+c6dfa2ddaf,g3ddfee87b4+a60788ef87,g50ff169b8f+2eb0e556e8,g52b1c1532d+90ebb246c7,g555ede804d+a60788ef87,g591dd9f2cf+ba8caea58f,g5ec818987f+864ee9cddb,g858d7b2824+9ee1ab4172,g876c692160+a40945ebb7,g8a8a8dda67+90ebb246c7,g8cdfe0ae6a+4fd9e222a8,g99cad8db69+5e309b7bc6,g9ddcbc5298+a1346535a5,ga1e77700b3+df8f93165b,ga8c6da7877+aa12a14d27,gae46bcf261+c6dfa2ddaf,gb0e22166c9+8634eb87fb,gb3f2274832+d0da15e3be,gba4ed39666+1ac82b564f,gbb8dafda3b+5dfd9c994b,gbeb006f7da+97157f9740,gc28159a63d+c6dfa2ddaf,gc86a011abf+9ee1ab4172,gcf0d15dbbd+a60788ef87,gdaeeff99f8+1cafcb7cd4,gdc0c513512+9ee1ab4172,ge79ae78c31+c6dfa2ddaf,geb67518f79+ba1859f325,geb961e4c1e+f9439d1e6f,gee10cc3b42+90ebb246c7,gf1cff7945b+9ee1ab4172,w.2024.12
LSST Data Management Base Package
Loading...
Searching...
No Matches
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.
 
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))
 
static void medianFilter (ImageT const &img, ImageT &outimg, int halfsize)
 Run a spatial median filter over the given input img, writing the results to out.
 
static void 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.
 
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.
 
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.
 
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.
 
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 , typename VariancePixelT >
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
408
409 // when doing stray flux: the footprints and pixels, which we'll
410 // combine into the return 'strays' HeavyFootprint at the end.
412 std::vector<std::vector<afwGeom::Span>> straySpans(tfoots.size());
416
417 int ix0 = img.getX0();
418 int iy0 = img.getY0();
419 geom::Box2I sumbb = tsum->getBBox();
420 int sumx0 = sumbb.getMinX();
421 int sumy0 = sumbb.getMinY();
422
423 for (size_t i=0; i<tfoots.size(); ++i) {
428 }
429
430 bool always = (strayFluxOptions & STRAYFLUX_TO_POINT_SOURCES_ALWAYS);
431
432 typedef std::uint16_t itype;
434
435 if (strayFluxOptions & STRAYFLUX_NEAREST_FOOTPRINT) {
436 // Compute the map of which footprint is closest to each
437 // pixel in the bbox.
438 typedef std::uint16_t dtype;
441
444
445 if (!always && ispsf.size()) {
446 // create a temp list that has empty footprints in place
447 // of all the point sources.
448 auto empty = std::make_shared<det::Footprint>();
449 empty->setPeakSchema(foot.getPeaks().getSchema());
450 for (size_t i=0; i<tfoots.size(); ++i) {
451 if (ispsf[i]) {
452 templist.push_back(empty);
453 } else {
454 templist.push_back(tfoots[i]);
455 }
456 }
457 footlist = &templist;
458 }
459 nearestFootprint(*footlist, nearest, dist);
460 }
461
462 // Go through the (parent) Footprint looking for stray flux:
463 // pixels that are not claimed by any template, and positive.
464 for (afwGeom::Span const & s : *foot.getSpans()) {
465 int y = s.getY();
466 int x0 = s.getX0();
467 int x1 = s.getX1();
468 typename ImageT::x_iterator tsum_it =
469 tsum->row_begin(y - sumy0) + (x0 - sumx0);
470 typename MaskedImageT::x_iterator in_it =
471 img.row_begin(y - iy0) + (x0 - ix0);
472 double contrib[tfoots.size()];
473
474 for (int x = x0; x <= x1; ++x, ++tsum_it, ++in_it) {
475 // Skip pixels that are covered by at least one
476 // template (*tsum_it > 0) or the input is not
477 // positive (*in_it <= 0).
478 if ((*tsum_it > 0) || (*in_it).image() <= 0) {
479 continue;
480 }
481
482 if (strayFluxOptions & STRAYFLUX_R_TO_FOOTPRINT) {
483 // we'll compute these just-in-time
484 for (size_t i=0; i<tfoots.size(); ++i) {
485 contrib[i] = -1.0;
486 }
487 } else if (strayFluxOptions & STRAYFLUX_NEAREST_FOOTPRINT) {
488 for (size_t i=0; i<tfoots.size(); ++i) {
489 contrib[i] = 0.0;
490 }
491 int i = (*nearest)[geom::Point2I(x, y)];
492 contrib[i] = 1.0;
493 } else {
494 // R_TO_PEAK
495 for (size_t i=0; i<tfoots.size(); ++i) {
496 // Split the stray flux by 1/(1+r^2) to peaks
497 int dx, dy;
498 dx = pkx[i] - x;
499 dy = pky[i] - y;
500 contrib[i] = 1. / (1. + dx*dx + dy*dy);
501 }
502 }
503
504 // Round 1: skip point sources unless STRAYFLUX_TO_POINT_SOURCES_ALWAYS
505 // are we going to assign stray flux to ptsrcs?
506 bool ptsrcs = always;
507 double csum = 0.;
508 for (size_t i=0; i<tfoots.size(); ++i) {
509 // if we're skipping point sources and this is a point source...
510 if ((!ptsrcs) && ispsf.size() && ispsf[i]) {
511 continue;
512 }
513 if (contrib[i] == -1.0) {
514 contrib[i] = _get_contrib_r_to_footprint(x, y, tfoots[i]);
515 }
516 csum += contrib[i];
517 }
518 if ((csum == 0.) &&
519 (strayFluxOptions &
521 // No extended sources -- assign to pt sources
522 //LOGL_DEBUG(_log, "necessary to assign stray flux to point sources");
523 ptsrcs = true;
524 for (size_t i=0; i<tfoots.size(); ++i) {
525 if (contrib[i] == -1.0) {
526 contrib[i] = _get_contrib_r_to_footprint(x, y, tfoots[i]);
527 }
528 csum += contrib[i];
529 }
530 }
531
532 // Drop small contributions...
533 double strayclip = (clipStrayFluxFraction * csum);
534 csum = 0.;
535 for (size_t i=0; i<tfoots.size(); ++i) {
536 // skip ptsrcs?
537 if ((!ptsrcs) && ispsf.size() && ispsf[i]) {
538 contrib[i] = 0.;
539 continue;
540 }
541 // skip small contributions
542 if (contrib[i] < strayclip) {
543 contrib[i] = 0.;
544 continue;
545 }
546 csum += contrib[i];
547 }
548
549 for (size_t i=0; i<tfoots.size(); ++i) {
550 if (contrib[i] == 0.) {
551 continue;
552 }
553 // the stray flux to give to template i
554 double p = (contrib[i] / csum) * (*in_it).image();
555
556 if (!strayfoot[i]) {
557 strayfoot[i] = std::make_shared<det::Footprint>();
558 strayfoot[i]->setPeakSchema(foot.getPeaks().getSchema());
559 }
560 straySpans[i].push_back(afwGeom::Span(y, x, x));
561 straypix[i].push_back(p);
562 straymask[i].push_back((*in_it).mask());
563 strayvar[i].push_back((*in_it).variance());
564 }
565 }
566 }
567
568 // Store the stray flux in HeavyFootprints
569 for (size_t i=0; i<tfoots.size(); ++i) {
570 if (strayfoot[i]) {
571 strayfoot[i]->setSpans(std::make_shared<afwGeom::SpanSet>(straySpans[i]));
572 }
573 if (!strayfoot[i]) {
575 } else {
579 HeavyFootprintPtrT heavy(new HeavyFootprint(*strayfoot[i]));
580 ndarray::Array<ImagePixelT,1,1> himg = heavy->getImageArray();
581 typename std::vector<ImagePixelT>::const_iterator spix;
582 typename std::vector<MaskPixelT>::const_iterator smask;
583 typename std::vector<VariancePixelT>::const_iterator svar;
584 typename ndarray::Array<ImagePixelT,1,1>::Iterator hpix;
585 typename ndarray::Array<MaskPixelT,1,1>::Iterator mpix;
586 typename ndarray::Array<VariancePixelT,1,1>::Iterator vpix;
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}
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.
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
std::shared_ptr< lsst::afw::detection::HeavyFootprint< ImagePixelT, MaskPixelT, VariancePixelT > > HeavyFootprintPtrT
static const int STRAYFLUX_TO_POINT_SOURCES_ALWAYS
T end(T... args)
Point< int, 2 > Point2I
Definition Point.h:321
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
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 <em>strayFluxOptions</em> includes <em>STRAYFLUX_TO_POINT_SOURCES_ALWAYS</em>, then point sources are always included in the 1/(1+r^2) splitting. Otherwise, if <em>STRAYFLUX_TO_POINT_SOURCES_WHEN_NECESSARY</em>, point sources are included only if there are no extended sources nearby. If any stray-flux portion is less than <em>clipStrayFluxFraction</em>, 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
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
lsst::afw::image::Image< ImagePixelT > ImageT
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

◆ 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
1255 afwGeom::SpanSet::const_iterator fwd = spans.begin();
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
const_MaskedImageLocator< typename Image::xy_locator, typename Mask::xy_locator, typename Variance::xy_locator > const_xy_locator
A const_locator for a MaskedImage.
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
std::shared_ptr< lsst::afw::image::Mask< MaskPixelT > > MaskPtrT
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();
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:117

◆ 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}

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