LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
LSST Data Management Base Package
_saturated.cc
Go to the documentation of this file.
1 /*
2  * Handle saturated pixels when making colour images
3  */
4 #include <cmath>
5 #include "boost/format.hpp"
6 #include "lsst/afw/detection.h"
8 #include "Rgb.h"
9 
10 namespace lsst {
11 namespace afw {
12 namespace display {
13 
14 namespace {
15 template <typename ImageT>
16 class SetPixels {
17 public:
18  explicit SetPixels() {}
19 
20  void setValue(float value) { _value = value; }
21 
22  void operator()(lsst::geom::Point2I const& point, typename ImageT::Pixel& arrayInput) {
23  arrayInput = _value;
24  }
25 
26 private:
27  float _value{0};
28 };
29 } // namespace
30 
31 template <typename ImageT>
32 void replaceSaturatedPixels(ImageT& rim, // R image (e.g. i)
33  ImageT& gim, // G image (e.g. r)
34  ImageT& bim, // B image (e.g. g)
35  int borderWidth, // width of border used to estimate colour of saturated regions
36  float saturatedPixelValue // the brightness of a saturated pixel, once fixed
37 ) {
38  int const width = rim.getWidth(), height = rim.getHeight();
39  int const x0 = rim.getX0(), y0 = rim.getY0();
40 
41  if (width != gim.getWidth() || height != gim.getHeight() || x0 != gim.getX0() || y0 != gim.getY0()) {
42  throw LSST_EXCEPT(
44  str(boost::format("R image has different size/origin from G image "
45  "(%dx%d+%d+%d v. %dx%d+%d+%d") %
46  width % height % x0 % y0 % gim.getWidth() % gim.getHeight() % gim.getX0() % gim.getY0()));
47  }
48  if (width != bim.getWidth() || height != bim.getHeight() || x0 != bim.getX0() || y0 != bim.getY0()) {
49  throw LSST_EXCEPT(
51  str(boost::format("R image has different size/origin from B image "
52  "(%dx%d+%d+%d v. %dx%d+%d+%d") %
53  width % height % x0 % y0 % bim.getWidth() % bim.getHeight() % bim.getX0() % bim.getY0()));
54  }
55 
56  bool const useMaxPixel = !std::isfinite(saturatedPixelValue);
57 
58  SetPixels<typename ImageT::Image> setR, setG, setB; // functors used to set pixel values
59 
60  // Find all the saturated pixels in any of the three image
61  int const npixMin = 1; // minimum number of pixels in an object
62  afw::image::MaskPixel const SAT = rim.getMask()->getPlaneBitMask("SAT");
64 
65  detection::FootprintSet sat(*rim.getMask(), satThresh, npixMin);
66  sat.merge(detection::FootprintSet(*gim.getMask(), satThresh, npixMin));
67  sat.merge(detection::FootprintSet(*bim.getMask(), satThresh, npixMin));
68  // go through the list of saturated regions, determining the mean colour of the surrounding pixels
69  using FootprintList = detection::FootprintSet::FootprintList;
71  for (const auto& foot : *feet) {
72  auto const bigFoot = std::make_shared<detection::Footprint>(foot->getSpans()->dilated(borderWidth),
73  foot->getRegion());
74 
75  double sumR = 0, sumG = 0, sumB = 0; // sum of all non-saturated adjoining pixels
76  double maxR = 0, maxG = 0, maxB = 0; // maximum of non-saturated adjoining pixels
77 
78  for (auto span : *bigFoot->getSpans()) {
79  int const y = span.getY() - y0;
80  if (y < 0 || y >= height) {
81  continue;
82  }
83  int sx0 = span.getX0() - x0;
84  if (sx0 < 0) {
85  sx0 = 0;
86  }
87  int sx1 = span.getX1() - x0;
88  if (sx1 >= width) {
89  sx1 = width - 1;
90  }
91 
92  for (typename ImageT::iterator rptr = rim.at(sx0, y), rend = rim.at(sx1 + 1, y),
93  gptr = gim.at(sx0, y), bptr = bim.at(sx0, y);
94  rptr != rend; ++rptr, ++gptr, ++bptr) {
95  if (!((rptr.mask() | gptr.mask() | bptr.mask()) & SAT)) {
96  float val = rptr.image();
97  sumR += val;
98  if (val > maxR) {
99  maxR = val;
100  }
101 
102  val = gptr.image();
103  sumG += val;
104  if (val > maxG) {
105  maxG = val;
106  }
107 
108  val = bptr.image();
109  sumB += val;
110  if (val > maxB) {
111  maxB = val;
112  }
113  }
114  }
115  }
116  // OK, we have the mean fluxes for the pixels surrounding this set of saturated pixels
117  // so we can figure out the proper values to use for the saturated ones
118  float R = 0, G = 0, B = 0; // mean intensities
119  if (sumR + sumB + sumG > 0) {
120  if (sumR > sumG) {
121  if (sumR > sumB) {
122  R = useMaxPixel ? maxR : saturatedPixelValue;
123 
124  G = (R * sumG) / sumR;
125  B = (R * sumB) / sumR;
126  } else {
127  B = useMaxPixel ? maxB : saturatedPixelValue;
128  R = (B * sumR) / sumB;
129  G = (B * sumG) / sumB;
130  }
131  } else {
132  if (sumG > sumB) {
133  G = useMaxPixel ? maxG : saturatedPixelValue;
134  R = (G * sumR) / sumG;
135  B = (G * sumB) / sumG;
136  } else {
137  B = useMaxPixel ? maxB : saturatedPixelValue;
138  R = (B * sumR) / sumB;
139  G = (B * sumG) / sumB;
140  }
141  }
142  }
143  // Now that we know R, G, and B we can fix the values
144  setR.setValue(R);
145  foot->getSpans()->applyFunctor(setR, *rim.getImage());
146  setG.setValue(G);
147  foot->getSpans()->applyFunctor(setG, *gim.getImage());
148  setB.setValue(B);
149  foot->getSpans()->applyFunctor(setB, *bim.getImage());
150  }
151 }
152 
154  image::MaskedImage<float>& bim, int borderWidth,
155  float saturatedPixelValue);
156 } // namespace display
157 } // namespace afw
158 } // namespace lsst
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
int y
Definition: SpanSet.cc:48
A set of Footprints, associated with a MaskedImage.
Definition: FootprintSet.h:53
std::shared_ptr< FootprintList > getFootprints()
: Return the Footprints of detected objects
Definition: FootprintSet.h:156
void merge(FootprintSet const &rhs, int tGrow=0, int rGrow=0, bool isotropic=true)
Merge a FootprintSet into *this.
std::vector< std::shared_ptr< Footprint > > FootprintList
The FootprintSet's set of Footprints.
Definition: FootprintSet.h:56
A Threshold is used to pass a threshold value to detection algorithms.
Definition: Threshold.h:43
@ BITMASK
Use (pixels & (given mask))
Definition: Threshold.h:48
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:73
Reports invalid arguments.
Definition: Runtime.h:66
T isfinite(T... args)
void replaceSaturatedPixels(ImageT &rim, ImageT &gim, ImageT &bim, int borderWidth, float saturatedPixelValue)
Definition: _saturated.cc:32
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.
ImageT val
Definition: CR.cc:146