LSSTApplications  17.0+11,17.0+34,17.0+56,17.0+57,17.0+59,17.0+7,17.0-1-g377950a+33,17.0.1-1-g114240f+2,17.0.1-1-g4d4fbc4+28,17.0.1-1-g55520dc+49,17.0.1-1-g5f4ed7e+52,17.0.1-1-g6dd7d69+17,17.0.1-1-g8de6c91+11,17.0.1-1-gb9095d2+7,17.0.1-1-ge9fec5e+5,17.0.1-1-gf4e0155+55,17.0.1-1-gfc65f5f+50,17.0.1-1-gfc6fb1f+20,17.0.1-10-g87f9f3f+1,17.0.1-11-ge9de802+16,17.0.1-16-ga14f7d5c+4,17.0.1-17-gc79d625+1,17.0.1-17-gdae4c4a+8,17.0.1-2-g26618f5+29,17.0.1-2-g54f2ebc+9,17.0.1-2-gf403422+1,17.0.1-20-g2ca2f74+6,17.0.1-23-gf3eadeb7+1,17.0.1-3-g7e86b59+39,17.0.1-3-gb5ca14a,17.0.1-3-gd08d533+40,17.0.1-30-g596af8797,17.0.1-4-g59d126d+4,17.0.1-4-gc69c472+5,17.0.1-6-g5afd9b9+4,17.0.1-7-g35889ee+1,17.0.1-7-gc7c8782+18,17.0.1-9-gc4bbfb2+3,w.2019.22
LSSTDataManagementBasePackage
Namespaces | Functions
lsst.afw.display Namespace Reference

Namespaces

 displayLib
 
 ds9
 
 ds9Regions
 
 interface
 
 rgb
 
 utils
 
 virtualDevice
 

Functions

 PYBIND11_MODULE (_simpleFits, mod)
 
 PYBIND11_MODULE (rgb, mod)
 
template<typename ImageT >
void replaceSaturatedPixels (ImageT &rim, ImageT &gim, ImageT &bim, int borderWidth=2, float saturatedPixelValue=65535)
 
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 void replaceSaturatedPixels (image::MaskedImage< float > &rim, image::MaskedImage< float > &gim, image::MaskedImage< float > &bim, int borderWidth, float saturatedPixelValue)
 
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)
 

Function Documentation

◆ getZScale() [1/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 }
T empty(T... args)
T front(T... args)
T end(T... args)
T min(T... args)
T make_pair(T... args)
T max(T... args)
T size(T... args)
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
STL class.
T begin(T... args)
T back(T... args)
T sort(T... args)

◆ 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 std::pair<double, double> lsst::afw::display::getZScale ( image::Image< float > const &  image,
int const  nSamples,
double const  contrast 
)

◆ PYBIND11_MODULE() [1/2]

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

Definition at line 37 of file rgb.cc.

37  {
38  /* Module level */
39  mod.def("replaceSaturatedPixels", replaceSaturatedPixels<lsst::afw::image::MaskedImage<float>>, "rim"_a,
40  "gim"_a, "bim"_a, "borderWidth"_a = 2, "saturatedPixelValue"_a = 65535);
41  mod.def("getZScale", getZScale<std::uint16_t>, "image"_a, "nsamples"_a = 1000, "contrast"_a = 0.25);
42  mod.def("getZScale", getZScale<float>, "image"_a, "nsamples"_a = 1000, "contrast"_a = 0.25);
43 }
void replaceSaturatedPixels(ImageT &rim, ImageT &gim, ImageT &bim, int borderWidth=2, float saturatedPixelValue=65535)
Definition: saturated.cc:30
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:74

◆ PYBIND11_MODULE() [2/2]

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

Definition at line 53 of file _simpleFits.cc.

53  {
54  declareAll<image::Image<std::uint16_t>>(mod);
55  declareAll<image::Image<std::uint64_t>>(mod);
56  declareAll<image::Image<int>>(mod);
57  declareAll<image::Image<float>>(mod);
58  declareAll<image::Image<double>>(mod);
59  declareAll<image::Mask<std::uint16_t>>(mod);
60  declareAll<image::Mask<image::MaskPixel>>(mod);
61 }

