19 template <
typename ImagePixelT,
typename MaskPixelT,
typename VariancePixelT>
22 template <
typename ImagePixelT,
typename MaskPixelT,
typename VariancePixelT>
26 template <
typename ImagePixelT,
typename MaskPixelT,
typename VariancePixelT>
29 template <
typename ImagePixelT,
typename MaskPixelT,
typename VariancePixelT>
32 template <
typename ImagePixelT,
typename MaskPixelT,
typename VariancePixelT>
35 template <
typename ImagePixelT,
typename MaskPixelT,
typename VariancePixelT>
55 const itype nil = 0xffff;
63 foot->getSpans()->setImage(*argmin,
static_cast<itype
>(i));
64 foot->getSpans()->setImage(*dist,
static_cast<dtype
>(1));
69 int const height = dist->getHeight();
70 int const width = dist->getWidth();
73 for (
int y = 0;
y != height; ++
y) {
76 for (
int x = 0;
x != width; ++
x, ++dim.x(), ++aim.x()) {
84 dim(0, 0) = width + height;
88 dtype ndist = dim(0,-1) + 1;
89 if (ndist < dim(0,0)) {
96 dtype ndist = dim(-1,0) + 1;
97 if (ndist < dim(0,0)) {
106 for (
int y = height - 1;
y >= 0; --
y) {
109 for (
int x = width - 1;
x >= 0; --
x, --dim.x(), --aim.x()) {
112 if (
y + 1 < height) {
113 dtype ndist = dim(0,1) + 1;
114 if (ndist < dim(0,0)) {
121 dtype ndist = dim(1,0) + 1;
122 if (ndist < dim(0,0)) {
145 template<
typename ImagePixelT,
typename MaskPixelT,
typename VariancePixelT>
151 int S = halfsize*2 + 1;
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));
163 ImagePixelT vals[S*S];
164 for (
int y=halfsize;
y<H-halfsize; ++
y) {
167 inpix !=
end; ++inpix.x(), ++optr) {
168 for (
int i=0; i<SS; ++i)
169 vals[i] = inpix[locs[i]];
176 for (
int y=0;
y<2*halfsize; ++
y) {
179 iy = H - 1 - (
y-halfsize);
182 for (; iptr !=
end; ++iptr,++optr)
185 for (
int y=halfsize;
y<H-halfsize; ++
y) {
188 for (; iptr !=
end; ++iptr,++optr)
193 for (; iptr !=
end; ++iptr,++optr)
223 template<
typename ImagePixelT,
typename MaskPixelT,
typename VariancePixelT>
230 int cx =
peak.getIx();
231 int cy =
peak.getIy();
232 int ix0 = img.
getX0();
233 int iy0 = img.
getY0();
246 for (s = 0; s <
std::max(DW,DH); s += S) {
248 for (p=0; p<S; p++) {
288 for (
int i=0; i<(8*L); i++,
x += dx,
y += dy) {
292 if (i % (2*L) == 0) {
295 dx = ( leg % 2) * (-1 + 2*(leg/2));
297 dy = ((leg+1) % 2) * ( 1 - 2*(leg/2));
300 int px = cx +
x - ix0;
301 int py = cy +
y - iy0;
303 if (px < 0 || px >= iW || py < 0 || py >= iH)
306 ImagePixelT pix = (*shadowingImg)(px,
py);
313 const double A = 0.3;
319 ds0 = (double(
y) / double(
x)) - A;
322 for (shx=1; shx<=S; shx++) {
323 int xsign = (
x>0?1:-1);
325 psx = cx +
x + (xsign*shx) - ix0;
326 if (psx < 0 || psx >= iW)
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)
334 img(psx, psy) =
std::min(img(psx, psy), pix);
340 ds0 = (double(
x) / double(
y)) - A;
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)
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)
354 img(psx, psy) =
std::min(img(psx, psy), pix);
360 shadowingImg->assign(img);
364 static double _get_contrib_r_to_footprint(
int x,
int y,
370 int dx = sp.getX0() -
x;
383 int dy = sp.getY() -
y;
384 minr2 =
std::min(minr2, (
double)(mindx*mindx + dy*dy));
387 return 1. / (1. + minr2);
391 template<
typename ImagePixelT,
typename MaskPixelT,
typename VariancePixelT>
396 MaskedImageT
const& img,
397 int strayFluxOptions,
402 double clipStrayFluxFraction,
417 int ix0 = img.getX0();
418 int iy0 = img.getY0();
423 for (
size_t i=0; i<tfoots.size(); ++i) {
430 bool always = (strayFluxOptions & STRAYFLUX_TO_POINT_SOURCES_ALWAYS);
435 if (strayFluxOptions & STRAYFLUX_NEAREST_FOOTPRINT) {
445 if (!always && ispsf.
size()) {
448 auto empty = std::make_shared<det::Footprint>();
450 for (
size_t i=0; i<tfoots.size(); ++i) {
457 footlist = &templist;
459 nearestFootprint(*footlist, nearest, dist);
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()];
474 for (
int x = x0;
x <= x1; ++
x, ++tsum_it, ++in_it) {
478 if ((*tsum_it > 0) || (*in_it).image() <= 0) {
482 if (strayFluxOptions & STRAYFLUX_R_TO_FOOTPRINT) {
484 for (
size_t i=0; i<tfoots.size(); ++i) {
487 }
else if (strayFluxOptions & STRAYFLUX_NEAREST_FOOTPRINT) {
488 for (
size_t i=0; i<tfoots.size(); ++i) {
491 int i = nearest->get0(
x,
y);
495 for (
size_t i=0; i<tfoots.size(); ++i) {
500 contrib[i] = 1. / (1. + dx*dx + dy*dy);
506 bool ptsrcs = always;
508 for (
size_t i=0; i<tfoots.size(); ++i) {
510 if ((!ptsrcs) && ispsf.
size() && ispsf[i]) {
513 if (contrib[i] == -1.0) {
514 contrib[i] = _get_contrib_r_to_footprint(
x,
y, tfoots[i]);
520 STRAYFLUX_TO_POINT_SOURCES_WHEN_NECESSARY)) {
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]);
533 double strayclip = (clipStrayFluxFraction * csum);
535 for (
size_t i=0; i<tfoots.size(); ++i) {
537 if ((!ptsrcs) && ispsf.
size() && ispsf[i]) {
542 if (contrib[i] < strayclip) {
549 for (
size_t i=0; i<tfoots.size(); ++i) {
550 if (contrib[i] == 0.) {
554 double p = (contrib[i] / csum) * (*in_it).image();
557 strayfoot[i] = std::make_shared<det::Footprint>();
563 strayvar[i].
push_back((*in_it).variance());
569 for (
size_t i=0; i<tfoots.size(); ++i) {
571 strayfoot[i]->setSpans(std::make_shared<afwGeom::SpanSet>(straySpans[i]));
579 HeavyFootprintPtrT heavy(
new HeavyFootprint(*strayfoot[i]));
580 ndarray::Array<ImagePixelT,1,1> himg = heavy->getImageArray();
588 assert((
size_t)strayfoot[i]->getArea() == straypix[i].size());
590 for (spix = straypix[i].
begin(),
591 smask = straymask[i].
begin(),
592 svar = strayvar [i].
begin(),
594 mpix = heavy->getMaskArray().begin(),
595 vpix = heavy->getVarianceArray().begin();
596 spix != straypix[i].end();
597 ++spix, ++smask, ++svar, ++hpix, ++mpix, ++vpix) {
607 template<
typename ImagePixelT,
typename MaskPixelT,
typename VariancePixelT>
617 for (
size_t i=0; i<timgs.
size(); ++i) {
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));
690 template<
typename ImagePixelT,
typename MaskPixelT,
typename VariancePixelT>
702 int strayFluxOptions,
703 double clipStrayFluxFraction
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());
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");
717 if (!img.getBBox().contains(foot.
getBBox())) {
719 "Image bbox MUST contain parent footprint");
730 bool findStrayFlux = (strayFluxOptions & ASSIGN_STRAYFLUX);
732 int ix0 = img.getX0();
733 int iy0 = img.getY0();
741 if (!tsum->getBBox().contains(foot.
getBBox())) {
743 "Template sum image MUST contain parent footprint");
750 _sum_templates(timgs, tsum);
753 for (
size_t i=0; i<timgs.
size(); ++i) {
754 ImagePtrT timg = timgs[i];
756 MaskedImagePtrT port(
new MaskedImageT(timg->getDimensions()));
757 port->setXY0(timg->getXY0());
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) {
781 double frac =
std::max((ImagePixelT)0.,
static_cast<ImagePixelT
>(*tptr)) / (*tsum_it);
785 out_it.mask() = (*in_it).mask();
786 out_it.variance() = (*in_it).variance();
787 out_it.image() = (*in_it).image() * frac;
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());
800 (
boost::format(
"'pkx' and 'pky' must be the same length as templates (%d,%d vs %d)")
803 _find_stray_flux(foot, tsum, img, strayFluxOptions, tfoots,
804 ispsf, pkx, pky, clipStrayFluxFraction, strays);
826 int cx,
int cy,
bool forward=
true)
827 : _real(
real), _cx(cx), _cy(cy), _forward(forward)
837 return _real ==
other;
840 return _real !=
other;
843 return _real <=
other;
846 return _real <
other;
849 return _real >=
other;
852 return _real >
other;
856 return (_real ==
other._real) &&
858 (_forward ==
other._forward);
861 return !(*
this ==
other);
882 return _real >= _end;
888 return _real->getX0() - _cx;
890 return _cx - _real->getX1();
895 return _real->getX1() - _cx;
897 return _cx - _real->getX0();
901 return std::abs(_real->getY() - _cy);
904 return _real->getX0();
907 return _real->getX1();
910 return _real->getY();
972 template<
typename ImagePixelT,
typename MaskPixelT,
typename VariancePixelT>
979 auto sfoot = std::make_shared<det::Footprint>();
980 sfoot->setPeakSchema(foot.getPeaks().getSchema());
994 if (peakspan == spans.
begin()) {
998 "Failed to find span containing (%i,%i): before the beginning of this footprint", cx, cy);
1010 LOGL_WARN(_log,
"Failed to find span containing (%i,%i): nearest is %i, [%i,%i]. "
1011 "Footprint bbox is [%i,%i],[%i,%i]",
1018 LOGL_DEBUG(_log,
"Span containing (%i,%i): (x=[%i,%i], y=%i)",
1059 int fdxlo = fwd.
dxlo();
1060 int bdxlo = back.
dxlo();
1067 for (fend = fwd; fend.
notDone(); ++fend) {
1068 if (fend.
dy() != dy)
1071 for (bend = back; bend.
notDone(); ++bend) {
1072 if (bend.
dy() != dy)
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(),
1081 if (bdxlo > fdxlo) {
1086 while ((fwd != fend) && (fwd.
dxhi() < bdxlo)) {
1091 LOGL_DEBUG(_log,
"Advanced to forward span %i, [%i, %i]",
1092 fy, fwd.
x0(), fwd.
x1());
1095 }
else if (fdxlo > bdxlo) {
1100 while ((back != bend) && (back.
dxhi() < fdxlo)) {
1105 LOGL_DEBUG(_log,
"Advanced to backward span %i, [%i, %i]",
1106 by, back.
x0(), back.
x1());
1111 if ((back == bend) || (fwd == fend)) {
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);
1142 LOGL_DEBUG(_log,
"Stepped forward to span %i, [%i, %i]",
1143 fwd.
y(), fwd.
x0(), fwd.
x1());
1150 LOGL_DEBUG(_log,
"Stepped backward to span %i, [%i, %i]",
1151 back.
y(), back.
x0(), back.
x1());
1155 if ((back == bend) || (fwd == fend)) {
1170 sfoot->setSpans(std::make_shared<afwGeom::SpanSet>(
std::move(tmpSpans)));
1186 template<
typename ImagePixelT,
typename MaskPixelT,
typename VariancePixelT>
1197 bool* patchedEdges) {
1201 *patchedEdges =
false;
1203 int cx =
peak.getIx();
1204 int cy =
peak.getIy();
1220 "Image too small for symmetrized footprint");
1225 bool touchesEdge =
false;
1227 LOGL_DEBUG(_log,
"Checking footprint for EDGE bits");
1230 MaskPixelT edgebit =
mask->getPlaneBitMask(
"EDGE");
1232 fwd != spans.
end(); ++fwd) {
1233 int x0 = fwd->getX0();
1234 int x1 = fwd->getX1();
1236 mask->x_at(x0 -
mask->getX0(), fwd->getY() -
mask->getY0());
1237 for (
int x=x0;
x<=x1; ++
x, ++xiter) {
1238 if ((*xiter) & edgebit) {
1247 LOGL_DEBUG(_log,
"Footprint includes an EDGE pixel.");
1260 for (; fwd <= back; fwd++, back--) {
1261 int fy = fwd->getY();
1262 int by = back->getY();
1264 for (
int fx=fwd->getX0(), bx=back->getX1();
1279 ImagePixelT pixf = theimg->get0(fx, fy);
1280 ImagePixelT pixb = theimg->get0(bx, by);
1281 ImagePixelT pix =
std::min(pixf, pixb);
1283 pix =
std::max(pix,
static_cast<ImagePixelT
>(0));
1285 targetimg->set0(fx, fy, pix);
1286 targetimg->set0(bx, by, pix);
1302 LOGL_DEBUG(_log,
"Footprint touches EDGE: start bbox [%i,%i],[%i,%i]",
1306 for (fwd = ospans.
begin(); fwd != ospans.
end(); ++fwd) {
1307 int y = fwd->getY();
1308 int x = fwd->getX0();
1310 int ym = cy + (cy -
y);
1311 int xm = cx + (cx -
x);
1321 LOGL_DEBUG(_log,
"Footprint touches EDGE: grown bbox [%i,%i],[%i,%i]",
1326 sfoot->getSpans()->copyImage(*targetimg, *targetimg2);
1328 LOGL_DEBUG(_log,
"Symmetric footprint spans:");
1330 for (fwd = sspans.
begin(); fwd != sspans.
end(); ++fwd) {
1331 LOGL_DEBUG(_log,
" %s", fwd->toString().c_str());
1337 for (fwd = ospans.
begin(); fwd != ospans.
end(); ++fwd) {
1338 int y = fwd->getY();
1339 int x0 = fwd->getX0();
1340 int x1 = fwd->getX1();
1342 int ym = cy + (cy -
y);
1343 int xm0 = cx + (cx - x0);
1344 int xm1 = cx + (cx - x1);
1355 x0 = cx + (cx - (imbb.
getMinX() - 1));
1358 x1 = cx + (cx - (imbb.
getMaxX() + 1));
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);
1365 targetimg2->x_at(x0 - targetimg2->getX0(),
y - targetimg2->getY0());
1366 for (
int x=x0;
x<=x1; ++
x, ++outiter, ++initer) {
1367 *outiter = initer.image();
1371 sfoot->setSpans(std::make_shared<afwGeom::SpanSet>(
std::move(newSpans)));
1372 targetimg = targetimg2;
1375 *patchedEdges = touchesEdge;
1383 template<
typename ImagePixelT,
typename MaskPixelT,
typename VariancePixelT>
1388 ImagePixelT thresh) {
1398 int const y = sp->getY();
1399 int const x0 = sp->getX0();
1400 int const x1 = sp->getX1();
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) {
1416 template<
typename ImagePixelT,
typename MaskPixelT,
typename VariancePixelT>
1421 ImagePixelT thresh) {
1425 auto significant = std::make_shared<det::Footprint>();
1428 int const x0 = img->getX0(), y0 = img->getY0();
1433 int const y = span.
getY();
1435 typename ImageT::const_x_iterator
iter = img->x_at(
x - x0,
y - y0);
1436 bool onSpan =
false;
1439 if (*
iter >= thresh) {
1442 }
else if (onSpan) {
1451 significant->setSpans(std::make_shared<afwGeom::SpanSet>(
std::move(tmpSpans)));