#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.

                    {
 
     typedef typename det::Footprint::SpanList SpanList;
     typedef typename det::HeavyFootprint<ImagePixelT, MaskPixelT, VariancePixelT> HeavyFootprint;
     typedef typename boost::shared_ptr< HeavyFootprint > HeavyFootprintPtrT;
 
     // when doing stray flux: the footprints and pixels, which we'll
     // combine into the return 'strays' HeavyFootprint at the end.
     std::vector<PTR(det::Footprint) > strayfoot;
     std::vector<std::vector<ImagePixelT> > straypix;
     std::vector<std::vector<MaskPixelT> > straymask;
     std::vector<std::vector<VariancePixelT> > strayvar;
 
     int ix0 = img.getX0();
     int iy0 = img.getY0();
     geom::Box2I sumbb = tsum->getBBox();
     int sumx0 = sumbb.getMinX();
     int sumy0 = sumbb.getMinY();
 
     for (size_t i=0; i<tfoots.size(); ++i) {
         strayfoot.push_back(PTR(det::Footprint)());
         straypix.push_back(std::vector<ImagePixelT>());
         straymask.push_back(std::vector<MaskPixelT>());
         strayvar.push_back(std::vector<VariancePixelT>());
     }
 
     bool always = (strayFluxOptions & STRAYFLUX_TO_POINT_SOURCES_ALWAYS);
 
     typedef boost::uint16_t itype;
     PTR(image::Image<itype>) nearest;
 
     if (strayFluxOptions & STRAYFLUX_NEAREST_FOOTPRINT) {
         // Compute the map of which footprint is closest to each
         // pixel in the bbox.
         typedef boost::uint16_t dtype;
         PTR(image::Image<dtype>) dist(new image::Image<dtype>(sumbb));
         nearest = PTR(image::Image<itype>)(new image::Image<itype>(sumbb));
 
         std::vector<PTR(det::Footprint)> templist;
         std::vector<PTR(det::Footprint)>* footlist = &tfoots;
 
         if (!always && ispsf.size()) {
             // create a temp list that has empty footprints in place
             // of all the point sources.
             PTR(det::Footprint) empty(new det::Footprint(foot.getPeaks().getSchema()));
             for (size_t i=0; i<tfoots.size(); ++i) {
                 if (ispsf[i]) {
                     templist.push_back(empty);
                 } else {
                     templist.push_back(tfoots[i]);
                 }
             }
             footlist = &templist;
         }
         nearestFootprint(*footlist, nearest, dist);
     }
 
     // Go through the (parent) Footprint looking for stray flux:
     // pixels that are not claimed by any template, and positive.
     const SpanList& spans = foot.getSpans();
     for (SpanList::const_iterator s = spans.begin(); s < spans.end(); s++) {
         int y = (*s)->getY();
         int x0 = (*s)->getX0();
         int x1 = (*s)->getX1();
         typename ImageT::x_iterator tsum_it =
             tsum->row_begin(y - sumy0) + (x0 - sumx0);
         typename MaskedImageT::x_iterator in_it =
             img.row_begin(y - iy0) + (x0 - ix0);
         double contrib[tfoots.size()];
 
         for (int x = x0; x <= x1; ++x, ++tsum_it, ++in_it) {
             // Skip pixels that are covered by at least one
             // template (*tsum_it > 0) or the input is not
             // positive (*in_it <= 0).
             if ((*tsum_it > 0) || (*in_it).image() <= 0) {
                 continue;
             }
 
             if (strayFluxOptions & STRAYFLUX_R_TO_FOOTPRINT) {
                 // we'll compute these just-in-time
                 for (size_t i=0; i<tfoots.size(); ++i) {
                     contrib[i] = -1.0;
                 }
             } else if (strayFluxOptions & STRAYFLUX_NEAREST_FOOTPRINT) {
                 for (size_t i=0; i<tfoots.size(); ++i) {
                     contrib[i] = 0.0;
                 }
                 int i = nearest->get0(x, y);
                 contrib[i] = 1.0;
             } else {
                 // R_TO_PEAK
                 for (size_t i=0; i<tfoots.size(); ++i) {
                     // Split the stray flux by 1/(1+r^2) to peaks
                     int dx, dy;
                     dx = pkx[i] - x;
                     dy = pky[i] - y;
                     contrib[i] = 1. / (1. + dx*dx + dy*dy);
                 }
             }
 
             // Round 1: skip point sources unless STRAYFLUX_TO_POINT_SOURCES_ALWAYS
             // are we going to assign stray flux to ptsrcs?
             bool ptsrcs = always;
             double csum = 0.;
             for (size_t i=0; i<tfoots.size(); ++i) {
                 // if we're skipping point sources and this is a point source...
                 if ((!ptsrcs) && ispsf.size() && ispsf[i]) {
                     continue;
                 }
                 if (contrib[i] == -1.0) {
                     contrib[i] = _get_contrib_r_to_footprint(x, y, tfoots[i]);
                 }
                 csum += contrib[i];
             }
             if ((csum == 0.) &&
                 (strayFluxOptions &
                  STRAYFLUX_TO_POINT_SOURCES_WHEN_NECESSARY)) {
                 // No extended sources -- assign to pt sources
                 //log.debugf("necessary to assign stray flux to point sources");
                 ptsrcs = true;
                 for (size_t i=0; i<tfoots.size(); ++i) {
                     if (contrib[i] == -1.0) {
                         contrib[i] = _get_contrib_r_to_footprint(x, y, tfoots[i]);
                     }
                     csum += contrib[i];
                 }
             }
 
             // Drop small contributions...
             double strayclip = (clipStrayFluxFraction * csum);
             csum = 0.;
             for (size_t i=0; i<tfoots.size(); ++i) {
                 // skip ptsrcs?
                 if ((!ptsrcs) && ispsf.size() && ispsf[i]) {
                     contrib[i] = 0.;
                     continue;
                 }
                 // skip small contributions
                 if (contrib[i] < strayclip) {
                     contrib[i] = 0.;
                     continue;
                 }
                 csum += contrib[i];
             }
 
             for (size_t i=0; i<tfoots.size(); ++i) {
                 if (contrib[i] == 0.) {
                     continue;
                 }
                 // the stray flux to give to template i
                 double p = (contrib[i] / csum) * (*in_it).image();
                 if (!strayfoot[i]) {
                     strayfoot[i] = PTR(det::Footprint)(new det::Footprint(foot.getPeaks().getSchema()));
                 }
                 strayfoot[i]->addSpanInSeries(y, x, x); // XXX this may be rather heavy in memory use
                 straypix[i].push_back(p);
                 straymask[i].push_back((*in_it).mask());
                 strayvar[i].push_back((*in_it).variance());
             }
         }
     }
 
     // Store the stray flux in HeavyFootprints
     for (size_t i=0; i<tfoots.size(); ++i) {
         if (!strayfoot[i]) {
             strays.push_back(HeavyFootprintPtrT());
         } else {
             HeavyFootprintPtrT heavy(new HeavyFootprint(*strayfoot[i]));
             ndarray::Array<ImagePixelT,1,1> himg = heavy->getImageArray();
             typename std::vector<ImagePixelT>::const_iterator spix;
             typename std::vector<MaskPixelT>::const_iterator smask;
             typename std::vector<VariancePixelT>::const_iterator svar;
             typename ndarray::Array<ImagePixelT,1,1>::Iterator hpix;
             typename ndarray::Array<MaskPixelT,1,1>::Iterator mpix;
             typename ndarray::Array<VariancePixelT,1,1>::Iterator vpix;
 
             assert((size_t)strayfoot[i]->getNpix() == straypix[i].size());
 
             for (spix = straypix[i].begin(),
                      smask = straymask[i].begin(),
                      svar  = strayvar [i].begin(),
                      hpix = himg.begin(),
                      mpix = heavy->getMaskArray().begin(),
                      vpix = heavy->getVarianceArray().begin();
                  spix != straypix[i].end();
                  ++spix, ++smask, ++svar, ++hpix, ++mpix, ++vpix) {
                 *hpix = *spix;
                 *mpix = *smask;
                 *vpix = *svar;
             }
             strays.push_back(heavy);
         }
     }
 }

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.

                                {
     geom::Box2I sumbb = tsum->getBBox();
     int sumx0 = sumbb.getMinX();
     int sumy0 = sumbb.getMinY();
 
     // Compute  tsum = the sum of templates
     for (size_t i=0; i<timgs.size(); ++i) {
         ImagePtrT timg = timgs[i];
         geom::Box2I tbb = timg->getBBox();
         int tx0 = tbb.getMinX();
         int ty0 = tbb.getMinY();
         // To handle "ramped" templates that can extend outside the
         // parent, clip the bbox.  Note that we saved tx0,ty0 BEFORE
         // doing this!
         tbb.clip(sumbb);
         int copyx0 = tbb.getMinX();
         // Here we iterate over the template bbox -- we could instead
         // iterate over the "tfoot"s.
         for (int y=tbb.getMinY(); y<=tbb.getMaxY(); ++y) {
             typename ImageT::x_iterator in_it = timg->row_begin(y - ty0) + (copyx0 - tx0);
             typename ImageT::x_iterator inend = in_it + tbb.getWidth();
             typename ImageT::x_iterator tsum_it =
                 tsum->row_begin(y - sumy0) + (copyx0 - sumx0);
             for (; in_it != inend; ++in_it, ++tsum_it) {
                 *tsum_it += std::max((ImagePixelT)0., static_cast<ImagePixelT>(*in_it));
             }
         }
     }
 
 }

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.

       {
     typedef typename det::Footprint::SpanList SpanList;
 
     if (timgs.size() != tfoots.size()) {
         throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
             (boost::format("Template images must be the same length as template footprints (%d vs %d)")
                 % timgs.size() % tfoots.size()).str());
     }
 
     for (size_t i=0; i<timgs.size(); ++i) {
         if (!timgs[i]->getBBox().contains(tfoots[i]->getBBox())) {
             throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                               "Template image MUST contain template footprint");
         }
         if (!img.getBBox().contains(foot.getBBox())) {
             throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                               "Image bbox MUST contain parent footprint");
         }
         // template bounding-boxes *can* extend outside the parent
         // footprint if we are ramping templates with significant flux
         // at the edges.  We handle this below.
     }
 
     // the apportioned flux return value
     std::vector<MaskedImagePtrT> portions;
 
     pexLog::Log log(pexLog::Log::getDefaultLog(),
                     "lsst.meas.deblender.apportionFlux");
     bool findStrayFlux = (strayFluxOptions & ASSIGN_STRAYFLUX);
 
     int ix0 = img.getX0();
     int iy0 = img.getY0();
     geom::Box2I fbb = foot.getBBox();
 
     if (!tsum) {
         tsum = ImagePtrT(new ImageT(fbb.getDimensions()));
         tsum->setXY0(fbb.getMinX(), fbb.getMinY());
     }
 
     if (!tsum->getBBox().contains(foot.getBBox())) {
         throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
                           "Template sum image MUST contain parent footprint");
     }
 
     geom::Box2I sumbb = tsum->getBBox();
     int sumx0 = sumbb.getMinX();
     int sumy0 = sumbb.getMinY();
 
     _sum_templates(timgs, tsum);
 
     // Compute flux portions
     for (size_t i=0; i<timgs.size(); ++i) {
         ImagePtrT timg = timgs[i];
         // Initialize return value:
         MaskedImagePtrT port(new MaskedImageT(timg->getDimensions()));
         port->setXY0(timg->getXY0());
         portions.push_back(port);
 
         // Split flux = image * template / tsum
         geom::Box2I tbb = timg->getBBox();
         int tx0 = tbb.getMinX();
         int ty0 = tbb.getMinY();
         // As above
         tbb.clip(sumbb);
         int copyx0 = tbb.getMinX();
         for (int y=tbb.getMinY(); y<=tbb.getMaxY(); ++y) {
             typename MaskedImageT::x_iterator in_it =
                 img.row_begin(y - iy0) + (copyx0 - ix0);
             typename ImageT::x_iterator tptr =
                 timg->row_begin(y - ty0) + (copyx0 - tx0);
             typename ImageT::x_iterator tend = tptr + tbb.getWidth();
             typename ImageT::x_iterator tsum_it =
                 tsum->row_begin(y - sumy0) + (copyx0 - sumx0);
             typename MaskedImageT::x_iterator out_it =
                 port->row_begin(y - ty0) + (copyx0 - tx0);
             for (; tptr != tend; ++tptr, ++in_it, ++out_it, ++tsum_it) {
                 if (*tsum_it == 0) {
                     continue;
                 }
                 double frac = std::max((ImagePixelT)0., static_cast<ImagePixelT>(*tptr)) / (*tsum_it);
                 //if (frac == 0) {
                 // treat mask planes differently?
                 // }
                 out_it.mask()     = (*in_it).mask();
                 out_it.variance() = (*in_it).variance()*frac;
                 out_it.image()    = (*in_it).image() * frac;
             }
         }
     }
 
     if (findStrayFlux) {
         if ((ispsf.size() > 0) && (ispsf.size() != timgs.size())) {
             throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
                 (boost::format("'ispsf' must be the same length as templates (%d vs %d)")
                      % ispsf.size() % timgs.size()).str());
         }
         if ((pkx.size() != timgs.size()) || (pky.size() != timgs.size())) {
             throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
                 (boost::format("'pkx' and 'pky' must be the same length as templates (%d,%d vs %d)")
                     % pkx.size() % pky.size() % timgs.size()).str());
         }
         _find_stray_flux(foot, tsum, img, strayFluxOptions, tfoots,
                          ispsf, pkx, pky, clipStrayFluxFraction, strays);
     }
     return portions;
 }

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.

                         {
 
     typedef typename MaskedImageT::const_xy_locator xy_loc;
     typedef typename det::Footprint::SpanList SpanList;
 
     *patchedEdges = false;
 
     int cx = peak.getIx();
     int cy = peak.getIy();
 
     pexLog::Log log(pexLog::Log::getDefaultLog(),
                     "lsst.meas.deblender.symmetricFootprint");
 
     if (!img.getBBox(image::PARENT).contains(foot.getBBox())) {
         throw LSST_EXCEPT(lsst::pex::exceptions::LengthError, "Image too small for footprint");
     }
 
     FootprintPtrT sfoot = symmetrizeFootprint(foot, cx, cy);
     if (!sfoot) {
         return std::pair<ImagePtrT, FootprintPtrT>(ImagePtrT(), sfoot);
     }
 
     if (!img.getBBox(image::PARENT).contains(sfoot->getBBox())) {
         throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
                           "Image too small for symmetrized footprint");
     }
     const SpanList spans = sfoot->getSpans();
 
     // does this footprint touch an EDGE?
     bool touchesEdge = false;
     if (patchEdge) {
         log.debugf("Checking footprint for EDGE bits");
         MaskPtrT mask = img.getMask();
         bool edge = false;
         MaskPixelT edgebit = mask->getPlaneBitMask("EDGE");
         for (SpanList::const_iterator fwd=spans.begin();
              fwd != spans.end(); ++fwd) {
             int x0 = (*fwd)->getX0();
             int x1 = (*fwd)->getX1();
             typename MaskT::x_iterator xiter =
                 mask->x_at(x0 - mask->getX0(), (*fwd)->getY() - mask->getY0());
             for (int x=x0; x<=x1; ++x, ++xiter) {
                 if ((*xiter) & edgebit) {
                     edge = true;
                     break;
                 }
             }
             if (edge)
                 break;
         }
         if (edge) {
             log.debugf("Footprint includes an EDGE pixel.");
             touchesEdge = true;
         }
     }
 
     // The result image:
     ImagePtrT targetimg(new ImageT(sfoot->getBBox()));
 
     SpanList::const_iterator fwd  = spans.begin();
     SpanList::const_iterator back = spans.end()-1;
 
     ImagePtrT theimg = img.getImage();
 
     for (; fwd <= back; fwd++, back--) {
         int fy = (*fwd)->getY();
         int by = (*back)->getY();
 
         for (int fx=(*fwd)->getX0(), bx=(*back)->getX1();
              fx <= (*fwd)->getX1();
              fx++, bx--) {
             // FIXME -- CURRENTLY WE IGNORE THE MASK PLANE!  options
             // include ORing the mask bits, or being clever about
             // ignoring some masked pixels, or copying the mask bits
             // of the min pixel
 
             // We have already checked the bounding box, so this should always be satisfied
             assert(theimg->getBBox(image::PARENT).contains(geom::Point2I(fx, fy)));
             assert(theimg->getBBox(image::PARENT).contains(geom::Point2I(bx, by)));
 
             // FIXME -- we could do this with image iterators instead.
             // But first profile to show that it's necessary and an
             // improvement.
             ImagePixelT pixf = theimg->get0(fx, fy);
             ImagePixelT pixb = theimg->get0(bx, by);
             ImagePixelT pix = std::min(pixf, pixb);
             if (minZero) {
                 pix = std::max(pix, static_cast<ImagePixelT>(0));
             }
             targetimg->set0(fx, fy, pix);
             targetimg->set0(bx, by, pix);
 
         }
     }
 
 
     if (touchesEdge) {
         // Find spans whose mirrors fall outside the image bounds,
         // grow the footprint to include those spans, and plug in
         // their pixel values.
         geom::Box2I bb = sfoot->getBBox();
 
         // Actually, it's not necessarily the IMAGE bounds that count
         //-- the footprint may not go right to the image edge.
         //geom::Box2I imbb = img.getBBox();
         geom::Box2I imbb = foot.getBBox();
 
         log.debugf("Footprint touches EDGE: start bbox [%i,%i],[%i,%i]",
                    bb.getMinX(), bb.getMaxX(), bb.getMinY(), bb.getMaxY());
         // original footprint spans
         const SpanList ospans = foot.getSpans();
         for (fwd = ospans.begin(); fwd != ospans.end(); ++fwd) {
             int y = (*fwd)->getY();
             int x = (*fwd)->getX0();
             // mirrored coords
             int ym = cy + (cy - y);
             int xm = cx + (cx - x);
             if (!imbb.contains(geom::Point2I(xm, ym))) {
                 bb.include(geom::Point2I(x, y));
             }
             x = (*fwd)->getX1();
             xm = cx + (cx - x);
             if (!imbb.contains(geom::Point2I(xm, ym))) {
                 bb.include(geom::Point2I(x, y));
             }
         }
         log.debugf("Footprint touches EDGE: grown bbox [%i,%i],[%i,%i]",
                    bb.getMinX(), bb.getMaxX(), bb.getMinY(), bb.getMaxY());
 
         // New template image
         ImagePtrT targetimg2(new ImageT(bb));
         copyWithinFootprint(*sfoot, targetimg, targetimg2);
 
         log.debugf("Symmetric footprint spans:");
         const SpanList sspans = sfoot->getSpans();
         for (fwd = sspans.begin(); fwd != sspans.end(); ++fwd) {
             log.debugf("  %s", (*fwd)->toString().c_str());
         }
 
         // copy original 'img' pixels for the portion of spans whose
         // mirrors are out of bounds.
         for (fwd = ospans.begin(); fwd != ospans.end(); ++fwd) {
             int y   = (*fwd)->getY();
             int x0  = (*fwd)->getX0();
             int x1 = (*fwd)->getX1();
             // mirrored coords
             int ym  = cy + (cy - y);
             int xm0 = cx + (cx - x0);
             int xm1 = cx + (cx - x1);
             bool in0 = imbb.contains(geom::Point2I(xm0, ym));
             bool in1 = imbb.contains(geom::Point2I(xm1, ym));
             if (in0 && in1) {
                 // both endpoints of the symmetric span are in bounds; nothing to do
                 continue;
             }
             // clip to the part of the span where the mirror is out of bounds
             if (in0) {
                 // the mirror of x0 is in-bounds; move x0 to be the first pixel
                 // whose mirror would be out-of-bounds
                 x0 = cx + (cx - (imbb.getMinX() - 1));
             }
             if (in1) {
                 x1 = cx + (cx - (imbb.getMaxX() + 1));
             }
             log.debugf("Span y=%i, x=[%i,%i] has mirror (%i,[%i,%i]) out-of-bounds; clipped to %i,[%i,%i]",
                        y, (*fwd)->getX0(), (*fwd)->getX1(), ym, xm1, xm0, y, x0, x1);
             typename MaskedImageT::x_iterator initer =
                 img.x_at(x0 - img.getX0(), y - img.getY0());
             typename ImageT::x_iterator outiter =
                 targetimg2->x_at(x0 - targetimg2->getX0(), y - targetimg2->getY0());
             for (int x=x0; x<=x1; ++x, ++outiter, ++initer) {
                 *outiter = initer.image();
             }
             sfoot->addSpan(y, x0, x1);
         }
         sfoot->normalize();
         targetimg = targetimg2;
     }
 
     *patchedEdges = touchesEdge;
     return std::pair<ImagePtrT, FootprintPtrT>(targetimg, sfoot);
 }

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.

                                              {
     pexLog::Log log(pexLog::Log::getDefaultLog(),
                     "lsst.meas.deblender.getSignificantEdgePixels");
 
     sfoot->normalize();                 // Algorithms require spans sorted by y
 
     PTR(det::Footprint) significant(new det::Footprint(sfoot->getPeaks().getSchema()));
 
     int const x0 = img->getX0(), y0 = img->getY0();
     CONST_PTR(det::Footprint) const edges = sfoot->findEdgePixels();
     for (det::Footprint::SpanList::const_iterator ss = edges->getSpans().begin();
          ss != edges->getSpans().end(); ++ss) {
         geom::Span const& span = **ss;
         int const y = span.getY();
         int x = span.getX0();
         typename ImageT::const_x_iterator iter = img->x_at(x - x0, y - y0);
         bool onSpan = false;            // Are we in a span of interest
         int xSpan;                      // Starting x of span
         for (; x <= span.getX1(); ++x, ++iter) {
             if (*iter >= thresh) {
                 onSpan = true;
                 xSpan = x;
             } else if (onSpan) {
                 onSpan = false;
                 significant->addSpanInSeries(y, xSpan, x - 1);
             }
         }
         if (onSpan) {
             significant->addSpanInSeries(y, xSpan, span.getX1());
         }
     }
 
     return significant;
 }

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.

                                              {
     typedef typename det::Footprint::SpanList SpanList;
 
     pexLog::Log log(pexLog::Log::getDefaultLog(),
                     "lsst.meas.deblender.hasSignificantFluxAtEdge");
 
     // Find edge template pixels with significant flux -- perhaps
     // because their symmetric pixels were outside the footprint?
     // (clipped by an image edge, etc)
     PTR(det::Footprint) edge = sfoot->findEdgePixels();
 
     const SpanList spans = edge->getSpans();
     for (SpanList::const_iterator sp = spans.begin(); sp != spans.end(); ++sp) {
         int const y  = (*sp)->getY();
         int const x0 = (*sp)->getX0();
         int const x1 = (*sp)->getX1();
         int x;
         typename ImageT::const_x_iterator xiter;
         for (xiter = img->x_at(x0 - img->getX0(), y - img->getY0()), x=x0; x<=x1; ++x, ++xiter) {
             if (*xiter >= thresh) {
                 return true;
             }
         }
     }
     return false;
 }

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.

                                {
 
     int cx = peak.getIx();
     int cy = peak.getIy();
     int ix0 = img.getX0();
     int iy0 = img.getY0();
     int iW = img.getWidth();
     int iH = img.getHeight();
 
     ImagePtrT shadowingImg = ImagePtrT(new ImageT(img, true));
 
     int DW = std::max(cx - img.getX0(), img.getX0() + img.getWidth() - cx);
     int DH = std::max(cy - img.getY0(), img.getY0() + img.getHeight() - cy);
 
     const int S = 5;
 
     // Work out from the peak in chunks of "S" pixels.
     int s;
     for (s = 0; s < std::max(DW,DH); s += S) {
         int p;
         for (p=0; p<S; p++) {
             // visit pixels with L_inf distance = s + p from the
             // center (ie, the s+p'th square ring of pixels)
             // L is the half-length of the ring (box).
             int L = s+p;
             int x = L, y = -L;
             int dx = 0, dy = 0; // initialized here to satisfy the
                                 // compiler; initialized for real
                                 // below (first time through loop)
             /*
             int i;
             int leg;
             for (i=0; i<(8*L); i++, x += dx, y += dy) {
                 if (i % (2*L) == 0) {
                     leg = (i/(2*L));
                     dx = ( leg    % 2) * (-1 + 2*(leg/2));
                     dy = ((leg+1) % 2) * ( 1 - 2*(leg/2));
                 }
              */
 
             /*
              We visit pixels in a box of "radius" L, in this order:
 
              L=1:
 
              4 3 2
              5   1
              6 7 0
 
              L=2:
 
              8  7  6  5  4
              9           3
              10          2
              11          1
              12 13 14 15 0
 
              Note that the number of pixel visited is 8*L, and that we
              change "dx" or "dy" each "2*L" steps.
              */
             for (int i=0; i<(8*L); i++, x += dx, y += dy) {
                 // time to change directions?  (Note that this runs
                 // the first time through the loop, initializing dx,dy
                 // appropriately.)
                 if (i % (2*L) == 0) {
                     int leg = (i/(2*L));
                     // dx = [ 0, -1,  0, 1 ][leg]
                     dx = ( leg    % 2) * (-1 + 2*(leg/2));
                     // dy = [ 1,  0, -1, 0 ][leg]
                     dy = ((leg+1) % 2) * ( 1 - 2*(leg/2));
                 }
                 //printf("  i=%i, leg=%i, dx=%i, dy=%i, x=%i, y=%i\n", i, leg, dx, dy, x, y);
                 int px = cx + x - ix0;
                 int py = cy + y - iy0;
                 // If the shadowing pixel is out of bounds, nothing to do.
                 if (px < 0 || px >= iW || py < 0 || py >= iH)
                     continue;
                 // The pixel casting the shadow
                 ImagePixelT pix = (*shadowingImg)(px,py);
 
                 // Cast this pixel's shadow S pixels long in a cone.
                 // We compute the range of slopes (or inverse-slopes)
                 // shadowed by the pixel, [ds0,ds1]
                 double ds0, ds1;
                 // Range of slopes shadowed
                 const double A = 0.3;
                 int shx, shy;
                 int psx, psy;
                 // Are we traversing a vertical edge of the box?
                 if (dx == 0) {
                     // (if so, then "x" is +- L, so no div-by-zero)
                     ds0 = (double(y) / double(x)) - A;
                     ds1 = ds0 + 2.0 * A;
                     // cast the shadow on column x + sign(x)*shx
                     for (shx=1; shx<=S; shx++) {
                         int xsign = (x>0?1:-1);
                         // the column being shadowed
                         psx = cx + x + (xsign*shx) - ix0;
                         if (psx < 0 || psx >= iW)
                             continue;
                         // shadow covers a range of y values based on slope
                         for (shy  = iround(shx * ds0);
                              shy <= iround(shx * ds1); shy++) {
                             psy = cy + y + xsign*shy - iy0;
                             if (psy < 0 || psy >= iH)
                                 continue;
                             img(psx, psy) = std::min(img(psx, psy), pix);
                         }
                     }
 
                 } else {
                     // We're traversing a horizontal edge of the box; y = +-L
                     ds0 = (double(x) / double(y)) - A;
                     ds1 = ds0 + 2.0 * A;
                     // Cast shadow on row y + sign(y)*shy
                     for (shy=1; shy<=S; shy++) {
                         int ysign = (y>0?1:-1);
                         psy = cy + y + (ysign*shy) - iy0;
                         if (psy < 0 || psy >= iH)
                             continue;
                         // shadow covers a range of x vals based on slope
                         for (shx  = iround(shy * ds0);
                              shx <= iround(shy * ds1); shx++) {
                             psx = cx + x + ysign*shx - ix0;
                             if (psx < 0 || psx >= iW)
                                 continue;
                             img(psx, psy) = std::min(img(psx, psy), pix);
                         }
                     }
                 }
             }
         }
         //*shadowingImg <<= *img;
         shadowingImg->operator<<=(img);
     }
 }

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.

                            {
     int S = halfsize*2 + 1;
     int SS = S*S;
     typedef typename ImageT::xy_locator xy_loc;
     xy_loc pix = img.xy_at(halfsize,halfsize);
     std::vector<typename xy_loc::cached_location_t> locs;
     for (int i=0; i<S; ++i) {
         for (int j=0; j<S; ++j) {
             locs.push_back(pix.cache_location(j-halfsize, i-halfsize));
         }
     }
     int W = img.getWidth();
     int H = img.getHeight();
     ImagePixelT vals[S*S];
     for (int y=halfsize; y<H-halfsize; ++y) {
         xy_loc inpix = img.xy_at(halfsize, y), end = img.xy_at(W-halfsize, y);
         for (typename ImageT::x_iterator optr = out.row_begin(y) + halfsize;
              inpix != end; ++inpix.x(), ++optr) {
             for (int i=0; i<SS; ++i)
                 vals[i] = inpix[locs[i]];
             std::nth_element(vals, vals+SS/2, vals+SS);
             *optr = vals[SS/2];
         }
     }
 
     // grumble grumble margins
     for (int y=0; y<2*halfsize; ++y) {
         int iy = y;
         if (y >= halfsize)
             iy = H - 1 - (y-halfsize);
         typename ImageT::x_iterator optr = out.row_begin(iy);
         typename ImageT::x_iterator iptr = img.row_begin(iy), end=img.row_end(iy);
         for (; iptr != end; ++iptr,++optr)
             *optr = *iptr;
     }
     for (int y=halfsize; y<H-halfsize; ++y) {
         typename ImageT::x_iterator optr = out.row_begin(y);
         typename ImageT::x_iterator iptr = img.row_begin(y), end=img.row_begin(y)+halfsize;
         for (; iptr != end; ++iptr,++optr)
             *optr = *iptr;
         iptr = img.row_begin(y) + ((W-1) - halfsize);
         end  = img.row_begin(y) + (W-1);
         optr = out.row_begin(y) + ((W-1) - halfsize);
         for (; iptr != end; ++iptr,++optr)
             *optr = *iptr;
     }
 
 }

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.

                     {
 
     typedef typename det::Footprint::SpanList SpanList;
 
     PTR(det::Footprint) sfoot(new det::Footprint(foot.getPeaks().getSchema()));
     const SpanList spans = foot.getSpans();
     assert(foot.isNormalized());
 
     pexLog::Log log(pexLog::Log::getDefaultLog(),
                     "lsst.meas.deblender.symmetrizeFootprint");
 
     // Find the Span containing the peak.
     PTR(det::Span) target(new det::Span(cy, cx, cx));
     SpanList::const_iterator peakspan =
         std::upper_bound(spans.begin(), spans.end(), target, span_ptr_compare);
     // upper_bound returns the last position where "target" could be inserted;
     // ie, the first Span larger than "target".  The Span containing "target"
     // should be peakspan-1 or (if the peak is on the first pixel in the span,
     // peakspan.
     PTR(det::Span) sp;
     if (peakspan == spans.begin()) {
         sp = *peakspan;
         if (!sp->contains(cx, cy)) {
             log.warnf(
                 "Failed to find span containing (%i,%i): before the beginning of this footprint", cx, cy);
             return PTR(det::Footprint)();
         }
     } else {
         --peakspan;
         sp = *peakspan;
 
         if (!sp->contains(cx, cy)) {
             ++peakspan;
             sp = *peakspan;
             if (!sp->contains(cx, cy)) {
                 geom::Box2I fbb = foot.getBBox();
                 log.warnf("Failed to find span containing (%i,%i): nearest is %i, [%i,%i].  "
                           "Footprint bbox is [%i,%i],[%i,%i]",
                           cx, cy, sp->getY(), sp->getX0(), sp->getX1(),
                           fbb.getMinX(), fbb.getMaxX(), fbb.getMinY(), fbb.getMaxY());
                 return PTR(det::Footprint)();
             }
         }
     }
     log.debugf("Span containing (%i,%i): (x=[%i,%i], y=%i)",
                cx, cy, sp->getX0(), sp->getX1(), sp->getY());
 
     // The symmetric templates are essentially an AND of the footprint
     // pixels and its 180-degree-rotated self, rotated around the
     // peak (cx,cy).
     //
     // We iterate forward and backward simultaneously, starting from
     // the span containing the peak and moving out, row by row.
     // 
     // In the loop below, we search for the next pair of Spans that
     // overlap (in "dx" from the center), output the overlapping
     // portion of the Spans, and advance either the "fwd" or "back"
     // iterator.  When we fail to find an overlapping pair of Spans,
     // we move on to the next row.
     //
     // [The following paragraph is somewhat obsoleted by the
     // RelativeSpanIterator class, which performs some of the renaming
     // and the dx,dy coords.]
     //
     // '''In reading the code, "forward", "advancing", etc, are all
     // from the perspective of the "fwd" iterator (the one going
     // forward through the Span list, from low to high Y and then low
     // to high X).  It will help to imagine making a copy of the
     // footprint and rotating it around the center pixel by 180
     // degrees, so that "fwd" and "back" are both iterating the same
     // direction; we're then just finding the AND of those two
     // iterators, except we have to work in dx,dy coordinates rather
     // than original x,y coords, and the accessors for "back" are
     // opposite.'''
 
     RelativeSpanIterator fwd (peakspan, spans, cx, cy, true);
     RelativeSpanIterator back(peakspan, spans, cx, cy, false);
 
     int dy = 0;
     while (fwd.notDone() && back.notDone()) {
         // forward and backward "y"; just symmetric around cy
         int fy = cy + dy;
         int by = cy - dy;
         // delta-x of the beginnings of the spans, for "fwd" and "back"
         int fdxlo =  fwd.dxlo();
         int bdxlo = back.dxlo();
 
         // First find:
         //    fend -- first span in the next row, or end(); ie,
         //            the end of this row in the forward direction
         //    bend -- the end of this row in the backward direction
         RelativeSpanIterator fend, bend;
         for (fend = fwd; fend.notDone(); ++fend) {
             if (fend.dy() != dy)
                 break;
         }
         for (bend = back; bend.notDone(); ++bend) {
             if (bend.dy() != dy)
                 break;
         }
 
         log.debugf("dy=%i, fy=%i, fx=[%i, %i],   by=%i, fx=[%i, %i],  fdx=%i, bdx=%i",
                    dy, fy, fwd.x0(), fwd.x1(), by, back.x0(), back.x1(),
                    fdxlo, bdxlo);
 
         // Find possibly-overlapping span
         if (bdxlo > fdxlo) {
             log.debugf("Advancing forward.");
             // While the "forward" span is entirely to the "left" of the "backward" span,
             // (in dx coords), ie, |---fwd---X   X---back---|
             // and we are comparing the edges marked X
             while ((fwd != fend) && (fwd.dxhi() < bdxlo)) {
                 fwd++;
                 if (fwd == fend) {
                     log.debugf("Reached fend");
                 } else {
                     log.debugf("Advanced to forward span %i, [%i, %i]",
                                fy, fwd.x0(), fwd.x1());
                 }
             }
         } else if (fdxlo > bdxlo) {
             log.debugf("Advancing backward.");
             // While the "backward" span is entirely to the "left" of the "foreward" span,
             // (in dx coords), ie, |---back---X   X---fwd---|
             // and we are comparing the edges marked X
             while ((back != bend) && (back.dxhi() < fdxlo)) {
                 back++;
                 if (back == bend) {
                     log.debugf("Reached bend");
                 } else {
                     log.debugf("Advanced to backward span %i, [%i, %i]",
                                by, back.x0(), back.x1());
                 }
             }
         }
 
         if ((back == bend) || (fwd == fend)) {
             // We reached the end of the row without finding spans that could
             // overlap.  Move onto the next dy.
             if (back == bend) {
                 log.debugf("Reached bend");
             }
             if (fwd == fend) {
                 log.debugf("Reached fend");
             }
             back = bend;
             fwd  = fend;
             dy++;
             continue;
         }
 
         // Spans may overlap -- find the overlapping part.
         int dxlo = std::max(fwd.dxlo(), back.dxlo());
         int dxhi = std::min(fwd.dxhi(), back.dxhi());
         if (dxlo <= dxhi) {
             log.debugf("Adding span fwd %i, [%i, %i],  back %i, [%i, %i]",
                        fy, cx+dxlo, cx+dxhi, by, cx-dxhi, cx-dxlo);
             sfoot->addSpan(fy, cx + dxlo, cx + dxhi);
             sfoot->addSpan(by, cx - dxhi, cx - dxlo);
         }
 
         // Advance the one whose "hi" edge is smallest
         if (fwd.dxhi() < back.dxhi()) {
             fwd++;
             if (fwd == fend) {
                 log.debugf("Stepped to fend");
             } else {
                 log.debugf("Stepped forward to span %i, [%i, %i]",
                            fwd.y(), fwd.x0(), fwd.x1());
             }
         } else {
             back++;
             if (back == bend) {
                 log.debugf("Stepped to bend");
             } else {
                 log.debugf("Stepped backward to span %i, [%i, %i]",
                            back.y(), back.x0(), back.x1());
             }
         }
 
         if ((back == bend) || (fwd == fend)) {
             // Reached the end of the row.  On to the next dy!
             if (back == bend) {
                 log.debugf("Reached bend");
             }
             if (fwd == fend) {
                 log.debugf("Reached fend");
             }
             back = bend;
             fwd  = fend;
             dy++;
             continue;
         }
 
     }
     sfoot->normalize();
     return sfoot;
 }

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:

/home/lsstsw/stack/Linux64/meas_deblender/11.0.rc2+12/include/lsst/meas/deblender/Baseline.h
/home/lsstsw/stack/Linux64/meas_deblender/11.0.rc2+12/src/Baseline.cc

Public Types

Static Public Member Functions

Static Public Attributes

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 >

Member Typedef Documentation

Member Function Documentation

Member Data Documentation

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