LSST Applications g0b6bd0c080+a72a5dd7e6,g1182afd7b4+2a019aa3bb,g17e5ecfddb+2b8207f7de,g1d67935e3f+06cf436103,g38293774b4+ac198e9f13,g396055baef+6a2097e274,g3b44f30a73+6611e0205b,g480783c3b1+98f8679e14,g48ccf36440+89c08d0516,g4b93dc025c+98f8679e14,g5c4744a4d9+a302e8c7f0,g613e996a0d+e1c447f2e0,g6c8d09e9e7+25247a063c,g7271f0639c+98f8679e14,g7a9cd813b8+124095ede6,g9d27549199+a302e8c7f0,ga1cf026fa3+ac198e9f13,ga32aa97882+7403ac30ac,ga786bb30fb+7a139211af,gaa63f70f4e+9994eb9896,gabf319e997+ade567573c,gba47b54d5d+94dc90c3ea,gbec6a3398f+06cf436103,gc6308e37c7+07dd123edb,gc655b1545f+ade567573c,gcc9029db3c+ab229f5caf,gd01420fc67+06cf436103,gd877ba84e5+06cf436103,gdb4cecd868+6f279b5b48,ge2d134c3d5+cc4dbb2e3f,ge448b5faa6+86d1ceac1d,gecc7e12556+98f8679e14,gf3ee170dca+25247a063c,gf4ac96e456+ade567573c,gf9f5ea5b4d+ac198e9f13,gff490e6085+8c2580be5c,w.2022.27
LSST Data Management Base Package
Namespaces | Functions
lsst.afw.display Namespace Reference

Namespaces

namespace  ds9
 
namespace  ds9Regions
 
namespace  interface
 
namespace  rgb
 
namespace  utils
 
namespace  virtualDevice
 

Functions

 PYBIND11_MODULE (_simpleFits, mod)
 
 PYBIND11_MODULE (_rgb, mod)
 
template<typename ImageT >
void replaceSaturatedPixels (ImageT &rim, ImageT &gim, ImageT &bim, int borderWidth, float saturatedPixelValue)
 
template void replaceSaturatedPixels (image::MaskedImage< float > &rim, image::MaskedImage< float > &gim, image::MaskedImage< float > &bim, int borderWidth, float saturatedPixelValue)
 
template<class T >
std::pair< double, double > getZScale (image::Image< T > const &image, int const nSamples=1000, double const contrast=0.25)
 Calculate an IRAF/ds9-style zscaling. More...
 
template std::pair< double, double > getZScale (image::Image< std::uint16_t > const &image, int const nSamples, double const contrast)
 
template std::pair< double, double > getZScale (image::Image< float > const &image, int const nSamples, double const contrast)
 
template<typename ImageT >
void writeBasicFits (int fd, ImageT const &data, geom::SkyWcs const *Wcs, char const *title)
 
template<typename ImageT >
void writeBasicFits (std::string const &filename, ImageT const &data, geom::SkyWcs const *Wcs, char const *title)
 
template<typename ImageT >
void writeBasicFits (int fd, ImageT const &data, lsst::afw::geom::SkyWcs const *Wcs=nullptr, char const *title=nullptr)
 
template<typename ImageT >
void writeBasicFits (std::string const &filename, ImageT const &data, lsst::afw::geom::SkyWcs const *Wcs=nullptr, const char *title=nullptr)
 

Function Documentation

◆ getZScale() [1/3]

template std::pair< double, double > lsst::afw::display::getZScale ( image::Image< float > const &  image,
int const  nSamples,
double const  contrast 
)

◆ getZScale() [2/3]

template std::pair< double, double > lsst::afw::display::getZScale ( image::Image< std::uint16_t > const &  image,
int const  nSamples,
double const  contrast 
)

◆ getZScale() [3/3]

template<class T >
std::pair< double, double > lsst::afw::display::getZScale ( image::Image< T > const &  image,
int const  nSamples = 1000,
double const  contrast = 0.25 
)

Calculate an IRAF/ds9-style zscaling.

