LSST Applications g0b6bd0c080+a72a5dd7e6,g1182afd7b4+2a019aa3bb,g17e5ecfddb+2b8207f7de,g1d67935e3f+06cf436103,g38293774b4+ac198e9f13,g396055baef+6a2097e274,g3b44f30a73+6611e0205b,g480783c3b1+98f8679e14,g48ccf36440+89c08d0516,g4b93dc025c+98f8679e14,g5c4744a4d9+a302e8c7f0,g613e996a0d+e1c447f2e0,g6c8d09e9e7+25247a063c,g7271f0639c+98f8679e14,g7a9cd813b8+124095ede6,g9d27549199+a302e8c7f0,ga1cf026fa3+ac198e9f13,ga32aa97882+7403ac30ac,ga786bb30fb+7a139211af,gaa63f70f4e+9994eb9896,gabf319e997+ade567573c,gba47b54d5d+94dc90c3ea,gbec6a3398f+06cf436103,gc6308e37c7+07dd123edb,gc655b1545f+ade567573c,gcc9029db3c+ab229f5caf,gd01420fc67+06cf436103,gd877ba84e5+06cf436103,gdb4cecd868+6f279b5b48,ge2d134c3d5+cc4dbb2e3f,ge448b5faa6+86d1ceac1d,gecc7e12556+98f8679e14,gf3ee170dca+25247a063c,gf4ac96e456+ade567573c,gf9f5ea5b4d+ac198e9f13,gff490e6085+8c2580be5c,w.2022.27
LSST Data Management Base Package
CR.cc
Go to the documentation of this file.
1// -*- LSST-C++ -*-
2
3/*
4 * LSST Data Management System
5 * Copyright 2008-2016 AURA/LSST.
6 *
7 * This product includes software developed by the
8 * LSST Project (http://www.lsst.org/).
9 *
10 * This program is free software: you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation, either version 3 of the License, or
13 * (at your option) any later version.
14 *
15 * This program is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the LSST License Statement and
21 * the GNU General Public License along with this program. If not,
22 * see <https://www.lsstcorp.org/LegalNotices/>.
23 */
24
32#include <stdexcept>
33#include <algorithm>
34#include <cassert>
35#include <string>
36#include <typeinfo>
37
38#include <iostream>
39
40#include "boost/format.hpp"
41
42#include "lsst/pex/exceptions.h"
43#include "lsst/log/Log.h"
44#include "lsst/pex/exceptions.h"
46#include "lsst/afw/geom.h"
52
57namespace lsst {
58namespace afw {
59namespace detection {
63class IdSpan {
64public:
67
68 explicit IdSpan(int id, int y, int x0, int x1) : id(id), y(y), x0(x0), x1(x1) {}
69 int id; /* ID for object */
70 int y; /* Row wherein IdSpan dwells */
71 int x0, x1; /* inclusive range of columns */
72};
76 struct IdSpanCompar {
78 if (a->id < b->id) {
79 return true;
80 } else if (a->id > b->id) {
81 return false;
82 } else {
83 if (a->y < b->y) {
84 return true;
85 } else if (a->y > b->y) {
86 return false;
87 } else {
88 return (a->x0 < b->x0) ? true : false;
89 }
90 }
91 }
92};
96int resolve_alias(const std::vector<int> &aliases, /* list of aliases */
97 int id) { /* alias to look up */
98 int resolved = id; /* resolved alias */
99
100 while (id != aliases[id]) {
101 resolved = id = aliases[id];
102 }
103
104 return (resolved);
105}
106} // namespace detection
107} // namespace afw
108} // namespace lsst
109
110namespace lsst {
111namespace meas {
112namespace algorithms {
113
114namespace {
115
116template <typename ImageT, typename MaskT>
117void removeCR(afw::image::MaskedImage<ImageT, MaskT> &mi,
118 std::vector<std::shared_ptr<afw::detection::Footprint>> &CRs, double const bkgd, MaskT const,
119 MaskT const saturBit, MaskT const badMask, bool const debias, bool const grow);
120
121template <typename ImageT>
122bool condition_3(ImageT *estimate, double const peak, double const mean_ns, double const mean_we,
123 double const mean_swne, double const mean_nwse, double const dpeak, double const dmean_ns,
124 double const dmean_we, double const dmean_swne, double const dmean_nwse, double const thresH,
125 double const thresV, double const thresD, double const cond3Fac);
126
127/************************************************************************************************************/
128//
129// A class to hold a detected pixel
130template <typename ImageT>
131struct CRPixel {
132 typedef typename std::shared_ptr<CRPixel> Ptr;
133
134 CRPixel(int _col, int _row, ImageT _val, int _id = -1) : id(_id), col(_col), row(_row), val(_val) {
135 _i = ++i;
136 }
137 ~CRPixel() {}
138
139 bool operator<(const CRPixel &a) const { return _i < a._i; }
140
141 int get_i() const { return _i; }
142
143 int id; // Unique ID for cosmic ray (not cosmic ray pixel)
144 int col; // position
145 int row; // of pixel
146 ImageT val; // initial value of pixel
147private:
148 static int i; // current value of running index
149 int mutable _i; // running index
150};
151
152template <typename ImageT>
153int CRPixel<ImageT>::i = 0;
154
155template <typename ImageT>
156struct Sort_CRPixel_by_id {
157 bool operator()(CRPixel<ImageT> const &a, CRPixel<ImageT> const &b) const { return a.id < b.id; }
158};
159
160/*****************************************************************************/
161/*
162 * This is the code to see if a given pixel is bad
163 *
164 * Note that the pixel in question is at index 0, so its value is pt_0[0]
165 */
166template <typename MaskedImageT>
167bool is_cr_pixel(typename MaskedImageT::Image::Pixel *corr, // corrected value
168 typename MaskedImageT::xy_locator loc, // locator for this pixel
169 double const minSigma, // minSigma, or -threshold if negative
170 double const thresH, double const thresV, double const thresD, // for condition #3
171 double const bkgd, // unsubtracted background level
172 double const cond3Fac // fiddle factor for condition #3
173) {
174 typedef typename MaskedImageT::Image::Pixel ImagePixel;
175 //
176 // Unpack some values
177 //
178 ImagePixel const v_00 = loc.image(0, 0);
179
180 if (v_00 < 0) {
181 return false;
182 }
183 /*
184 * condition #1 is not applied on a pixel-by-pixel basis
185 */
186 ;
187 /*
188 * condition #2
189 */
190 ImagePixel const mean_we = (loc.image(-1, 0) + loc.image(1, 0)) / 2; // avgs of surrounding 8 pixels
191 ImagePixel const mean_ns = (loc.image(0, 1) + loc.image(0, -1)) / 2;
192 ImagePixel const mean_swne = (loc.image(-1, -1) + loc.image(1, 1)) / 2;
193 ImagePixel const mean_nwse = (loc.image(-1, 1) + loc.image(1, -1)) / 2;
194
195 if (minSigma < 0) { /* |thres_sky_sigma| is threshold */
196 if (v_00 < -minSigma) {
197 return false;
198 }
199 } else {
200 double const thres_sky_sigma = minSigma * sqrt(loc.variance(0, 0));
201
202 if (v_00 < mean_ns + thres_sky_sigma && v_00 < mean_we + thres_sky_sigma &&
203 v_00 < mean_swne + thres_sky_sigma && v_00 < mean_nwse + thres_sky_sigma) {
204 return false;
205 }
206 }
207 /*
208 * condition #3
209 *
210 * Note that this uses mean_ns etc. even if minSigma is negative
211 */
212 double const dv_00 = sqrt(loc.variance(0, 0));
213 // standard deviation of means of surrounding pixels
214 double const dmean_we = sqrt(loc.variance(-1, 0) + loc.variance(1, 0)) / 2;
215 double const dmean_ns = sqrt(loc.variance(0, 1) + loc.variance(0, -1)) / 2;
216 double const dmean_swne = sqrt(loc.variance(-1, -1) + loc.variance(1, 1)) / 2;
217 double const dmean_nwse = sqrt(loc.variance(-1, 1) + loc.variance(1, -1)) / 2;
218
219 if (!condition_3(corr, v_00 - bkgd, mean_ns - bkgd, mean_we - bkgd, mean_swne - bkgd, mean_nwse - bkgd,
220 dv_00, dmean_ns, dmean_we, dmean_swne, dmean_nwse, thresH, thresV, thresD, cond3Fac)) {
221 return false;
222 }
223 /*
224 * OK, it's a contaminated pixel
225 */
226 *corr += static_cast<ImagePixel>(bkgd);
227
228 return true;
229}
230
231/************************************************************************************************************/
232//
233// Worker routine to process the pixels adjacent to a span (including the points just
234// to the left and just to the right)
235//
236template <typename MaskedImageT>
237void checkSpanForCRs(afw::detection::Footprint *extras, // Extra spans get added to this Footprint
238 std::vector<CRPixel<typename MaskedImageT::Image::Pixel>> &crpixels,
239 // a list of pixels containing CRs
240 int const y, // the row to process
241 int const x0, int const x1, // range of pixels in the span (inclusive)
242 MaskedImageT &image,
243 double const minSigma, // minSigma
244 double const thresH, double const thresV, double const thresD, // for cond. #3
245 double const bkgd, // unsubtracted background level
246 double const cond3Fac, // fiddle factor for condition #3
247 bool const keep // if true, don't remove the CRs
248) {
249 typedef typename MaskedImageT::Image::Pixel MImagePixel;
250 typename MaskedImageT::xy_locator loc = image.xy_at(x0 - 1, y); // locator for data
251
252 int const imageX0 = image.getX0();
253 int const imageY0 = image.getY0();
254
255 for (int x = x0 - 1; x <= x1 + 1; ++x) {
256 MImagePixel corr = 0; // new value for pixel
257 if (is_cr_pixel<MaskedImageT>(&corr, loc, minSigma, thresH, thresV, thresD, bkgd, cond3Fac)) {
258 if (keep) {
259 crpixels.push_back(CRPixel<MImagePixel>(x + imageX0, y + imageY0, loc.image()));
260 }
261 loc.image() = corr;
262
263 std::vector<afw::geom::Span> spanList(extras->getSpans()->begin(), extras->getSpans()->end());
264 spanList.push_back(afw::geom::Span(y + imageY0, x + imageX0, x + imageX0));
265 extras->setSpans(std::make_shared<afw::geom::SpanSet>(std::move(spanList)));
266 }
267 ++loc.x();
268 }
269}
270
271/************************************************************************************************************/
272/*
273 * Find the sum of the pixels in a Footprint
274 */
275template <typename ImageT>
276class CountsInCR {
277public:
278 CountsInCR(double const bkgd) : _bkgd(bkgd), _sum(0.0) {}
279
280 void operator()(geom::Point2I const &point, ImageT const &val) { _sum += val - _bkgd; }
281
282 virtual void reset(afw::detection::Footprint const &) {}
283 virtual void reset() { _sum = 0.0; }
284
285 double getCounts() const { return _sum; }
286
287private:
288 double const _bkgd; // the Image's background level
289 ImageT _sum; // the sum of all DN in the Footprint, corrected for bkgd
290};
291} // namespace
292
293/*
294 * Restore all the pixels in crpixels to their pristine state
295 */
296template <typename ImageT>
297static void reinstateCrPixels(
298 ImageT *image, // the image in question
299 std::vector<CRPixel<typename ImageT::Pixel>> const &crpixels // a list of pixels with CRs
300) {
301 if (crpixels.empty()) return;
302
303 typedef typename std::vector<CRPixel<typename ImageT::Pixel>>::const_iterator crpixel_iter;
304 for (crpixel_iter crp = crpixels.begin(), end = crpixels.end(); crp < end - 1; ++crp) {
305 *image->at(crp->col - image->getX0(), crp->row - image->getY0()) = crp->val;
306 }
307}
308
314template <typename MaskedImageT>
316 MaskedImageT &mimage,
317 afw::detection::Psf const &psf,
318 double const bkgd,
319 daf::base::PropertySet const &ps,
320 bool const keep
321) {
322 typedef typename MaskedImageT::Image ImageT;
323 typedef typename ImageT::Pixel ImagePixel;
324 typedef typename MaskedImageT::Mask::Pixel MaskPixel;
325
326 // Parse the PropertySet
327 double const minSigma = ps.getAsDouble("minSigma"); // min sigma over sky in pixel for CR candidate
328 double const minDn = ps.getAsDouble("min_DN"); // min number of DN in an CRs
329 double const cond3Fac = ps.getAsDouble("cond3_fac"); // fiddle factor for condition #3
330 double const cond3Fac2 = ps.getAsDouble("cond3_fac2"); // 2nd fiddle factor for condition #3
331 int const niteration = ps.getAsInt("niteration"); // Number of times to look for contaminated
332 // pixels near CRs
333 int const nCrPixelMax = ps.getAsInt("nCrPixelMax"); // maximum number of contaminated pixels
334 /*
335 * thresholds for 3rd condition
336 *
337 * Realise PSF at center of image
338 */
339 std::shared_ptr<afw::math::Kernel const> kernel = psf.getLocalKernel();
340 if (!kernel) {
341 throw LSST_EXCEPT(pexExcept::NotFoundError, "Psf is unable to return a kernel");
342 }
344 afw::detection::Psf::Image(geom::ExtentI(kernel->getWidth(), kernel->getHeight()));
345 kernel->computeImage(psfImage, true);
346
347 int const xc = kernel->getCtr().getX(); // center of PSF
348 int const yc = kernel->getCtr().getY();
349
350 double const I0 = psfImage(xc, yc);
351 double const thresH =
352 cond3Fac2 * (0.5 * (psfImage(xc - 1, yc) + psfImage(xc + 1, yc))) / I0; // horizontal
353 double const thresV = cond3Fac2 * (0.5 * (psfImage(xc, yc - 1) + psfImage(xc, yc + 1))) / I0; // vertical
354 double const thresD = cond3Fac2 *
355 (0.25 * (psfImage(xc - 1, yc - 1) + psfImage(xc + 1, yc + 1) +
356 psfImage(xc - 1, yc + 1) + psfImage(xc + 1, yc - 1))) /
357 I0; // diag
358 /*
359 * Setup desired mask planes
360 */
361 MaskPixel const badBit = mimage.getMask()->getPlaneBitMask("BAD"); // Generic bad pixels
362 MaskPixel const crBit = mimage.getMask()->getPlaneBitMask("CR"); // CR-contaminated pixels
363 MaskPixel const interpBit = mimage.getMask()->getPlaneBitMask("INTRP"); // Interpolated pixels
364 MaskPixel const saturBit = mimage.getMask()->getPlaneBitMask("SAT"); // Saturated pixels
365 MaskPixel const nodataBit = mimage.getMask()->getPlaneBitMask("NO_DATA"); // Non data pixels
366
367 MaskPixel const badMask = (badBit | interpBit | saturBit | nodataBit); // naughty pixels
368 /*
369 * Go through the frame looking at each pixel (except the edge ones which we ignore)
370 */
371 int const ncol = mimage.getWidth();
372 int const nrow = mimage.getHeight();
373
374 std::vector<CRPixel<ImagePixel>> crpixels; // storage for detected CR-contaminated pixels
375 typedef typename std::vector<CRPixel<ImagePixel>>::iterator crpixel_iter;
376 typedef typename std::vector<CRPixel<ImagePixel>>::reverse_iterator crpixel_riter;
377
378 for (int j = 1; j < nrow - 1; ++j) {
379 typename MaskedImageT::xy_locator loc = mimage.xy_at(1, j); // locator for data
380
381 for (int i = 1; i < ncol - 1; ++i, ++loc.x()) {
382 ImagePixel corr = 0;
383 if (!is_cr_pixel<MaskedImageT>(&corr, loc, minSigma, thresH, thresV, thresD, bkgd, cond3Fac)) {
384 continue;
385 }
386 /*
387 * condition #4
388 */
389 if (loc.mask() & badMask) {
390 continue;
391 }
392 if ((loc.mask(-1, 1) | loc.mask(0, 1) | loc.mask(1, 1) | loc.mask(-1, 0) | loc.mask(1, 0) |
393 loc.mask(-1, -1) | loc.mask(0, -1) | loc.mask(1, -1)) &
394 interpBit) {
395 continue;
396 }
397 /*
398 * OK, it's a CR
399 *
400 * replace CR-contaminated pixels with reasonable values as we go through
401 * image, which increases the detection rate
402 */
403 crpixels.push_back(CRPixel<ImagePixel>(i + mimage.getX0(), j + mimage.getY0(), loc.image()));
404 loc.image() = corr; /* just a preliminary estimate */
405
406 if (static_cast<int>(crpixels.size()) > nCrPixelMax) {
407 reinstateCrPixels(mimage.getImage().get(), crpixels);
408
410 (boost::format("Too many CR pixels (max %d)") % nCrPixelMax).str());
411 }
412 }
413 }
414 /*
415 * We've found them on a pixel-by-pixel basis, now merge those pixels
416 * into cosmic rays
417 */
418 std::vector<int> aliases; // aliases for initially disjoint parts of CRs
419 aliases.reserve(1 + crpixels.size() / 2); // initial size of aliases
420
421 std::vector<afw::detection::IdSpan::Ptr> spans; // y:x0,x1 for objects
422 spans.reserve(aliases.capacity()); // initial size of spans
423
424 aliases.push_back(0); // 0 --> 0
425
431 int ncr = 0; // number of detected cosmic rays
432 if (!crpixels.empty()) {
433 int id; // id number for a CR
434 int x0 = -1, x1 = -1, y = -1; // the beginning and end column, and row of this span in a CR
435
436 // I am dummy
437 CRPixel<ImagePixel> dummy(0, -1, 0, -1);
438 crpixels.push_back(dummy);
439 // printf("Created dummy CR: i %i, id %i, col %i, row %i, val %g\n", dummy.get_i(), dummy.id,
440 // dummy.col, dummy.row, (double)dummy.val);
441 for (crpixel_iter crp = crpixels.begin(); crp < crpixels.end() - 1; ++crp) {
442 // printf("Looking at CR: i %i, id %i, col %i, row %i, val %g\n", crp->get_i(), crp->id, crp->col,
443 // crp->row, (double)crp->val);
444
445 if (crp->id < 0) { // not already assigned
446 crp->id = ++ncr; // a new CR
447 aliases.push_back(crp->id);
448 y = crp->row;
449 x0 = x1 = crp->col;
450 // printf(" Assigned ID %i; looking at row %i, start col %i\n", crp->id, crp->row, crp->col);
451 }
452 id = crp->id;
453 // printf(" Next CRpix has i=%i, id=%i, row %i, col %i\n", crp[1].get_i(), crp[1].id, crp[1].row,
454 // crp[1].col);
455
456 if (crp[1].row == crp[0].row && crp[1].col == crp[0].col + 1) {
457 // printf(" Adjoining! Set next CRpix id = %i; x1=%i\n", crp[1].id, x1);
458 crp[1].id = id;
459 ++x1;
460 } else {
461 assert(y >= 0 && x0 >= 0 && x1 >= 0);
463 // printf(" Not adjoining; adding span id=%i, y=%i, x = [%i, %i]\n", id, y, x0, x1);
464 }
465 }
466 }
467
468 // At the end of this loop, all crpixel entries have been assigned an ID,
469 // except for the "dummy" entry at the end of the array.
470 if (crpixels.size() > 0) {
471 for (crpixel_iter cp = crpixels.begin(); cp != crpixels.end() - 1; cp++) {
472 assert(cp->id >= 0);
473 assert(cp->col >= 0);
474 assert(cp->row >= 0);
475 }
476 // dummy:
477 assert(crpixels[crpixels.size() - 1].id == -1);
478 assert(crpixels[crpixels.size() - 1].col == 0);
479 assert(crpixels[crpixels.size() - 1].row == -1);
480 }
481
482 for (std::vector<afw::detection::IdSpan::Ptr>::iterator sp = spans.begin(), end = spans.end(); sp != end;
483 sp++) {
484 assert((*sp)->id >= 0);
485 assert((*sp)->y >= 0);
486 assert((*sp)->x0 >= 0);
487 assert((*sp)->x1 >= (*sp)->x0);
488 for (std::vector<afw::detection::IdSpan::Ptr>::iterator sp2 = sp + 1; sp2 != end; sp2++) {
489 assert((*sp2)->y >= (*sp)->y);
490 if ((*sp2)->y == (*sp)->y) {
491 assert((*sp2)->x0 > (*sp)->x1);
492 }
493 }
494 }
495
496 /*
497 * See if spans touch each other
498 */
499 for (std::vector<afw::detection::IdSpan::Ptr>::iterator sp = spans.begin(), end = spans.end(); sp != end;
500 ++sp) {
501 int const y = (*sp)->y;
502 int const x0 = (*sp)->x0;
503 int const x1 = (*sp)->x1;
504
505 // this loop will probably run for only a few steps
506 for (std::vector<afw::detection::IdSpan::Ptr>::iterator sp2 = sp + 1; sp2 != end; ++sp2) {
507 if ((*sp2)->y == y) {
508 // on this row (but not adjoining columns, since it would have been merged into this span);
509 // keep looking.
510 continue;
511 } else if ((*sp2)->y != (y + 1)) {
512 // sp2 is more than one row below; can't be connected.
513 break;
514 } else if ((*sp2)->x0 > (x1 + 1)) {
515 // sp2 is more than one column away to the right; can't be connected
516 break;
517 } else if ((*sp2)->x1 >= (x0 - 1)) {
518 // touches
519 int r1 = afw::detection::resolve_alias(aliases, (*sp)->id);
520 int r2 = afw::detection::resolve_alias(aliases, (*sp2)->id);
521 aliases[r1] = r2;
522 }
523 }
524 }
525
526 /*
527 * Resolve aliases; first alias chains, then the IDs in the spans
528 */
529 for (unsigned int i = 0; i != spans.size(); ++i) {
530 spans[i]->id = afw::detection::resolve_alias(aliases, spans[i]->id);
531 }
532
533 /*
534 * Sort spans by ID, so we can sweep through them once
535 */
536 if (spans.size() > 0) {
537 std::sort(spans.begin(), spans.end(), afw::detection::IdSpanCompar());
538 }
539
540 /*
541 * Build Footprints from spans
542 */
544
545 if (spans.size() > 0) {
546 int id = spans[0]->id;
547 unsigned int i0 = 0; // initial value of i
548 for (unsigned int i = i0; i <= spans.size(); ++i) { // <= size to catch the last object
549 if (i == spans.size() || spans[i]->id != id) {
551
553 spanList.reserve(i - i0);
554 for (; i0 < i; ++i0) {
555 spanList.push_back(afw::geom::Span(spans[i0]->y, spans[i0]->x0, spans[i0]->x1));
556 }
557 cr->setSpans(std::make_shared<afw::geom::SpanSet>(std::move(spanList)));
558 CRs.push_back(cr);
559 }
560
561 if (i < spans.size()) {
562 id = spans[i]->id;
563 }
564 }
565 }
566
567 reinstateCrPixels(mimage.getImage().get(), crpixels);
568 /*
569 * apply condition #1
570 */
571 CountsInCR<typename ImageT::Pixel> CountDN(bkgd);
572 for (std::vector<std::shared_ptr<afw::detection::Footprint>>::iterator cr = CRs.begin(), end = CRs.end();
573 cr != end; ++cr) {
574 // find the sum of pixel values within the CR
575 (*cr)->getSpans()->applyFunctor(CountDN, *mimage.getImage());
576
577 LOGL_DEBUG("TRACE4.lsst.algorithms.CR", "CR at (%d, %d) has %g DN", (*cr)->getBBox().getMinX(),
578 (*cr)->getBBox().getMinY(), CountDN.getCounts());
579 if (CountDN.getCounts() < minDn) { /* not bright enough */
580 LOGL_DEBUG("TRACE5.lsst.algorithms.CR", "Erasing CR");
581
582 cr = CRs.erase(cr);
583 --cr; // back up to previous CR (we're going to increment it)
584 --end;
585 }
586 CountDN.reset();
587 }
588 ncr = CRs.size(); /* some may have been too faint */
589 /*
590 * We've found them all, time to kill them all
591 */
592 bool const debias_values = true;
593 bool grow = false;
594 LOGL_DEBUG("TRACE2.lsst.algorithms.CR", "Removing initial list of CRs");
595 removeCR(mimage, CRs, bkgd, crBit, saturBit, badMask, debias_values, grow);
596#if 0 // Useful to see phase 2 in display; debugging only
597 (void)setMaskFromFootprintList(mimage.getMask().get(), CRs,
598 mimage.getMask()->getPlaneBitMask("DETECTED"));
599#endif
600 /*
601 * Now that we've removed them, go through image again, examining area around
602 * each CR for extra bad pixels. Note that we set cond3Fac = 0 for this pass
603 *
604 * We iterate niteration times; niter==1 was sufficient for SDSS data, but megacam
605 * CCDs are different -- who knows for other devices?
606 */
607 bool too_many_crs = false; // we've seen too many CR pixels
608 int nextra = 0; // number of pixels added to list of CRs
609 for (int i = 0; i != niteration && !too_many_crs; ++i) {
610 LOGL_DEBUG("TRACE1.lsst.algorithms.CR", "Starting iteration %d", i);
612 fiter != CRs.end(); fiter++) {
614 /*
615 * Are all those `CR' pixels interpolated? If so, don't grow it
616 */
617 {
618 // this work should be taken on in DM-9538
619 // std::shared_ptr<afw::detection::Footprint> om = footprintAndMask(cr, mimage.getMask(),
620 // interpBit);
621 // Placeholder until DM-9538 then remove empty footprint
622 auto om = std::make_shared<afw::detection::Footprint>();
623 int const npix = (om) ? om->getArea() : 0;
624
625 if (static_cast<std::size_t>(npix) == cr->getArea()) {
626 continue;
627 }
628 }
629 /*
630 * No; some of the suspect pixels aren't interpolated
631 */
632 afw::detection::Footprint extra; // extra pixels added to cr
633 for (auto siter = cr->getSpans()->begin(); siter != cr->getSpans()->end(); siter++) {
634 auto const span = siter;
635
636 /*
637 * Check the lines above and below the span. We're going to check a 3x3 region around
638 * the pixels, so we need a buffer around the edge. We check the pixels just to the
639 * left/right of the span, so the buffer needs to be 2 pixels (not just 1) in the
640 * column direction, but only 1 in the row direction.
641 */
642 int const y = span->getY() - mimage.getY0();
643 if (y < 2 || y >= nrow - 2) {
644 continue;
645 }
646 int x0 = span->getX0() - mimage.getX0();
647 int x1 = span->getX1() - mimage.getX0();
648 x0 = (x0 < 2) ? 2 : (x0 > ncol - 3) ? ncol - 3 : x0;
649 x1 = (x1 < 2) ? 2 : (x1 > ncol - 3) ? ncol - 3 : x1;
650
651 checkSpanForCRs(&extra, crpixels, y - 1, x0, x1, mimage, minSigma / 2, thresH, thresV, thresD,
652 bkgd, 0, keep);
653 checkSpanForCRs(&extra, crpixels, y, x0, x1, mimage, minSigma / 2, thresH, thresV, thresD,
654 bkgd, 0, keep);
655 checkSpanForCRs(&extra, crpixels, y + 1, x0, x1, mimage, minSigma / 2, thresH, thresV, thresD,
656 bkgd, 0, keep);
657 }
658
659 if (extra.getSpans()->size() > 0) { // we added some pixels
660 if (nextra + static_cast<int>(crpixels.size()) > nCrPixelMax) {
661 too_many_crs = true;
662 break;
663 }
664
665 nextra += extra.getArea();
666
667 std::vector<afw::geom::Span> tmpSpanList(cr->getSpans()->begin(), cr->getSpans()->end());
668 for (auto const &spn : (*extra.getSpans())) {
669 tmpSpanList.push_back(spn);
670 }
671 cr->setSpans(std::make_shared<afw::geom::SpanSet>(std::move(tmpSpanList)));
672 }
673 }
674
675 if (nextra == 0) {
676 break;
677 }
678 }
679 /*
680 * mark those pixels as CRs
681 */
682 if (!too_many_crs) {
683 for (auto const &foot : CRs) {
684 foot->getSpans()->setMask(*mimage.getMask().get(), crBit);
685 }
686 }
687 /*
688 * Maybe reinstate initial values; n.b. the same pixel may appear twice, so we want the
689 * first value stored (hence the uses of rbegin/rend)
690 *
691 * We have to do this if we decide _not_ to remove certain CRs,
692 * for example those which lie next to saturated pixels
693 */
694 if (keep || too_many_crs) {
695 if (crpixels.size() > 0) {
696 int const imageX0 = mimage.getX0();
697 int const imageY0 = mimage.getY0();
698
699 std::sort(crpixels.begin(), crpixels.end()); // sort into birth order
700
701 crpixel_riter rend = crpixels.rend();
702 for (crpixel_riter crp = crpixels.rbegin(); crp != rend; ++crp) {
703 if (crp->row == -1)
704 // dummy; skip it.
705 continue;
706 mimage.at(crp->col - imageX0, crp->row - imageY0).image() = crp->val;
707 }
708 }
709 } else {
710 if (true || nextra > 0) {
711 grow = true;
712 LOGL_DEBUG("TRACE2.lsst.algorithms.CR", "Removing final list of CRs, grow = %d", grow);
713 removeCR(mimage, CRs, bkgd, crBit, saturBit, badMask, debias_values, grow);
714 }
715 /*
716 * we interpolated over all CR pixels, so set the interp bits too
717 */
718 for (auto const &foot : CRs) {
719 foot->getSpans()->setMask(*mimage.getMask().get(), static_cast<MaskPixel>(crBit | interpBit));
720 }
721 }
722
723 if (too_many_crs) { // we've cleaned up, so we can throw the exception
725 (boost::format("Too many CR pixels (max %d)") % nCrPixelMax).str());
726 }
727
728 return CRs;
729}
730
731/*****************************************************************************/
732namespace {
733/*
734 * Is condition 3 true?
735 */
736template <typename ImageT>
737bool condition_3(ImageT *estimate, // estimate of true value of pixel
738 double const peak, // counts in central pixel (no sky)
739 double const mean_ns, // mean in NS direction (no sky)
740 double const mean_we, // " " WE " " " "
741 double const mean_swne, // " " SW-NE " " " "
742 double const mean_nwse, // " " NW-SE " " " "
743 double const dpeak, // standard deviation of peak value
744 double const dmean_ns, // s.d. of mean in NS direction
745 double const dmean_we, // " " " " WE " "
746 double const dmean_swne, // " " " " SW-NE " "
747 double const dmean_nwse, // " " " " NW-SE " "
748 double const thresH, // horizontal threshold
749 double const thresV, // vertical threshold
750 double const thresD, // diagonal threshold
751 double const cond3Fac // fiddle factor for noise
752) {
753 if (thresV * (peak - cond3Fac * dpeak) > mean_ns + cond3Fac * dmean_ns) {
754 *estimate = (ImageT)mean_ns;
755 return true;
756 }
757
758 if (thresH * (peak - cond3Fac * dpeak) > mean_we + cond3Fac * dmean_we) {
759 *estimate = mean_we;
760 return true;
761 }
762
763 if (thresD * (peak - cond3Fac * dpeak) > mean_swne + cond3Fac * dmean_swne) {
764 *estimate = mean_swne;
765 return true;
766 }
767
768 if (thresD * (peak - cond3Fac * dpeak) > mean_nwse + cond3Fac * dmean_nwse) {
769 *estimate = mean_nwse;
770 return true;
771 }
772
773 return false;
774}
775
776/************************************************************************************************************/
777/*
778 * Interpolate over a CR's pixels
779 */
780template <typename MaskedImageT>
781class RemoveCR {
782public:
783 RemoveCR(MaskedImageT const &mimage, double const bkgd, typename MaskedImageT::Mask::Pixel badMask,
784 bool const debias, afw::math::Random &rand)
785 : _image(mimage),
786 _bkgd(bkgd),
787 _ncol(mimage.getWidth()),
788 _nrow(mimage.getHeight()),
789 _badMask(badMask),
790 _debias(debias),
791 _rand(rand) {}
792
793 void operator()(geom::Point2I const &point, int) {
794 int const &x = point.getX();
795 int const &y = point.getY();
796 typename MaskedImageT::xy_locator loc = _image.xy_at(x - _image.getX0(), y - _image.getY0());
797 typedef typename MaskedImageT::Image::Pixel MImagePixel;
799 int ngood = 0; // number of good values on min
800
801 MImagePixel const minval =
802 _bkgd - 2 * sqrt(loc.variance()); // min. acceptable pixel value after interp
803 /*
804 * W-E row
805 */
806 if (x - 2 >= 0 && x + 2 < _ncol) {
807 if ((loc.mask(-2, 0) | _badMask) || (loc.mask(-1, 0) | _badMask) || (loc.mask(1, 0) | _badMask) ||
808 (loc.mask(2, 0) | _badMask)) {
809 ; // estimate is contaminated
810 } else {
811 MImagePixel const v_m2 = loc.image(-2, 0);
812 MImagePixel const v_m1 = loc.image(-1, 0);
813 MImagePixel const v_p1 = loc.image(1, 0);
814 MImagePixel const v_p2 = loc.image(2, 0);
815
816 MImagePixel const tmp = interp::lpc_1_c1 * (v_m1 + v_p1) + interp::lpc_1_c2 * (v_m2 + v_p2);
817
818 if (tmp > minval && tmp < min) {
819 min = tmp;
820 ngood++;
821 }
822 }
823 }
824 /*
825 * N-S column
826 */
827 if (y - 2 >= 0 && y + 2 < _nrow) {
828 if ((loc.mask(0, -2) | _badMask) || (loc.mask(0, -1) | _badMask) || (loc.mask(0, 1) | _badMask) ||
829 (loc.mask(0, 2) | _badMask)) {
830 ; /* estimate is contaminated */
831 } else {
832 MImagePixel const v_m2 = loc.image(0, -2);
833 MImagePixel const v_m1 = loc.image(0, -1);
834 MImagePixel const v_p1 = loc.image(0, 1);
835 MImagePixel const v_p2 = loc.image(0, 2);
836
837 MImagePixel const tmp = interp::lpc_1_c1 * (v_m1 + v_p1) + interp::lpc_1_c2 * (v_m2 + v_p2);
838
839 if (tmp > minval && tmp < min) {
840 min = tmp;
841 ngood++;
842 }
843 }
844 }
845 /*
846 * SW--NE diagonal
847 */
848 if (x - 2 >= 0 && x + 2 < _ncol && y - 2 >= 0 && y + 2 < _nrow) {
849 if ((loc.mask(-2, -2) | _badMask) || (loc.mask(-1, -1) | _badMask) ||
850 (loc.mask(1, 1) | _badMask) || (loc.mask(2, 2) | _badMask)) {
851 ; /* estimate is contaminated */
852 } else {
853 MImagePixel const v_m2 = loc.image(-2, -2);
854 MImagePixel const v_m1 = loc.image(-1, -1);
855 MImagePixel const v_p1 = loc.image(1, 1);
856 MImagePixel const v_p2 = loc.image(2, 2);
857
858 MImagePixel const tmp =
859 interp::lpc_1s2_c1 * (v_m1 + v_p1) + interp::lpc_1s2_c2 * (v_m2 + v_p2);
860
861 if (tmp > minval && tmp < min) {
862 min = tmp;
863 ngood++;
864 }
865 }
866 }
867 /*
868 * SE--NW diagonal
869 */
870 if (x - 2 >= 0 && x + 2 < _ncol && y - 2 >= 0 && y + 2 < _nrow) {
871 if ((loc.mask(2, -2) | _badMask) || (loc.mask(1, -1) | _badMask) ||
872 (loc.mask(-1, 1) | _badMask) || (loc.mask(-2, 2) | _badMask)) {
873 ; /* estimate is contaminated */
874 } else {
875 MImagePixel const v_m2 = loc.image(2, -2);
876 MImagePixel const v_m1 = loc.image(1, -1);
877 MImagePixel const v_p1 = loc.image(-1, 1);
878 MImagePixel const v_p2 = loc.image(-2, 2);
879
880 MImagePixel const tmp =
881 interp::lpc_1s2_c1 * (v_m1 + v_p1) + interp::lpc_1s2_c2 * (v_m2 + v_p2);
882
883 if (tmp > minval && tmp < min) {
884 min = tmp;
885 ngood++;
886 }
887 }
888 }
889 /*
890 * Have we altogether failed to find an acceptable value? If so interpolate
891 * using the full-up interpolation code both vertically and horizontally
892 * and take the average. This can fail for large enough defects (e.g. CRs
893 * lying in bad columns), in which case the interpolator returns -1. If
894 * both directions fail, use the background value.
895 */
896 if (ngood == 0) {
897 std::pair<bool, MImagePixel const> val_h = interp::singlePixel(x, y, _image, true, minval);
898 std::pair<bool, MImagePixel const> val_v = interp::singlePixel(x, y, _image, false, minval);
899
900 if (!val_h.first) {
901 if (!val_v.first) { // Still no good value. Guess wildly
902 min = _bkgd + sqrt(loc.variance()) * _rand.gaussian();
903 } else {
904 min = val_v.second;
905 }
906 } else {
907 if (val_v.first) {
908 min = val_h.second;
909 } else {
910 min = (val_v.second + val_h.second) / 2;
911 }
912 }
913 }
914 /*
915 * debias the minimum; If more than one uncontaminated estimate was
916 * available, estimate the bias to be simply that due to choosing the
917 * minimum of two Gaussians. In fact, even some of the "good" pixels
918 * may have some extra charge, so even if ngood > 2, still use this
919 * estimate
920 */
921 if (ngood > 0) {
922 LOGL_DEBUG("TRACE3.lsst.algorithms.CR", "Adopted min==%g at (%d, %d) (ngood=%d)",
923 static_cast<double>(min), x, y, ngood);
924 }
925
926 if (_debias && ngood > 1) {
927 min -= interp::min2GaussianBias * sqrt(loc.variance()) * _rand.gaussian();
928 }
929
930 loc.image() = min;
931 }
932
933private:
934 MaskedImageT const &_image;
935 double _bkgd;
936 int _ncol, _nrow;
937 typename MaskedImageT::Mask::Pixel _badMask;
938 bool _debias;
939 afw::math::Random &_rand;
940};
941
942/************************************************************************************************************/
943/*
944 * actually remove CRs from the frame
945 */
946template <typename ImageT, typename MaskT>
947void removeCR(afw::image::MaskedImage<ImageT, MaskT> &mi, // image to search
948 std::vector<std::shared_ptr<afw::detection::Footprint>> &CRs, // list of cosmic rays
949 double const bkgd, // non-subtracted background
950 MaskT const, // Bit value used to label CRs
951 MaskT const saturBit, // Bit value used to label saturated pixels
952 MaskT const badMask, // Bit mask for bad pixels
953 bool const debias, // statistically debias values?
954 bool const grow // Grow CRs?
955) {
956 afw::math::Random rand; // a random number generator
957 /*
958 * replace the values of cosmic-ray contaminated pixels with 1-dim 2nd-order weighted means Cosmic-ray
959 * contaminated pixels have already been given a mask value, crBit
960 *
961 * If there are no good options (i.e. all estimates are contaminated), try using just pixels that are not
962 * CRs; failing that, interpolate in the row- or column direction over as large a distance as is required
963 *
964 * XXX SDSS (and we) go through this list backwards; why?
965 */
966
967 // a functor to remove a CR
968 RemoveCR<afw::image::MaskedImage<ImageT, MaskT>> removeCR(mi, bkgd, badMask, debias, rand);
969
970 for (std::vector<std::shared_ptr<afw::detection::Footprint>>::reverse_iterator fiter = CRs.rbegin();
971 fiter != CRs.rend(); ++fiter) {
973 /*
974 * If I grow this CR does it touch saturated pixels? If so, don't
975 * interpolate and add CR pixels to saturated mask
976 */
977 /* This is commented out pending work from DM-9538
978 if (grow && cr->getArea() < 100) {
979 try {
980 bool const isotropic = false; // use a slow isotropic grow?
981 auto gcr = std::make_shared<afw::detection::Footprint>(cr->getSpans()->dilate(1),
982 cr.getRegion()); std::shared_ptr<afw::detection::Footprint> const saturPixels =
983 footprintAndMask(gcr, mi.getMask(), saturBit); auto const saturPixels = footprintAndMask(gcr,
984 mi.getMask(), saturBit);
985
986 if (saturPixels->getArea() > 0) { // pixel is adjacent to a saturation trail
987 setMaskFromFootprint(mi.getMask().get(), *saturPixels, saturBit);
988
989 continue;
990 }
991 } catch(pex::exceptions::LengthError &) {
992 continue;
993 }
994 }
995 */
996 /*
997 * OK, fix it
998 */
999 // applyFunctor must have an argument other than the functor, however in this case
1000 // additional arguments are unneeded, so pass a constant 1 to be iterated over
1001 // and ignore
1002 cr->getSpans()->applyFunctor(removeCR, 1);
1003 }
1004}
1005} // namespace
1006
1007/************************************************************************************************************/
1008//
1009// Explicit instantiations
1010// \cond
1011#define INSTANTIATE(TYPE) \
1012 template std::vector<std::shared_ptr<afw::detection::Footprint>> findCosmicRays( \
1013 afw::image::MaskedImage<TYPE> &image, afw::detection::Psf const &psf, double const bkgd, \
1014 daf::base::PropertySet const &ps, bool const keep)
1015
1016INSTANTIATE(float);
1017INSTANTIATE(double); // Why do we need double images?
1018// \endcond
1019} // namespace algorithms
1020} // namespace meas
1021} // namespace lsst
int min
int end
double x
#define INSTANTIATE(FROMSYS, TOSYS)
Definition: Detector.cc:484
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
LSST DM logging module built on log4cxx.
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
Definition: Log.h:515
Random & _rand
Definition: RandomImage.cc:43
int y
Definition: SpanSet.cc:48
table::Key< int > b
table::Key< int > a
T begin(T... args)
T capacity(T... args)
Class to describe the properties of a detected object from an image.
Definition: Footprint.h:63
std::shared_ptr< geom::SpanSet > getSpans() const
Return a shared pointer to the SpanSet.
Definition: Footprint.h:115
std::size_t getArea() const
Return the number of pixels in this Footprint.
Definition: Footprint.h:173
run-length code for part of object
Definition: CR.cc:63
std::shared_ptr< IdSpan > Ptr
Definition: CR.cc:65
std::shared_ptr< const IdSpan > ConstPtr
Definition: CR.cc:66
IdSpan(int id, int y, int x0, int x1)
Definition: CR.cc:68
A polymorphic base class for representing an image's Point Spread Function.
Definition: Psf.h:76
image::Image< Pixel > Image
Image type returned by computeImage.
Definition: Psf.h:83
A range of pixels within one row of an Image.
Definition: Span.h:47
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:51
Class for storing generic metadata.
Definition: PropertySet.h:66
int getAsInt(std::string const &name) const
Get the last value for a bool/char/short/int property name (possibly hierarchical).
double getAsDouble(std::string const &name) const
Get the last value for any arithmetic property name (possibly hierarchical).
Reports attempts to exceed implementation-defined length limits for some classes.
Definition: Runtime.h:76
Reports attempts to access elements using an invalid key.
Definition: Runtime.h:151
T empty(T... args)
T end(T... args)
T erase(T... args)
T max(T... args)
T move(T... args)
int resolve_alias(const std::vector< int > &aliases, int id)
Follow a chain of aliases, returning the final resolved value.
Definition: CR.cc:96
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
std::int32_t MaskPixel
default type for Masks and MaskedImage Masks
lsst::afw::detection::Footprint Footprint
Definition: Source.h:57
double const min2GaussianBias
Mean value of the minimum of two N(0,1) variates.
Definition: Interp.h:62
double const lpc_1_c1
LPC coefficients for sigma = 1, S/N = infty.
Definition: Interp.h:51
double const lpc_1s2_c1
LPC coefficients for sigma = 1/sqrt(2), S/N = infty.
Definition: Interp.h:57
std::pair< bool, typename MaskedImageT::Image::Pixel > singlePixel(int x, int y, MaskedImageT const &image, bool horizontal, double minval)
Return a boolean status (true: interpolation is OK) and the interpolated value for a pixel,...
Definition: Interp.cc:2125
std::vector< std::shared_ptr< afw::detection::Footprint > > findCosmicRays(MaskedImageT &image, afw::detection::Psf const &psf, double const bkgd, daf::base::PropertySet const &ps, bool const keep=false)
Find cosmic rays in an Image, and mask and remove them.
Definition: CR.cc:315
float Pixel
Typedefs to be used for pixel values.
Definition: common.h:37
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
A base class for image defects.
T push_back(T... args)
T rand(T... args)
T rbegin(T... args)
T rend(T... args)
T reserve(T... args)
T size(T... args)
T sort(T... args)
T sqrt(T... args)
ImageT val
Definition: CR.cc:146
int id
Definition: CR.cc:143
int row
Definition: CR.cc:145
int col
Definition: CR.cc:144
comparison functor; sort by ID, then by row (y), then by column range start (x0)
Definition: CR.cc:76
bool operator()(const IdSpan::ConstPtr a, const IdSpan::ConstPtr b)
Definition: CR.cc:77
Key< int > psf
Definition: Exposure.cc:65