LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
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"
47 #include "lsst/afw/detection/Psf.h"
49 #include "lsst/afw/math/Random.h"
52 
57 namespace lsst {
58 namespace afw {
59 namespace detection {
63 class IdSpan {
64 public:
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 };
96 int 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 
110 namespace lsst {
111 namespace meas {
112 namespace algorithms {
113 
114 namespace {
115 
116 template <typename ImageT, typename MaskT>
117 void 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 
121 template <typename ImageT>
122 bool 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
130 template <typename ImageT>
131 struct 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
147 private:
148  static int i; // current value of running index
149  int mutable _i; // running index
150 };
151 
152 template <typename ImageT>
153 int CRPixel<ImageT>::i = 0;
154 
155 template <typename ImageT>
156 struct 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  */
166 template <typename MaskedImageT>
167 bool 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 //
236 template <typename MaskedImageT>
237 void 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  */
275 template <typename ImageT>
276 class CountsInCR {
277 public:
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 
287 private:
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  */
296 template <typename ImageT>
297 static 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 
314 template <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  }
343  afw::detection::Psf::Image psfImage =
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);
611  for (std::vector<std::shared_ptr<afw::detection::Footprint>>::iterator fiter = CRs.begin();
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 /*****************************************************************************/
732 namespace {
733 /*
734  * Is condition 3 true?
735  */
736 template <typename ImageT>
737 bool 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  */
780 template <typename MaskedImageT>
781 class RemoveCR {
782 public:
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 
933 private:
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  */
946 template <typename ImageT, typename MaskT>
947 void 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 
1016 INSTANTIATE(float);
1017 INSTANTIATE(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
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