LSST Applications g04e9c324dd+8c5ae1fdc5,g134cb467dc+1b3060144d,g18429d2f64+f642bf4753,g199a45376c+0ba108daf9,g1fd858c14a+2dcf163641,g262e1987ae+7b8c96d2ca,g29ae962dfc+3bd6ecb08a,g2cef7863aa+aef1011c0b,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+53e1a9e7c5,g4595892280+fef73a337f,g47891489e3+2efcf17695,g4d44eb3520+642b70b07e,g53246c7159+8c5ae1fdc5,g67b6fd64d1+2efcf17695,g67fd3c3899+b70e05ef52,g74acd417e5+317eb4c7d4,g786e29fd12+668abc6043,g87389fa792+8856018cbb,g89139ef638+2efcf17695,g8d7436a09f+3be3c13596,g8ea07a8fe4+9f5ccc88ac,g90f42f885a+a4e7b16d9b,g97be763408+ad77d7208f,g9dd6db0277+b70e05ef52,ga681d05dcb+a3f46e7fff,gabf8522325+735880ea63,gac2eed3f23+2efcf17695,gb89ab40317+2efcf17695,gbf99507273+8c5ae1fdc5,gd8ff7fe66e+b70e05ef52,gdab6d2f7ff+317eb4c7d4,gdc713202bf+b70e05ef52,gdfd2d52018+b10e285e0f,ge365c994fd+310e8507c4,ge410e46f29+2efcf17695,geaed405ab2+562b3308c0,gffca2db377+8c5ae1fdc5,w.2025.35
LSST Data Management Base Package
Loading...
Searching...
No Matches
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
53/*
54 * TODO: These should go into afw --- actually, there're already there, but
55 * in an anon namespace (DM-35750).
56 */
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};
73
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};
93
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 */
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
430
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 /*
469 * See if spans touch each other
470 */
471 for (std::vector<afw::detection::IdSpan::Ptr>::iterator sp = spans.begin(), end = spans.end(); sp != end;
472 ++sp) {
473 int const y = (*sp)->y;
474 int const x0 = (*sp)->x0;
475 int const x1 = (*sp)->x1;
476
477 // this loop will probably run for only a few steps
478 for (std::vector<afw::detection::IdSpan::Ptr>::iterator sp2 = sp + 1; sp2 != end; ++sp2) {
479 if ((*sp2)->y == y) {
480 // on this row (but not adjoining columns, since it would have been merged into this span);
481 // keep looking.
482 continue;
483 } else if ((*sp2)->y != (y + 1)) {
484 // sp2 is more than one row below; can't be connected.
485 break;
486 } else if ((*sp2)->x0 > (x1 + 1)) {
487 // sp2 is more than one column away to the right; can't be connected
488 break;
489 } else if ((*sp2)->x1 >= (x0 - 1)) {
490 // touches
491 int r1 = afw::detection::resolve_alias(aliases, (*sp)->id);
492 int r2 = afw::detection::resolve_alias(aliases, (*sp2)->id);
493 aliases[r1] = r2;
494 }
495 }
496 }
497
498 /*
499 * Resolve aliases; first alias chains, then the IDs in the spans
500 */
501 for (unsigned int i = 0; i != spans.size(); ++i) {
502 spans[i]->id = afw::detection::resolve_alias(aliases, spans[i]->id);
503 }
504
505 /*
506 * Sort spans by ID, so we can sweep through them once
507 */
508 if (spans.size() > 0) {
509 std::sort(spans.begin(), spans.end(), afw::detection::IdSpanCompar());
510 }
511
512 /*
513 * Build Footprints from spans
514 */
516
517 if (spans.size() > 0) {
518 int id = spans[0]->id;
519 unsigned int i0 = 0; // initial value of i
520 for (unsigned int i = i0; i <= spans.size(); ++i) { // <= size to catch the last object
521 if (i == spans.size() || spans[i]->id != id) {
523
525 spanList.reserve(i - i0);
526 for (; i0 < i; ++i0) {
527 spanList.push_back(afw::geom::Span(spans[i0]->y, spans[i0]->x0, spans[i0]->x1));
528 }
529 cr->setSpans(std::make_shared<afw::geom::SpanSet>(std::move(spanList)));
530 CRs.push_back(cr);
531 }
532
533 if (i < spans.size()) {
534 id = spans[i]->id;
535 }
536 }
537 }
538
539 reinstateCrPixels(mimage.getImage().get(), crpixels);
540 /*
541 * apply condition #1
542 */
543 CountsInCR<typename ImageT::Pixel> CountDN(bkgd);
544 for (std::vector<std::shared_ptr<afw::detection::Footprint>>::iterator cr = CRs.begin(), end = CRs.end();
545 cr != end; ++cr) {
546 // find the sum of pixel values within the CR
547 (*cr)->getSpans()->applyFunctor(CountDN, *mimage.getImage());
548
549 LOGL_DEBUG("TRACE4.lsst.algorithms.CR", "CR at (%d, %d) has %g DN", (*cr)->getBBox().getMinX(),
550 (*cr)->getBBox().getMinY(), CountDN.getCounts());
551 if (CountDN.getCounts() < minDn) { /* not bright enough */
552 LOGL_DEBUG("TRACE5.lsst.algorithms.CR", "Erasing CR");
553
554 cr = CRs.erase(cr);
555 --cr; // back up to previous CR (we're going to increment it)
556 --end;
557 }
558 CountDN.reset();
559 }
560 ncr = CRs.size(); /* some may have been too faint */
561 /*
562 * We've found them all, time to kill them all
563 */
564 bool const debias_values = true;
565 bool grow = false;
566 LOGL_DEBUG("TRACE2.lsst.algorithms.CR", "Removing initial list of CRs");
567 removeCR(mimage, CRs, bkgd, crBit, saturBit, badMask, debias_values, grow);
568#if 0 // Useful to see phase 2 in display; debugging only
569 (void)setMaskFromFootprintList(mimage.getMask().get(), CRs,
570 mimage.getMask()->getPlaneBitMask("DETECTED"));
571#endif
572 /*
573 * Now that we've removed them, go through image again, examining area around
574 * each CR for extra bad pixels. Note that we set cond3Fac = 0 for this pass
575 *
576 * We iterate niteration times; niter==1 was sufficient for SDSS data, but megacam
577 * CCDs are different -- who knows for other devices?
578 */
579 bool too_many_crs = false; // we've seen too many CR pixels
580 int nextra = 0; // number of pixels added to list of CRs
581 for (int i = 0; i != niteration && !too_many_crs; ++i) {
582 LOGL_DEBUG("TRACE1.lsst.algorithms.CR", "Starting iteration %d", i);
584 fiter != CRs.end(); fiter++) {
586 /*
587 * Are all those `CR' pixels interpolated? If so, don't grow it
588 */
589 {
590 // this work should be taken on in DM-9538
591 // std::shared_ptr<afw::detection::Footprint> om = footprintAndMask(cr, mimage.getMask(),
592 // interpBit);
593 // Placeholder until DM-9538 then remove empty footprint
595 int const npix = (om) ? om->getArea() : 0;
596
597 if (static_cast<std::size_t>(npix) == cr->getArea()) {
598 continue;
599 }
600 }
601 /*
602 * No; some of the suspect pixels aren't interpolated
603 */
604 afw::detection::Footprint extra; // extra pixels added to cr
605 for (auto siter = cr->getSpans()->begin(); siter != cr->getSpans()->end(); siter++) {
606 auto const span = siter;
607
608 /*
609 * Check the lines above and below the span. We're going to check a 3x3 region around
610 * the pixels, so we need a buffer around the edge. We check the pixels just to the
611 * left/right of the span, so the buffer needs to be 2 pixels (not just 1) in the
612 * column direction, but only 1 in the row direction.
613 */
614 int const y = span->getY() - mimage.getY0();
615 if (y < 2 || y >= nrow - 2) {
616 continue;
617 }
618 int x0 = span->getX0() - mimage.getX0();
619 int x1 = span->getX1() - mimage.getX0();
620 x0 = (x0 < 2) ? 2 : (x0 > ncol - 3) ? ncol - 3 : x0;
621 x1 = (x1 < 2) ? 2 : (x1 > ncol - 3) ? ncol - 3 : x1;
622
623 checkSpanForCRs(&extra, crpixels, y - 1, x0, x1, mimage, minSigma / 2, thresH, thresV, thresD,
624 bkgd, 0, keep);
625 checkSpanForCRs(&extra, crpixels, y, x0, x1, mimage, minSigma / 2, thresH, thresV, thresD,
626 bkgd, 0, keep);
627 checkSpanForCRs(&extra, crpixels, y + 1, x0, x1, mimage, minSigma / 2, thresH, thresV, thresD,
628 bkgd, 0, keep);
629 }
630
631 if (extra.getSpans()->size() > 0) { // we added some pixels
632 if (nextra + static_cast<int>(crpixels.size()) > nCrPixelMax) {
633 too_many_crs = true;
634 break;
635 }
636
637 nextra += extra.getArea();
638
639 std::vector<afw::geom::Span> tmpSpanList(cr->getSpans()->begin(), cr->getSpans()->end());
640 for (auto const &spn : (*extra.getSpans())) {
641 tmpSpanList.push_back(spn);
642 }
643 cr->setSpans(std::make_shared<afw::geom::SpanSet>(std::move(tmpSpanList)));
644 }
645 }
646
647 if (nextra == 0) {
648 break;
649 }
650 }
651 /*
652 * mark those pixels as CRs
653 */
654 if (!too_many_crs) {
655 for (auto const &foot : CRs) {
656 foot->getSpans()->setMask(*mimage.getMask().get(), crBit);
657 }
658 }
659 /*
660 * Maybe reinstate initial values; n.b. the same pixel may appear twice, so we want the
661 * first value stored (hence the uses of rbegin/rend)
662 *
663 * We have to do this if we decide _not_ to remove certain CRs,
664 * for example those which lie next to saturated pixels
665 */
666 if (keep || too_many_crs) {
667 if (crpixels.size() > 0) {
668 int const imageX0 = mimage.getX0();
669 int const imageY0 = mimage.getY0();
670
671 std::sort(crpixels.begin(), crpixels.end()); // sort into birth order
672
673 crpixel_riter rend = crpixels.rend();
674 for (crpixel_riter crp = crpixels.rbegin(); crp != rend; ++crp) {
675 if (crp->row == -1)
676 // dummy; skip it.
677 continue;
678 mimage.at(crp->col - imageX0, crp->row - imageY0).image() = crp->val;
679 }
680 }
681 } else {
682 if (true || nextra > 0) {
683 grow = true;
684 LOGL_DEBUG("TRACE2.lsst.algorithms.CR", "Removing final list of CRs, grow = %d", grow);
685 removeCR(mimage, CRs, bkgd, crBit, saturBit, badMask, debias_values, grow);
686 }
687 /*
688 * we interpolated over all CR pixels, so set the interp bits too
689 */
690 for (auto const &foot : CRs) {
691 foot->getSpans()->setMask(*mimage.getMask().get(), static_cast<MaskPixel>(crBit | interpBit));
692 }
693 }
694
695 if (too_many_crs) { // we've cleaned up, so we can throw the exception
697 (boost::format("Too many CR pixels (max %d)") % nCrPixelMax).str());
698 }
699
700 return CRs;
701}
702
703/*****************************************************************************/
704namespace {
705/*
706 * Is condition 3 true?
707 */
708template <typename ImageT>
709bool condition_3(ImageT *estimate, // estimate of true value of pixel
710 double const peak, // counts in central pixel (no sky)
711 double const mean_ns, // mean in NS direction (no sky)
712 double const mean_we, // " " WE " " " "
713 double const mean_swne, // " " SW-NE " " " "
714 double const mean_nwse, // " " NW-SE " " " "
715 double const dpeak, // standard deviation of peak value
716 double const dmean_ns, // s.d. of mean in NS direction
717 double const dmean_we, // " " " " WE " "
718 double const dmean_swne, // " " " " SW-NE " "
719 double const dmean_nwse, // " " " " NW-SE " "
720 double const thresH, // horizontal threshold
721 double const thresV, // vertical threshold
722 double const thresD, // diagonal threshold
723 double const cond3Fac // fiddle factor for noise
724) {
725 if (thresV * (peak - cond3Fac * dpeak) > mean_ns + cond3Fac * dmean_ns) {
726 *estimate = (ImageT)mean_ns;
727 return true;
728 }
729
730 if (thresH * (peak - cond3Fac * dpeak) > mean_we + cond3Fac * dmean_we) {
731 *estimate = mean_we;
732 return true;
733 }
734
735 if (thresD * (peak - cond3Fac * dpeak) > mean_swne + cond3Fac * dmean_swne) {
736 *estimate = mean_swne;
737 return true;
738 }
739
740 if (thresD * (peak - cond3Fac * dpeak) > mean_nwse + cond3Fac * dmean_nwse) {
741 *estimate = mean_nwse;
742 return true;
743 }
744
745 return false;
746}
747
748/************************************************************************************************************/
749/*
750 * Interpolate over a CR's pixels
751 */
752template <typename MaskedImageT>
753class RemoveCR {
754public:
755 RemoveCR(MaskedImageT const &mimage, double const bkgd, typename MaskedImageT::Mask::Pixel badMask,
756 bool const debias, afw::math::Random &rand)
757 : _image(mimage),
758 _bkgd(bkgd),
759 _ncol(mimage.getWidth()),
760 _nrow(mimage.getHeight()),
761 _badMask(badMask),
762 _debias(debias),
763 _rand(rand) {}
764
765 void operator()(geom::Point2I const &point, int) {
766 int const &x = point.getX();
767 int const &y = point.getY();
768 typename MaskedImageT::xy_locator loc = _image.xy_at(x - _image.getX0(), y - _image.getY0());
769 typedef typename MaskedImageT::Image::Pixel MImagePixel;
771 int ngood = 0; // number of good values on min
772
773 MImagePixel const minval =
774 _bkgd - 2 * sqrt(loc.variance()); // min. acceptable pixel value after interp
775 /*
776 * W-E row
777 */
778 if (x - 2 >= 0 && x + 2 < _ncol) {
779 if ((loc.mask(-2, 0) | _badMask) || (loc.mask(-1, 0) | _badMask) || (loc.mask(1, 0) | _badMask) ||
780 (loc.mask(2, 0) | _badMask)) {
781 ; // estimate is contaminated
782 } else {
783 MImagePixel const v_m2 = loc.image(-2, 0);
784 MImagePixel const v_m1 = loc.image(-1, 0);
785 MImagePixel const v_p1 = loc.image(1, 0);
786 MImagePixel const v_p2 = loc.image(2, 0);
787
788 MImagePixel const tmp = interp::lpc_1_c1 * (v_m1 + v_p1) + interp::lpc_1_c2 * (v_m2 + v_p2);
789
790 if (tmp > minval && tmp < min) {
791 min = tmp;
792 ngood++;
793 }
794 }
795 }
796 /*
797 * N-S column
798 */
799 if (y - 2 >= 0 && y + 2 < _nrow) {
800 if ((loc.mask(0, -2) | _badMask) || (loc.mask(0, -1) | _badMask) || (loc.mask(0, 1) | _badMask) ||
801 (loc.mask(0, 2) | _badMask)) {
802 ; /* estimate is contaminated */
803 } else {
804 MImagePixel const v_m2 = loc.image(0, -2);
805 MImagePixel const v_m1 = loc.image(0, -1);
806 MImagePixel const v_p1 = loc.image(0, 1);
807 MImagePixel const v_p2 = loc.image(0, 2);
808
809 MImagePixel const tmp = interp::lpc_1_c1 * (v_m1 + v_p1) + interp::lpc_1_c2 * (v_m2 + v_p2);
810
811 if (tmp > minval && tmp < min) {
812 min = tmp;
813 ngood++;
814 }
815 }
816 }
817 /*
818 * SW--NE diagonal
819 */
820 if (x - 2 >= 0 && x + 2 < _ncol && y - 2 >= 0 && y + 2 < _nrow) {
821 if ((loc.mask(-2, -2) | _badMask) || (loc.mask(-1, -1) | _badMask) ||
822 (loc.mask(1, 1) | _badMask) || (loc.mask(2, 2) | _badMask)) {
823 ; /* estimate is contaminated */
824 } else {
825 MImagePixel const v_m2 = loc.image(-2, -2);
826 MImagePixel const v_m1 = loc.image(-1, -1);
827 MImagePixel const v_p1 = loc.image(1, 1);
828 MImagePixel const v_p2 = loc.image(2, 2);
829
830 MImagePixel const tmp =
831 interp::lpc_1s2_c1 * (v_m1 + v_p1) + interp::lpc_1s2_c2 * (v_m2 + v_p2);
832
833 if (tmp > minval && tmp < min) {
834 min = tmp;
835 ngood++;
836 }
837 }
838 }
839 /*
840 * SE--NW diagonal
841 */
842 if (x - 2 >= 0 && x + 2 < _ncol && y - 2 >= 0 && y + 2 < _nrow) {
843 if ((loc.mask(2, -2) | _badMask) || (loc.mask(1, -1) | _badMask) ||
844 (loc.mask(-1, 1) | _badMask) || (loc.mask(-2, 2) | _badMask)) {
845 ; /* estimate is contaminated */
846 } else {
847 MImagePixel const v_m2 = loc.image(2, -2);
848 MImagePixel const v_m1 = loc.image(1, -1);
849 MImagePixel const v_p1 = loc.image(-1, 1);
850 MImagePixel const v_p2 = loc.image(-2, 2);
851
852 MImagePixel const tmp =
853 interp::lpc_1s2_c1 * (v_m1 + v_p1) + interp::lpc_1s2_c2 * (v_m2 + v_p2);
854
855 if (tmp > minval && tmp < min) {
856 min = tmp;
857 ngood++;
858 }
859 }
860 }
861 /*
862 * Have we altogether failed to find an acceptable value? If so interpolate
863 * using the full-up interpolation code both vertically and horizontally
864 * and take the average. This can fail for large enough defects (e.g. CRs
865 * lying in bad columns), in which case the interpolator returns -1. If
866 * both directions fail, use the background value.
867 */
868 if (ngood == 0) {
869 std::pair<bool, MImagePixel const> val_h = interp::singlePixel(x, y, _image, true, minval);
870 std::pair<bool, MImagePixel const> val_v = interp::singlePixel(x, y, _image, false, minval);
871
872 if (!val_h.first) {
873 if (!val_v.first) { // Still no good value. Guess wildly
874 min = _bkgd + sqrt(loc.variance()) * _rand.gaussian();
875 } else {
876 min = val_v.second;
877 }
878 } else {
879 if (val_v.first) {
880 min = val_h.second;
881 } else {
882 min = (val_v.second + val_h.second) / 2;
883 }
884 }
885 }
886 /*
887 * debias the minimum; If more than one uncontaminated estimate was
888 * available, estimate the bias to be simply that due to choosing the
889 * minimum of two Gaussians. In fact, even some of the "good" pixels
890 * may have some extra charge, so even if ngood > 2, still use this
891 * estimate
892 */
893 if (ngood > 0) {
894 LOGL_DEBUG("TRACE3.lsst.algorithms.CR", "Adopted min==%g at (%d, %d) (ngood=%d)",
895 static_cast<double>(min), x, y, ngood);
896 }
897
898 if (_debias && ngood > 1) {
899 min -= interp::min2GaussianBias * sqrt(loc.variance()) * _rand.gaussian();
900 }
901
902 loc.image() = min;
903 }
904
905private:
906 MaskedImageT const &_image;
907 double _bkgd;
908 int _ncol, _nrow;
909 typename MaskedImageT::Mask::Pixel _badMask;
910 bool _debias;
911 afw::math::Random &_rand;
912};
913
914/************************************************************************************************************/
915/*
916 * actually remove CRs from the frame
917 */
918template <typename ImageT, typename MaskT>
919void removeCR(afw::image::MaskedImage<ImageT, MaskT> &mi, // image to search
920 std::vector<std::shared_ptr<afw::detection::Footprint>> &CRs, // list of cosmic rays
921 double const bkgd, // non-subtracted background
922 MaskT const, // Bit value used to label CRs
923 MaskT const saturBit, // Bit value used to label saturated pixels
924 MaskT const badMask, // Bit mask for bad pixels
925 bool const debias, // statistically debias values?
926 bool const grow // Grow CRs?
927) {
928 afw::math::Random rand; // a random number generator
929 /*
930 * replace the values of cosmic-ray contaminated pixels with 1-dim 2nd-order weighted means Cosmic-ray
931 * contaminated pixels have already been given a mask value, crBit
932 *
933 * If there are no good options (i.e. all estimates are contaminated), try using just pixels that are not
934 * CRs; failing that, interpolate in the row- or column direction over as large a distance as is required
935 *
936 * XXX SDSS (and we) go through this list backwards; why?
937 */
938
939 // a functor to remove a CR
940 RemoveCR<afw::image::MaskedImage<ImageT, MaskT>> removeCR(mi, bkgd, badMask, debias, rand);
941
942 for (std::vector<std::shared_ptr<afw::detection::Footprint>>::reverse_iterator fiter = CRs.rbegin();
943 fiter != CRs.rend(); ++fiter) {
944 std::shared_ptr<afw::detection::Footprint> cr = *fiter;
945 /*
946 * If I grow this CR does it touch saturated pixels? If so, don't
947 * interpolate and add CR pixels to saturated mask
948 */
949 /* This is commented out pending work from DM-9538
950 if (grow && cr->getArea() < 100) {
951 try {
952 bool const isotropic = false; // use a slow isotropic grow?
953 auto gcr = std::make_shared<afw::detection::Footprint>(cr->getSpans()->dilate(1),
954 cr.getRegion()); std::shared_ptr<afw::detection::Footprint> const saturPixels =
955 footprintAndMask(gcr, mi.getMask(), saturBit); auto const saturPixels = footprintAndMask(gcr,
956 mi.getMask(), saturBit);
957
958 if (saturPixels->getArea() > 0) { // pixel is adjacent to a saturation trail
959 setMaskFromFootprint(mi.getMask().get(), *saturPixels, saturBit);
960
961 continue;
962 }
963 } catch(pex::exceptions::LengthError &) {
964 continue;
965 }
966 }
967 */
968 /*
969 * OK, fix it
970 */
971 // applyFunctor must have an argument other than the functor, however in this case
972 // additional arguments are unneeded, so pass a constant 1 to be iterated over
973 // and ignore
974 cr->getSpans()->applyFunctor(removeCR, 1);
975 }
976}
977} // namespace
978
979/************************************************************************************************************/
980//
981// Explicit instantiations
982// \cond
983#define INSTANTIATE(TYPE) \
984 template std::vector<std::shared_ptr<afw::detection::Footprint>> findCosmicRays( \
985 afw::image::MaskedImage<TYPE> &image, afw::detection::Psf const &psf, double const bkgd, \
986 daf::base::PropertySet const &ps, bool const keep)
987
988INSTANTIATE(float);
989INSTANTIATE(double); // Why do we need double images?
990// \endcond
991} // namespace algorithms
992} // namespace meas
993} // namespace lsst
#define INSTANTIATE(FROMSYS, TOSYS)
Definition Detector.cc:509
#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
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:82
std::shared_ptr< math::Kernel const > getLocalKernel(lsst::geom::Point2D position, image::Color color=image::Color()) const
Return a FixedKernel corresponding to the Psf image at the given point.
Definition Psf.cc:149
virtual lsst::geom::Point2D getAveragePosition() const
Return the average position of the stars used to construct the Psf.
Definition Psf.cc:189
image::Image< Pixel > Image
Image type returned by computeImage.
Definition Psf.h:89
A range of pixels within one row of an Image.
Definition Span.h:47
Class for storing generic metadata.
Definition PropertySet.h:67
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
T empty(T... args)
T end(T... args)
T erase(T... args)
T make_shared(T... args)
T max(T... args)
T min(T... args)
T move(T... args)
bool operator<(const String &lhs, const String &rhs)
int resolve_alias(const std::vector< int > &aliases, int id)
Follow a chain of aliases, returning the final resolved value.
Definition CR.cc:96
Extent< int, 2 > ExtentI
Definition Extent.h:396
Point< int, 2 > Point2I
Definition Point.h:321
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
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)
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