To quote Frank Valdes (http://iraf.net/forum/viewtopic.php?showtopic=134139)

ZSCALE ALGORITHM

The zscale algorithm is designed to display the  image  values  near
the  median  image  value  without  the  time  consuming  process of
computing a full image histogram.  This is particularly  useful  for
astronomical  images  which  generally  have a very peaked histogram
corresponding to  the  background  sky  in  direct  imaging  or  the
continuum in a two dimensional spectrum.

The  sample  of pixels, specified by values greater than zero in the
sample mask zmask or by an  image  section,  is  selected  up  to  a
maximum  of nsample pixels.  If a bad pixel mask is specified by the
bpmask parameter then any pixels with mask values which are  greater
than  zero  are not counted in the sample.  Only the first pixels up
to the limit are selected where the order is by line beginning  from
the  first line.  If no mask is specified then a grid of pixels with
even spacing along lines and columns that  make  up  a  number  less
than or equal to the maximum sample size is used.

If  a  contrast of zero is specified (or the zrange flag is used and
the image does not have a  valid  minimum/maximum  value)  then  the
minimum  and maximum of the sample is used for the intensity mapping
range.

If the contrast  is  not  zero  the  sample  pixels  are  ranked  in
brightness  to  form  the  function  I(i) where i is the rank of the
pixel and I is its value.  Generally the midpoint of  this  function
(the  median) is very near the peak of the image histogram and there
is a well defined slope about the midpoint which is related  to  the
width  of the histogram.  At the ends of the I(i) function there are
a few very bright and dark pixels due to objects and defects in  the
field.   To  determine  the  slope  a  linear  function  is fit with
iterative rejection;

<code>
        I(i) = intercept + slope * (i - midpoint)
</code>

If more than half of the points are rejected then there is  no  well
defined  slope  and  the full range of the sample defines z1 and z2.
Otherwise the endpoints of the linear function  are  used  (provided
they are within the original range of the sample):

<code>
        z1 = I(midpoint) + (slope / contrast) * (1 - midpoint)
        z2 = I(midpoint) + (slope / contrast) * (npoints - midpoint)
</code>

As  can  be  seen,  the parameter contrast may be used to adjust the
contrast produced by this algorithm.
Parameters
imageThe image we wish to stretch
nSamplesNumber of samples to use
contrastStretch parameter; see description

Definition at line 167 of file _scaling.cc.

167 {
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}
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
T back(T... args)
T begin(T... args)
Reports errors that are due to events beyond the control of the program.
Definition: Runtime.h:104
T empty(T... args)
T end(T... args)
T front(T... args)
T make_pair(T... args)
T max(T... args)
T min(T... args)
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
T size(T... args)
T sort(T... args)

◆ PYBIND11_MODULE() [1/2]

lsst::afw::display::PYBIND11_MODULE ( _rgb  ,
mod   
)

Definition at line 38 of file _rgb.cc.

38 {
39 /* Module level */
40 lsst::utils::python::WrapperCollection wrappers(mod, "lsst.afw.display.rbg");
41 wrappers.wrap([](auto &mod) {
42 mod.def("replaceSaturatedPixels", replaceSaturatedPixels<lsst::afw::image::MaskedImage<float>>,
43 "rim"_a, "gim"_a, "bim"_a, "borderWidth"_a = 2, "saturatedPixelValue"_a = 65535);
44 mod.def("getZScale", getZScale<std::uint16_t>, "image"_a, "nsamples"_a = 1000, "contrast"_a = 0.25);
45 mod.def("getZScale", getZScale<float>, "image"_a, "nsamples"_a = 1000, "contrast"_a = 0.25);
46 });
47 wrappers.finish();
48}
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:73
void replaceSaturatedPixels(ImageT &rim, ImageT &gim, ImageT &bim, int borderWidth, float saturatedPixelValue)
Definition: _saturated.cc:32

◆ PYBIND11_MODULE() [2/2]

lsst::afw::display::PYBIND11_MODULE ( _simpleFits  ,
mod   
)

Definition at line 58 of file _simpleFits.cc.

58 {
59 lsst::utils::python::WrapperCollection wrappers(mod, "lsst.afw.display");
60 declareAll<image::Image<std::uint16_t>>(wrappers);
61 declareAll<image::Image<std::uint64_t>>(wrappers);
62 declareAll<image::Image<int>>(wrappers);
63 declareAll<image::Image<float>>(wrappers);
64 declareAll<image::Image<double>>(wrappers);
65 declareAll<image::Mask<std::uint16_t>>(wrappers);
66 declareAll<image::Mask<image::MaskPixel>>(wrappers);
67 wrappers.finish();
68}

◆ replaceSaturatedPixels() [1/2]

template void lsst::afw::display::replaceSaturatedPixels ( image::MaskedImage< float > &  rim,
image::MaskedImage< float > &  gim,
image::MaskedImage< float > &  bim,
int  borderWidth,
float  saturatedPixelValue 
)

◆ replaceSaturatedPixels() [2/2]

template<typename ImageT >
void lsst::afw::display::replaceSaturatedPixels ( ImageT &  rim,
ImageT &  gim,
ImageT &  bim,
int  borderWidth,
float  saturatedPixelValue 
)

Definition at line 32 of file _saturated.cc.

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");
63 detection::Threshold const satThresh(SAT, detection::Threshold::BITMASK);
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;
70 std::shared_ptr<FootprintList> feet = sat.getFootprints();
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}
int y
Definition: SpanSet.cc:48
Reports invalid arguments.
Definition: Runtime.h:66
T isfinite(T... args)
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
ImageT val
Definition: CR.cc:146

◆ writeBasicFits() [1/4]

template<typename ImageT >
void lsst::afw::display::writeBasicFits ( int  fd,
ImageT const &  data,
geom::SkyWcs const *  Wcs,
char const *  title 
)

Definition at line 356 of file simpleFits.cc.

