27#ifndef ioutils_fits_fitsFile_hpp
28#define ioutils_fits_fitsFile_hpp
30#include "../../mxlib.hpp"
32#include "../../improc/eigenImage.hpp"
51template <
typename dataT,
class verboseT = verbose::d>
55 friend class fitsFile_test;
162 long naxes(
int dim );
215 long *fpix{
nullptr };
216 long *lpix{
nullptr };
217 long *inc{
nullptr };
246 fpix =
new long[
naxis];
247 lpix =
new long[
naxis];
248 inc =
new long[
naxis];
250 catch(
const std::bad_alloc &e )
253 std::string(
"allocating pixel read arrays: " ) + e.what() );
254#ifdef MXLIB_TRAP_ALLOC_ERRORS
260 catch(
const std::exception &e )
263 std::string(
"allocating pixel read arrays: " ) + e.what() );
264#ifdef MXLIB_TRAP_ALLOC_ERRORS
272 internal::mxlib_error_report<verboseT>(
error_t::exception,
"allocating pixel read arrays" );
273#ifdef MXLIB_TRAP_ALLOC_ERRORS
281 if( fpix ==
nullptr )
283 return internal::mxlib_error_report<verboseT>(
error_t::allocerr,
"fpix is nullptr" );
286 if( lpix ==
nullptr )
288 return internal::mxlib_error_report<verboseT>(
error_t::allocerr,
"lpix is nullptr" );
293 return internal::mxlib_error_report<verboseT>(
error_t::allocerr,
"inc is nullptr" );
339 fitsHeader<verboseT> &head
349 const std::string &fname
359 fitsHeader<verboseT> &head,
360 const std::string &fname
365 const std::vector<std::string> &flist
370 std::vector<fitsHeader<verboseT>> &heads,
372 const std::vector<std::string> &flist
393 template <
typename arrT>
408 template <
typename arrT>
411 fitsHeader<verboseT> &head
425 template <
typename arrT>
428 const std::string &fname
442 template <
typename arrT>
445 fitsHeader<verboseT> &head,
446 const std::string &fname
462 template <
typename cubeT>
464 const std::vector<std::string> &flist,
465 std::vector<fitsHeader<verboseT>> *heads =
480 template <
typename cubeT>
483 std::vector<fitsHeader<verboseT>> &heads,
485 const std::vector<std::string> &flist
511 const std::string &fname
522 std::vector<fitsHeader<verboseT>> &heads,
523 const std::vector<std::string> &flist
542 fitsHeader<verboseT> *head
567 fitsHeader<verboseT> &head
592 fitsHeader<verboseT> &head
616 template <
typename arrT>
635 template <
typename arrT>
638 fitsHeader<verboseT> &head
680template <
typename dataT,
class verboseT>
685template <
typename dataT,
class verboseT>
691template <
typename dataT,
class verboseT>
695 errc = internal::mxlib_error_report<verboseT>( fileName( fname,
false ) );
698template <
typename dataT,
class verboseT>
703 internal::mxlib_error_report<verboseT>( fileName( fname, doopen ) );
706template <
typename dataT,
class verboseT>
711 errc = internal::mxlib_error_report<verboseT>( fileName( fname, doopen ) );
714template <
typename dataT,
class verboseT>
719 internal::mxlib_error_report<verboseT>( close() );
724template <
typename dataT,
class verboseT>
730template <
typename dataT,
class verboseT>
748template <
typename dataT,
class verboseT>
754template <
typename dataT,
class verboseT>
757 if( dim >= m_naxis || dim > 2)
765template <
typename dataT,
class verboseT>
773 if( m_fileName ==
"" )
780 fits_open_file( &m_fptr, m_fileName.c_str(), READONLY, &fstatus );
784 return internal::mxlib_error_report<verboseT>(
fits_status2error_t( fstatus ),
"Opening file " + m_fileName );
788 fits_get_img_dim( m_fptr, &m_naxis, &fstatus );
792 "Getting number of axes in file " + m_fileName );
796 if( m_naxes ==
nullptr )
798 return internal::mxlib_error_report<verboseT>(
error_t::allocerr,
"m_naxes is nullptr" );
802 fits_get_img_size( m_fptr, m_naxis, m_naxes, &fstatus );
806 "Getting dimensions in file " + m_fileName );
814template <
typename dataT,
class verboseT>
820template <
typename dataT,
class verboseT>
829 fits_close_file( m_fptr, &fstatus );
833 return internal::mxlib_error_report<verboseT>(
fits_status2error_t( fstatus ),
"Closing file " + m_fileName );
843template <
typename dataT,
class verboseT>
855template <
typename dataT,
class verboseT>
868 if( m_x0 > -1 && m_y0 > -1 && m_xpix > -1 && m_ypix > -1 && m_naxis == 2 )
870 return m_xpix * m_ypix;
874 for(
int i = 0; i < m_naxis && i < 3; ++i )
883template <
typename dataT,
class verboseT>
892 if( axis >= m_naxis || axis > 2 )
900 if( m_x0 > -1 && m_y0 > -1 && m_xpix > -1 && m_ypix > -1 && m_naxis == 2 )
910 return m_naxes[axis];
914template <
typename dataT,
class verboseT>
917 error_t errc = pixarr.allocate( m_naxis );
920 if( m_naxes ==
nullptr )
925 if( m_x0 > -1 && m_y0 > -1 && m_xpix > -1 && m_ypix > -1 && m_naxis == 2 )
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;
937 if( m_x0 < 0 && m_y0 < 0 && m_xpix < 0 && m_ypix < 0 && m_z0 < 0 && m_zframes < 0 )
939 for(
int i = 0; i < m_naxis && i < 3; i++ )
942 pixarr.lpix[i] = m_naxes[i];
948 if( m_x0 > -1 && m_y0 > -1 && m_xpix > -1 && m_ypix > -1 )
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;
961 pixarr.lpix[0] = m_naxes[0];
963 pixarr.lpix[1] = m_naxes[1];
969 if( m_z0 > -1 && m_zframes > -1 )
971 pixarr.fpix[2] = m_z0 + 1;
972 pixarr.lpix[2] = pixarr.fpix[2] + m_zframes - 1;
978 pixarr.lpix[2] = m_naxes[2];
991template <
typename dataT,
class verboseT>
1007 fits_read_subset( m_fptr,
1017 if( fstatus && fstatus != END_OF_FILE )
1020 "Reading data from " + m_fileName );
1026template <
typename dataT,
class verboseT>
1034template <
typename dataT,
class verboseT>
1042template <
typename dataT,
class verboseT>
1052template <
typename dataT,
class verboseT>
1055 if( flist.size() == 0 )
1060 long sz0 = 0, sz1 = 0;
1062 for(
int i = 0; i < flist.size(); ++i )
1069 sz0 = getSize( 0, errc );
1072 sz1 = getSize( 1, errc );
1079template <
typename dataT,
class verboseT>
1082 const std::vector<std::string> &flist )
1084 if( flist.size() == 0 )
1089 long sz0 = 0, sz1 = 0;
1091 for(
size_t i = 0; i < flist.size(); ++i )
1101 sz0 = getSize( 0, errc );
1104 sz1 = getSize( 1, errc );
1115template <typename arrT, class verboseT, bool isCube = improc::is_eigenCube<arrT>::value>
1116struct eigenArrResize
1119 error_t resize( arrT &arr,
int xsz,
int ysz,
int zsz )
1123 arr.resize( xsz, ysz, zsz );
1125 catch(
const std::bad_alloc &e )
1128 std::string(
"resizing array: " ) + e.what() );
1129#ifdef MXLIB_TRAP_ALLOC_ERRORS
1135 catch(
const std::exception &e )
1138 std::string(
"resizing array: " ) + e.what() );
1139#ifdef MXLIB_TRAP_ALLOC_ERRORS
1148#ifdef MXLIB_TRAP_ALLOC_ERRORS
1159template <
typename arrT,
class verboseT>
1160struct eigenArrResize<arrT, verboseT, false>
1163 error_t resize( arrT &arr,
int xsz,
int ysz, [[maybe_unused]]
int zsz )
1167 arr.resize( xsz, ysz );
1169 catch(
const std::bad_alloc &e )
1172 std::string(
"resizing array: " ) + e.what() );
1173#ifdef MXLIB_TRAP_ALLOC_ERRORS
1179 catch(
const std::exception &e )
1182 std::string(
"resizing array: " ) + e.what() );
1183#ifdef MXLIB_TRAP_ALLOC_ERRORS
1192#ifdef MXLIB_TRAP_ALLOC_ERRORS
1203template <
typename dataT,
class verboseT>
1204template <
typename arrT>
1219 eigenArrResize<arrT, verboseT> arrresz;
1223 pixarrs.lpix[0] - pixarrs.fpix[0] + 1,
1224 pixarrs.lpix[1] - pixarrs.fpix[1] + 1,
1225 pixarrs.lpix[2] - pixarrs.fpix[2] + 1 ) );
1227 else if( m_naxis > 1 )
1230 arrresz.resize( im, pixarrs.lpix[0] - pixarrs.fpix[0] + 1, pixarrs.lpix[1] - pixarrs.fpix[1] + 1, 1 ) );
1234 mxlib_error_check( arrresz.resize( im, pixarrs.lpix[0] - pixarrs.fpix[0] + 1, 1, 1 ) );
1237 fits_read_subset( m_fptr,
1238 fitsType<typename arrT::Scalar>(),
1247 if( fstatus && fstatus != END_OF_FILE )
1250 "Reading data from " + m_fileName );
1256template <
typename dataT,
class verboseT>
1257template <
typename arrT>
1261 errc = read( data );
1264 return internal::mxlib_error_report<verboseT>( errc );
1267 errc = readHeader( head );
1270 return internal::mxlib_error_report<verboseT>( errc );
1276template <
typename dataT,
class verboseT>
1277template <
typename arrT>
1281 errc = fileName( fname );
1284 return internal::mxlib_error_report<verboseT>( errc );
1287 errc = read( data );
1290 return internal::mxlib_error_report<verboseT>( errc );
1295template <
typename dataT,
class verboseT>
1296template <
typename arrT>
1300 errc = fileName( fname );
1303 return internal::mxlib_error_report<verboseT>( errc );
1306 errc = read( data );
1309 return internal::mxlib_error_report<verboseT>( errc );
1312 errc = readHeader( head );
1315 return internal::mxlib_error_report<verboseT>( errc );
1321template <
typename dataT,
class verboseT>
1322template <
typename cubeT>
1324 const std::vector<std::string> &flist,
1330 if( flist.size() == 0 )
1336 errc = fileName( flist[0], 1 );
1339 return internal::mxlib_error_report<verboseT>( errc );
1343 errc = calcPixarrs( pixarrs );
1347 return internal::mxlib_error_report<verboseT>( errc );
1350 cube.resize( pixarrs.lpix[0] - pixarrs.fpix[0] + 1, pixarrs.lpix[1] - pixarrs.fpix[1] + 1, flist.size() );
1353 fits_read_subset( m_fptr,
1354 fitsType<typename cubeT::Scalar>(),
1359 (
void *)cube.image( 0 ).data(),
1363 if( fstatus && fstatus != END_OF_FILE )
1366 "Reading data from " + m_fileName );
1371 errc = readHeader( ( *heads )[0] );
1374 return internal::mxlib_error_report<verboseT>( errc );
1379 for(
int i = 1; i < flist.size(); ++i )
1381 errc = fileName( flist[i], 1 );
1384 return internal::mxlib_error_report<verboseT>( errc );
1388 fits_read_subset( m_fptr,
1389 fitsType<typename cubeT::Scalar>(),
1394 (
void *)cube.image( i ).data(),
1398 if( fstatus && fstatus != END_OF_FILE )
1401 "Reading data from " + m_fileName );
1406 errc = readHeader( ( *heads )[i] );
1409 return internal::mxlib_error_report<verboseT>( errc );
1417template <
typename dataT,
class verboseT>
1418template <
typename cubeT>
1421 const std::vector<std::string> &flist )
1426template <
typename dataT,
class verboseT>
1431 char keyword[FLEN_KEYWORD];
1432 char value[FLEN_VALUE];
1436 typename std::list<headerIteratorT> head_keys;
1437 typename std::list<headerIteratorT>::iterator head_keys_it;
1440 bool head_keys_only =
false;
1441 if( head.
size() > 0 )
1443 head_keys_only =
true;
1444 headerIteratorT headIt = head.
begin();
1445 while( headIt != head.
end() )
1447 head_keys.push_back( headIt );
1460 comment =
new char[FLEN_COMMENT];
1472 fits_get_hdrspace( m_fptr, &keysexist, &morekeys, &fstatus );
1482 "Reading header from " + m_fileName );
1485 for(
int i = 0; i < keysexist; i++ )
1487 fits_read_keyn( m_fptr, i + 1, keyword, value, comment, &fstatus );
1497 "Reading header from " + m_fileName );
1500 if( !head_keys_only )
1502 if( strcmp( keyword,
"COMMENT" ) == 0 )
1504 head.template append<fitsCommentType>( keyword, fitsCommentType( value ), comment );
1506 else if( strcmp( keyword,
"HISTORY" ) == 0 )
1508 head.template append<fitsHistoryType>( keyword, fitsHistoryType( value ), comment );
1513 head.
append( keyword, value, comment );
1518 head_keys_it = head_keys.begin();
1519 while( head_keys_it != head_keys.end() )
1521 if( ( *( *head_keys_it ) ).keyword() == keyword )
1523 head[keyword].value( (
const char *)value );
1526 head[keyword].comment( comment );
1529 head_keys.
erase( head_keys_it );
1537 if( head_keys.empty() )
1552template <
typename dataT,
class verboseT>
1561template <
typename dataT,
class verboseT>
1563 const std::vector<std::string> &flist )
1565 if(heads.size() != 0 && heads.size() != flist.size())
1567 return internal::mxlib_error_report<verboseT>(
error_t::invalidarg,
"head vector is not empty and not same size as file list");
1570 if(heads.size() == 0)
1572 heads.resize(flist.size());
1576 for(
int i = 0; i < flist.size(); ++i )
1586template <
typename dataT,
class verboseT>
1610 if( m_naxes ==
nullptr )
1612 return internal::mxlib_error_report<verboseT>(
error_t::allocerr,
"m_naxes is nullptr" );
1625 std::string forceFileName =
"!" + m_fileName;
1627 fits_create_file( &m_fptr, forceFileName.c_str(), &fstatus );
1630 return internal::mxlib_error_report<verboseT>(
fits_status2error_t( fstatus ),
"Creating " + m_fileName );
1635 fits_create_img( m_fptr, fitsBITPIX<dataT>(), m_naxis, m_naxes, &fstatus );
1639 "Creating image in" + m_fileName );
1647 LONGLONG nelements = 1;
1649 for(
int i = 0; i < m_naxis && i < 3; ++i )
1651 nelements *= m_naxes[i];
1655 fits_write_pix( m_fptr, fitsType<dataT>(), fpixel, nelements, (
void *)im, &fstatus );
1658 return internal::mxlib_error_report<verboseT>(
fits_status2error_t( fstatus ),
"Writing data " + m_fileName );
1665 for( it = head->
begin(); it != head->
end(); ++it )
1667 error_t errc = it->write( m_fptr );
1671 "Writing keyword " + m_fileName );
1679template <
typename dataT,
class verboseT>
1685template <
typename dataT,
class verboseT>
1691template <
typename dataT,
class verboseT>
1699template <
typename dataT,
class verboseT>
1701 const std::string &fname,
const dataT *im,
int d1,
int d2,
int d3,
fitsHeader<verboseT> &head )
1708template <
typename dataT,
class verboseT>
1709template <
typename arrT>
1717template <
typename dataT,
class verboseT>
1718template <
typename arrT>
1723 return write( fname, im.data(), im.rows(), im.cols(), planes( im ), head );
1726template <
typename dataT,
class verboseT>
1735template <
typename dataT,
class verboseT>
1744template <
typename dataT,
class verboseT>
1751template <
typename dataT,
class verboseT>
1755 m_zframes = zframes;
Class to manage interactions with a FITS file.
long m_y0
The starting y-pixel to read from.
long m_zframes
The number of frames to read from a cube.
long m_xpix
The number of x-pixels to read.
int naxis()
Get the current value of m_naxis.
long getSize(error_t &errc)
Get the total size.
void setCubeReadSize()
Set to read all frames from a cube.
long naxes(int dim)
Get the current value of m_naxes for the specified dimension.
error_t open()
Open the file and gets its dimensions.
std::string fileName()
Get the current value of m_fileName.
std::string m_fileName
The path to the file.
fitsfile * m_fptr
The cfitsio data structure.
int getDimensions(error_t &errc)
Get the number of dimensions (i.e. m_naxis)
long m_z0
The starting frame to read from a cube.
void setReadSize()
Set to read all the pixels in the file.
long m_ypix
The number of y-pixels to read.
int m_noComment
Flag to control whether the comment string is read.
fitsFile()
Default constructor.
int m_naxis
The dimensions of the image (1D, 2D, or 3D)
dataT m_nulval
The value to replace null values with.
error_t close()
Close the file.
error_t readHeader(fitsHeader< verboseT > &head)
Read the header from the fits file.
error_t calcPixarrs(pixarrT &pixarr)
Fill in the read-size arrays for reading a subset (always used)
long m_naxes[3]
The size of each dimension.
long m_x0
The starting x-pixel to read from.
error_t read(dataT *data)
Read the contents of the FITS file into an array.
int m_anynul
Records whether any null values were replaced.
void construct()
One time initialization common to all constructors.
error_t write(const dataT *im, int d1, int d2, int d3, fitsHeader< verboseT > *head)
Write the contents of a raw array to the FITS file.
bool m_isOpen
Flag indicating whether the file is open or not.
Declares and defines utilities to work with FITS files.
error_t
The mxlib error codes.
@ noerror
No error has occurred.
@ std_exception
An exception was thrown.
@ exception
An exception was thrown.
@ std_bad_alloc
A bad allocation exception was thrown.
@ allocerr
An error occurred during memory allocation.
@ paramnotset
A parameter was not set.
@ invalidconfig
A config setting was invalid.
@ invalidarg
An argument was invalid.
#define mxlib_error_check(fxn)
Perform an error check, if an error occurs report it and return the error. Does not return on no erro...
#define mxlib_error_return(fxn)
Perform an error check, if an error occurs report it, and return the error code even if no error.
fitsFile< long > fitsFilel
A fitsFile to work in signed long integers.
fitsFile< unsigned short > fitsFileus
A fitsFile to work in unsigned short integers.
fitsFile< short > fitsFiles
A fitsFile to work in signed short integers.
fitsFile< char > fitsFilec
A fitsFile to work in signed characters.
fitsFile< int > fitsFilei
A fitsFile to work in signed integers.
fitsFile< unsigned char > fitsFileuc
A fitsFile to work in unsigned characters.
fitsFile< unsigned int > fitsFileui
A fitsFile to work in unsigned integers.
fitsFile< double > fitsFiled
A fitsFile to work in double precision.
fitsFile< float > fitsFilef
A fitsFile to work in single precision floats.
static constexpr error_t fits_status2error_t(const int &err)
Convert a FITS status code to error_t.
Function object to return the number of planes for any Eigen like object, whether 2D or a 3D cube.