LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
_scaling.cc
Go to the documentation of this file.
1 /*
2  * Calculate some scalings. Legacy code from NAOJ, which could probably be rewritten in python
3  * if the need arose. It could certainly use more LSST primitives (e.g. line fitting)
4  */
5 #include <cstdint>
6 #include <cmath>
7 #include <algorithm>
8 #include <cstdio>
9 #include <vector>
10 #include <stdexcept>
11 
12 #include "lsst/pex/exceptions.h"
13 #include "lsst/afw/image/Image.h"
14 
15 namespace lsst {
16 namespace afw {
17 namespace display {
18 
19 template <class T>
20 static void getSample(image::Image<T> const& image, std::size_t const nSamples, std::vector<T>& vSample) {
21  int const width = image.getWidth();
22  int const height = image.getHeight();
23 
24  // extract, from image, about nSample samples
25  // such that they form a grid.
26  vSample.reserve(nSamples);
27 
28  int const initialStride =
29  std::max(1, static_cast<int>(std::sqrt(width * height / static_cast<double>(nSamples))));
30 
31  for (int stride = initialStride; stride >= 1; --stride) {
32  vSample.clear();
33 
34  for (int y = 0; y < height; y += stride) {
35  for (int x = 0; x < width; x += stride) {
36  T const elem = image(x, y);
37  if (std::isfinite(elem)) {
38  vSample.push_back(elem);
39  }
40  }
41  }
42 
43  // if more than 80% of nSamples were sampled, OK.
44  if (5 * vSample.size() > 4 * nSamples) {
45  break;
46  }
47  }
48 }
49 
50 static inline double computeSigma(std::vector<double> const& vFlat, std::vector<int> const& vBadPix,
51  std::size_t const nGoodPix) {
52  if (nGoodPix <= 1) {
53  return 0;
54  }
55 
56  double sumz = 0, sumsq = 0;
57  for (unsigned i = 0; i < vFlat.size(); ++i) {
58  if (!vBadPix[i]) {
59  double const z = vFlat[i];
60 
61  sumz += z;
62  sumsq += z * z;
63  }
64  }
65 
66  double const goodPix = nGoodPix;
67  double const tmp = sumsq / (goodPix - 1) - sumz * sumz / (goodPix * (goodPix - 1));
68 
69  return (tmp > 0) ? std::sqrt(tmp) : 0;
70 }
71 
72 template <class T>
73 static std::pair<double, double> fitLine(
74  int* nGoodPixOut, // returned; it'd be nice to use std::tuple from C++11
75  std::vector<T> const& vSample, double const nSigmaClip, int const nGrow, int const minpix,
76  int const nIter) {
77  // map the indices of vSample to [-1.0, 1.0]
78  double const xscale = 2.0 / (vSample.size() - 1);
79  std::vector<double> xnorm;
80  xnorm.reserve(vSample.size());
81  for (std::size_t i = 0; i < vSample.size(); ++i) {
82  xnorm.push_back(i * xscale - 1.0);
83  }
84 
85  // Mask that is used in k-sigma clipping
86  std::vector<int> vBadPix(vSample.size(), 0);
87 
88  int nGoodPix = vSample.size();
89  int nGoodPixOld = nGoodPix + 1;
90 
91  // values to be obtained
92  double intercept = 0;
93  double slope = 0;
94 
95  for (int iteration = 0; iteration < nIter; ++iteration) {
96  if (nGoodPix < minpix || nGoodPix >= nGoodPixOld) {
97  break;
98  }
99 
100  double sum = nGoodPix;
101  double sumx = 0, sumy = 0, sumxx = 0, sumxy = 0;
102  for (std::size_t i = 0; i < vSample.size(); ++i) {
103  if (!vBadPix[i]) {
104  double const x = xnorm[i];
105  double const y = vSample[i];
106 
107  sumx += x;
108  sumy += y;
109  sumxx += x * x;
110  sumxy += x * y;
111  }
112  }
113 
114  double delta = sum * sumxx - sumx * sumx;
115 
116  // slope and intercept
117  intercept = (sumxx * sumy - sumx * sumxy) / delta;
118  slope = (sum * sumxy - sumx * sumy) / delta;
119 
120  // residue
121  std::vector<double> vFlat;
122  vFlat.reserve(vSample.size());
123  for (unsigned i = 0; i < vSample.size(); ++i) {
124  vFlat.push_back(vSample[i] - (xnorm[i] * slope + intercept));
125  }
126 
127  // Threshold of k-sigma clipping
128  double const sigma = computeSigma(vFlat, vBadPix, nGoodPix);
129  double const hcut = sigma * nSigmaClip;
130  double const lcut = -hcut;
131 
132  // revise vBadPix
133  nGoodPixOld = nGoodPix;
134  nGoodPix = 0;
135 
136  for (unsigned i = 0; i < vSample.size(); ++i) {
137  double val = vFlat[i];
138  if (val < lcut || hcut < val) {
139  vBadPix[i] = 1;
140  }
141  }
142 
143  // blurr vBadPix
144  std::vector<int> vBadPix_new;
145  vBadPix_new.reserve(vSample.size());
146  for (unsigned x = 0; x < vSample.size(); ++x) {
147  int imin = (static_cast<int>(x) > nGrow) ? x - nGrow : -1;
148  int val = 0;
149  for (int i = x; i > imin; --i) {
150  val += vBadPix[i];
151  }
152  vBadPix_new.push_back(val ? 1 : 0);
153  if (!val) {
154  ++nGoodPix;
155  }
156  }
157  vBadPix = vBadPix_new;
158  }
159 
160  // return the scale of x-axis
161  *nGoodPixOut = nGoodPix;
162 
163  return std::make_pair(intercept - slope, slope * xscale);
164 }
165 
166 template <class T>
167 std::pair<double, double> getZScale(image::Image<T> const& image, int const nSamples, double const contrast) {
168  // extract samples
169  std::vector<T> vSample;
170  getSample(image, nSamples, vSample);
171  int nPix = vSample.size();
172 
173  if (vSample.empty()) {
174  throw LSST_EXCEPT(pex::exceptions::RuntimeError, "ZScale: No pixel in image is finite");
175  }
176 
177  std::sort(vSample.begin(), vSample.end());
178 
179  // max, min, median
180  // N.b. you can get a median in linear time, but we need the sorted array for fitLine()
181  // If we wanted to speed this up, the best option would be to quantize
182  // the pixel values and build a histogram
183  double const zmin = vSample.front();
184  double const zmax = vSample.back();
185  int const iCenter = nPix / 2;
186  T median = (nPix & 1) ? vSample[iCenter] : (vSample[iCenter] + vSample[iCenter + 1]) / 2;
187 
188  // fit a line to the sorted sample
189  const int maxRejectionRatio = 2;
190  const int npixelsMin = 5;
191 
192  int minpix = std::max(npixelsMin, nPix / maxRejectionRatio);
193  int nGrow = std::max(1, nPix / 100);
194 
195  const double nSigmaClip = 2.5;
196  const int nIterations = 5;
197 
198  int nGoodPix = 0;
199  std::pair<double, double> ret = fitLine(&nGoodPix, vSample, nSigmaClip, nGrow, minpix, nIterations);
200 #if 0 // unused, but calculated and potentially useful
201  double const zstart = ret.first;
202 #endif
203  double const zslope = ret.second;
204 
205  double z1, z2;
206  if (nGoodPix < minpix) {
207  z1 = zmin;
208  z2 = zmax;
209  } else {
210  double const slope = zslope / contrast;
211 
212  z1 = std::max(zmin, median - iCenter * slope);
213  z2 = std::min(zmax, median + (nPix - iCenter - 1) * slope);
214  }
215 
216  return std::make_pair(z1, z2);
217 }
218 //
219 // Explicit instantiations
220 #define INSTANTIATE_GETZSCALE(T) \
221  template std::pair<double, double> getZScale(image::Image<T> const& image, int const nSamples, \
222  double const contrast)
223 
226 } // namespace display
227 } // namespace afw
228 } // namespace lsst
#define INSTANTIATE_GETZSCALE(T)
Definition: _scaling.cc:220
double x
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:50
afw::table::Key< afw::table::Array< ImagePixelT > > image
double z
Definition: Match.cc:44
int y
Definition: SpanSet.cc:49
T back(T... args)
T begin(T... args)
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:58
Reports errors that are due to events beyond the control of the program.
Definition: Runtime.h:104
T clear(T... args)
T empty(T... args)
T end(T... args)
T front(T... args)
T isfinite(T... args)
T make_pair(T... args)
T max(T... args)
T min(T... args)
std::pair< double, double > getZScale(image::Image< T > const &image, int const nSamples, double const contrast)
Calculate an IRAF/ds9-style zscaling.
Definition: _scaling.cc:167
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
A base class for image defects.
T push_back(T... args)
T reserve(T... args)
T size(T... args)
T sort(T... args)
T sqrt(T... args)
ImageT val
Definition: CR.cc:146