360 {
361 /*
362 * Allocate cards for FITS headers
363 */
364 std::list<Card> cards;
365 /*
366 * What sort if image is it?
367 */
368 int bitpix = lsst::afw::fits::getBitPix<typename ImageT::Pixel>();
369 if (bitpix == 20) { // cfitsio for "Unsigned short"
370 cards.emplace_back("BZERO", 32768.0, "");
371 cards.emplace_back("BSCALE", 1.0, "");
372 bitpix = 16;
373 } else if (bitpix == 0) {
374 throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError, "Unsupported image type");
375 }
376 /*
377 * Generate WcsA, pixel coordinates, allowing for X0 and Y0
378 */
379 addWcs("A", cards, data.getX0(), data.getY0());
380 /*
381 * Now WcsB, so that pixel (0,0) is correctly labelled (but ignoring XY0)
382 */
383 addWcs("B", cards);
384
385 if (title) {
386 cards.emplace_back("OBJECT", title, "Image being displayed");
387 }
388 /*
389 * Was there something else?
390 */
391 if (Wcs == nullptr) {
392 addWcs("", cards); // works around a ds9 bug that WCSA/B is ignored if no Wcs is present
393 } else {
394 using NameList = std::vector<std::string>;
395
396 auto shift = lsst::geom::Extent2D(-data.getX0(), -data.getY0());
397 auto newWcs = Wcs->copyAtShiftedPixelOrigin(shift);
398
399 std::shared_ptr<lsst::daf::base::PropertySet> metadata = newWcs->getFitsMetadata();
400
401 NameList paramNames = metadata->paramNames();
402
403 for (auto const &paramName : paramNames) {
404 if (paramName == "SIMPLE" || paramName == "BITPIX" || paramName == "NAXIS" || paramName == "NAXIS1" || paramName == "NAXIS2" ||
405 paramName == "XTENSION" || paramName == "PCOUNT" || paramName == "GCOUNT") {
406 continue;
407 }
408 std::type_info const &type = metadata->typeOf(paramName);
409 if (type == typeid(bool)) {
410 cards.emplace_back(paramName, metadata->get<bool>(paramName));
411 } else if (type == typeid(int)) {
412 cards.emplace_back(paramName, metadata->get<int>(paramName));
413 } else if (type == typeid(float)) {
414 cards.emplace_back(paramName, metadata->get<float>(paramName));
415 } else if (type == typeid(double)) {
416 cards.emplace_back(paramName, metadata->get<double>(paramName));
417 } else {
418 cards.emplace_back(paramName, metadata->get<std::string>(paramName));
419 }
420 }
421 }
422 /*
423 * Basic FITS stuff
424 */
425 const int naxis = 2; // == NAXIS
426 int naxes[naxis]; /* values of NAXIS1 etc */
427 naxes[0] = data.getWidth();
428 naxes[1] = data.getHeight();
429
430 write_fits_hdr(fd, bitpix, naxis, naxes, cards, 1);
431 for (int y = 0; y != data.getHeight(); ++y) {
432 if (write_fits_data(fd, bitpix, (char *)(data.row_begin(y)), (char *)(data.row_end(y))) < 0) {
434 (boost::format("Error writing data for row %d") % y).str());
435 }
436 }
437
438 pad_to_fits_record(fd, data.getWidth() * data.getHeight(), bitpix);
439}
char * data
Definition: BaseRecord.cc:61
table::Key< int > type
Definition: Detector.cc:163
T emplace_back(T... args)
T get(T... args)
Extent< double, 2 > Extent2D
Definition: Extent.h:400

◆ writeBasicFits() [2/4]

template<typename ImageT >
void lsst::afw::display::writeBasicFits ( int  fd,
ImageT const &  data,
lsst::afw::geom::SkyWcs const *  Wcs = nullptr,
char const *  title = nullptr 
)

◆ writeBasicFits() [3/4]

template<typename ImageT >
void lsst::afw::display::writeBasicFits ( std::string const &  filename,
ImageT const &  data,
geom::SkyWcs const *  Wcs,
char const *  title 
)

Definition at line 442 of file simpleFits.cc.

446 {
447 int fd;
448 if ((filename.c_str())[0] == '|') { // a command
449 const char *cmd = filename.c_str() + 1;
450 while (isspace(*cmd)) {
451 cmd++;
452 }
453
454 fd = fileno(popen(cmd, "w"));
455 } else {
456 fd = creat(filename.c_str(), 777);
457 }
458
459 if (fd < 0) {
461 (boost::format("Cannot open \"%s\"") % filename).str());
462 }
463
464 try {
465 writeBasicFits(fd, data, Wcs, title);
467 (void)close(fd);
468 throw;
469 }
470
471 (void)close(fd);
472}
Provides consistent interface for LSST exceptions.
Definition: Exception.h:107
void writeBasicFits(std::string const &filename, ImageT const &data, geom::SkyWcs const *Wcs, char const *title)
Definition: simpleFits.cc:442

◆ writeBasicFits() [4/4]

template<typename ImageT >
void lsst::afw::display::writeBasicFits ( std::string const &  filename,
ImageT const &  data,
lsst::afw::geom::SkyWcs const *  Wcs = nullptr,
const char *  title = nullptr 
)