◆ replaceSaturatedPixels() [1/2]

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

Definition at line 30 of file saturated.cc.

35  {
36  int const width = rim.getWidth(), height = rim.getHeight();
37  int const x0 = rim.getX0(), y0 = rim.getY0();
38 
39  if (width != gim.getWidth() || height != gim.getHeight() || x0 != gim.getX0() || y0 != gim.getY0()) {
40  throw LSST_EXCEPT(
41  pex::exceptions::InvalidParameterError,
42  str(boost::format("R image has different size/origin from G image "
43  "(%dx%d+%d+%d v. %dx%d+%d+%d") %
44  width % height % x0 % y0 % gim.getWidth() % gim.getHeight() % gim.getX0() % gim.getY0()));
45  }
46  if (width != bim.getWidth() || height != bim.getHeight() || x0 != bim.getX0() || y0 != bim.getY0()) {
47  throw LSST_EXCEPT(
48  pex::exceptions::InvalidParameterError,
49  str(boost::format("R image has different size/origin from B image "
50  "(%dx%d+%d+%d v. %dx%d+%d+%d") %
51  width % height % x0 % y0 % bim.getWidth() % bim.getHeight() % bim.getX0() % bim.getY0()));
52  }
53 
54  bool const useMaxPixel = !std::isfinite(saturatedPixelValue);
55 
56  SetPixels<typename ImageT::Image> setR, setG, setB; // functors used to set pixel values
57 
58  // Find all the saturated pixels in any of the three image
59  int const npixMin = 1; // minimum number of pixels in an object
60  afw::image::MaskPixel const SAT = rim.getMask()->getPlaneBitMask("SAT");
61  detection::Threshold const satThresh(SAT, detection::Threshold::BITMASK);
62 
63  detection::FootprintSet sat(*rim.getMask(), satThresh, npixMin);
64  sat.merge(detection::FootprintSet(*gim.getMask(), satThresh, npixMin));
65  sat.merge(detection::FootprintSet(*bim.getMask(), satThresh, npixMin));
66  // go through the list of saturated regions, determining the mean colour of the surrounding pixels
67  typedef detection::FootprintSet::FootprintList FootprintList;
68  std::shared_ptr<FootprintList> feet = sat.getFootprints();
69  for (FootprintList::const_iterator ptr = feet->begin(), end = feet->end(); ptr != end; ++ptr) {
71  auto const bigFoot = std::make_shared<detection::Footprint>(foot->getSpans()->dilated(borderWidth),
72  foot->getRegion());
73 
74  double sumR = 0, sumG = 0, sumB = 0; // sum of all non-saturated adjoining pixels
75  double maxR = 0, maxG = 0, maxB = 0; // maximum of non-saturated adjoining pixels
76 
77  for (auto span = bigFoot->getSpans()->begin(), send = bigFoot->getSpans()->end(); span != send;
78  ++span) {
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 }
uint64_t * ptr
Definition: RangeSet.cc:88
int y
Definition: SpanSet.cc:49
ImageT val
Definition: CR.cc:146
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:168
T isfinite(T... args)
std::int32_t MaskPixel
default type for Masks and MaskedImage Masks
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
int end

◆ replaceSaturatedPixels() [2/2]

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

◆ writeBasicFits() [1/2]

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

Definition at line 363 of file simpleFits.cc.

367  {
368  /*
369  * Allocate cards for FITS headers
370  */
371  std::list<Card> cards;
372  /*
373  * What sort if image is it?
374  */
375  int bitpix = lsst::afw::fits::getBitPix<typename ImageT::Pixel>();
376  if (bitpix == 20) { // cfitsio for "Unsigned short"
377  cards.push_back(Card("BZERO", 32768.0, ""));
378  cards.push_back(Card("BSCALE", 1.0, ""));
379  bitpix = 16;
380  } else if (bitpix == 0) {
381  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError, "Unsupported image type");
382  }
383  /*
384  * Generate WcsA, pixel coordinates, allowing for X0 and Y0
385  */
386  addWcs("A", cards, data.getX0(), data.getY0());
387  /*
388  * Now WcsB, so that pixel (0,0) is correctly labelled (but ignoring XY0)
389  */
390  addWcs("B", cards);
391 
392  if (title) {
393  cards.push_back(Card("OBJECT", title, "Image being displayed"));
394  }
395  /*
396  * Was there something else?
397  */
398  if (Wcs == NULL) {
399  addWcs("", cards); // works around a ds9 bug that WCSA/B is ignored if no Wcs is present
400  } else {
401  typedef std::vector<std::string> NameList;
402 
403  auto shift = lsst::geom::Extent2D(-data.getX0(), -data.getY0());
404  auto newWcs = Wcs->copyAtShiftedPixelOrigin(shift);
405 
406  std::shared_ptr<lsst::daf::base::PropertySet> metadata = newWcs->getFitsMetadata();
407 
408  NameList paramNames = metadata->paramNames();
409 
410  for (NameList::const_iterator i = paramNames.begin(), end = paramNames.end(); i != end; ++i) {
411  if (*i == "SIMPLE" || *i == "BITPIX" || *i == "NAXIS" || *i == "NAXIS1" || *i == "NAXIS2" ||
412  *i == "XTENSION" || *i == "PCOUNT" || *i == "GCOUNT") {
413  continue;
414  }
415  std::type_info const &type = metadata->typeOf(*i);
416  if (type == typeid(bool)) {
417  cards.push_back(Card(*i, metadata->get<bool>(*i)));
418  } else if (type == typeid(int)) {
419  cards.push_back(Card(*i, metadata->get<int>(*i)));
420  } else if (type == typeid(float)) {
421  cards.push_back(Card(*i, metadata->get<float>(*i)));
422  } else if (type == typeid(double)) {
423  cards.push_back(Card(*i, metadata->get<double>(*i)));
424  } else {
425  cards.push_back(Card(*i, metadata->get<std::string>(*i)));
426  }
427  }
428  }
429  /*
430  * Basic FITS stuff
431  */
432  const int naxis = 2; // == NAXIS
433  int naxes[naxis]; /* values of NAXIS1 etc */
434  naxes[0] = data.getWidth();
435  naxes[1] = data.getHeight();
436 
437  write_fits_hdr(fd, bitpix, naxis, naxes, cards, 1);
438  for (int y = 0; y != data.getHeight(); ++y) {
439  if (write_fits_data(fd, bitpix, (char *)(data.row_begin(y)), (char *)(data.row_end(y))) < 0) {
441  (boost::format("Error writing data for row %d") % y).str());
442  }
443  }
444 
445  pad_to_fits_record(fd, data.getWidth() * data.getHeight(), bitpix);
446 }
int y
Definition: SpanSet.cc:49
STL class.
T push_back(T... args)
char * data
Definition: BaseRecord.cc:62
table::Key< int > type
Definition: Detector.cc:167
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:168
STL class.
T get(T... args)
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
Extent< double, 2 > Extent2D
Definition: Extent.h:400
int end
Reports errors that are due to events beyond the control of the program.
Definition: Runtime.h:104

◆ writeBasicFits() [2/2]

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 449 of file simpleFits.cc.

453  {
454  int fd;
455  if ((filename.c_str())[0] == '|') { // a command
456  const char *cmd = filename.c_str() + 1;
457  while (isspace(*cmd)) {
458  cmd++;
459  }
460 
461  fd = fileno(popen(cmd, "w"));
462  } else {
463  fd = creat(filename.c_str(), 777);
464  }
465 
466  if (fd < 0) {
468  (boost::format("Cannot open \"%s\"") % filename).str());
469  }
470 
471  try {
472  writeBasicFits(fd, data, Wcs, title);
474  (void)close(fd);
475  throw;
476  }
477 
478  (void)close(fd);
479 }
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:449
char * data
Definition: BaseRecord.cc:62
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:168
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
T isspace(T... args)
Reports errors that are due to events beyond the control of the program.
Definition: Runtime.h:104