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
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
30namespace 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
35using 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
48namespace image = lsst::afw::image;
50
51#define FITS_SIZE 2880
52
54class Card {
55public:
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 */
81int 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 */
126namespace {
127void 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 */
144namespace {
145void 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 */
161void 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 */
181void 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
204int 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 */
259void 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
278int 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
339void 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
351namespace lsst {
352namespace afw {
353namespace display {
354
355template <typename ImageT>
356void 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
441template <typename ImageT>
442void 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
484INSTANTIATE_IMAGE(float);
485INSTANTIATE_IMAGE(double);
487
488INSTANTIATE_MASK(std::uint16_t);
489INSTANTIATE_MASK(image::MaskPixel);
491} // namespace display
492} // namespace afw
493} // namespace lsst
char * data
Definition BaseRecord.cc:61
int end
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
uint64_t * ptr
Definition RangeSet.cc:95
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
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)
T get(T... args)
#define INSTANTIATE_IMAGE(IMAGE)
Definition ImagePca.cc:125
T memcpy(T... args)
T memset(T... args)
void writeBasicFits(int fd, ImageT const &data, geom::SkyWcs const *Wcs, char const *title)
Extent< double, 2 > Extent2D
Definition Extent.h:400
T perror(T... args)
#define FITS_SIZE
Definition simpleFits.cc:51
ImageT val
Definition CR.cc:146