LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
simpleFits.cc
Go to the documentation of this file.
1 /*
2  * LSST Data Management System
3  * Copyright 2008, 2009, 2010 LSST Corporation.
4  *
5  * This product includes software developed by the
6  * LSST Project (http://www.lsst.org/).
7  *
8  * This program is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the LSST License Statement and
19  * the GNU General Public License along with this program. If not,
20  * see <http://www.lsstcorp.org/LegalNotices/>.
21  */
22 
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 }
34 #endif
35 using namespace posix;
36 #include <cstdlib>
37 #include <cstdio>
38 #include <cstring>
39 #include <cassert>
40 #include <cctype>
41 
42 #include "lsst/pex/exceptions.h"
43 #include "boost/any.hpp"
44 #include "lsst/afw/fits.h"
45 
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() {}
70 
71  int write(int fd, int ncard, char *record) const;
72 
73  std::string keyword;
74  boost::any value;
75  std::string comment;
76 };
77 
78 /*****************************************************************************/
79 /*
80  * Write a Card
81  */
82 int Card::write(int fd,
83  int ncard,
84  char *record
85  ) const {
86  char *card = &record[80*ncard];
87 
88  if (value.type() == typeid(std::string)) {
89  const char *str = boost::any_cast<std::string>(value).c_str();
90  if (keyword == "" ||
91  keyword == "COMMENT" || keyword == "END" || keyword == "HISTORY") {
92  sprintf(card, "%-8.8s%-72s", keyword.c_str(), str);
93  } else {
94  sprintf(card,"%-8.8s= '%s' %c%-*s",
95  keyword.c_str(), str,
96  (comment == "" ? ' ' : '/'),
97  (int)(80 - 14 - strlen(str)), comment.c_str());
98  }
99  } else {
100  sprintf(card, "%-8.8s= ", keyword.c_str());
101  card += 10;
102  if (value.type() == typeid(bool)) {
103  sprintf(card, "%20s", boost::any_cast<bool>(value) ? "T" : "F");
104  } else if (value.type() == typeid(int)) {
105  sprintf(card, "%20d", boost::any_cast<int>(value));
106  } else if (value.type() == typeid(double)) {
107  sprintf(card, "%20.10f", boost::any_cast<double>(value));
108  } else if (value.type() == typeid(float)) {
109  sprintf(card, "%20.7f", boost::any_cast<float>(value));
110  }
111  card += 20;
112  sprintf(card, " %c%-48s", (comment == "" ? ' ' : '/'), comment.c_str());
113  }
114 /*
115  * Write record if full
116  */
117  if (++ncard == 36) {
118  if (posix::write(fd, record, FITS_SIZE) != FITS_SIZE) {
119  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError, "Cannot write header record");
120  }
121  ncard = 0;
122  }
123 
124  return ncard;
125 }
127 
128 /*****************************************************************************/
129 /*
130  * Utilities
131  *
132  * Flip high-order bit so as to write unsigned short to FITS. Grrr.
133  */
134 namespace {
135  void flip_high_bit(char *arr, // array that needs bits swapped
136  const int n) { // number of bytes in arr
137  if (n%2 != 0) {
138  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
139  (boost::format("Attempt to bit flip odd number of bytes: %d") % n).str());
140  }
141 
142  unsigned short* uarr = reinterpret_cast<unsigned short *>(arr);
143  for(unsigned short *end = uarr + n/2; uarr < end; ++uarr) {
144  *uarr ^= 0x8000;
145  }
146  }
147 }
148 
149 /*
150  * Byte swap ABABAB -> BABABAB in place
151  */
152 namespace {
153  void swap_2(char *arr, // array to swap
154  const int n) { // number of bytes
155  if (n%2 != 0) {
156  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
157  (boost::format("Attempt to byte swap odd number of bytes: %d") % n).str());
158  }
159 
160  for(char *end = arr + n;arr < end;arr += 2) {
161  char t = arr[0];
162  arr[0] = arr[1];
163  arr[1] = t;
164  }
165  }
166  /*
167  * Byte swap ABCDABCD -> DCBADCBA in place (e.g. sun <--> vax)
168  */
169  void swap_4(char *arr, // array to swap
170  const int n) { // number of bytes
171  if (n%4 != 0) {
172  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
173  (boost::format("Attempt to byte swap non-multiple of 4 bytes: %d") % n).str());
174  }
175 
176  for(char *end = arr + n;arr < end;arr += 4) {
177  char t = arr[0];
178  arr[0] = arr[3];
179  arr[3] = t;
180  t = arr[1];
181  arr[1] = arr[2];
182  arr[2] = t;
183  }
184  }
185 
186  /*
187  * Byte swap ABCDEFGH -> HGFEDCBA in place (e.g. sun <--> vax)
188  */
189  void swap_8(char *arr, // array to swap
190  const int n) { // number of bytes
191  if (n%8 != 0) {
192  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
193  (boost::format("Attempt to byte swap non-multiple of 8 bytes: %d") % n).str());
194  }
195 
196  for(char *end = arr + n;arr < end;arr += 8) {
197  char t = arr[0];
198  arr[0] = arr[7];
199  arr[7] = t;
200  t = arr[1];
201  arr[1] = arr[6];
202  arr[6] = t;
203  t = arr[2];
204  arr[2] = arr[5];
205  arr[5] = t;
206  t = arr[3];
207  arr[3] = arr[4];
208  arr[4] = t;
209  }
210  }
211 
212 /*****************************************************************************/
213 
214  int write_fits_hdr(int fd,
215  int bitpix,
216  int naxis,
217  int *naxes,
218  std::list<Card>& cards, /* extra header cards */
219  int primary) /* is this the primary HDU? */
220  {
221  int i;
222  char record[FITS_SIZE + 1]; /* write buffer */
223 
224  int ncard = 0;
225  if (primary) {
226  Card card("SIMPLE", true);
227  ncard = card.write(fd, ncard, record);
228  } else {
229  Card card("XTENSION", "IMAGE");
230  ncard = card.write(fd, ncard, record);
231  }
232 
233  {
234  Card card("BITPIX", bitpix);
235  ncard = card.write(fd, ncard, record);
236  }
237  {
238  Card card("NAXIS", naxis);
239  ncard = card.write(fd, ncard, record);
240  }
241  for(i = 0; i < naxis; i++) {
242  char key[] = "NAXIS.";
243  sprintf(key, "NAXIS%d", i + 1);
244  Card card(key, naxes[i]);
245  ncard = card.write(fd, ncard, record);
246  }
247  if (primary) {
248  Card card("EXTEND", true, "There may be extensions");
249  ncard = card.write(fd,ncard,record);
250  }
251 /*
252  * Write extra header cards
253  */
254  for (std::list<Card>::const_iterator card = cards.begin(); card != cards.end(); card++) {
255  ncard = card->write(fd,ncard,record);
256  }
257 
258  {
259  Card card("END", "");
260  ncard = card.write(fd,ncard,record);
261  }
262  while(ncard != 0) {
263  Card card("", "");
264  ncard = card.write(fd,ncard,record);
265  }
266 
267  return 0;
268  }
269 
270 /*
271  * Pad out to a FITS record boundary
272  */
273  void pad_to_fits_record(int fd, // output file descriptor
274  int npixel, // number of pixels already written to HDU
275  int bitpix // bitpix for this datatype
276  ) {
277  const int bytes_per_pixel = (bitpix > 0 ? bitpix : -bitpix)/8;
278  int nbyte = npixel*bytes_per_pixel;
279 
280  if (nbyte%FITS_SIZE != 0) {
281  char record[FITS_SIZE + 1]; /* write buffer */
282 
283  nbyte = FITS_SIZE - nbyte%FITS_SIZE;
284  memset(record, ' ', nbyte);
285  if (write(fd, record, nbyte) != nbyte) {
286  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
287  "error padding file to multiple of fits block size");
288  }
289  }
290  }
291 
292  int write_fits_data(int fd,
293  int bitpix,
294  char *begin,
295  char *end
296  ) {
297  const int bytes_per_pixel = (bitpix > 0 ? bitpix : -bitpix)/8;
298  int swap_bytes = 0; // the default
299 #if defined(LSST_LITTLE_ENDIAN) // we'll need to byte swap FITS
300  if (bytes_per_pixel > 1) {
301  swap_bytes = 1;
302  }
303 #endif
304 
305  char *buff = NULL; // I/O buffer
306  bool allocated = false; // do I need to free it?
307  if (swap_bytes || bitpix == 16) {
308  buff = new char[FITS_SIZE*bytes_per_pixel];
309  allocated = true;
310  }
311 
312  int nbyte = end - begin;
313  int nwrite = (nbyte > FITS_SIZE) ? FITS_SIZE : nbyte;
314  for (char *ptr = begin; ptr != end; nbyte -= nwrite, ptr += nwrite) {
315  if (end - ptr < nwrite) {
316  nwrite = end - ptr;
317  }
318 
319  if (swap_bytes) {
320  memcpy(buff, ptr, nwrite);
321  if (bitpix == 16) { // flip high-order bit
322  flip_high_bit(buff, nwrite);
323  }
324 
325  if (bytes_per_pixel == 2) {
326  swap_2(buff, nwrite);
327  } else if (bytes_per_pixel == 4) {
328  swap_4(buff, nwrite);
329  } else if (bytes_per_pixel == 8) {
330  swap_8(buff, nwrite);
331  } else {
332  fprintf(stderr,"You cannot get here\n");
333  abort();
334  }
335  } else {
336  if (bitpix == 16) { // flip high-order bit
337  memcpy(buff, ptr, nwrite);
338  flip_high_bit(buff, nwrite);
339  } else {
340  buff = ptr;
341  }
342  }
343 
344  if (write(fd, buff, nwrite) != nwrite) {
345  perror("Error writing image: ");
346  break;
347  }
348  }
349 
350  if (allocated) {
351  delete buff;
352  }
353 
354  return (nbyte == 0 ? 0 : -1);
355  }
356 }
357 
358 namespace lsst { namespace afw { namespace display {
359 
360 template<typename ImageT>
361 void writeBasicFits(int fd, // file descriptor to write to
362  ImageT const& data, // The data to write
363  image::Wcs const* Wcs, // which Wcs to use for pixel
364  char const *title // title to write to DS9
365  ) {
366  /*
367  * Allocate cards for FITS headers
368  */
369  std::list<Card> cards;
370  /*
371  * What sort if image is it?
372  */
373  int bitpix = lsst::afw::fits::getBitPix<typename ImageT::Pixel>();
374  if (bitpix == 20) { // cfitsio for "Unsigned short"
375  cards.push_back(Card("BZERO", 32768.0, ""));
376  cards.push_back(Card("BSCALE", 1.0, ""));
377  bitpix = 16;
378  } else if (bitpix == 0) {
379  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError, "Unsupported image type");
380  }
381  /*
382  * Generate WcsA, pixel coordinates, allowing for X0 and Y0
383  */
384  std::string wcsName = "A";
385  cards.push_back(Card(str(boost::format("CRVAL1%s") % wcsName),
386  data.getX0(), "(output) Column pixel of Reference Pixel"));
387  cards.push_back(Card(str(boost::format("CRVAL2%s") % wcsName),
388  data.getY0(), "(output) Row pixel of Reference Pixel"));
389  cards.push_back(Card(str(boost::format("CRPIX1%s") % wcsName), 1.0,
390  "Column Pixel Coordinate of Reference"));
391  cards.push_back(Card(str(boost::format("CRPIX2%s") % wcsName), 1.0, "Row Pixel Coordinate of Reference"));
392  cards.push_back(Card(str(boost::format("CTYPE1%s") % wcsName), "LINEAR", "Type of projection"));
393  cards.push_back(Card(str(boost::format("CTYPE1%s") % wcsName), "LINEAR", "Type of projection"));
394  cards.push_back(Card(str(boost::format("CUNIT1%s") % wcsName), "PIXEL", "Column unit"));
395  cards.push_back(Card(str(boost::format("CUNIT2%s") % wcsName), "PIXEL", "Row unit"));
396  /*
397  * Now WcsB, so that pixel (0,0) is correctly labelled (but ignoring XY0)
398  */
399  wcsName = "B";
400  cards.push_back(Card(str(boost::format("CRVAL1%s") % wcsName), 0,
401  "(output) Column pixel of Reference Pixel"));
402  cards.push_back(Card(str(boost::format("CRVAL2%s") % wcsName), 0,
403  "(output) Row pixel of Reference Pixel"));
404  cards.push_back(Card(str(boost::format("CRPIX1%s") % wcsName), 1.0,
405  "Column Pixel Coordinate of Reference"));
406  cards.push_back(Card(str(boost::format("CRPIX2%s") % wcsName), 1.0, "Row Pixel Coordinate of Reference"));
407  cards.push_back(Card(str(boost::format("CTYPE1%s") % wcsName), "LINEAR", "Type of projection"));
408  cards.push_back(Card(str(boost::format("CTYPE1%s") % wcsName), "LINEAR", "Type of projection"));
409  cards.push_back(Card(str(boost::format("CUNIT1%s") % wcsName), "PIXEL", "Column unit"));
410  cards.push_back(Card(str(boost::format("CUNIT2%s") % wcsName), "PIXEL", "Row unit"));
411 
412  if (title) {
413  cards.push_back(Card("OBJECT", title, "Image being displayed"));
414  }
415  /*
416  * Was there something else?
417  */
418  if (Wcs != NULL) {
419  typedef std::vector<std::string> NameList;
420 
421  image::Wcs::Ptr newWcs = Wcs->clone(); //Create a copy
422  newWcs->shiftReferencePixel(-data.getX0(), -data.getY0());
423 
424  lsst::daf::base::PropertySet::Ptr metadata = newWcs->getFitsMetadata();
425 
426  NameList paramNames = metadata->paramNames();
427 
428  for (NameList::const_iterator i = paramNames.begin(), end = paramNames.end(); i != end; ++i) {
429  if (*i == "SIMPLE" ||
430  *i == "BITPIX" ||
431  *i == "NAXIS" ||
432  *i == "NAXIS1" ||
433  *i == "NAXIS2" ||
434  *i == "XTENSION" ||
435  *i == "PCOUNT" ||
436  *i == "GCOUNT"
437  ) {
438  continue;
439  }
440  std::type_info const &type = metadata->typeOf(*i);
441  if (type == typeid(bool)) {
442  cards.push_back(Card(*i, metadata->get<bool>(*i)));
443  } else if (type == typeid(int)) {
444  cards.push_back(Card(*i, metadata->get<int>(*i)));
445  } else if (type == typeid(float)) {
446  cards.push_back(Card(*i, metadata->get<float>(*i)));
447  } else if (type == typeid(double)) {
448  cards.push_back(Card(*i, metadata->get<double>(*i)));
449  } else {
450  cards.push_back(Card(*i, metadata->get<std::string>(*i)));
451  }
452  }
453  }
454  /*
455  * Basic FITS stuff
456  */
457  const int naxis = 2; // == NAXIS
458  int naxes[naxis]; /* values of NAXIS1 etc */
459  naxes[0] = data.getWidth();
460  naxes[1] = data.getHeight();
461 
462  write_fits_hdr(fd, bitpix, naxis, naxes, cards, 1);
463  for (int y = 0; y != data.getHeight(); ++y) {
464  if (write_fits_data(fd, bitpix, (char *)(data.row_begin(y)), (char *)(data.row_end(y))) < 0){
465  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
466  (boost::format("Error writing data for row %d") % y).str());
467  }
468  }
469 
470  pad_to_fits_record(fd, data.getWidth()*data.getHeight(), bitpix);
471 }
472 
473 /******************************************************************************/
474 
475 template<typename ImageT>
476 void writeBasicFits(std::string const& filename, // file to write, or "| cmd"
477  ImageT const& data, // The data to write
478  image::Wcs const* Wcs, // which Wcs to use for pixel
479  char const* title // title to write to DS9
480  ) {
481  int fd;
482  if ((filename.c_str())[0] == '|') { // a command
483  const char *cmd = filename.c_str() + 1;
484  while (isspace(*cmd)) {
485  cmd++;
486  }
487 
488  fd = fileno(popen(cmd, "w"));
489  } else {
490  fd = creat(filename.c_str(), 777);
491  }
492 
493  if (fd < 0) {
494  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError,
495  (boost::format("Cannot open \"%s\"") % filename).str());
496  }
497 
498  try {
499  writeBasicFits(fd, data, Wcs, title);
501  (void)close(fd);
502  throw;
503  }
504 
505  (void)close(fd);
506 }
507 
509 #define INSTANTIATE(IMAGET) \
510  template void writeBasicFits(int, IMAGET const&, image::Wcs const *, char const *); \
511  template void writeBasicFits(std::string const&, IMAGET const&, image::Wcs const *, char const *)
512 
513 #define INSTANTIATE_IMAGE(T) INSTANTIATE(lsst::afw::image::Image<T>)
514 #define INSTANTIATE_MASK(T) INSTANTIATE(lsst::afw::image::Mask<T>)
515 
516 INSTANTIATE_IMAGE(boost::uint16_t);
517 INSTANTIATE_IMAGE(int);
518 INSTANTIATE_IMAGE(float);
519 INSTANTIATE_IMAGE(double);
520 INSTANTIATE_IMAGE(boost::uint64_t);
521 
522 INSTANTIATE_MASK(boost::uint16_t);
524 
525 }}}
int y
table::Key< std::string > name
Definition: ApCorrMap.cc:71
void writeBasicFits(int fd, ImageT const &data, image::Wcs const *Wcs, char const *title)
Definition: simpleFits.cc:361
boost::shared_ptr< Wcs > Ptr
Definition: Wcs.h:113
Definitions to write a FITS image.
boost::shared_ptr< PropertySet > Ptr
Definition: PropertySet.h:90
virtual Ptr clone(void) const
Definition: Wcs.cc:552
Implementation of the WCS standard for a any projection.
Definition: Wcs.h:107
#define INSTANTIATE_IMAGE(IMAGE)
Definition: SdssShape.cc:884
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
lsst::daf::base::PropertySet PropertySet
Definition: Wcs.cc:58
bool any(CoordinateExpr< N > const &expr)
Return true if any elements are true.
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
#define FITS_SIZE
Definition: simpleFits.cc:51
ImageT val
Definition: CR.cc:154
Include files required for standard LSST Exception handling.