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
simpleFits.cc
Go to the documentation of this file.
1 /*
2  * This file is part of afw.
3  *
4  * Developed for the LSST Data Management System.
5  * This product includes software developed by the LSST Project
6  * (https://www.lsst.org).
7  * See the COPYRIGHT file at the top-level directory of this distribution
8  * for details of code ownership.
9  *
10  * This program is free software: you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation, either version 3 of the License, or
13  * (at your option) any later version.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the GNU General Public License
21  * along with this program. If not, see <https://www.gnu.org/licenses/>.
22  */
23 
24 /*
25  * Write a FITS image to a file descriptor; useful for talking to DS9
26  *
27  * This version knows about LSST data structures
28  */
29 #ifndef DOXYGEN // Doxygen doesn't like includes inside namespaces
30 namespace posix { // here so no-one includes them first outside namespace posix {}
31 #include <unistd.h>
32 #include <fcntl.h>
33 } // namespace posix
34 #endif // !DOXYGEN
35 using namespace posix;
36 #include <cstdint>
37 #include <cstdlib>
38 #include <cstdio>
39 #include <cstring>
40 #include <cassert>
41 #include <cctype>
42 #include <any>
43 
44 #include "lsst/pex/exceptions.h"
45 #include "lsst/afw/fits.h"
46 #include "simpleFits.h"
47 
48 namespace image = lsst::afw::image;
50 
51 #define FITS_SIZE 2880
52 
54 class Card {
55 public:
56  Card(const std::string &name, bool val, const char *commnt = "")
57  : keyword(name), value(val), comment(commnt) {}
58  Card(const std::string &name, int val, const char *commnt = "")
59  : keyword(name), value(val), comment(commnt) {}
60  Card(const std::string &name, double val, const char *commnt = "")
61  : keyword(name), value(val), comment(commnt) {}
62  Card(const std::string &name, float val, const char *commnt = "")
63  : keyword(name), value(val), comment(commnt) {}
64  Card(const std::string &name, const std::string &val, const char *commnt = "")
65  : keyword(name), value(val), comment(commnt) {}
66  Card(const std::string &name, const char *val, const char *commnt = "")
67  : keyword(name), value(std::string(val)), comment(commnt) {}
68 
69  ~Card() = default;
70 
71  int write(int fd, int ncard, char *record) const;
72 
73  std::string keyword;
74  std::any value;
75  std::string comment;
76 };
77 
78 /*
79  * Write a Card
80  */
81 int Card::write(int fd, int ncard, char *record) const {
82  char *card = &record[80 * ncard];
83 
84  if (value.type() == typeid(std::string)) {
85  std::string const &str = std::any_cast<std::string>(value);
86  if (keyword == "" || keyword == "COMMENT" || keyword == "END" || keyword == "HISTORY") {
87  sprintf(card, "%-8.8s%-72s", keyword.c_str(), str.c_str());
88  } else {
89  sprintf(card, "%-8.8s= '%s' %c%-*s", keyword.c_str(), str.c_str(), (comment == "" ? ' ' : '/'),
90  (int)(80 - 14 - str.size()), comment.c_str());
91  }
92  } else {
93  sprintf(card, "%-8.8s= ", keyword.c_str());
94  card += 10;
95  if (value.type() == typeid(bool)) {
96  sprintf(card, "%20s", std::any_cast<bool>(value) ? "T" : "F");
97  } else if (value.type() == typeid(int)) {
98  sprintf(card, "%20d", std::any_cast<int>(value));
99  } else if (value.type() == typeid(double)) {
100  sprintf(card, "%20.10f", std::any_cast<double>(value));
101  } else if (value.type() == typeid(float)) {
102  sprintf(card, "%20.7f", std::any_cast<float>(value));
103  }
104  card += 20;
105  sprintf(card, " %c%-48s", (comment == "" ? ' ' : '/'), comment.c_str());
106  }
107  /*
108  * Write record if full
109  */
110  if (++ncard == 36) {
111  if (posix::write(fd, record, FITS_SIZE) != FITS_SIZE) {
112  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError, "Cannot write header record");
113  }
114  ncard = 0;
115  }
116 
117  return ncard;
118 }
120 
121 /*
122  * Utilities
123  *
124  * Flip high-order bit so as to write unsigned short to FITS. Grrr.
125  */
126 namespace {
127 void flip_high_bit(char *arr, // array that needs bits swapped
128  const int n) { // number of bytes in arr
129  if (n % 2 != 0) {
131  (boost::format("Attempt to bit flip odd number of bytes: %d") % n).str());
132  }
133 
134  unsigned short *uarr = reinterpret_cast<unsigned short *>(arr);
135  for (unsigned short *end = uarr + n / 2; uarr < end; ++uarr) {
136  *uarr ^= 0x8000;
137  }
138 }
139 } // namespace
140 
141 /*
142  * Byte swap ABABAB -> BABABAB in place
143  */
144 namespace {
145 void swap_2(char *arr, // array to swap
146  const int n) { // number of bytes
147  if (n % 2 != 0) {
149  (boost::format("Attempt to byte swap odd number of bytes: %d") % n).str());
150  }
151 
152  for (char *end = arr + n; arr < end; arr += 2) {
153  char t = arr[0];
154  arr[0] = arr[1];
155  arr[1] = t;
156  }
157 }
158 /*
159  * Byte swap ABCDABCD -> DCBADCBA in place (e.g. sun <--> vax)
160  */
161 void swap_4(char *arr, // array to swap
162  const int n) { // number of bytes
163  if (n % 4 != 0) {
165  (boost::format("Attempt to byte swap non-multiple of 4 bytes: %d") % n).str());
166  }
167 
168  for (char *end = arr + n; arr < end; arr += 4) {
169  char t = arr[0];
170  arr[0] = arr[3];
171  arr[3] = t;
172  t = arr[1];
173  arr[1] = arr[2];
174  arr[2] = t;
175  }
176 }
177 
178 /*
179  * Byte swap ABCDEFGH -> HGFEDCBA in place (e.g. sun <--> vax)
180  */
181 void swap_8(char *arr, // array to swap
182  const int n) { // number of bytes
183  if (n % 8 != 0) {
185  (boost::format("Attempt to byte swap non-multiple of 8 bytes: %d") % n).str());
186  }
187 
188  for (char *end = arr + n; arr < end; arr += 8) {
189  char t = arr[0];
190  arr[0] = arr[7];
191  arr[7] = t;
192  t = arr[1];
193  arr[1] = arr[6];
194  arr[6] = t;
195  t = arr[2];
196  arr[2] = arr[5];
197  arr[5] = t;
198  t = arr[3];
199  arr[3] = arr[4];
200  arr[4] = t;
201  }
202 }
203 
204 int write_fits_hdr(int fd, int bitpix, int naxis, int *naxes, std::list<Card> &cards, /* extra header cards */
205  int primary) /* is this the primary HDU? */
206 {
207  int i;
208  char record[FITS_SIZE + 1]; /* write buffer */
209 
210  int ncard = 0;
211  if (primary) {
212  Card card("SIMPLE", true);
213  ncard = card.write(fd, ncard, record);
214  } else {
215  Card card("XTENSION", "IMAGE");
216  ncard = card.write(fd, ncard, record);
217  }
218 
219  {
220  Card card("BITPIX", bitpix);
221  ncard = card.write(fd, ncard, record);
222  }
223  {
224  Card card("NAXIS", naxis);
225  ncard = card.write(fd, ncard, record);
226  }
227  for (i = 0; i < naxis; i++) {
228  char key[] = "NAXIS.";
229  sprintf(key, "NAXIS%d", i + 1);
230  Card card(key, naxes[i]);
231  ncard = card.write(fd, ncard, record);
232  }
233  if (primary) {
234  Card card("EXTEND", true, "There may be extensions");
235  ncard = card.write(fd, ncard, record);
236  }
237  /*
238  * Write extra header cards
239  */
240  for (std::list<Card>::const_iterator card = cards.begin(); card != cards.end(); card++) {
241  ncard = card->write(fd, ncard, record);
242  }
243 
244  {
245  Card card("END", "");
246  ncard = card.write(fd, ncard, record);
247  }
248  while (ncard != 0) {
249  Card card("", "");
250  ncard = card.write(fd, ncard, record);
251  }
252 
253  return 0;
254 }
255 
256 /*
257  * Pad out to a FITS record boundary
258  */
259 void pad_to_fits_record(int fd, // output file descriptor
260  int npixel, // number of pixels already written to HDU
261  int bitpix // bitpix for this datatype
262 ) {
263  const int bytes_per_pixel = (bitpix > 0 ? bitpix : -bitpix) / 8;
264  int nbyte = npixel * bytes_per_pixel;
265 
266  if (nbyte % FITS_SIZE != 0) {
267  char record[FITS_SIZE + 1]; /* write buffer */
268 
269  nbyte = FITS_SIZE - nbyte % FITS_SIZE;
270  memset(record, ' ', nbyte);
271  if (write(fd, record, nbyte) != nbyte) {
273  "error padding file to multiple of fits block size");
274  }
275  }
276 }
277 
278 int write_fits_data(int fd, int bitpix, char *begin, char *end) {
279  const int bytes_per_pixel = (bitpix > 0 ? bitpix : -bitpix) / 8;
280  int swap_bytes = 0; // the default
281 #if defined(LSST_LITTLE_ENDIAN) // we'll need to byte swap FITS
282  if (bytes_per_pixel > 1) {
283  swap_bytes = 1;
284  }
285 #endif
286 
287  char *buff = nullptr; // I/O buffer
288  bool allocated = false; // do I need to free it?
289  if (swap_bytes || bitpix == 16) {
290  buff = new char[FITS_SIZE * bytes_per_pixel];
291  allocated = true;
292  }
293 
294  int nbyte = end - begin;
295  int nwrite = (nbyte > FITS_SIZE) ? FITS_SIZE : nbyte;
296  for (char *ptr = begin; ptr != end; nbyte -= nwrite, ptr += nwrite) {
297  if (end - ptr < nwrite) {
298  nwrite = end - ptr;
299  }
300 
301  if (swap_bytes) {
302  memcpy(buff, ptr, nwrite);
303  if (bitpix == 16) { // flip high-order bit
304  flip_high_bit(buff, nwrite);
305  }
306 
307  if (bytes_per_pixel == 2) {
308  swap_2(buff, nwrite);
309  } else if (bytes_per_pixel == 4) {
310  swap_4(buff, nwrite);
311  } else if (bytes_per_pixel == 8) {
312  swap_8(buff, nwrite);
313  } else {
314  fprintf(stderr, "You cannot get here\n");
315  abort();
316  }
317  } else {
318  if (bitpix == 16) { // flip high-order bit
319  memcpy(buff, ptr, nwrite);
320  flip_high_bit(buff, nwrite);
321  } else {
322  buff = ptr;
323  }
324  }
325 
326  if (write(fd, buff, nwrite) != nwrite) {
327  perror("Error writing image: ");
328  break;
329  }
330  }
331 
332  if (allocated) {
333  delete buff;
334  }
335 
336  return (nbyte == 0 ? 0 : -1);
337 }
338 
339 void addWcs(std::string const &wcsName, std::list<Card> &cards, int x0 = 0.0, int y0 = 0.0) {
340  cards.emplace_back(str(boost::format("CRVAL1%s") % wcsName), x0, "(output) Column pixel of Reference Pixel");
341  cards.emplace_back(str(boost::format("CRVAL2%s") % wcsName), y0, "(output) Row pixel of Reference Pixel");
342  cards.emplace_back(str(boost::format("CRPIX1%s") % wcsName), 1.0, "Column Pixel Coordinate of Reference");
343  cards.emplace_back(str(boost::format("CRPIX2%s") % wcsName), 1.0, "Row Pixel Coordinate of Reference");
344  cards.emplace_back(str(boost::format("CTYPE1%s") % wcsName), "LINEAR", "Type of projection");
345  cards.emplace_back(str(boost::format("CTYPE1%s") % wcsName), "LINEAR", "Type of projection");
346  cards.emplace_back(str(boost::format("CUNIT1%s") % wcsName), "PIXEL", "Column unit");
347  cards.emplace_back(str(boost::format("CUNIT2%s") % wcsName), "PIXEL", "Row unit");
348 }
349 } // namespace
350 
351 namespace lsst {
352 namespace afw {
353 namespace display {
354 
355 template <typename ImageT>
356 void writeBasicFits(int fd, // file descriptor to write to
357  ImageT const &data, // The data to write
358  geom::SkyWcs const *Wcs, // which Wcs to use for pixel
359  char const *title // title to write to DS9
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 }
440 
441 template <typename ImageT>
442 void writeBasicFits(std::string const &filename, // file to write, or "| cmd"
443  ImageT const &data, // The data to write
444  geom::SkyWcs const *Wcs, // which Wcs to use for pixel
445  char const *title // title to write to DS9
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 }
473 
475 #define INSTANTIATE(IMAGET) \
476  template void writeBasicFits(int, IMAGET const &, geom::SkyWcs const *, char const *); \
477  template void writeBasicFits(std::string const &, IMAGET const &, geom::SkyWcs const *, char const *)
478 
479 #define INSTANTIATE_IMAGE(T) INSTANTIATE(lsst::afw::image::Image<T>)
480 #define INSTANTIATE_MASK(T) INSTANTIATE(lsst::afw::image::Mask<T>)
481 
483 INSTANTIATE_IMAGE(int);
484 INSTANTIATE_IMAGE(float);
485 INSTANTIATE_IMAGE(double);
487 
488 INSTANTIATE_MASK(std::uint16_t);
489 INSTANTIATE_MASK(image::MaskPixel);
491 } // namespace display
492 } // namespace afw
493 } // namespace lsst
table::Key< std::string > name
Definition: Amplifier.cc:116
char * data
Definition: BaseRecord.cc:61
int end
table::Key< int > type
Definition: Detector.cc:163
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
afw::table::Key< afw::table::Array< ImagePixelT > > image
uint64_t * ptr
Definition: RangeSet.cc:88
int y
Definition: SpanSet.cc:48
T abort(T... args)
T begin(T... args)
T c_str(T... args)
Class for storing generic metadata.
Definition: PropertySet.h:66
std::vector< std::string > paramNames(bool topLevelOnly=true) const
A variant of names that excludes the names of subproperties.
std::type_info const & typeOf(std::string const &name) const
Get the type of values for a property name (possibly hierarchical).
T get(std::string const &name) const
Get the last value for a property name (possibly hierarchical).
Provides consistent interface for LSST exceptions.
Definition: Exception.h:107
Reports errors that are due to events beyond the control of the program.
Definition: Runtime.h:104
T emplace_back(T... args)
T end(T... args)
T sprintf(T... args)
#define INSTANTIATE_IMAGE(IMAGE)
Definition: ImagePca.cc:125
T memcpy(T... args)
T memset(T... args)
void writeBasicFits(std::string const &filename, ImageT const &data, geom::SkyWcs const *Wcs, char const *title)
Definition: simpleFits.cc:442
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
void write(OutputArchiveHandle &handle) const override
bool any(CoordinateExpr< N > const &expr) noexcept
Return true if any elements are true.
Extent< double, 2 > Extent2D
Definition: Extent.h:400
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
A base class for image defects.
T perror(T... args)
#define FITS_SIZE
Definition: simpleFits.cc:51
T size(T... args)
ImageT val
Definition: CR.cc:146