LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
CR.cc
Go to the documentation of this file.
1 // -*- LSST-C++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2008, 2009, 2010 LSST Corporation.
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 <http://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 
41 #include "boost/format.hpp"
42 
43 #include "lsst/pex/exceptions.h"
44 #include "lsst/pex/logging/Trace.h"
45 #include "lsst/pex/exceptions.h"
48 #include "lsst/afw/geom.h"
49 #include "lsst/afw/detection/Psf.h"
51 #include "lsst/afw/math/Random.h"
54 
59 namespace lsst {
60 namespace afw {
61 namespace detection {
65 class IdSpan {
66 public:
67  typedef boost::shared_ptr<IdSpan> Ptr;
68  typedef boost::shared_ptr<const IdSpan> ConstPtr;
69 
70  explicit IdSpan(int id, int y, int x0, int x1) : id(id), y(y), x0(x0), x1(x1) {}
71  int id; /* ID for object */
72  int y; /* Row wherein IdSpan dwells */
73  int x0, x1; /* inclusive range of columns */
74 };
78 struct IdSpanCompar : public std::binary_function<const IdSpan::ConstPtr, const IdSpan::ConstPtr, bool> {
80  if (a->id < b->id) {
81  return true;
82  } else if(a->id > b->id) {
83  return false;
84  } else {
85  return (a->y < b->y) ? true : false;
86  }
87  }
88 };
92 int resolve_alias(const std::vector<int>& aliases, /* list of aliases */
93  int id) { /* alias to look up */
94  int resolved = id; /* resolved alias */
95 
96  while (id != aliases[id]) {
97  resolved = id = aliases[id];
98  }
99 
100  return(resolved);
101 }
102 }}}
103 
104 namespace lsst {
105 namespace meas {
106 namespace algorithms {
107 
108 namespace geom = lsst::afw::geom;
109 namespace math = lsst::afw::math;
110 namespace image = lsst::afw::image;
111 namespace detection = lsst::afw::detection;
112 namespace pexLogging = lsst::pex::logging;
113 
114 namespace {
115 
116 template<typename ImageT, typename MaskT>
117 void removeCR(image::MaskedImage<ImageT, MaskT> & mi, std::vector<detection::Footprint::Ptr> & CRs,
118  double const bkgd, MaskT const , MaskT const saturBit, MaskT const badMask,
119  bool const debias, bool const grow);
120 
121 template<typename ImageT>
122 bool condition_3(ImageT *estimate, double const peak,
123  double const mean_ns, double const mean_we, double const mean_swne, double const mean_nwse,
124  double const dpeak,
125  double const dmean_ns, double const dmean_we,double const dmean_swne,double const dmean_nwse,
126  double const thresH, double const thresV, double const thresD,
127  double const cond3Fac
128  );
129 
130 /************************************************************************************************************/
131 //
132 // A class to hold a detected pixel
133 template<typename ImageT>
134 struct CRPixel {
135  typedef typename boost::shared_ptr<CRPixel> Ptr;
136 
137  CRPixel(int _col, int _row, ImageT _val, int _id = -1) :
138  id(_id), col(_col), row(_row), val(_val) {
139  _i = ++i;
140  }
141  ~CRPixel() {}
142 
143  bool operator< (const CRPixel& a) const {
144  return _i < a._i;
145  }
146 
147  int get_i() const {
148  return _i;
149  }
150 
151  int id; // Unique ID for cosmic ray (not cosmic ray pixel)
152  int col; // position
153  int row; // of pixel
154  ImageT val; // initial value of pixel
155 private:
156  static int i; // current value of running index
157  int mutable _i; // running index
158 };
159 
160 template<typename ImageT>
161 int CRPixel<ImageT>::i = 0;
162 
163 template<typename ImageT>
164 struct Sort_CRPixel_by_id {
165  bool operator() (CRPixel<ImageT> const & a, CRPixel<ImageT> const & b) const {
166  return a.id < b.id;
167  }
168 };
169 
170 /*****************************************************************************/
171 /*
172  * This is the code to see if a given pixel is bad
173  *
174  * Note that the pixel in question is at index 0, so its value is pt_0[0]
175  */
176 template <typename MaskedImageT>
177 bool is_cr_pixel(typename MaskedImageT::Image::Pixel *corr, // corrected value
178  typename MaskedImageT::xy_locator loc, // locator for this pixel
179  double const minSigma, // minSigma, or -threshold if negative
180  double const thresH, double const thresV, double const thresD, // for condition #3
181  double const bkgd, // unsubtracted background level
182  double const cond3Fac // fiddle factor for condition #3
183  )
184 {
185  typedef typename MaskedImageT::Image::Pixel ImagePixel;
186  //
187  // Unpack some values
188  //
189  ImagePixel const v_00 = loc.image(0, 0);
190 
191  if (v_00 < 0) {
192  return false;
193  }
194  /*
195  * condition #1 is not applied on a pixel-by-pixel basis
196  */
197  ;
198  /*
199  * condition #2
200  */
201  ImagePixel const mean_we = (loc.image(-1, 0) + loc.image( 1, 0))/2; // avgs of surrounding 8 pixels
202  ImagePixel const mean_ns = (loc.image( 0, 1) + loc.image( 0, -1))/2;
203  ImagePixel const mean_swne = (loc.image(-1, -1) + loc.image( 1, 1))/2;
204  ImagePixel const mean_nwse = (loc.image(-1, 1) + loc.image( 1, -1))/2;
205 
206  if (minSigma < 0) { /* |thres_sky_sigma| is threshold */
207  if (v_00 < -minSigma) {
208  return false;
209  }
210  } else {
211  double const thres_sky_sigma = minSigma*sqrt(loc.variance(0, 0));
212 
213  if (v_00 < mean_ns + thres_sky_sigma &&
214  v_00 < mean_we + thres_sky_sigma &&
215  v_00 < mean_swne + thres_sky_sigma &&
216  v_00 < mean_nwse + thres_sky_sigma) {
217  return false;
218  }
219  }
220 /*
221  * condition #3
222  *
223  * Note that this uses mean_ns etc. even if minSigma is negative
224  */
225  double const dv_00 = sqrt(loc.variance( 0, 0));
226  // standard deviation of means of surrounding pixels
227  double const dmean_we = sqrt(loc.variance(-1, 0) + loc.variance( 1, 0))/2;
228  double const dmean_ns = sqrt(loc.variance( 0, 1) + loc.variance( 0, -1))/2;
229  double const dmean_swne = sqrt(loc.variance(-1, -1) + loc.variance( 1, 1))/2;
230  double const dmean_nwse = sqrt(loc.variance(-1, 1) + loc.variance( 1, -1))/2;
231 
232  if (!condition_3(corr,
233  v_00 - bkgd, mean_ns - bkgd, mean_we - bkgd, mean_swne - bkgd, mean_nwse - bkgd,
234  dv_00, dmean_ns, dmean_we, dmean_swne, dmean_nwse,
235  thresH, thresV, thresD, cond3Fac)){
236  return false;
237  }
238 /*
239  * OK, it's a contaminated pixel
240  */
241  *corr += static_cast<ImagePixel>(bkgd);
242 
243  return true;
244 }
245 
246 /************************************************************************************************************/
247 //
248 // Worker routine to process the pixels adjacent to a span (including the points just
249 // to the left and just to the right)
250 //
251 template <typename MaskedImageT>
252 void checkSpanForCRs(detection::Footprint *extras, // Extra spans get added to this Footprint
253  std::vector<CRPixel<typename MaskedImageT::Image::Pixel> >& crpixels,
254  // a list of pixels containing CRs
255  int const y, // the row to process
256  int const x0, int const x1, // range of pixels in the span (inclusive)
257  MaskedImageT& image,
258  double const minSigma, // minSigma
259  double const thresH, double const thresV, double const thresD, // for cond. #3
260  double const bkgd, // unsubtracted background level
261  double const cond3Fac, // fiddle factor for condition #3
262  bool const keep // if true, don't remove the CRs
263  )
264 {
265  typedef typename MaskedImageT::Image::Pixel MImagePixel;
266  typename MaskedImageT::xy_locator loc = image.xy_at(x0 - 1, y); // locator for data
267 
268  int const imageX0 = image.getX0();
269  int const imageY0 = image.getY0();
270 
271  for (int x = x0 - 1; x <= x1 + 1; ++x) {
272  MImagePixel corr = 0; // new value for pixel
273  if (is_cr_pixel<MaskedImageT>(&corr, loc, minSigma, thresH, thresV, thresD,
274  bkgd, cond3Fac)) {
275  if (keep) {
276  crpixels.push_back(CRPixel<MImagePixel>(x + imageX0, y + imageY0, loc.image()));
277  }
278  loc.image() = corr;
279 
280  extras->addSpan(y + imageY0, x + imageX0, x + imageX0);
281  }
282  ++loc.x();
283  }
284 }
285 
286 /************************************************************************************************************/
287 /*
288  * Find the sum of the pixels in a Footprint
289  */
290 template <typename ImageT>
291 class CountsInCR : public detection::FootprintFunctor<ImageT> {
292 public:
293  CountsInCR(ImageT const& mimage, double const bkgd) :
294  detection::FootprintFunctor<ImageT>(mimage),
295  _bkgd(bkgd),
296  _sum(0.0) {}
297 
298  // method called for each pixel by apply()
299  void operator()(typename ImageT::xy_locator loc, // locator pointing at the pixel
300  int, int
301  ) {
302  _sum += *loc - _bkgd;
303  }
304 
305  virtual void reset(detection::Footprint const&) {}
306  virtual void reset() {
307  _sum = 0.0;
308  }
309 
310  double getCounts() const { return _sum; }
311 private:
312  double const _bkgd; // the Image's background level
313  typename ImageT::Pixel _sum; // the sum of all DN in the Footprint, corrected for bkgd
314 };
315 }
316 
317 /*
318  * Restore all the pixels in crpixels to their pristine state
319  */
320 template <typename ImageT>
321 static void reinstateCrPixels(
322  ImageT *image, // the image in question
323  std::vector<CRPixel<typename ImageT::Pixel> > const& crpixels // a list of pixels with CRs
324  )
325 {
326  if (crpixels.empty()) return;
327 
328  typedef typename std::vector<CRPixel<typename ImageT::Pixel> >::const_iterator crpixel_iter;
329  for (crpixel_iter crp = crpixels.begin(), end = crpixels.end(); crp < end - 1 ; ++crp) {
330  *image->at(crp->col - image->getX0(), crp->row - image->getY0()) = crp->val;
331  }
332 }
333 
339 template <typename MaskedImageT>
340 std::vector<detection::Footprint::Ptr>
341 findCosmicRays(MaskedImageT &mimage,
342  detection::Psf const &psf,
343  double const bkgd,
344  lsst::pex::policy::Policy const &policy,
345  bool const keep
346  ) {
347  typedef typename MaskedImageT::Image ImageT;
348  typedef typename ImageT::Pixel ImagePixel;
349  typedef typename MaskedImageT::Mask::Pixel MaskPixel;
350 
351  // Parse the Policy
352  double const minSigma = policy.getDouble("minSigma"); // min sigma over sky in pixel for CR candidate
353  double const minDn = policy.getDouble("min_DN"); // min number of DN in an CRs
354  double const cond3Fac = policy.getDouble("cond3_fac"); // fiddle factor for condition #3
355  double const cond3Fac2 = policy.getDouble("cond3_fac2"); // 2nd fiddle factor for condition #3
356  int const niteration = policy.getInt("niteration"); // Number of times to look for contaminated
357  // pixels near CRs
358  int const nCrPixelMax = policy.getInt("nCrPixelMax"); // maximum number of contaminated pixels
359 /*
360  * thresholds for 3rd condition
361  *
362  * Realise PSF at center of image
363  */
365  afw::geom::Point2D(mimage.getWidth() / 2.0, mimage.getHeight() / 2.0)
366  );
367  if (!kernel) {
368  throw LSST_EXCEPT(pexExcept::NotFoundError, "Psf is unable to return a kernel");
369  }
370  detection::Psf::Image psfImage = detection::Psf::Image(geom::ExtentI(kernel->getWidth(), kernel->getHeight()));
371  kernel->computeImage(psfImage, true);
372 
373  int const xc = kernel->getCtrX(); // center of PSF
374  int const yc = kernel->getCtrY();
375 
376  double const I0 = psfImage(xc, yc);
377  double const thresH = cond3Fac2*(0.5*(psfImage(xc - 1, yc) + psfImage(xc + 1, yc)))/I0; // horizontal
378  double const thresV = cond3Fac2*(0.5*(psfImage(xc, yc - 1) + psfImage(xc, yc + 1)))/I0; // vertical
379  double const thresD = cond3Fac2*(0.25*(psfImage(xc - 1, yc - 1) + psfImage(xc + 1, yc + 1) +
380  psfImage(xc - 1, yc + 1) + psfImage(xc + 1, yc - 1)))/I0; // diag
381 /*
382  * Setup desired mask planes
383  */
384  MaskPixel const badBit = mimage.getMask()->getPlaneBitMask("BAD"); // Generic bad pixels
385  MaskPixel const crBit = mimage.getMask()->getPlaneBitMask("CR"); // CR-contaminated pixels
386  MaskPixel const interpBit = mimage.getMask()->getPlaneBitMask("INTRP"); // Interpolated pixels
387  MaskPixel const saturBit = mimage.getMask()->getPlaneBitMask("SAT"); // Saturated pixels
388 
389  MaskPixel const badMask = (badBit | interpBit | saturBit); // naughty pixels
390 /*
391  * Go through the frame looking at each pixel (except the edge ones which we ignore)
392  */
393  int const ncol = mimage.getWidth();
394  int const nrow = mimage.getHeight();
395 
396  std::vector<CRPixel<ImagePixel> > crpixels; // storage for detected CR-contaminated pixels
397  typedef typename std::vector<CRPixel<ImagePixel> >::iterator crpixel_iter;
398  typedef typename std::vector<CRPixel<ImagePixel> >::reverse_iterator crpixel_riter;
399 
400  for (int j = 1; j < nrow - 1; ++j) {
401  typename MaskedImageT::xy_locator loc = mimage.xy_at(1, j); // locator for data
402 
403  for (int i = 1; i < ncol - 1; ++i, ++loc.x()) {
404  ImagePixel corr = 0;
405  if (!is_cr_pixel<MaskedImageT>(&corr, loc, minSigma,
406  thresH, thresV, thresD, bkgd, cond3Fac)) {
407  continue;
408  }
409 /*
410  * condition #4
411  */
412  if (loc.mask() & badMask) {
413  continue;
414  }
415  if ((loc.mask(-1, 1) | loc.mask(0, 1) | loc.mask(1, 1) |
416  loc.mask(-1, 0) | loc.mask(1, 0) |
417  loc.mask(-1, -1) | loc.mask(0, -1) | loc.mask(1, -1)) & interpBit) {
418  continue;
419  }
420 /*
421  * OK, it's a CR
422  *
423  * replace CR-contaminated pixels with reasonable values as we go through
424  * image, which increases the detection rate
425  */
426  crpixels.push_back(CRPixel<ImagePixel>(i + mimage.getX0(), j + mimage.getY0(), loc.image()));
427  loc.image() = corr; /* just a preliminary estimate */
428 
429  if (static_cast<int>(crpixels.size()) > nCrPixelMax) {
430  reinstateCrPixels(mimage.getImage().get(), crpixels);
431 
432  throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
433  (boost::format("Too many CR pixels (max %d)") % nCrPixelMax).str());
434  }
435  }
436  }
437 /*
438  * We've found them on a pixel-by-pixel basis, now merge those pixels
439  * into cosmic rays
440  */
441  std::vector<int> aliases; // aliases for initially disjoint parts of CRs
442  aliases.reserve(1 + crpixels.size()/2); // initial size of aliases
443 
444  std::vector<detection::IdSpan::Ptr> spans; // y:x0,x1 for objects
445  spans.reserve(aliases.capacity()); // initial size of spans
446 
447  aliases.push_back(0); // 0 --> 0
448 
454  int ncr = 0; // number of detected cosmic rays
455  if (!crpixels.empty()) {
456  int id; // id number for a CR
457  int x0 = -1, x1 = -1, y = -1; // the beginning and end column, and row of this span in a CR
458 
459  // I am dummy
460  CRPixel<ImagePixel> dummy(0, -1, 0, -1);
461  crpixels.push_back(dummy);
462  //printf("Created dummy CR: i %i, id %i, col %i, row %i, val %g\n", dummy.get_i(), dummy.id, dummy.col, dummy.row, (double)dummy.val);
463  for (crpixel_iter crp = crpixels.begin(); crp < crpixels.end() - 1 ; ++crp) {
464  //printf("Looking at CR: i %i, id %i, col %i, row %i, val %g\n", crp->get_i(), crp->id, crp->col, crp->row, (double)crp->val);
465 
466  if (crp->id < 0) { // not already assigned
467  crp->id = ++ncr; // a new CR
468  aliases.push_back(crp->id);
469  y = crp->row;
470  x0 = x1 = crp->col;
471  //printf(" Assigned ID %i; looking at row %i, start col %i\n", crp->id, crp->row, crp->col);
472  }
473  id = crp->id;
474  //printf(" Next CRpix has i=%i, id=%i, row %i, col %i\n", crp[1].get_i(), crp[1].id, crp[1].row, crp[1].col);
475 
476  if (crp[1].row == crp[0].row && crp[1].col == crp[0].col + 1) {
477  //printf(" Adjoining! Set next CRpix id = %i; x1=%i\n", crp[1].id, x1);
478  crp[1].id = id;
479  ++x1;
480  } else {
481  assert (y >= 0 && x0 >= 0 && x1 >= 0);
482  spans.push_back(detection::IdSpan::Ptr(new detection::IdSpan(id, y, x0, x1)));
483  //printf(" Not adjoining; adding span id=%i, y=%i, x = [%i, %i]\n", id, y, x0, x1);
484  }
485  }
486  }
487 
488  // At the end of this loop, all crpixel entries have been assigned an ID,
489  // except for the "dummy" entry at the end of the array.
490  if (crpixels.size() > 0) {
491  for (crpixel_iter cp = crpixels.begin(); cp != crpixels.end() - 1; cp++) {
492  assert(cp->id >= 0);
493  assert(cp->col >= 0);
494  assert(cp->row >= 0);
495  }
496  // dummy:
497  assert(crpixels[crpixels.size()-1].id == -1);
498  assert(crpixels[crpixels.size()-1].col == 0);
499  assert(crpixels[crpixels.size()-1].row == -1);
500  }
501 
502  for (std::vector<detection::IdSpan::Ptr>::iterator sp = spans.begin(), end = spans.end(); sp != end; sp++) {
503  assert((*sp)->id >= 0);
504  assert((*sp)->y >= 0);
505  assert((*sp)->x0 >= 0);
506  assert((*sp)->x1 >= (*sp)->x0);
507  for (std::vector<detection::IdSpan::Ptr>::iterator sp2 = sp + 1; sp2 != end; sp2++) {
508  assert((*sp2)->y >= (*sp)->y);
509  if ((*sp2)->y == (*sp)->y) {
510  assert((*sp2)->x0 > (*sp)->x1);
511  }
512  }
513  }
514 
515 /*
516  * See if spans touch each other
517  */
518  for (std::vector<detection::IdSpan::Ptr>::iterator sp = spans.begin(), end = spans.end();
519  sp != end; ++sp) {
520  int const y = (*sp)->y;
521  int const x0 = (*sp)->x0;
522  int const x1 = (*sp)->x1;
523 
524  // this loop will probably run for only a few steps
525  for (std::vector<detection::IdSpan::Ptr>::iterator sp2 = sp + 1; sp2 != end; ++sp2) {
526  if ((*sp2)->y == y) {
527  // on this row (but not adjoining columns, since it would have been merged into this span);
528  // keep looking.
529  continue;
530  } else if ((*sp2)->y != (y + 1)) {
531  // sp2 is more than one row below; can't be connected.
532  break;
533  } else if ((*sp2)->x0 > (x1 + 1)) {
534  // sp2 is more than one column away to the right; can't be connected
535  break;
536  } else if ((*sp2)->x1 >= (x0 - 1)) {
537  // touches
538  int r1 = detection::resolve_alias(aliases, (*sp)->id);
539  int r2 = detection::resolve_alias(aliases, (*sp2)->id);
540  aliases[r1] = r2;
541  }
542  }
543  }
544 
545 
546 
547 
548 
549 /*
550  * Resolve aliases; first alias chains, then the IDs in the spans
551  */
552  for (unsigned int i = 0; i != spans.size(); ++i) {
553  spans[i]->id = detection::resolve_alias(aliases, spans[i]->id);
554  }
555 
556 /*
557  * Sort spans by ID, so we can sweep through them once
558  */
559  if (spans.size() > 0) {
560  std::sort(spans.begin(), spans.end(), detection::IdSpanCompar());
561  }
562 
563 /*
564  * Build Footprints from spans
565  */
566  std::vector<detection::Footprint::Ptr> CRs; // our cosmic rays
567 
568  if (spans.size() > 0) {
569  int id = spans[0]->id;
570  unsigned int i0 = 0; // initial value of i
571  for (unsigned int i = i0; i <= spans.size(); ++i) { // <= size to catch the last object
572  if (i == spans.size() || spans[i]->id != id) {
574 
575  for (; i0 < i; ++i0) {
576  cr->addSpan(spans[i0]->y, spans[i0]->x0, spans[i0]->x1);
577  }
578  CRs.push_back(cr);
579  }
580 
581  if (i < spans.size()) {
582  id = spans[i]->id;
583  }
584  }
585  }
586 
587  reinstateCrPixels(mimage.getImage().get(), crpixels);
588 /*
589  * apply condition #1
590  */
591  CountsInCR<ImageT> CountDN(*mimage.getImage(), bkgd);
592  for (std::vector<detection::Footprint::Ptr>::iterator cr = CRs.begin(), end = CRs.end();
593  cr != end; ++cr) {
594  CountDN.apply(**cr); // find the sum of pixel values within the CR
595 
596  pexLogging::TTrace<10>("algorithms.CR", "CR at (%d, %d) has %g DN",
597  (*cr)->getBBox().getMinX(), (*cr)->getBBox().getMinY(), CountDN.getCounts());
598  if (CountDN.getCounts() < minDn) { /* not bright enough */
599  pexLogging::TTrace<11>("algorithms.CR", "Erasing CR");
600 
601  cr = CRs.erase(cr);
602  --cr; // back up to previous CR (we're going to increment it)
603  --end;
604  }
605  }
606  ncr = CRs.size(); /* some may have been too faint */
607 /*
608  * We've found them all, time to kill them all
609  */
610  bool const debias_values = true;
611  bool grow = false;
612  pexLogging::TTrace<2>("algorithms.CR", "Removing initial list of CRs");
613  removeCR(mimage, CRs, bkgd, crBit, saturBit, badMask, debias_values, grow);
614 #if 0 // Useful to see phase 2 in ds9; debugging only
615  (void)setMaskFromFootprintList(mimage.getMask().get(), CRs,
616  mimage.getMask()->getPlaneBitMask("DETECTED"));
617 #endif
618 /*
619  * Now that we've removed them, go through image again, examining area around
620  * each CR for extra bad pixels. Note that we set cond3Fac = 0 for this pass
621  *
622  * We iterate niteration times; niter==1 was sufficient for SDSS data, but megacam
623  * CCDs are different -- who knows for other devices?
624  */
625  bool too_many_crs = false; // we've seen too many CR pixels
626  int nextra = 0; // number of pixels added to list of CRs
627  for (int i = 0; i != niteration && !too_many_crs; ++i) {
628  pexLogging::TTrace<1>("algorithms.CR", "Starting iteration %d", i);
629  for (std::vector<detection::Footprint::Ptr>::iterator fiter = CRs.begin();
630  fiter != CRs.end(); fiter++) {
631  detection::Footprint::Ptr cr = *fiter;
632 /*
633  * Are all those `CR' pixels interpolated? If so, don't grow it
634  */
635  {
636  detection::Footprint::Ptr om = footprintAndMask(cr, mimage.getMask(), interpBit);
637  int const npix = (om) ? om->getNpix() : 0;
638 
639  if (npix == cr->getNpix()) {
640  continue;
641  }
642  }
643 /*
644  * No; some of the suspect pixels aren't interpolated
645  */
646  detection::Footprint extra; // extra pixels added to cr
647  for (detection::Footprint::SpanList::const_iterator siter = cr->getSpans().begin();
648  siter != cr->getSpans().end(); siter++) {
649  detection::Span::Ptr const span = *siter;
650 
651  /*
652  * Check the lines above and below the span. We're going to check a 3x3 region around
653  * the pixels, so we need a buffer around the edge. We check the pixels just to the
654  * left/right of the span, so the buffer needs to be 2 pixels (not just 1) in the
655  * column direction, but only 1 in the row direction.
656  */
657  int const y = span->getY() - mimage.getY0();
658  if (y < 2 || y >= nrow - 2) {
659  continue;
660  }
661  int x0 = span->getX0() - mimage.getX0();
662  int x1 = span->getX1() - mimage.getX0();
663  x0 = (x0 < 2) ? 2 : (x0 > ncol - 3) ? ncol - 3 : x0;
664  x1 = (x1 < 2) ? 2 : (x1 > ncol - 3) ? ncol - 3 : x1;
665 
666  checkSpanForCRs(&extra, crpixels, y - 1, x0, x1, mimage,
667  minSigma/2, thresH, thresV, thresD, bkgd, 0, keep);
668  checkSpanForCRs(&extra, crpixels, y, x0, x1, mimage,
669  minSigma/2, thresH, thresV, thresD, bkgd, 0, keep);
670  checkSpanForCRs(&extra, crpixels, y + 1, x0, x1, mimage,
671  minSigma/2, thresH, thresV, thresD, bkgd, 0, keep);
672  }
673 
674  if (extra.getSpans().size() > 0) { // we added some pixels
675  if (nextra + static_cast<int>(crpixels.size()) > nCrPixelMax) {
676  too_many_crs = true;
677  break;
678  }
679 
680  nextra += extra.getNpix();
681 
682  detection::Footprint::SpanList &espans = extra.getSpans();
683  for (detection::Footprint::SpanList::const_iterator siter = espans.begin();
684  siter != espans.end(); siter++) {
685  cr->addSpan(**siter);
686  }
687  cr->normalize();
688  }
689  }
690 
691  if (nextra == 0) {
692  break;
693  }
694  }
695 /*
696  * mark those pixels as CRs
697  */
698  if (!too_many_crs) {
699  (void)setMaskFromFootprintList(mimage.getMask().get(), CRs, crBit);
700  }
701 /*
702  * Maybe reinstate initial values; n.b. the same pixel may appear twice, so we want the
703  * first value stored (hence the uses of rbegin/rend)
704  *
705  * We have to do this if we decide _not_ to remove certain CRs,
706  * for example those which lie next to saturated pixels
707  */
708  if (keep || too_many_crs) {
709  if (crpixels.size() > 0) {
710  int const imageX0 = mimage.getX0();
711  int const imageY0 = mimage.getY0();
712 
713  std::sort(crpixels.begin(), crpixels.end()); // sort into birth order
714 
715  crpixel_riter rend = crpixels.rend();
716  for (crpixel_riter crp = crpixels.rbegin(); crp != rend; ++crp) {
717  if (crp->row == -1)
718  // dummy; skip it.
719  continue;
720  mimage.at(crp->col - imageX0, crp->row - imageY0).image() = crp->val;
721  }
722  }
723  } else {
724  if (true || nextra > 0) {
725  grow = true;
726  pexLogging::TTrace<2>("algorithms.CR", "Removing final list of CRs, grow = %d", grow);
727  removeCR(mimage, CRs, bkgd, crBit, saturBit, badMask, debias_values, grow);
728  }
729 /*
730  * we interpolated over all CR pixels, so set the interp bits too
731  */
732  (void)setMaskFromFootprintList(mimage.getMask().get(), CRs,
733  static_cast<MaskPixel>(crBit | interpBit));
734  }
735 
736  if (too_many_crs) { // we've cleaned up, so we can throw the exception
737  throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
738  (boost::format("Too many CR pixels (max %d)") % nCrPixelMax).str());
739  }
740 
741  return CRs;
742 }
743 
744 /*****************************************************************************/
745 namespace {
746 /*
747  * Is condition 3 true?
748  */
749 template<typename ImageT>
750 bool condition_3(ImageT *estimate, // estimate of true value of pixel
751  double const peak, // counts in central pixel (no sky)
752  double const mean_ns, // mean in NS direction (no sky)
753  double const mean_we, // " " WE " " " "
754  double const mean_swne, // " " SW-NE " " " "
755  double const mean_nwse, // " " NW-SE " " " "
756  double const dpeak, // standard deviation of peak value
757  double const dmean_ns, // s.d. of mean in NS direction
758  double const dmean_we, // " " " " WE " "
759  double const dmean_swne, // " " " " SW-NE " "
760  double const dmean_nwse, // " " " " NW-SE " "
761  double const thresH, // horizontal threshold
762  double const thresV, // vertical threshold
763  double const thresD, // diagonal threshold
764  double const cond3Fac // fiddle factor for noise
765  )
766 {
767  if (thresV*(peak - cond3Fac*dpeak) > mean_ns + cond3Fac*dmean_ns) {
768  *estimate = (ImageT)mean_ns;
769  return true;
770  }
771 
772  if (thresH*(peak - cond3Fac*dpeak) > mean_we + cond3Fac*dmean_we) {
773  *estimate = mean_we;
774  return true;
775  }
776 
777  if (thresD*(peak - cond3Fac*dpeak) > mean_swne + cond3Fac*dmean_swne) {
778  *estimate = mean_swne;
779  return true;
780  }
781 
782  if (thresD*(peak - cond3Fac*dpeak) > mean_nwse + cond3Fac*dmean_nwse) {
783  *estimate = mean_nwse;
784  return true;
785  }
786 
787  return false;
788 }
789 
790 /************************************************************************************************************/
791 /*
792  * Interpolate over a CR's pixels
793  */
794 template <typename MaskedImageT>
795 class RemoveCR : public detection::FootprintFunctor<MaskedImageT> {
796 public:
797  RemoveCR(MaskedImageT const& mimage,
798  double const bkgd,
799  typename MaskedImageT::Mask::Pixel badMask,
800  bool const debias,
802  ) : detection::FootprintFunctor<MaskedImageT>(mimage),
803  _bkgd(bkgd),
804  _ncol(mimage.getWidth()),
805  _nrow(mimage.getHeight()),
806  _badMask(badMask),
807  _debias(debias),
808  _rand(rand) {}
809 
810  // method called for each pixel by apply()
811  void operator()(typename MaskedImageT::xy_locator loc, // locator pointing at the pixel
812  int x, // column-position of pixel
813  int y // row-position of pixel
814  ) {
815  typedef typename MaskedImageT::Image::Pixel MImagePixel;
816  MImagePixel min = std::numeric_limits<MImagePixel>::max();
817  int ngood = 0; // number of good values on min
818 
819  MImagePixel const minval = _bkgd - 2*sqrt(loc.variance()); // min. acceptable pixel value after interp
820 /*
821  * W-E row
822  */
823  if (x - 2 >= 0 && x + 2 < _ncol) {
824  if ((loc.mask(-2, 0) | _badMask) || (loc.mask(-1, 0) | _badMask) ||
825  (loc.mask( 1, 0) | _badMask) || (loc.mask( 2, 0) | _badMask)) {
826  ; // estimate is contaminated
827  } else {
828  MImagePixel const v_m2 = loc.image(-2, 0);
829  MImagePixel const v_m1 = loc.image(-1, 0);
830  MImagePixel const v_p1 = loc.image( 1, 0);
831  MImagePixel const v_p2 = loc.image( 2, 0);
832 
833  MImagePixel const tmp =
834  interp::lpc_1_c1*(v_m1 + v_p1) + interp::lpc_1_c2*(v_m2 + v_p2);
835 
836  if (tmp > minval && tmp < min) {
837  min = tmp;
838  ngood++;
839  }
840  }
841  }
842 /*
843  * N-S column
844  */
845  if (y - 2 >= 0 && y + 2 < _nrow) {
846  if ((loc.mask(0, -2) | _badMask) || (loc.mask(0, -1) | _badMask) ||
847  (loc.mask(0, 1) | _badMask) || (loc.mask(0, 2) | _badMask)) {
848  ; /* estimate is contaminated */
849  } else {
850  MImagePixel const v_m2 = loc.image(0, -2);
851  MImagePixel const v_m1 = loc.image(0, -1);
852  MImagePixel const v_p1 = loc.image(0, 1);
853  MImagePixel const v_p2 = loc.image(0, 2);
854 
855  MImagePixel const tmp =
856  interp::lpc_1_c1*(v_m1 + v_p1) + interp::lpc_1_c2*(v_m2 + v_p2);
857 
858  if (tmp > minval && tmp < min) {
859  min = tmp;
860  ngood++;
861  }
862  }
863  }
864 /*
865  * SW--NE diagonal
866  */
867  if (x - 2 >= 0 && x + 2 < _ncol && y - 2 >= 0 && y + 2 < _nrow) {
868  if ((loc.mask(-2, -2) | _badMask) || (loc.mask(-1, -1) | _badMask) ||
869  (loc.mask( 1, 1) | _badMask) || (loc.mask( 2, 2) | _badMask)) {
870  ; /* estimate is contaminated */
871  } else {
872  MImagePixel const v_m2 = loc.image(-2, -2);
873  MImagePixel const v_m1 = loc.image(-1, -1);
874  MImagePixel const v_p1 = loc.image( 1, 1);
875  MImagePixel const v_p2 = loc.image( 2, 2);
876 
877  MImagePixel const tmp =
878  interp::lpc_1s2_c1*(v_m1 + v_p1) + interp::lpc_1s2_c2*(v_m2 + v_p2);
879 
880  if (tmp > minval && tmp < min) {
881  min = tmp;
882  ngood++;
883  }
884  }
885  }
886 /*
887  * SE--NW diagonal
888  */
889  if (x - 2 >= 0 && x + 2 < _ncol && y - 2 >= 0 && y + 2 < _nrow) {
890  if ((loc.mask( 2, -2) | _badMask) || (loc.mask( 1, -1) | _badMask) ||
891  (loc.mask(-1, 1) | _badMask) || (loc.mask(-2, 2) | _badMask)) {
892  ; /* estimate is contaminated */
893  } else {
894  MImagePixel const v_m2 = loc.image( 2, -2);
895  MImagePixel const v_m1 = loc.image( 1, -1);
896  MImagePixel const v_p1 = loc.image(-1, 1);
897  MImagePixel const v_p2 = loc.image(-2, 2);
898 
899  MImagePixel const tmp =
900  interp::lpc_1s2_c1*(v_m1 + v_p1) + interp::lpc_1s2_c2*(v_m2 + v_p2);
901 
902  if (tmp > minval && tmp < min) {
903  min = tmp;
904  ngood++;
905  }
906  }
907  }
908 /*
909  * Have we altogether failed to find an acceptable value? If so interpolate
910  * using the full-up interpolation code both vertically and horizontally
911  * and take the average. This can fail for large enough defects (e.g. CRs
912  * lying in bad columns), in which case the interpolator returns -1. If
913  * both directions fail, use the background value.
914  */
915  if (ngood == 0) {
916  std::pair<bool, MImagePixel const> val_h =
917  interp::singlePixel(x, y, this->getImage(), true, minval);
918  std::pair<bool, MImagePixel const> val_v =
919  interp::singlePixel(x, y, this->getImage(), false, minval);
920 
921  if (!val_h.first) {
922  if (!val_v.first) { // Still no good value. Guess wildly
923  min = _bkgd + sqrt(loc.variance())*_rand.gaussian();
924  } else {
925  min = val_v.second;
926  }
927  } else {
928  if (val_v.first) {
929  min = val_h.second;
930  } else {
931  min = (val_v.second + val_h.second)/2;
932  }
933  }
934  }
935 /*
936  * debias the minimum; If more than one uncontaminated estimate was
937  * available, estimate the bias to be simply that due to choosing the
938  * minimum of two Gaussians. In fact, even some of the "good" pixels
939  * may have some extra charge, so even if ngood > 2, still use this
940  * estimate
941  */
942  if (ngood > 0) {
943  pexLogging::TTrace<5>("algorithms.CR", "Adopted min==%g at (%d, %d) (ngood=%d)",
944  static_cast<double>(min), x, y, ngood);
945  }
946 
947  if (_debias && ngood > 1) {
948  min -= interp::min2GaussianBias*sqrt(loc.variance())*_rand.gaussian();
949  }
950 
951  loc.image() = min;
952  }
953 private:
954  double _bkgd;
955  int _ncol, _nrow;
956  typename MaskedImageT::Mask::Pixel _badMask;
957  bool _debias;
959 };
960 
961 /************************************************************************************************************/
962 /*
963  * actually remove CRs from the frame
964  */
965 template<typename ImageT, typename MaskT>
966 void removeCR(image::MaskedImage<ImageT, MaskT> & mi, // image to search
967  std::vector<detection::Footprint::Ptr> & CRs, // list of cosmic rays
968  double const bkgd, // non-subtracted background
969  MaskT const , // Bit value used to label CRs
970  MaskT const saturBit, // Bit value used to label saturated pixels
971  MaskT const badMask, // Bit mask for bad pixels
972  bool const debias, // statistically debias values?
973  bool const grow // Grow CRs?
974  )
975 {
976  lsst::afw::math::Random rand; // a random number generator
977  /*
978  * replace the values of cosmic-ray contaminated pixels with 1-dim 2nd-order weighted means Cosmic-ray
979  * contaminated pixels have already been given a mask value, crBit
980  *
981  * If there are no good options (i.e. all estimates are contaminated), try using just pixels that are not
982  * CRs; failing that, interpolate in the row- or column direction over as large a distance as is required
983  *
984  * XXX SDSS (and we) go through this list backwards; why?
985  */
986 
987  // a functor to remove a CR
988  RemoveCR<image::MaskedImage<ImageT, MaskT> > removeCR(mi, bkgd, badMask, debias, rand);
989 
990  for (std::vector<detection::Footprint::Ptr>::reverse_iterator fiter = CRs.rbegin();
991  fiter != CRs.rend(); ++fiter) {
992  detection::Footprint::Ptr cr = *fiter;
993 /*
994  * If I grow this CR does it touch saturated pixels? If so, don't
995  * interpolate and add CR pixels to saturated mask
996  */
997  if (grow && cr->getNpix() < 100) {
998  try {
999  bool const isotropic = false; // use a slow isotropic grow?
1000  detection::Footprint::Ptr gcr = growFootprint(cr, 1, isotropic);
1001  detection::Footprint::Ptr const saturPixels = footprintAndMask(gcr, mi.getMask(), saturBit);
1002 
1003  if (saturPixels->getNpix() > 0) { // pixel is adjacent to a saturation trail
1004  setMaskFromFootprint(mi.getMask().get(), *saturPixels, saturBit);
1005 
1006  continue;
1007  }
1008  } catch(lsst::pex::exceptions::LengthError &) {
1009  continue;
1010  }
1011  }
1012 /*
1013  * OK, fix it
1014  */
1015  removeCR.apply(*cr);
1016  }
1017 }
1018 }
1019 
1020 /************************************************************************************************************/
1021 //
1022 // Explicit instantiations
1023 // \cond
1024 #define INSTANTIATE(TYPE) \
1025  template \
1026  std::vector<detection::Footprint::Ptr> \
1027  findCosmicRays(lsst::afw::image::MaskedImage<TYPE> &image, \
1028  detection::Psf const &psf, \
1029  double const bkgd, \
1030  lsst::pex::policy::Policy const& policy, \
1031  bool const keep \
1032  )
1033 
1034 INSTANTIATE(float);
1035 INSTANTIATE(double); // Why do we need double images?
1036 // \endcond
1037 }}}
int y
boost::shared_ptr< const IdSpan > ConstPtr
Definition: CR.cc:68
boost::uint16_t MaskPixel
An include file to include the header files for lsst::afw::geom.
boost::shared_ptr< math::Kernel const > getLocalKernel(geom::Point2D position=makeNullPoint(), image::Color color=image::Color()) const
Return a FixedKernel corresponding to the Psf image at the given point.
std::vector< boost::shared_ptr< lsst::afw::detection::Footprint > > findCosmicRays(MaskedImageT &image, lsst::afw::detection::Psf const &psf, double const bkgd, lsst::pex::policy::Policy const &policy, bool const keep=false)
Find cosmic rays in an Image, and mask and remove them.
Definition: CR.cc:341
double const min2GaussianBias
Mean value of the minimum of two N(0,1) variates.
Definition: Interp.h:60
Represent a set of pixels of an arbitrary shape and size.
boost::shared_ptr< Kernel const > ConstPtr
Definition: Kernel.h:142
int _i
Definition: CR.cc:157
a container for holding hierarchical configuration data in memory.
Definition: Policy.h:169
boost::shared_ptr< Footprint > footprintAndMask(boost::shared_ptr< Footprint > const &foot, typename image::Mask< MaskT >::Ptr const &mask, MaskT const bitmask)
Return a Footprint that&#39;s the intersection of a Footprint with a Mask.
definition of the Trace messaging facilities
double getDouble(const std::string &name) const
Definition: Policy.h:617
boost::shared_ptr< Footprint > Ptr
Definition: Footprint.h:67
int getInt(const std::string &name) const
Definition: Policy.h:603
int const x0
Definition: saturated.cc:45
boost::shared_ptr< IdSpan > Ptr
Definition: CR.cc:67
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
boost::shared_ptr< Footprint > growFootprint(Footprint const &foot, int nGrow, bool isotropic=true)
ImageT::Pixel _sum
Definition: CR.cc:313
bool operator()(const IdSpan::ConstPtr a, const IdSpan::ConstPtr b)
Definition: CR.cc:79
std::vector< Span::Ptr > SpanList
The Footprint&#39;s Span list.
Definition: Footprint.h:71
int _id
Definition: Mask.cc:679
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:77
IdSpan(int id, int y, int x0, int x1)
Definition: CR.cc:70
int _nrow
Definition: CR.cc:955
A set of pixels in an Image.
Definition: Footprint.h:62
int x
int row
Definition: CR.cc:153
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
#define INSTANTIATE(T)
lsst::afw::math::Random & _rand
Definition: CR.cc:958
MaskPtr getMask(bool const noThrow=false) const
Return a (Ptr to) the MaskedImage&#39;s mask.
Definition: MaskedImage.h:879
const Span & addSpan(const int y, const int x0, const int x1)
std::pair< bool, typename MaskedImageT::Image::Pixel > singlePixel(int x, int y, MaskedImageT const &image, bool horizontal, double minval)
Definition: Interp.cc:2134
Random number generator class.
MaskT setMaskFromFootprint(lsst::afw::image::Mask< MaskT > *mask, Footprint const &footprint, MaskT const bitmask)
OR bitmask into all the Mask&#39;s pixels that are in the Footprint.
int id
Definition: CR.cc:151
int resolve_alias(const std::vector< int > &aliases, int id)
Definition: CR.cc:92
afw::table::Key< double > b
bool _debias
Definition: CR.cc:957
int _ncol
Definition: CR.cc:955
Implementation of the Class MaskedImage.
A polymorphic base class for representing an image&#39;s Point Spread Function.
Definition: Psf.h:68
A functor class to allow users to process all the pixels in a Footprint.
ImageT val
Definition: CR.cc:154
image::Image< Pixel > Image
Image type returned by computeImage.
Definition: Psf.h:79
int col
Definition: CR.cc:152
MaskT setMaskFromFootprintList(lsst::afw::image::Mask< MaskT > *mask, std::vector< boost::shared_ptr< Footprint >> const &footprints, MaskT const bitmask)
OR bitmask into all the Mask&#39;s pixels which are in the set of Footprints.
double const _bkgd
Definition: CR.cc:312
Include files required for standard LSST Exception handling.
MaskedImageT::Mask::Pixel _badMask
Definition: CR.cc:956