LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
_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"
14
15namespace lsst {
16namespace afw {
17namespace display {
18
19template <class T>
20static 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
50static 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
72template <class T>
73static 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);
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
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
166template <class T>
167std::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
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
afw::table::Key< double > sigma
double z
Definition Match.cc:44
int y
Definition SpanSet.cc:48
T back(T... args)
T begin(T... args)
A class to represent a 2-dimensional array of pixels.
Definition Image.h:51
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
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