Line data Source code
1 : /** \file fitsFile.hpp
2 : * \brief Declares and defines a class to work with a FITS file
3 : * \ingroup fits_processing_files
4 : * \author Jared R. Males (jaredmales@gmail.com)
5 : *
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 2015-2022 Jared R. Males (jaredmales@gmail.com)
10 : //
11 : // This file is part of mxlib.
12 : //
13 : // mxlib is free software: you can redistribute it and/or modify
14 : // it under the terms of the GNU General Public License as published by
15 : // the Free Software Foundation, either version 3 of the License, or
16 : // (at your option) any later version.
17 : //
18 : // mxlib is distributed in the hope that it will be useful,
19 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
20 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 : // GNU General Public License for more details.
22 : //
23 : // You should have received a copy of the GNU General Public License
24 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
25 : //***********************************************************************//
26 :
27 : #ifndef ioutils_fits_fitsFile_hpp
28 : #define ioutils_fits_fitsFile_hpp
29 :
30 : #include <memory>
31 :
32 : #include "../../mxlib.hpp"
33 :
34 : #include "../../improc/eigenImage.hpp"
35 :
36 : #include "fitsUtils.hpp"
37 : #include "fitsHeader.hpp"
38 :
39 : namespace mx
40 : {
41 : namespace fits
42 : {
43 :
44 : namespace fitsFileDetail
45 : {
46 :
47 : /** \cond */
48 : /// Resettable CFITSIO operation table used by fitsFile and its deterministic failure tests.
49 : struct fitsFileCfitsioOps
50 : {
51 : int ( *openFile )( fitsfile **, const char *, int, int * ); ///< Open an existing FITS file.
52 : decltype( &ffgidm ) getImageDim; ///< Read the number of image axes.
53 : decltype( &ffgisz ) getImageSize; ///< Read image-axis sizes.
54 : decltype( &ffclos ) closeFile; ///< Close a FITS file.
55 : decltype( &ffgsv ) readSubset; ///< Read a FITS pixel subset.
56 : decltype( &ffghsp ) getHeaderSpace; ///< Read the number of FITS header cards.
57 : decltype( &ffgkyn ) readKey; ///< Read one FITS header card.
58 : decltype( &ffinit ) createFile; ///< Create a FITS file.
59 : decltype( &ffcrim ) createImage; ///< Create the primary FITS image.
60 : decltype( &ffppx ) writePixels; ///< Write FITS pixels.
61 : long *( *allocateLongs )( size_t ); ///< Allocate a pixel-coordinate array.
62 : };
63 :
64 : /// Access the process-wide CFITSIO operation table.
65 : fitsFileCfitsioOps &fitsFileCfitsioOpsInstance();
66 :
67 : /// Restore every CFITSIO operation to its production implementation.
68 : void resetFitsFileCfitsioOps();
69 : /** \endcond */
70 :
71 : } // namespace fitsFileDetail
72 :
73 : /// Class to manage interactions with a FITS file
74 : /** This class wraps the functionality of cfitsio.
75 : *
76 : * \tparam dataT is the datatype to use for in-memory storage of the image.
77 : * This does not have to match the data type
78 : * stored on disk for reading, but will be the type used for writing.
79 : *
80 : * \ingroup fits_processing
81 : */
82 : template <typename dataT, class verboseT = verbose::d>
83 : class fitsFile
84 : {
85 :
86 : friend class fitsFile_test;
87 :
88 : public:
89 : typedef typename fitsHeader<verboseT>::headerIteratorT headerIteratorT;
90 :
91 : protected:
92 : /// The path to the file
93 : std::string m_fileName;
94 :
95 : /// The cfitsio data structure
96 : fitsfile *m_fptr{ nullptr };
97 :
98 : /// The dimensions of the image (1D, 2D, or 3D)
99 : int m_naxis{ 0 };
100 :
101 : /// The size of each dimension
102 : long m_naxes[3];
103 :
104 : /// Flag indicating whether the file is open or not
105 : bool m_isOpen{ false };
106 :
107 : /// The value to replace null values with
108 : dataT m_nulval{ 0 };
109 :
110 : /// Records whether any null values were replaced
111 : int m_anynul{ 0 };
112 :
113 : /// Flag to control whether the comment string is read.
114 : int m_noComment{ 0 };
115 :
116 : /// The starting x-pixel to read from
117 : long m_x0{ -1 };
118 :
119 : /// The starting y-pixel to read from
120 : long m_y0{ -1 };
121 :
122 : /// The number of x-pixels to read
123 : long m_xpix{ -1 };
124 :
125 : /// The number of y-pixels to read
126 : long m_ypix{ -1 };
127 :
128 : /// The starting frame to read from a cube
129 : long m_z0{ -1 };
130 :
131 : /// The number of frames to read from a cube
132 : long m_zframes{ -1 };
133 :
134 : /// One time initialization common to all constructors
135 : void construct();
136 :
137 : public:
138 : /// Default constructor
139 : fitsFile();
140 :
141 : /// Default constructor with error code
142 : fitsFile( error_t &errc /**< [out] error_t code indicating success or error */ );
143 :
144 : /// Constructor with file name and error code
145 : /** The file is not opened.
146 : *
147 : */
148 : fitsFile( const std::string &fname, ///< [in] File name to set on construction
149 : error_t &errc /**< [out] error_t code indicating success or error */
150 : );
151 :
152 : /// Constructor with file name, and option to open.
153 : fitsFile( const std::string &fname, ///< [in] File name to set on construction
154 : bool doopen = true ///< [in] [optional] If true, then the file is opened (the default).
155 : );
156 :
157 : /// Constructor with file name, option to open, and error code
158 : fitsFile( const std::string &fname, ///< File name to set on construction
159 : bool doopen, ///< If true, then the file is opened (the default).
160 : error_t &errc /**< [out] error_t code indicating success or error */
161 : );
162 :
163 : /// Destructor
164 : ~fitsFile();
165 :
166 : /// Get the current value of m_fileName
167 : /**
168 : * \returns the current value of m_fileName.
169 : */
170 : std::string fileName();
171 :
172 : /// Set the file path, and optionally open the file.
173 : /**
174 : * \returns mx::error_t::noerror on success
175 : * \returns mx::error_t codes from \ref close or \ref open
176 : */
177 : error_t fileName( const std::string &fname, ///< The new file name.
178 : bool doopen = true ///< If true, then the file is opened (the default).
179 : );
180 :
181 : /// Get the current value of m_naxis
182 : /**
183 : * \returns the current value of m_naxis
184 : *
185 : */
186 : int naxis();
187 :
188 : /// Get the current value of m_naxes for the specified dimension
189 : /**
190 : * \returns the current value of m_naxes for the specified dimension. -1 if no such dimension
191 : *
192 : */
193 : long naxes( int dim /**< [in] the dimension */ );
194 :
195 : /// Open the file and gets its dimensions.
196 : /** File name needs to already have been set.
197 : * If the file has already been opened, this returns immediately with no re-open.
198 : *
199 : * \returns mx::error_t::noerror on success
200 : * \returns mx::invalidconfig if m_filename is empty.
201 : * \returns mx::error_t::bad_alloc on an allocation error
202 : * \returns mx::error_t::std_exception on other errors during allocations
203 : * \returns mx::error_t::exception on unexpected errors during allocations
204 : * \returns mx::error_t::allocerr if allocation fails without an exception
205 : * \returns mx::error_t::fits_* codes from cfitsio functions
206 : *
207 : */
208 : error_t open();
209 :
210 : /// Open the file, first setting the file path.
211 : /**
212 : * \returns mx::error_t::noerror on success
213 : * \returns mx::invalidconfig if m_filename is empty.
214 : * \returns mx::error_t::bad_alloc on an allocation error
215 : * \returns mx::error_t::std_exception on other errors during allocations
216 : * \returns mx::error_t::exception on unexpected errors during allocations
217 : * \returns mx::error_t::allocerr if allocation fails without an exception
218 : * \returns mx::error_t::fits_* codes from cfitsio functions
219 : */
220 : error_t open( const std::string &fname /**< The name of the file to open. */ );
221 :
222 : /// Close the file.
223 : /**
224 : * \returns mx::error_t::noerror on success
225 : * \returns mx::error_t::fits_* codes from cfitsio functions
226 : */
227 : error_t close();
228 :
229 : /// Get the number of dimensions (i.e. m_naxis)
230 : int getDimensions( error_t &errc );
231 :
232 : /// Get the total size
233 : long getSize( error_t &errc );
234 :
235 : /// Get the size of a specific dimension
236 : long getSize( size_t axis, error_t &errc );
237 :
238 : /** \name Reading Basic Arrays
239 : * These methods read FITS data into basic or raw arrays specified by a pointer.
240 : * @{
241 : */
242 :
243 : protected:
244 : struct pixarrT
245 : {
246 : long *fpix{ nullptr }; ///< Populated with the lower left pixel to read.
247 : long *lpix{ nullptr }; ///< Populated with the upper right pixel to read.
248 : long *inc{ nullptr }; ///< The increment.
249 :
250 53 : error_t allocate( int naxis )
251 : {
252 53 : if( naxis <= 0 )
253 : {
254 3 : return internal::mxlib_error_report<verboseT>( error_t::paramnotset, "naxis" );
255 : }
256 :
257 52 : if( fpix )
258 : {
259 1 : delete[] fpix;
260 1 : fpix = nullptr;
261 : }
262 :
263 52 : if( lpix )
264 : {
265 1 : delete[] lpix;
266 1 : lpix = nullptr;
267 : }
268 :
269 52 : if( inc )
270 : {
271 1 : delete[] inc;
272 1 : inc = nullptr;
273 : }
274 :
275 : try
276 : {
277 52 : fpix = fitsFileDetail::fitsFileCfitsioOpsInstance().allocateLongs( naxis );
278 49 : lpix = fitsFileDetail::fitsFileCfitsioOpsInstance().allocateLongs( naxis );
279 49 : inc = fitsFileDetail::fitsFileCfitsioOpsInstance().allocateLongs( naxis );
280 : }
281 4 : catch( const std::bad_alloc &e )
282 : {
283 1 : internal::mxlib_error_report<verboseT>( error_t::std_bad_alloc,
284 3 : std::string( "allocating pixel read arrays: " ) + e.what() );
285 : #ifdef MXLIB_TRAP_ALLOC_ERRORS
286 : return error_t::std_bad_alloc;
287 : #else
288 1 : throw;
289 : #endif
290 : }
291 2 : catch( const std::exception &e )
292 : {
293 1 : internal::mxlib_error_report<verboseT>( error_t::std_exception,
294 3 : std::string( "allocating pixel read arrays: " ) + e.what() );
295 : #ifdef MXLIB_TRAP_ALLOC_ERRORS
296 : return error_t::exception;
297 : #else
298 1 : throw;
299 : #endif
300 : }
301 2 : catch( ... )
302 : {
303 2 : internal::mxlib_error_report<verboseT>( error_t::exception, "allocating pixel read arrays" );
304 : #ifdef MXLIB_TRAP_ALLOC_ERRORS
305 : return error_t::exception;
306 : #else
307 1 : throw;
308 : #endif
309 : }
310 :
311 : // also check for null allocations (in case compiled with exceptions off)
312 49 : if( fpix == nullptr )
313 : {
314 9 : return internal::mxlib_error_report<verboseT>( error_t::allocerr, "fpix is nullptr" );
315 : }
316 :
317 46 : if( lpix == nullptr )
318 : {
319 3 : return internal::mxlib_error_report<verboseT>( error_t::allocerr, "lpix is nullptr" );
320 : }
321 :
322 45 : if( inc == nullptr )
323 : {
324 3 : return internal::mxlib_error_report<verboseT>( error_t::allocerr, "inc is nullptr" );
325 : }
326 :
327 44 : return error_t::noerror;
328 : }
329 :
330 51 : ~pixarrT()
331 : {
332 51 : if( fpix )
333 : {
334 45 : delete[] fpix;
335 : }
336 :
337 51 : if( lpix )
338 : {
339 47 : delete[] lpix;
340 : }
341 :
342 51 : if( inc )
343 : {
344 47 : delete[] inc;
345 : }
346 51 : }
347 : };
348 :
349 : /// Fill in the read-size arrays for reading a subset (always used)
350 : /**
351 : */
352 : error_t calcPixarrs( pixarrT &pixarr /**< [out] Populated with the allocated read-size arrays*/ );
353 :
354 : public:
355 : /// Read the contents of the FITS file into an array.
356 : /** The array pointed to by data must have been allocated.
357 : *
358 : * \returns 0 on success
359 : * \returns -1 on error
360 : */
361 : error_t read( dataT *data /**< [out] an allocated arrray large enough to hold the entire image */ );
362 :
363 : /// Read the contents of the FITS file into an array.
364 : /** The array pointed to by data must have been allocated.
365 : *
366 : * \returns 0 on success
367 : * \returns -1 on error
368 : */
369 : error_t read( dataT *data, ///< [out] an allocated arrray large enough to hold the entire image
370 : fitsHeader<verboseT> &head ///< [out] a fitsHeader object which is passed to \ref readHeader
371 : );
372 :
373 : /// Read the contents of the FITS file into an array.
374 : /** The array pointed to by data must have been allocated.
375 : *
376 : * \returns 0 on success
377 : * \returns -1 on error
378 : */
379 : error_t read( dataT *data, ///< [out] is an allocated arrray large enough to hold the entire image
380 : const std::string &fname ///< [in] is the file path, which is passed to \ref fileName
381 : );
382 :
383 : /// Read the contents of the FITS file into an array and read the header.
384 : /** The array pointed to by data must have been allocated.
385 : *
386 : * \returns 0 on success
387 : * \returns -1 on error
388 : */
389 : error_t read( dataT *data, ///< [out] an allocated arrray large enough to hold the entire image
390 : fitsHeader<verboseT> &head, ///< [out] a fitsHeader object which is passed to \ref readHeader
391 : const std::string &fname ///< [in] the file path, which is passed to \ref fileName
392 : );
393 :
394 : /// Read data from a vector list of files into an image cube
395 : error_t read( dataT *im, ///< [out] An allocated array large enough to hold all the images
396 : const std::vector<std::string> &flist ///< [in] The list of files to read.
397 : );
398 :
399 : /// Read data from a vector of files into an image cube with individual headers
400 : error_t read( dataT *im, /**< [out] An allocated array large enough to hold all the images */
401 : std::vector<fitsHeader<verboseT>> &heads, /**< [in/out] The vector of fits headers, allocated
402 : to contain one per image. */
403 : const std::vector<std::string> &flist /**< [in] The list of files to read. */
404 : );
405 :
406 : ///@}
407 :
408 : /** \name Reading Eigen Arrays
409 : * These methods read FITS data into array types with an Eigen-like interface.
410 : * @{
411 : */
412 :
413 : /// Read the contents of the FITS file into an Eigen array type (not a simple pointer).
414 : /** The type arrT can be any type with the following members defined:
415 : * - resize(int, int) (allocates memory)
416 : * - data() (returns a pointer to the underlying array)
417 : * - typedef arrT::Scalar (is the data type, does not have to be dataT)
418 : *
419 : * \tparam arrT is the type of array, see requirements above.
420 : *
421 : * \returns 0 on success
422 : * \returns -1 on error
423 : */
424 : template <typename arrT>
425 : error_t read( arrT &data /**< [out] is the array, which will be resized as
426 : necessary using its resize(int, int) member */ );
427 :
428 : /// Read the contents of the FITS file into an Eigen array type (not a simple pointer).
429 : /** The type arrT can be any type with the following members defined:
430 : * - resize(int, int) (allocates memory)
431 : * - data() (returns a pointer to the underlying array)
432 : * - typedef arrT::Scalar (is the data type, does not have to be dataT)
433 : *
434 : * \tparam arrT is the type of array, see requirements above.
435 : *
436 : * \returns 0 on success
437 : * \returns -1 on error
438 : */
439 : template <typename arrT>
440 : error_t read( arrT &data, /**< [out] is the array, which will be resized as necessary
441 : using its resize(int, int) member*/
442 : fitsHeader<verboseT> &head ///< [out] is a fitsHeader object which is passed to \ref readHeader
443 : );
444 :
445 : /// Read the contents of the FITS file into an Eigen array type (not a simple pointer).
446 : /** The type arrT can be any type with the following members defined:
447 : * - resize(int, int) (allocates memory)
448 : * - data() (returns a pointer to the underlying array)
449 : * - typedef arrT::Scalar (is the data type, does not have to be dataT)
450 : *
451 : * \tparam arrT is the type of array, see requirements above.
452 : *
453 : * \returns 0 on success
454 : * \returns -1 on error
455 : */
456 : template <typename arrT>
457 : error_t read( arrT &data, /**< [out] is the array, which will be resized as
458 : necessary using its resize(int, int) member*/
459 : const std::string &fname ///< [in] is the file path, which is passed to \ref fileName
460 : );
461 :
462 : /// Read the contents of the FITS file into an Eigen array type (not a simple pointer).
463 : /** The type arrT can be any type with the following members defined:
464 : * - resize(int, int) (allocates memory)
465 : * - data() (returns a pointer to the underlying array)
466 : * - typedef arrT::Scalar (is the data type, does not have to be dataT)
467 : *
468 : * \tparam arrT is the type of array, see requirements above.
469 : *
470 : * \returns 0 on success
471 : * \returns -1 on error
472 : */
473 : template <typename arrT>
474 : error_t read( arrT &data, /**< [out] the array, which will be resized as
475 : necessary using its resize(int, int) member*/
476 : fitsHeader<verboseT> &head, ///< [out] a fitsHeader object which is passed to \ref readHeader
477 : const std::string &fname ///< [in] the file path, which is passed to \ref fileName
478 : );
479 :
480 : /// Read data from a vector list of files into an image cube
481 : /** The type cubeT can be any type with the following members defined:
482 : * - resize(int, int, int) (allocates memory)
483 : * - data() (returns a pointer to the underlying array)
484 : * - typedef cubeT::Scalar (is the data type, does not have to be dataT)
485 : *
486 : * \note The images must all be the same size, or all be as large or larger than the subset specified.
487 : *
488 : * \tparam cubeT is the type of array, see requirements above.
489 : *
490 : * \returns 0 on success
491 : * \returns -1 on error
492 : */
493 : template <typename cubeT>
494 : error_t read( cubeT &cube, ///< [out] A cube which will be resized using its resize(int, int, int) member.
495 : const std::vector<std::string> &flist, ///< [in] The list of files to read.
496 : std::vector<fitsHeader<verboseT>> *heads =
497 : 0 ///< [out] [optional] A vector of fits headers, allocated to contain one per image.
498 : );
499 :
500 : /// Read data from a vector of files into an image cube with individual headers
501 : /** The type cubeT can be any type with the following members defined:
502 : * - resize(int, int, int) (allocates memory)
503 : * - data() (returns a pointer to the underlying array)
504 : * - typedef cubeT::Scalar (is the data type, does not have to be dataT)
505 : *
506 : * \tparam cubeT is the type of array, see requirements above.
507 : *
508 : * \returns 0 on success
509 : * \returns -1 on error
510 : */
511 : template <typename cubeT>
512 : error_t read( cubeT &cube, /**< [out] A cube which will be resized
513 : using its resize(int, int, int) member.*/
514 : std::vector<fitsHeader<verboseT>> &heads, /**< [out] The vector of fits headers, allocated to
515 : contain one per image.*/
516 : const std::vector<std::string> &flist ///< [in] The list of files to read.
517 : );
518 :
519 : ///@}
520 :
521 : /** \name Reading Headers
522 : * @{
523 : */
524 :
525 : /// Read the header from the fits file.
526 : /** If head is not empty, then only the keywords already in head are updated. Otherwise
527 : * the complete header is read.
528 : *
529 : * \returns 0 on success
530 : * \returns -1 on error
531 : */
532 : error_t readHeader( fitsHeader<verboseT> &head /**< [out] a fitsHeader object */ );
533 :
534 : /// Read the header from the fits file.
535 : /** If head is not empty, then only the keywords already in head are updated. Otherwise
536 : * the complete header is read.
537 : *
538 : * \returns 0 on success
539 : * \returns -1 on error
540 : */
541 : error_t readHeader( fitsHeader<verboseT> &head, ///< [out] a fitsHeader object
542 : const std::string &fname ///< [in] the file path, which is passed to \ref fileName
543 : );
544 :
545 : /// Read the headers from a list of FITS files.
546 : /** In each case, if the header is not empty then only the keywords already in head are updated. Otherwise
547 : * the complete header is read.
548 : *
549 : * \returns 0 on success
550 : * \returns -1 on error
551 : */
552 : error_t readHeader(
553 : std::vector<fitsHeader<verboseT>> &heads, /// A vector of fitsHeader objects to read into.
554 : const std::vector<std::string> &flist ///< [in] A list of files, each of which is passed to \ref fileName
555 : );
556 :
557 : ///@}
558 :
559 : /** \name Writing Basic Arrays
560 : * These methods write basic arrays specified by a pointer to FITS.
561 : * @{
562 : */
563 :
564 : /// Write the contents of a raw array to the FITS file.
565 : /**
566 : * \returns 0 on success
567 : * \returns -1 on error
568 : */
569 : error_t write( const dataT *im, ///< [in] is the array
570 : int d1, ///< [in] is the first dimension
571 : int d2, ///< [in] is the second dimension
572 : int d3, ///< [in] is the third dimenesion (minimum value is 1)
573 : fitsHeader<verboseT> *head ///< [in] a pointer to the header. Set to 0 if not used.
574 : );
575 :
576 : /// Write the contents of a raw array to the FITS file.
577 : /**
578 : * \returns 0 on success
579 : * \returns -1 on error
580 : */
581 : error_t write( const dataT *im, ///< [in] is the array
582 : int d1, ///< [in] is the first dimension
583 : int d2, ///< [in] is the second dimension
584 : int d3 ///< [in] is the third dimenesion (minimum value is 1)
585 : );
586 :
587 : /// Write the contents of a raw array to the FITS file.
588 : /**
589 : * Note: the type of the array must match dataT
590 : *
591 : * \returns 0 on success
592 : * \returns -1 on error
593 : */
594 : error_t write( const dataT *im, ///< [in] is the array
595 : int d1, ///< [in] is the first dimension
596 : int d2, ///< [in] is the second dimension
597 : int d3, ///< [in] is the third dimenesion (minimum value is 1)
598 : fitsHeader<verboseT> &head ///< [in] is the header
599 : );
600 :
601 : /// Write the contents of a raw array to the FITS file.
602 : /**
603 : * \returns 0 on success
604 : * \returns -1 on error
605 : */
606 : error_t write( const std::string &fname, ///< [in] is the name of the file.
607 : const dataT *im, ///< [in] is the array
608 : int d1, ///< [in] is the first dimension
609 : int d2, ///< [in] is the second dimension
610 : int d3 ///< [in] is the third dimenesion (minimum value is 1)
611 : );
612 :
613 : /// Write the contents of a raw array to the FITS file.
614 : /**
615 : * \returns 0 on success
616 : * \returns -1 on error
617 : */
618 : error_t write( const std::string &fname, ///< [in] is the name of the file.
619 : const dataT *im, ///< [in] is the array
620 : int d1, ///< [in] is the first dimension
621 : int d2, ///< [in] is the second dimension
622 : int d3, ///< [in] is the third dimenesion (minimum value is 1)
623 : fitsHeader<verboseT> &head ///< [in] is the header
624 : );
625 :
626 : ///@}
627 :
628 : /** \name Writing Eigen Arrays
629 : * These methods write array types with an Eigen-like interface.
630 : * @{
631 : */
632 :
633 : /// Write the contents of an Eigen-type array to a FITS file.
634 : /** The type arrT can be any type with the following members defined:
635 : * - data() (returns a pointer to the underlying array)
636 : * - rows() (returrns the number of rows)
637 : * - cols() (returns the number of columns)
638 : * - may have planes() defined
639 : *
640 : * Note: as with all write methods, the Scalar type of the array must match dataT
641 : *
642 : * \tparam arrT is the type of array, see requirements above.
643 : *
644 : * \returns 0 on success
645 : * \returns -1 on error
646 : */
647 : template <typename arrT>
648 : error_t write( const std::string &fname, ///< [in] is the name of the file.
649 : const arrT &im ///< [in] is the array
650 : );
651 :
652 : /// Write the contents of an Eigen-type array to a FITS file.
653 : /** The type arrT can be any type with the following members defined:
654 : * - data() (returns a pointer to the underlying array)
655 : * - rows() (returrns the number of rows)
656 : * - cols() (returns the number of columns)
657 : * - may have planes() defined.
658 : *
659 : * Note: as with all write methods, the Scalar type of the array must match dataT
660 : *
661 : * \tparam arrT is the type of array, see requirements above.
662 : *
663 : * \returns 0 on success
664 : * \returns -1 on error
665 : */
666 : template <typename arrT>
667 : error_t write( const std::string &fname, ///< [in] is the file path, which is passed to \ref fileName
668 : const arrT &im, ///< [in] is the array
669 : fitsHeader<verboseT> &head ///< [in] is a fitsHeader object which is passed to \ref readHeader
670 : );
671 :
672 : // int writeHeader( fitsHeader<verboseT> &head );
673 :
674 : ///@}
675 :
676 : /** \name Reading Subsets
677 : * It is often desirable to read only a subset of an image or images into memory. These methods
678 : * allow you to specify this.
679 : * @{
680 : */
681 :
682 : /// Set to read all the pixels in the file
683 : void setReadSize();
684 :
685 : /// Set to read only a subset of the pixels in the file
686 : /**
687 : */
688 : void setReadSize( long x0, ///< is the starting x-pixel to read
689 : long y0, ///< is the starting y-pixel to read
690 : long xpix, ///< is the number of x-pixels to read
691 : long ypix ///< is the number of y-pixels to read
692 : );
693 :
694 : /// Set to read all frames from a cube.
695 : /**
696 : */
697 : void setCubeReadSize();
698 :
699 : /// Set the number of frames to read from a cube.
700 : /**
701 : *
702 : */
703 : void setCubeReadSize( long z0, ///< is the starting frame to read
704 : long zframes ///< is the number of frames to read
705 : );
706 :
707 : ///@}
708 :
709 : }; // fitsFile
710 :
711 : template <typename dataT, class verboseT>
712 70 : fitsFile<dataT, verboseT>::fitsFile()
713 : {
714 70 : }
715 :
716 : template <typename dataT, class verboseT>
717 1 : fitsFile<dataT, verboseT>::fitsFile( error_t &errc )
718 : {
719 1 : errc = error_t::noerror;
720 1 : }
721 :
722 : template <typename dataT, class verboseT>
723 1 : fitsFile<dataT, verboseT>::fitsFile( const std::string &fname, error_t &errc )
724 : {
725 : // no errors are actually possible
726 1 : errc = internal::mxlib_error_report<verboseT>( fileName( fname, false ) ); // nothing printed if noerror
727 1 : }
728 :
729 : template <typename dataT, class verboseT>
730 4 : fitsFile<dataT, verboseT>::fitsFile( const std::string &fname, bool doopen )
731 : {
732 : // If an error happens on open(), then m_open will be false and this will persist
733 : // so no need to throw an exception
734 4 : internal::mxlib_error_report<verboseT>( fileName( fname, doopen ) ); // nothing printed if noerror
735 4 : }
736 :
737 : template <typename dataT, class verboseT>
738 1 : fitsFile<dataT, verboseT>::fitsFile( const std::string &fname, bool doopen, error_t &errc )
739 : {
740 : // If an error happens on open(), then m_open will be false and this will persist
741 : // so no need to throw an exception
742 1 : errc = internal::mxlib_error_report<verboseT>( fileName( fname, doopen ) ); // nothing printed if noerror
743 1 : }
744 :
745 : template <typename dataT, class verboseT>
746 77 : fitsFile<dataT, verboseT>::~fitsFile()
747 : {
748 77 : if( m_isOpen )
749 : {
750 19 : internal::mxlib_error_report<verboseT>( close() );
751 : }
752 77 : }
753 :
754 : template <typename dataT, class verboseT>
755 3 : std::string fitsFile<dataT, verboseT>::fileName()
756 : {
757 3 : return m_fileName;
758 : }
759 :
760 : template <typename dataT, class verboseT>
761 124 : error_t fitsFile<dataT, verboseT>::fileName( const std::string &fname, bool doopen )
762 : {
763 124 : if( m_isOpen )
764 : {
765 37 : mxlib_error_check( close() );
766 : }
767 :
768 124 : m_fileName = fname;
769 :
770 124 : if( doopen )
771 : {
772 62 : mxlib_error_check( open() );
773 : }
774 :
775 116 : return error_t::noerror;
776 : }
777 :
778 : template <typename dataT, class verboseT>
779 7 : int fitsFile<dataT, verboseT>::naxis()
780 : {
781 7 : return m_naxis;
782 : }
783 :
784 : template <typename dataT, class verboseT>
785 18 : long fitsFile<dataT, verboseT>::naxes( int dim )
786 : {
787 18 : if( dim >= m_naxis || dim > 2 )
788 : {
789 6 : return -1;
790 : }
791 :
792 12 : return m_naxes[dim];
793 : }
794 :
795 : template <typename dataT, class verboseT>
796 70 : error_t fitsFile<dataT, verboseT>::open()
797 : {
798 70 : if( m_isOpen ) // no-op
799 : {
800 1 : return error_t::noerror;
801 : }
802 :
803 69 : if( m_fileName == "" )
804 : {
805 3 : return internal::mxlib_error_report<verboseT>( error_t::invalidconfig, "File name is not set" );
806 : }
807 :
808 68 : int fstatus = 0;
809 :
810 68 : fitsFileDetail::fitsFileCfitsioOpsInstance().openFile( &m_fptr, m_fileName.c_str(), READONLY, &fstatus );
811 :
812 68 : if( fstatus )
813 : {
814 7 : return internal::mxlib_error_report<verboseT>( fits_status2error_t( fstatus ), "Opening file " + m_fileName );
815 : }
816 :
817 61 : fstatus = 0;
818 61 : fitsFileDetail::fitsFileCfitsioOpsInstance().getImageDim( m_fptr, &m_naxis, &fstatus );
819 61 : if( fstatus )
820 : {
821 1 : int closeStatus = 0;
822 1 : ffclos( m_fptr, &closeStatus );
823 1 : m_fptr = nullptr;
824 1 : return internal::mxlib_error_report<verboseT>( fits_status2error_t( fstatus ),
825 2 : "Getting number of axes in file " + m_fileName );
826 : }
827 :
828 60 : fstatus = 0;
829 60 : fitsFileDetail::fitsFileCfitsioOpsInstance().getImageSize( m_fptr, m_naxis, m_naxes, &fstatus );
830 60 : if( fstatus )
831 : {
832 1 : int closeStatus = 0;
833 1 : ffclos( m_fptr, &closeStatus );
834 1 : m_fptr = nullptr;
835 1 : return internal::mxlib_error_report<verboseT>( fits_status2error_t( fstatus ),
836 2 : "Getting dimensions in file " + m_fileName );
837 : }
838 :
839 59 : m_isOpen = true; // Only set this after opening is complete.
840 :
841 59 : return error_t::noerror;
842 : }
843 :
844 : template <typename dataT, class verboseT>
845 6 : error_t fitsFile<dataT, verboseT>::open( const std::string &fname )
846 : {
847 6 : mxlib_error_return( fileName( fname, true ) );
848 : }
849 :
850 : template <typename dataT, class verboseT>
851 115 : error_t fitsFile<dataT, verboseT>::close()
852 : {
853 115 : if( !m_isOpen )
854 : {
855 2 : return error_t::noerror; // No error.
856 : }
857 :
858 113 : int fstatus = 0;
859 113 : fitsFileDetail::fitsFileCfitsioOpsInstance().closeFile( m_fptr, &fstatus );
860 :
861 113 : if( fstatus )
862 : {
863 2 : return internal::mxlib_error_report<verboseT>( fits_status2error_t( fstatus ), "Closing file " + m_fileName );
864 : }
865 :
866 111 : m_fptr = nullptr;
867 111 : m_isOpen = false;
868 :
869 111 : return error_t::noerror;
870 : }
871 :
872 : template <typename dataT, class verboseT>
873 3 : int fitsFile<dataT, verboseT>::getDimensions( error_t &errc )
874 : {
875 3 : if( !m_isOpen )
876 : {
877 1 : errc = error_t::invalidconfig;
878 1 : return -1;
879 : }
880 :
881 2 : return m_naxis;
882 : }
883 :
884 : template <typename dataT, class verboseT>
885 6 : long fitsFile<dataT, verboseT>::getSize( error_t &errc )
886 : {
887 6 : if( !m_isOpen )
888 : {
889 1 : errc = error_t::invalidconfig;
890 1 : return -1;
891 : }
892 :
893 5 : long sz = 1;
894 :
895 5 : errc = error_t::noerror;
896 :
897 5 : if( m_naxis == 1 )
898 : {
899 2 : if( m_x0 > -1 && m_xpix > -1 )
900 : {
901 1 : return m_xpix;
902 : }
903 1 : return m_naxes[0];
904 : }
905 3 : else if( m_x0 > -1 && m_y0 > -1 && m_xpix > -1 && m_ypix > -1 && m_naxis == 2 )
906 : {
907 1 : return m_xpix * m_ypix;
908 : }
909 : else
910 : {
911 6 : for( int i = 0; i < m_naxis && i < 3; ++i )
912 : {
913 4 : sz *= m_naxes[i];
914 : }
915 : }
916 :
917 2 : return sz;
918 : }
919 :
920 : template <typename dataT, class verboseT>
921 21 : long fitsFile<dataT, verboseT>::getSize( size_t axis, error_t &errc )
922 : {
923 21 : if( !m_isOpen )
924 : {
925 1 : errc = error_t::invalidconfig;
926 1 : return -1;
927 : }
928 :
929 20 : if( axis >= m_naxis || axis > 2 )
930 : {
931 1 : errc = error_t::invalidarg;
932 1 : return -1;
933 : }
934 :
935 19 : errc = error_t::noerror;
936 :
937 19 : if( m_naxis == 1 && m_x0 > -1 && m_xpix > -1 )
938 : {
939 1 : return m_xpix;
940 : }
941 18 : else if( m_x0 > -1 && m_y0 > -1 && m_xpix > -1 && m_ypix > -1 && m_naxis == 2 )
942 : {
943 2 : if( axis == 0 )
944 : {
945 1 : return m_xpix;
946 : }
947 1 : return m_ypix;
948 : }
949 : else
950 : {
951 16 : return m_naxes[axis];
952 : }
953 : }
954 :
955 : template <typename dataT, class verboseT>
956 44 : error_t fitsFile<dataT, verboseT>::calcPixarrs( pixarrT &pixarr )
957 : {
958 44 : error_t errc = pixarr.allocate( m_naxis );
959 44 : if( errc != error_t::noerror )
960 : {
961 2 : return internal::mxlib_error_report<verboseT>( errc );
962 : }
963 :
964 42 : if( m_naxis == 1 )
965 : {
966 3 : if( m_x0 > -1 && m_xpix > -1 )
967 : {
968 1 : pixarr.fpix[0] = m_x0 + 1;
969 1 : pixarr.lpix[0] = pixarr.fpix[0] + m_xpix - 1;
970 : }
971 : else
972 : {
973 2 : pixarr.fpix[0] = 1;
974 2 : pixarr.lpix[0] = m_naxes[0];
975 : }
976 3 : pixarr.inc[0] = 1;
977 : }
978 39 : else if( m_x0 > -1 && m_y0 > -1 && m_xpix > -1 && m_ypix > -1 && m_naxis == 2 )
979 : {
980 3 : pixarr.fpix[0] = m_x0 + 1;
981 3 : pixarr.lpix[0] = pixarr.fpix[0] + m_xpix - 1;
982 3 : pixarr.fpix[1] = m_y0 + 1;
983 3 : pixarr.lpix[1] = pixarr.fpix[1] + m_ypix - 1;
984 :
985 3 : pixarr.inc[0] = 1;
986 3 : pixarr.inc[1] = 1;
987 : }
988 : else
989 : {
990 36 : if( m_x0 < 0 && m_y0 < 0 && m_xpix < 0 && m_ypix < 0 && m_z0 < 0 && m_zframes < 0 )
991 : {
992 102 : for( int i = 0; i < m_naxis && i < 3; i++ )
993 : {
994 69 : pixarr.fpix[i] = 1;
995 69 : pixarr.lpix[i] = m_naxes[i];
996 69 : pixarr.inc[i] = 1;
997 : }
998 33 : }
999 : else
1000 : {
1001 3 : if( m_x0 > -1 && m_y0 > -1 && m_xpix > -1 && m_ypix > -1 )
1002 : {
1003 2 : pixarr.fpix[0] = m_x0 + 1;
1004 2 : pixarr.lpix[0] = pixarr.fpix[0] + m_xpix - 1;
1005 2 : pixarr.fpix[1] = m_y0 + 1;
1006 2 : pixarr.lpix[1] = pixarr.fpix[1] + m_ypix - 1;
1007 :
1008 2 : pixarr.inc[0] = 1;
1009 2 : pixarr.inc[1] = 1;
1010 : }
1011 : else
1012 : {
1013 1 : pixarr.fpix[0] = 1;
1014 1 : pixarr.lpix[0] = m_naxes[0];
1015 1 : pixarr.fpix[1] = 1;
1016 1 : pixarr.lpix[1] = m_naxes[1];
1017 :
1018 1 : pixarr.inc[0] = 1;
1019 1 : pixarr.inc[1] = 1;
1020 : }
1021 :
1022 3 : if( m_z0 > -1 && m_zframes > -1 )
1023 : {
1024 2 : pixarr.fpix[2] = m_z0 + 1;
1025 2 : pixarr.lpix[2] = pixarr.fpix[2] + m_zframes - 1;
1026 2 : pixarr.inc[2] = 1;
1027 : }
1028 : else
1029 : {
1030 1 : pixarr.fpix[2] = 1;
1031 1 : pixarr.lpix[2] = m_naxes[2];
1032 1 : pixarr.inc[2] = 1;
1033 : }
1034 : }
1035 : }
1036 :
1037 42 : return error_t::noerror;
1038 : }
1039 :
1040 : /************************************************************/
1041 : /*** Basic Arrays ***/
1042 : /************************************************************/
1043 :
1044 : template <typename dataT, class verboseT>
1045 13 : error_t fitsFile<dataT, verboseT>::read( dataT *data )
1046 : {
1047 13 : if( !m_isOpen )
1048 : {
1049 2 : mxlib_error_check( open() );
1050 : }
1051 :
1052 13 : pixarrT pixarrs;
1053 :
1054 13 : mxlib_error_check( calcPixarrs( pixarrs ) );
1055 :
1056 : ///\todo test if there is a speed difference for full reads for fits_read_pix/subset
1057 :
1058 13 : int fstatus = 0;
1059 :
1060 26 : fitsFileDetail::fitsFileCfitsioOpsInstance().readSubset( m_fptr,
1061 : fitsType<dataT>(),
1062 : pixarrs.fpix,
1063 : pixarrs.lpix,
1064 : pixarrs.inc,
1065 13 : (void *)&m_nulval,
1066 : (void *)data,
1067 : &m_anynul,
1068 : &fstatus );
1069 :
1070 13 : if( fstatus && fstatus != END_OF_FILE )
1071 : {
1072 1 : return internal::mxlib_error_report<verboseT>( fits_status2error_t( fstatus ),
1073 2 : "Reading data from " + m_fileName );
1074 : }
1075 :
1076 12 : return error_t::noerror;
1077 13 : }
1078 :
1079 : template <typename dataT, class verboseT>
1080 1 : error_t fitsFile<dataT, verboseT>::read( dataT *data, fitsHeader<verboseT> &head )
1081 : {
1082 1 : mxlib_error_check( read( data ) );
1083 :
1084 1 : mxlib_error_return( readHeader( head ) );
1085 : }
1086 :
1087 : template <typename dataT, class verboseT>
1088 2 : error_t fitsFile<dataT, verboseT>::read( dataT *data, const std::string &fname )
1089 : {
1090 2 : mxlib_error_check( fileName( fname ) );
1091 :
1092 2 : mxlib_error_return( read( data ) );
1093 : }
1094 :
1095 : template <typename dataT, class verboseT>
1096 1 : error_t fitsFile<dataT, verboseT>::read( dataT *data, fitsHeader<verboseT> &head, const std::string &fname )
1097 : {
1098 1 : mxlib_error_check( fileName( fname ) );
1099 :
1100 1 : mxlib_error_check( read( data ) );
1101 :
1102 1 : mxlib_error_return( readHeader( head ) );
1103 : }
1104 :
1105 : template <typename dataT, class verboseT>
1106 2 : error_t fitsFile<dataT, verboseT>::read( dataT *im, const std::vector<std::string> &flist )
1107 : {
1108 2 : if( flist.size() == 0 )
1109 : {
1110 3 : return internal::mxlib_error_report<verboseT>( error_t::invalidarg, "Empty file list" );
1111 : }
1112 :
1113 1 : long sz0 = 0, sz1 = 0;
1114 :
1115 3 : for( int i = 0; i < flist.size(); ++i )
1116 : {
1117 2 : mxlib_error_check( fileName( flist[i], 1 ) );
1118 :
1119 2 : mxlib_error_check( read( im + i * sz0 * sz1 ) );
1120 :
1121 : error_t errc;
1122 2 : sz0 = getSize( 0, errc );
1123 2 : mxlib_error_check( errc );
1124 :
1125 2 : sz1 = getSize( 1, errc );
1126 2 : mxlib_error_check( errc );
1127 : }
1128 :
1129 1 : return error_t::noerror;
1130 : }
1131 :
1132 : template <typename dataT, class verboseT>
1133 3 : error_t fitsFile<dataT, verboseT>::read( dataT *im,
1134 : std::vector<fitsHeader<verboseT>> &heads,
1135 : const std::vector<std::string> &flist )
1136 : {
1137 3 : if( flist.size() == 0 )
1138 : {
1139 3 : return internal::mxlib_error_report<verboseT>( error_t::invalidarg, "Empty file list" );
1140 : }
1141 :
1142 2 : if( heads.size() != flist.size() )
1143 : {
1144 2 : return internal::mxlib_error_report<verboseT>( error_t::invalidarg,
1145 3 : "Header vector is not the same size as the file list" );
1146 : }
1147 :
1148 1 : long sz0 = 0, sz1 = 0;
1149 :
1150 6 : for( size_t i = 0; i < flist.size(); ++i )
1151 : {
1152 :
1153 5 : mxlib_error_check( fileName( flist[i], 1 ) );
1154 :
1155 5 : mxlib_error_check( read( im + i * sz0 * sz1 ) );
1156 :
1157 5 : mxlib_error_check( readHeader( heads[i] ) );
1158 :
1159 : error_t errc;
1160 5 : sz0 = getSize( 0, errc );
1161 5 : mxlib_error_check( errc );
1162 :
1163 5 : sz1 = getSize( 1, errc );
1164 5 : mxlib_error_check( errc );
1165 : }
1166 :
1167 1 : return error_t::noerror;
1168 : }
1169 :
1170 : /************************************************************/
1171 : /*** Eigen Arrays ***/
1172 : /************************************************************/
1173 :
1174 : template <typename arrT, class verboseT, bool isCube = improc::is_eigenCube<arrT>::value>
1175 : struct eigenArrResize
1176 : {
1177 : // If it's a cube, always pass zsz
1178 5 : error_t resize( arrT &arr, int xsz, int ysz, int zsz )
1179 : {
1180 : try
1181 : {
1182 5 : arr.resize( xsz, ysz, zsz );
1183 : }
1184 4 : catch( const std::bad_alloc &e )
1185 : {
1186 1 : internal::mxlib_error_report<verboseT>( error_t::std_bad_alloc,
1187 3 : std::string( "resizing array: " ) + e.what() );
1188 : #ifdef MXLIB_TRAP_ALLOC_ERRORS
1189 : return error_t::std_bad_alloc;
1190 : #else
1191 1 : throw;
1192 : #endif
1193 : }
1194 2 : catch( const std::exception &e )
1195 : {
1196 1 : internal::mxlib_error_report<verboseT>( error_t::std_exception,
1197 3 : std::string( "resizing array: " ) + e.what() );
1198 : #ifdef MXLIB_TRAP_ALLOC_ERRORS
1199 : return error_t::std_exception;
1200 : #else
1201 1 : throw;
1202 : #endif
1203 : }
1204 2 : catch( ... )
1205 : {
1206 2 : internal::mxlib_error_report<verboseT>( error_t::exception, "resizing array" );
1207 : #ifdef MXLIB_TRAP_ALLOC_ERRORS
1208 : return error_t::exception;
1209 : #else
1210 1 : throw;
1211 : #endif
1212 : }
1213 :
1214 2 : return error_t::noerror;
1215 : }
1216 : };
1217 :
1218 : template <typename arrT, class verboseT>
1219 : struct eigenArrResize<arrT, verboseT, false>
1220 : {
1221 : // If it's not a cube, never pass zsz
1222 14 : error_t resize( arrT &arr, int xsz, int ysz, [[maybe_unused]] int zsz )
1223 : {
1224 : try
1225 : {
1226 14 : arr.resize( xsz, ysz );
1227 : }
1228 4 : catch( const std::bad_alloc &e )
1229 : {
1230 1 : internal::mxlib_error_report<verboseT>( error_t::std_bad_alloc,
1231 3 : std::string( "resizing array: " ) + e.what() );
1232 : #ifdef MXLIB_TRAP_ALLOC_ERRORS
1233 : return error_t::std_bad_alloc;
1234 : #else
1235 1 : throw;
1236 : #endif
1237 : }
1238 2 : catch( const std::exception &e )
1239 : {
1240 1 : internal::mxlib_error_report<verboseT>( error_t::std_exception,
1241 3 : std::string( "resizing array: " ) + e.what() );
1242 : #ifdef MXLIB_TRAP_ALLOC_ERRORS
1243 : return error_t::std_exception;
1244 : #else
1245 1 : throw;
1246 : #endif
1247 : }
1248 2 : catch( ... )
1249 : {
1250 2 : internal::mxlib_error_report<verboseT>( error_t::exception, "resizing array" );
1251 : #ifdef MXLIB_TRAP_ALLOC_ERRORS
1252 : return error_t::exception;
1253 : #else
1254 1 : throw;
1255 : #endif
1256 : }
1257 :
1258 11 : return error_t::noerror;
1259 : }
1260 : };
1261 :
1262 : template <typename dataT, class verboseT>
1263 : template <typename arrT>
1264 13 : error_t fitsFile<dataT, verboseT>::read( arrT &im )
1265 : {
1266 : ///\todo this can probably be made part of one read function (or call read(data *)) with a call to resize with
1267 : /// SFINAE
1268 13 : int fstatus = 0;
1269 :
1270 13 : if( !m_isOpen )
1271 : {
1272 1 : mxlib_error_check( open() );
1273 : }
1274 :
1275 13 : pixarrT pixarrs;
1276 13 : mxlib_error_check( calcPixarrs( pixarrs ) );
1277 :
1278 : eigenArrResize<arrT, verboseT> arrresz;
1279 13 : if( m_naxis > 2 )
1280 : {
1281 2 : mxlib_error_check( arrresz.resize( im,
1282 : pixarrs.lpix[0] - pixarrs.fpix[0] + 1,
1283 : pixarrs.lpix[1] - pixarrs.fpix[1] + 1,
1284 : pixarrs.lpix[2] - pixarrs.fpix[2] + 1 ) );
1285 : }
1286 11 : else if( m_naxis > 1 )
1287 : {
1288 10 : mxlib_error_check(
1289 : arrresz.resize( im, pixarrs.lpix[0] - pixarrs.fpix[0] + 1, pixarrs.lpix[1] - pixarrs.fpix[1] + 1, 1 ) );
1290 : }
1291 : else
1292 : {
1293 1 : mxlib_error_check( arrresz.resize( im, pixarrs.lpix[0] - pixarrs.fpix[0] + 1, 1, 1 ) );
1294 : }
1295 :
1296 26 : fitsFileDetail::fitsFileCfitsioOpsInstance().readSubset( m_fptr,
1297 : fitsType<typename arrT::Scalar>(),
1298 : pixarrs.fpix,
1299 : pixarrs.lpix,
1300 : pixarrs.inc,
1301 13 : (void *)&m_nulval,
1302 13 : (void *)im.data(),
1303 : &m_anynul,
1304 : &fstatus );
1305 :
1306 13 : if( fstatus && fstatus != END_OF_FILE )
1307 : {
1308 2 : return internal::mxlib_error_report<verboseT>( fits_status2error_t( fstatus ),
1309 4 : "Reading data from " + m_fileName );
1310 : }
1311 :
1312 11 : return error_t::noerror;
1313 13 : }
1314 :
1315 : template <typename dataT, class verboseT>
1316 : template <typename arrT>
1317 : error_t fitsFile<dataT, verboseT>::read( arrT &data, fitsHeader<verboseT> &head )
1318 : {
1319 : error_t errc;
1320 : errc = read( data );
1321 : if( errc != error_t::noerror )
1322 : {
1323 : return internal::mxlib_error_report<verboseT>( errc );
1324 : }
1325 :
1326 : errc = readHeader( head );
1327 : if( errc != error_t::noerror )
1328 : {
1329 : return internal::mxlib_error_report<verboseT>( errc );
1330 : }
1331 :
1332 : return error_t::noerror;
1333 : }
1334 :
1335 : template <typename dataT, class verboseT>
1336 : template <typename arrT>
1337 7 : error_t fitsFile<dataT, verboseT>::read( arrT &data, const std::string &fname )
1338 : {
1339 : error_t errc;
1340 7 : errc = fileName( fname );
1341 7 : if( errc != error_t::noerror )
1342 : {
1343 1 : return internal::mxlib_error_report<verboseT>( errc );
1344 : }
1345 :
1346 6 : errc = read( data );
1347 6 : if( errc != error_t::noerror )
1348 : {
1349 1 : return internal::mxlib_error_report<verboseT>( errc );
1350 : }
1351 5 : return error_t::noerror;
1352 : }
1353 :
1354 : template <typename dataT, class verboseT>
1355 : template <typename arrT>
1356 7 : error_t fitsFile<dataT, verboseT>::read( arrT &data, fitsHeader<verboseT> &head, const std::string &fname )
1357 : {
1358 : error_t errc;
1359 7 : errc = fileName( fname );
1360 7 : if( errc != error_t::noerror )
1361 : {
1362 1 : return internal::mxlib_error_report<verboseT>( errc );
1363 : }
1364 :
1365 6 : errc = read( data );
1366 6 : if( errc != error_t::noerror )
1367 : {
1368 1 : return internal::mxlib_error_report<verboseT>( errc );
1369 : }
1370 :
1371 5 : errc = readHeader( head );
1372 5 : if( errc != error_t::noerror )
1373 : {
1374 1 : return internal::mxlib_error_report<verboseT>( errc );
1375 : }
1376 :
1377 4 : return error_t::noerror;
1378 : }
1379 :
1380 : template <typename dataT, class verboseT>
1381 : template <typename cubeT>
1382 12 : error_t fitsFile<dataT, verboseT>::read( cubeT &cube,
1383 : const std::vector<std::string> &flist,
1384 : std::vector<fitsHeader<verboseT>> *heads )
1385 : {
1386 : error_t errc;
1387 12 : int fstatus = 0;
1388 :
1389 12 : if( flist.size() == 0 )
1390 : {
1391 3 : return internal::mxlib_error_report<verboseT>( error_t::invalidarg, "Empty file list" );
1392 : }
1393 :
1394 11 : if( heads && heads->size() != flist.size() )
1395 : {
1396 2 : return internal::mxlib_error_report<verboseT>( error_t::invalidarg,
1397 3 : "Header vector is not the same size as the file list" );
1398 : }
1399 :
1400 : // Open the first file to get the dimensions.
1401 10 : errc = fileName( flist[0], 1 );
1402 10 : if( !!errc )
1403 : {
1404 1 : return internal::mxlib_error_report<verboseT>( errc );
1405 : }
1406 :
1407 9 : pixarrT pixarrs;
1408 9 : errc = calcPixarrs( pixarrs );
1409 :
1410 9 : if( !!errc )
1411 : {
1412 1 : return internal::mxlib_error_report<verboseT>( errc );
1413 : }
1414 :
1415 8 : cube.resize( pixarrs.lpix[0] - pixarrs.fpix[0] + 1, pixarrs.lpix[1] - pixarrs.fpix[1] + 1, flist.size() );
1416 :
1417 : // Now read first image.
1418 16 : fitsFileDetail::fitsFileCfitsioOpsInstance().readSubset( m_fptr,
1419 : fitsType<typename cubeT::Scalar>(),
1420 : pixarrs.fpix,
1421 : pixarrs.lpix,
1422 : pixarrs.inc,
1423 8 : (void *)&m_nulval,
1424 8 : (void *)cube.image( 0 ).data(),
1425 : &m_anynul,
1426 : &fstatus );
1427 :
1428 8 : if( fstatus && fstatus != END_OF_FILE )
1429 : {
1430 2 : return internal::mxlib_error_report<verboseT>( fits_status2error_t( fstatus ),
1431 4 : "Reading data from " + m_fileName );
1432 : }
1433 :
1434 6 : if( heads )
1435 : {
1436 4 : errc = readHeader( ( *heads )[0] );
1437 4 : if( errc != error_t::noerror )
1438 : {
1439 1 : return internal::mxlib_error_report<verboseT>( errc );
1440 : }
1441 : }
1442 :
1443 : // Now read in the rest.
1444 9 : for( int i = 1; i < flist.size(); ++i )
1445 : {
1446 8 : errc = fileName( flist[i], 1 );
1447 8 : if( errc != error_t::noerror )
1448 : {
1449 1 : return internal::mxlib_error_report<verboseT>( errc );
1450 : }
1451 :
1452 : // Now read image.
1453 14 : fitsFileDetail::fitsFileCfitsioOpsInstance().readSubset( m_fptr,
1454 : fitsType<typename cubeT::Scalar>(),
1455 : pixarrs.fpix,
1456 : pixarrs.lpix,
1457 : pixarrs.inc,
1458 7 : (void *)&m_nulval,
1459 7 : (void *)cube.image( i ).data(),
1460 : &m_anynul,
1461 : &fstatus );
1462 :
1463 7 : if( fstatus && fstatus != END_OF_FILE )
1464 : {
1465 2 : return internal::mxlib_error_report<verboseT>( fits_status2error_t( fstatus ),
1466 4 : "Reading data from " + m_fileName );
1467 : }
1468 :
1469 5 : if( heads )
1470 : {
1471 5 : errc = readHeader( ( *heads )[i] );
1472 5 : if( errc != error_t::noerror )
1473 : {
1474 1 : return internal::mxlib_error_report<verboseT>( errc );
1475 : }
1476 : }
1477 : }
1478 :
1479 1 : return error_t::noerror;
1480 9 : }
1481 :
1482 : template <typename dataT, class verboseT>
1483 : template <typename cubeT>
1484 6 : error_t fitsFile<dataT, verboseT>::read( cubeT &cube,
1485 : std::vector<fitsHeader<verboseT>> &heads,
1486 : const std::vector<std::string> &flist )
1487 : {
1488 6 : mxlib_error_return( read( cube, flist, &heads ) );
1489 : }
1490 :
1491 : template <typename dataT, class verboseT>
1492 34 : error_t fitsFile<dataT, verboseT>::readHeader( fitsHeader<verboseT> &head )
1493 : {
1494 34 : int fstatus = 0;
1495 :
1496 : char keyword[FLEN_KEYWORD];
1497 : char value[FLEN_VALUE];
1498 34 : std::unique_ptr<char[]> commentStorage;
1499 34 : char *comment = nullptr;
1500 :
1501 : // The keys to look for if head is already populated
1502 34 : typename std::list<headerIteratorT> head_keys;
1503 34 : typename std::list<headerIteratorT>::iterator head_keys_it;
1504 : // int num_head_keys;
1505 :
1506 34 : bool head_keys_only = false;
1507 34 : if( head.size() > 0 )
1508 : {
1509 18 : head_keys_only = true;
1510 18 : headerIteratorT headIt = head.begin();
1511 60 : while( headIt != head.end() )
1512 : {
1513 42 : head_keys.push_back( headIt );
1514 42 : ++headIt;
1515 : }
1516 : // num_head_keys = head.size();
1517 : }
1518 :
1519 : // If m_noComment is set, then we don't read in the comment
1520 34 : if( m_noComment )
1521 : {
1522 1 : comment = nullptr;
1523 : }
1524 : else
1525 : {
1526 33 : commentStorage = std::make_unique<char[]>( FLEN_COMMENT );
1527 33 : comment = commentStorage.get();
1528 : }
1529 :
1530 : int keysexist;
1531 : int morekeys;
1532 :
1533 34 : if( !m_isOpen )
1534 : {
1535 1 : mxlib_error_check( open() );
1536 : }
1537 :
1538 : // This gets the number of header keys to read
1539 33 : fitsFileDetail::fitsFileCfitsioOpsInstance().getHeaderSpace( m_fptr, &keysexist, &morekeys, &fstatus );
1540 :
1541 33 : if( fstatus )
1542 : {
1543 4 : return internal::mxlib_error_report<verboseT>( fits_status2error_t( fstatus ),
1544 8 : "Reading header from " + m_fileName );
1545 : }
1546 :
1547 289 : for( int i = 0; i < keysexist; i++ )
1548 : {
1549 279 : fitsFileDetail::fitsFileCfitsioOpsInstance().readKey( m_fptr, i + 1, keyword, value, comment, &fstatus );
1550 :
1551 279 : if( fstatus )
1552 : {
1553 1 : return internal::mxlib_error_report<verboseT>( fits_status2error_t( fstatus ),
1554 2 : "Reading header from " + m_fileName );
1555 : }
1556 :
1557 278 : if( !head_keys_only )
1558 : {
1559 102 : if( strcmp( keyword, "COMMENT" ) == 0 )
1560 : {
1561 84 : head.template append<fitsCommentType>( keyword, fitsCommentType( value ), comment );
1562 : }
1563 81 : else if( strcmp( keyword, "HISTORY" ) == 0 )
1564 : {
1565 4 : head.template append<fitsHistoryType>( keyword, fitsHistoryType( value ), comment );
1566 : }
1567 : else
1568 : {
1569 : // Otherwise we append it as an unknown type
1570 320 : head.append( keyword, value, comment );
1571 : }
1572 : }
1573 : else
1574 : {
1575 176 : head_keys_it = head_keys.begin();
1576 317 : while( head_keys_it != head_keys.end() )
1577 : {
1578 183 : if( ( *( *head_keys_it ) ).keyword() == keyword )
1579 : {
1580 42 : head[keyword].value( (const char *)value );
1581 42 : if( comment )
1582 : {
1583 164 : head[keyword].comment( comment );
1584 : }
1585 :
1586 42 : head_keys.erase( head_keys_it );
1587 :
1588 42 : break;
1589 : }
1590 141 : ++head_keys_it;
1591 : }
1592 :
1593 : // Quit if we're done.
1594 176 : if( head_keys.empty() )
1595 : {
1596 18 : break;
1597 : }
1598 : }
1599 : }
1600 :
1601 28 : return error_t::noerror;
1602 34 : }
1603 :
1604 : template <typename dataT, class verboseT>
1605 3 : error_t fitsFile<dataT, verboseT>::readHeader( fitsHeader<verboseT> &head, const std::string &fname )
1606 : {
1607 : error_t errc;
1608 3 : mxlib_error_check( fileName( fname ) );
1609 3 : mxlib_error_check( readHeader( head ) );
1610 1 : return error_t::noerror;
1611 : }
1612 :
1613 : template <typename dataT, class verboseT>
1614 4 : error_t fitsFile<dataT, verboseT>::readHeader( std::vector<fitsHeader<verboseT>> &heads,
1615 : const std::vector<std::string> &flist )
1616 : {
1617 4 : if( heads.size() != 0 && heads.size() != flist.size() )
1618 : {
1619 2 : return internal::mxlib_error_report<verboseT>( error_t::invalidarg,
1620 3 : "head vector is not empty and not same size as file list" );
1621 : }
1622 :
1623 3 : if( heads.size() == 0 )
1624 : {
1625 1 : heads.resize( flist.size() );
1626 : }
1627 :
1628 : error_t errc;
1629 12 : for( int i = 0; i < flist.size(); ++i )
1630 : {
1631 9 : mxlib_error_check( fileName( flist[i], 1 ) );
1632 :
1633 9 : mxlib_error_check( readHeader( heads[i] ) );
1634 : }
1635 :
1636 3 : return error_t::noerror;
1637 : }
1638 :
1639 : template <typename dataT, class verboseT>
1640 56 : error_t fitsFile<dataT, verboseT>::write( const dataT *im, int d1, int d2, int d3, fitsHeader<verboseT> *head )
1641 : {
1642 56 : int fstatus = 0;
1643 :
1644 56 : if( m_isOpen )
1645 : {
1646 1 : mxlib_error_check( close() );
1647 : }
1648 :
1649 56 : m_naxis = 1;
1650 56 : if( d2 > 0 )
1651 : {
1652 55 : if( d3 > 1 )
1653 : {
1654 3 : m_naxis = 3;
1655 : }
1656 : else
1657 : {
1658 52 : m_naxis = 2;
1659 : }
1660 : }
1661 :
1662 56 : m_naxes[0] = d1;
1663 56 : if( m_naxis > 1 )
1664 : {
1665 55 : m_naxes[1] = d2;
1666 : }
1667 56 : if( m_naxis > 2 )
1668 : {
1669 3 : m_naxes[2] = d3;
1670 : }
1671 :
1672 56 : std::string forceFileName = "!" + m_fileName;
1673 :
1674 56 : fitsFileDetail::fitsFileCfitsioOpsInstance().createFile( &m_fptr, forceFileName.c_str(), &fstatus );
1675 56 : if( fstatus )
1676 : {
1677 4 : return internal::mxlib_error_report<verboseT>( fits_status2error_t( fstatus ), "Creating " + m_fileName );
1678 : }
1679 52 : m_isOpen = true;
1680 :
1681 52 : fstatus = 0;
1682 52 : fitsFileDetail::fitsFileCfitsioOpsInstance().createImage( m_fptr, fitsBITPIX<dataT>(), m_naxis, m_naxes, &fstatus );
1683 52 : if( fstatus )
1684 : {
1685 1 : return internal::mxlib_error_report<verboseT>( fits_status2error_t( fstatus ),
1686 2 : "Creating image in" + m_fileName );
1687 : }
1688 :
1689 : long fpixel[3];
1690 51 : fpixel[0] = 1;
1691 51 : fpixel[1] = 1;
1692 51 : fpixel[2] = 1;
1693 :
1694 51 : LONGLONG nelements = 1;
1695 :
1696 154 : for( int i = 0; i < m_naxis && i < 3; ++i )
1697 : {
1698 103 : nelements *= m_naxes[i];
1699 : }
1700 :
1701 51 : fstatus = 0;
1702 51 : fitsFileDetail::fitsFileCfitsioOpsInstance()
1703 51 : .writePixels( m_fptr, fitsType<dataT>(), fpixel, nelements, (void *)im, &fstatus );
1704 51 : if( fstatus )
1705 : {
1706 1 : return internal::mxlib_error_report<verboseT>( fits_status2error_t( fstatus ), "Writing data " + m_fileName );
1707 : }
1708 :
1709 50 : if( head != 0 )
1710 : {
1711 35 : headerIteratorT it;
1712 :
1713 102 : for( it = head->begin(); it != head->end(); ++it )
1714 : {
1715 68 : error_t errc = it->write( m_fptr );
1716 68 : if( errc != error_t::noerror )
1717 : {
1718 1 : return internal::mxlib_error_report<verboseT>( errc, "Writing keyword " + m_fileName );
1719 : }
1720 : }
1721 : }
1722 :
1723 49 : mxlib_error_return( close() )
1724 56 : }
1725 :
1726 : template <typename dataT, class verboseT>
1727 1 : error_t fitsFile<dataT, verboseT>::write( const dataT *im, int d1, int d2, int d3 )
1728 : {
1729 1 : mxlib_error_return( write( im, d1, d2, d3, (fitsHeader<verboseT> *)0 ) );
1730 : }
1731 :
1732 : template <typename dataT, class verboseT>
1733 1 : error_t fitsFile<dataT, verboseT>::write( const dataT *im, int d1, int d2, int d3, fitsHeader<verboseT> &head )
1734 : {
1735 1 : mxlib_error_return( write( im, d1, d2, d3, &head ) );
1736 : }
1737 :
1738 : template <typename dataT, class verboseT>
1739 18 : error_t fitsFile<dataT, verboseT>::write( const std::string &fname, const dataT *im, int d1, int d2, int d3 )
1740 : {
1741 18 : mxlib_error_check( fileName( fname, false ) );
1742 :
1743 18 : mxlib_error_return( write( im, d1, d2, d3, (fitsHeader<verboseT> *)0 ) );
1744 : }
1745 :
1746 : template <typename dataT, class verboseT>
1747 36 : error_t fitsFile<dataT, verboseT>::write(
1748 : const std::string &fname, const dataT *im, int d1, int d2, int d3, fitsHeader<verboseT> &head )
1749 : {
1750 36 : mxlib_error_check( fileName( fname, false ) );
1751 :
1752 36 : mxlib_error_return( write( im, d1, d2, d3, &head ) );
1753 : }
1754 :
1755 : template <typename dataT, class verboseT>
1756 : template <typename arrT>
1757 16 : error_t fitsFile<dataT, verboseT>::write( const std::string &fname, const arrT &im )
1758 : {
1759 : improc::eigenArrPlanes<arrT> planes;
1760 :
1761 16 : mxlib_error_return( write( fname, im.data(), im.rows(), im.cols(), planes( im ) ) );
1762 : }
1763 :
1764 : template <typename dataT, class verboseT>
1765 : template <typename arrT>
1766 34 : error_t fitsFile<dataT, verboseT>::write( const std::string &fname, const arrT &im, fitsHeader<verboseT> &head )
1767 : {
1768 : improc::eigenArrPlanes<arrT> planes;
1769 :
1770 68 : return write( fname, im.data(), im.rows(), im.cols(), planes( im ), head );
1771 : }
1772 :
1773 : template <typename dataT, class verboseT>
1774 4 : void fitsFile<dataT, verboseT>::setReadSize()
1775 : {
1776 4 : m_x0 = -1;
1777 4 : m_y0 = -1;
1778 4 : m_xpix = -1;
1779 4 : m_ypix = -1;
1780 4 : }
1781 :
1782 : template <typename dataT, class verboseT>
1783 8 : void fitsFile<dataT, verboseT>::setReadSize( long x0, long y0, long xpix, long ypix )
1784 : {
1785 8 : m_x0 = x0;
1786 8 : m_y0 = y0;
1787 8 : m_xpix = xpix;
1788 8 : m_ypix = ypix;
1789 8 : }
1790 :
1791 : template <typename dataT, class verboseT>
1792 1 : void fitsFile<dataT, verboseT>::setCubeReadSize()
1793 : {
1794 1 : m_z0 = -1;
1795 1 : m_zframes = -1;
1796 1 : }
1797 :
1798 : template <typename dataT, class verboseT>
1799 3 : void fitsFile<dataT, verboseT>::setCubeReadSize( long z0, long zframes )
1800 : {
1801 3 : m_z0 = z0;
1802 3 : m_zframes = zframes;
1803 3 : }
1804 :
1805 : /** \ingroup fits_processing_typedefs
1806 : * @{
1807 : */
1808 :
1809 : /// A \ref fitsFile to work in signed characters
1810 : typedef fitsFile<char> fitsFilec;
1811 :
1812 : /// A \ref fitsFile to work in unsigned characters
1813 : typedef fitsFile<unsigned char> fitsFileuc;
1814 :
1815 : /// A \ref fitsFile to work in signed short integers
1816 : typedef fitsFile<short> fitsFiles;
1817 :
1818 : /// A \ref fitsFile to work in unsigned short integers
1819 : typedef fitsFile<unsigned short> fitsFileus;
1820 :
1821 : /// A \ref fitsFile to work in signed integers
1822 : typedef fitsFile<int> fitsFilei;
1823 :
1824 : /// A \ref fitsFile to work in unsigned integers
1825 : typedef fitsFile<unsigned int> fitsFileui;
1826 :
1827 : /// A \ref fitsFile to work in signed long integers
1828 : typedef fitsFile<long> fitsFilel;
1829 :
1830 : /// A \ref fitsFile to work in single precision floats
1831 : typedef fitsFile<float> fitsFilef;
1832 :
1833 : /// A \ref fitsFile to work in double precision
1834 : typedef fitsFile<double> fitsFiled;
1835 :
1836 : ///@}
1837 :
1838 : extern template class fitsFile<short, verbose::d>;
1839 : extern template class fitsFile<unsigned short, verbose::d>;
1840 : extern template class fitsFile<int, verbose::d>;
1841 : extern template class fitsFile<unsigned int, verbose::d>;
1842 : extern template class fitsFile<float, verbose::d>;
1843 : extern template class fitsFile<double, verbose::d>;
1844 :
1845 : } // namespace fits
1846 : } // namespace mx
1847 :
1848 : #endif // ioutils_fits_fitsFile_hpp
|
