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