fitsFile.hpp

Go to the documentation of this file.

/** \file fitsFile.hpp
 * \brief Declares and defines a class to work with a FITS file
 * \ingroup fits_processing_files
 * \author Jared R. Males (jaredmales@gmail.com)
 *
 */
 
//***********************************************************************//
// Copyright 2015-2022 Jared R. Males (jaredmales@gmail.com)
//
// This file is part of mxlib.
//
// mxlib is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// mxlib is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with mxlib.  If not, see <http://www.gnu.org/licenses/>.
//***********************************************************************//
 
#ifndef ioutils_fits_fitsFile_hpp
#define ioutils_fits_fitsFile_hpp
 
#include "../../mxError.hpp"
 
#include "../../improc/eigenImage.hpp"
 
#include "fitsUtils.hpp"
#include "fitsHeader.hpp"
 
namespace mx
{
namespace fits
{
 
/// Class to manage interactions with a FITS file
/** This class wraps the functionality of cfitsio.
 *
 * \tparam dataT is the datatype to use for in-memory storage of the image.  This does not have to match the data type
 * stored on disk for reading, but will be the type used for writing.
 *
 * \ingroup fits_processing
 */
template <typename dataT>
class fitsFile
{
 
    friend class fitsFile_test;
 
  protected:
    /// The path to the file
    std::string m_fileName;
 
    /// The cfitsio data structure
    fitsfile *m_fptr{ nullptr };
 
    /// The dimensions of the image (1D, 2D, 3D etc)
    int m_naxis;
 
    /// The size of each dimension
    long *m_naxes{ nullptr };
 
    /// Flag indicating whether the file is open or not
    bool m_isOpen{ false };
 
    /// The value to replace null values with
    dataT m_nulval{ 0 };
 
    /// Records whether any null values were replaced
    int m_anynul{ 0 };
 
    /// Flag to control whether the comment string is read.
    int m_noComment{ 0 };
 
    /// The starting x-pixel to read from
    long m_x0{ -1 };
 
    /// The starting y-pixel to read from
    long m_y0{ -1 };
 
    /// The number of x-pixels to read
    long m_xpix{ -1 };
 
    /// The number of y-pixels to read
    long m_ypix{ -1 };
 
    /// The starting frame to read from a cube
    long m_z0{ -1 };
 
    /// The number of frames to read from a cube
    long m_zframes{ -1 };
 
    /// One time initialization common to all constructors
    void construct();
 
  public:
    /// Default constructor
    fitsFile();
 
    /// Constructor with m_fileName, and option to open.
    fitsFile( const std::string &fname, ///< File name to set on construction
              bool doopen = true        ///< If true, then the file is opened (the default).
    );
 
    /// Destructor
    ~fitsFile();
 
    /// Get the current value of m_fileName
    /**
     * \returns the current value of m_fileName.
     */
    std::string fileName();
 
    /// Set the file path, and optionally open the file.
    /**
     * \returns 0 on success
     * \returns -1 on success
     */
    int fileName( const std::string &fname, ///< The new file name.
                  bool doopen = true        ///< If true, then the file is opened (the default).
    );
 
    /// Get the current value of m_naxis
    /**
     * \returns the current value of m_naxis
     *
     * \test Scenario: fitsFile calculating subimage sizes \ref tests_ioutils_fits_fitsFile_subimage_sizes "[test doc]"
     */
    int naxis();
 
    /// Get the current value of m_naxes for the specified dimension
    /**
     * \returns the current value of m_naxes for the specified dimension. -1 if no such dimension
     *
     * \test Scenario: fitsFile calculating subimage sizes \ref tests_ioutils_fits_fitsFile_subimage_sizes "[test doc]"
     */
    long naxes( int dim /**< [in] the dimension */ );
 
    /// Open the file and gets its dimensions.
    /** File name needs to already have been set.
     * If the file has already been opened, this returns immediately with no re-open.
     *
     * \returns 0 on success
     * \returns -1 on error
     */
    int open();
 
    /// Open the file, first setting the file path.
    /**
     * \returns 0 on success
     * \returns -1 on error
     */
    int open( const std::string &fname /**< The name of the file to open. */ );
 
    /// Close the file.
    /**
     * \returns 0 on success
     * \returns -1 on error
     */
    int close();
 
    /// Get the number of dimensions (i.e. m_naxis)
    int getDimensions();
 
    /// Get the total size
    long getSize();
 
    /// Get the size of a specific dimension
    long getSize( size_t axis );
 
    /** \name Reading Basic Arrays
     * These methods read FITS data into basic or raw arrays specified by a pointer.
     * @{
     */
 
  protected:
    /// Fill in the read-size arrays for reading a subset (always used)
    /** \note this allocates with new.  You are responsible for calling delete.
     *
     * \test Scenario: fitsFile calculating subimage sizes \ref tests_ioutils_fits_fitsFile_subimage_sizes "[test doc]"
     */
    void pixarrs( long **fpix, ///< Populated with the lower left pixel to read.  Is allocated.
                  long **lpix, ///< Populated with the upper right pixel to read. Is allocated.
                  long **inc   ///< The increment.  Is allocated.
    );
 
  public:
    /// Read the contents of the FITS file into an array.
    /** The array pointed to by data must have been allocated.
     *
     * \returns 0 on success
     * \returns -1 on error
     */
    int read( dataT *data /**< [out] an allocated arrray large enough to hold the entire image */ );
 
    /// Read the contents of the FITS file into an array.
    /** The array pointed to by data must have been allocated.
     *
     * \returns 0 on success
     * \returns -1 on error
     */
    int read( dataT *data,     ///< [out] an allocated arrray large enough to hold the entire image
              fitsHeader &head ///< [out] a fitsHeader object which is passed to \ref readHeader
    );
 
    /// Read the contents of the FITS file into an array.
    /** The array pointed to by data must have been allocated.
     *
     * \returns 0 on success
     * \returns -1 on error
     */
    int read( dataT *data,             ///< [out] is an allocated arrray large enough to hold the entire image
              const std::string &fname ///< [in] is the file path, which is passed to \ref fileName
    );
 
    /// Read the contents of the FITS file into an array and read the header.
    /** The array pointed to by data must have been allocated.
     *
     * \returns 0 on success
     * \returns -1 on error
     */
    int read( dataT *data,             ///< [out] an allocated arrray large enough to hold the entire image
              fitsHeader &head,        ///< [out] a fitsHeader object which is passed to \ref readHeader
              const std::string &fname ///< [in] the file path, which is passed to \ref fileName
    );
 
    /// Read data from a vector list of files into an image cube
    int read( dataT *im,                            ///< [out] An allocated array large enough to hold all the images
              const std::vector<std::string> &flist ///< [in] The list of files to read.
    );
 
    /// Read data from a vector of files into an image cube with individual headers
    int
    read( dataT *im,                      ///< [out] An allocated array large enough to hold all the images
          std::vector<fitsHeader> &heads, ///< [in.out] The vector of fits headers, allocated to contain one per image.
          const std::vector<std::string> &flist ///< [in] The list of files to read.
    );
 
    ///@}
 
    /** \name Reading Eigen Arrays
     * These methods read FITS data into array types with an Eigen-like interface.
     * @{
     */
 
    /// Read the contents of the FITS file into an Eigen array type (not a simple pointer).
    /** The type arrT can be any type with the following members defined:
     * - resize(int, int) (allocates memory)
     * - data() (returns a pointer to the underlying array)
     * - typedef arrT::Scalar (is the data type, does not have to be dataT)
     *
     * \tparam arrT is the type of array, see requirements above.
     *
     * \returns 0 on success
     * \returns -1 on error
     */
    template <typename arrT>
    int
    read( arrT &data /**< [out] is the array, which will be resized as necessary using its resize(int, int) member */ );
 
    /// Read the contents of the FITS file into an Eigen array type (not a simple pointer).
    /** The type arrT can be any type with the following members defined:
     * - resize(int, int) (allocates memory)
     * - data() (returns a pointer to the underlying array)
     * - typedef arrT::Scalar (is the data type, does not have to be dataT)
     *
     * \tparam arrT is the type of array, see requirements above.
     *
     * \returns 0 on success
     * \returns -1 on error
     */
    template <typename arrT>
    int read( arrT &data, ///< [out] is the array, which will be resized as necessary using its resize(int, int) member
              fitsHeader &head ///< [out] is a fitsHeader object which is passed to \ref readHeader
    );
 
    /// Read the contents of the FITS file into an Eigen array type (not a simple pointer).
    /** The type arrT can be any type with the following members defined:
     * - resize(int, int) (allocates memory)
     * - data() (returns a pointer to the underlying array)
     * - typedef arrT::Scalar (is the data type, does not have to be dataT)
     *
     * \tparam arrT is the type of array, see requirements above.
     *
     * \returns 0 on success
     * \returns -1 on error
     */
    template <typename arrT>
    int read( arrT &data, ///< [out] is the array, which will be resized as necessary using its resize(int, int) member
              const std::string &fname ///< [in] is the file path, which is passed to \ref fileName
    );
 
    /// Read the contents of the FITS file into an Eigen array type (not a simple pointer).
    /** The type arrT can be any type with the following members defined:
     * - resize(int, int) (allocates memory)
     * - data() (returns a pointer to the underlying array)
     * - typedef arrT::Scalar (is the data type, does not have to be dataT)
     *
     * \tparam arrT is the type of array, see requirements above.
     *
     * \returns 0 on success
     * \returns -1 on error
     */
    template <typename arrT>
    int read( arrT &data, ///< [out] the array, which will be resized as necessary using its resize(int, int) member
              fitsHeader &head,        ///< [out] a fitsHeader object which is passed to \ref readHeader
              const std::string &fname ///< [in] the file path, which is passed to \ref fileName
    );
 
    /// Read data from a vector list of files into an image cube
    /** The type cubeT can be any type with the following members defined:
     * - resize(int, int, int) (allocates memory)
     * - data() (returns a pointer to the underlying array)
     * - typedef cubeT::Scalar (is the data type, does not have to be dataT)
     *
     * \note The images must all be the same size, or all be as large or larger than the subset specified.
     *
     * \tparam cubeT is the type of array, see requirements above.
     *
     * \returns 0 on success
     * \returns -1 on error
     */
    template <typename cubeT>
    int read( cubeT &cube, ///< [out] A cube which will be resized using its resize(int, int, int) member.
              const std::vector<std::string> &flist, ///< [in] The list of files to read.
              std::vector<fitsHeader> *heads =
                  0 ///< [out] [optional] A vector of fits headers, allocated to contain one per image.
    );
 
    /// Read data from a vector of files into an image cube with individual headers
    /** The type cubeT can be any type with the following members defined:
     * - resize(int, int, int) (allocates memory)
     * - data() (returns a pointer to the underlying array)
     * - typedef cubeT::Scalar (is the data type, does not have to be dataT)
     *
     * \tparam cubeT is the type of array, see requirements above.
     *
     * \returns 0 on success
     * \returns -1 on error
     */
    template <typename cubeT>
    int read( cubeT &cube, ///< [out] A cube which will be resized using its resize(int, int, int) member.
              std::vector<fitsHeader> &heads, ///< [out] The vector of fits headers, allocated to contain one per image.
              const std::vector<std::string> &flist ///< [in] The list of files to read.
    );
 
    ///@}
 
    /** \name Reading Headers
     * @{
     */
 
    /// Read the header from the fits file.
    /** If head is not empty, then only the keywords already in head are updated.  Otherwise
     * the complete header is read.
     *
     * \returns 0 on success
     * \returns -1 on error
     */
    int readHeader( fitsHeader &head /**< [out] a fitsHeader object */ );
 
    /// Read the header from the fits file.
    /** If head is not empty, then only the keywords already in head are updated.  Otherwise
     * the complete header is read.
     *
     * \returns 0 on success
     * \returns -1 on error
     */
    int readHeader( fitsHeader &head,        ///< [out] a fitsHeader object
                    const std::string &fname ///< [in] the file path, which is passed to \ref fileName
    );
 
    /// Read the headers from a list of FITS files.
    /** In each case, if the header is not empty then only the keywords already in head are updated.  Otherwise
     * the complete header is read.
     *
     * \returns 0 on success
     * \returns -1 on error
     */
    int readHeader(
        std::vector<fitsHeader> &heads,       /// A vector of fitsHeader objects to read into.
        const std::vector<std::string> &flist ///< [in] A list of files, each of which is passed to \ref fileName
    );
 
    ///@}
 
    /** \name Writing Basic Arrays
     * These methods write basic arrays specified by a pointer to FITS.
     * @{
     */
 
    /// Write the contents of a raw array to the FITS file.
    /**
     * \returns 0 on success
     * \returns -1 on error
     */
    int write( const dataT *im, ///< [in] is the array
               int d1,          ///< [in] is the first dimension
               int d2,          ///< [in] is the second dimension
               int d3,          ///< [in] is the third dimenesion (minimum value is 1)
               fitsHeader *head ///< [in] a pointer to the header.  Set to 0 if not used.
    );
 
    /// Write the contents of a raw array to the FITS file.
    /**
     * \returns 0 on success
     * \returns -1 on error
     */
    int write( const dataT *im, ///< [in] is the array
               int d1,          ///< [in] is the first dimension
               int d2,          ///< [in] is the second dimension
               int d3           ///< [in] is the third dimenesion (minimum value is 1)
    );
 
    /// Write the contents of a raw array to the FITS file.
    /**
     * Note: the type of the array must match dataT
     *
     * \returns 0 on success
     * \returns -1 on error
     */
    int write( const dataT *im, ///< [in] is the array
               int d1,          ///< [in] is the first dimension
               int d2,          ///< [in] is the second dimension
               int d3,          ///< [in] is the third dimenesion (minimum value is 1)
               fitsHeader &head ///< [in] is the header
    );
 
    /// Write the contents of a raw array to the FITS file.
    /**
     * \returns 0 on success
     * \returns -1 on error
     */
    int write( const std::string &fname, ///< [in] is the name of the file.
               const dataT *im,          ///< [in] is the array
               int d1,                   ///< [in] is the first dimension
               int d2,                   ///< [in] is the second dimension
               int d3                    ///< [in] is the third dimenesion (minimum value is 1)
    );
 
    /// Write the contents of a raw array to the FITS file.
    /**
     * \returns 0 on success
     * \returns -1 on error
     */
    int write( const std::string &fname, ///< [in] is the name of the file.
               const dataT *im,          ///< [in] is the array
               int d1,                   ///< [in] is the first dimension
               int d2,                   ///< [in] is the second dimension
               int d3,                   ///< [in] is the third dimenesion (minimum value is 1)
               fitsHeader &head          ///< [in] is the header
    );
 
    ///@}
 
    /** \name Writing Eigen Arrays
     * These methods write array types with an Eigen-like interface.
     * @{
     */
 
    /// Write the contents of an Eigen-type array to a FITS file.
    /** The type arrT can be any type with the following members defined:
     * - data() (returns a pointer to the underlying array)
     * - rows() (returrns the number of rows)
     * - cols() (returns the number of columns)
     * - may have planes() defined
     *
     * Note: as with all write methods, the Scalar type of the array must match dataT
     *
     * \tparam arrT is the type of array, see requirements above.
     *
     * \returns 0 on success
     * \returns -1 on error
     */
    template <typename arrT>
    int write( const std::string &fname, ///< [in] is the name of the file.
               const arrT &im            ///< [in] is the array
    );
 
    /// Write the contents of an Eigen-type array to a FITS file.
    /** The type arrT can be any type with the following members defined:
     * - data() (returns a pointer to the underlying array)
     * - rows() (returrns the number of rows)
     * - cols() (returns the number of columns)
     * - may have planes() defined.
     *
     * Note: as with all write methods, the Scalar type of the array must match dataT
     *
     * \tparam arrT is the type of array, see requirements above.
     *
     * \returns 0 on success
     * \returns -1 on error
     */
    template <typename arrT>
    int write( const std::string &fname, ///< [in] is the file path, which is passed to \ref fileName
               const arrT &im,           ///< [in] is the array
               fitsHeader &head          ///< [in] is a fitsHeader object which is passed to \ref readHeader
    );
 
    // int writeHeader( fitsHeader &head );
 
    ///@}
 
    /** \name Reading Subsets
     * It is often desirable to read only a subset of an image or images into memory.  These methods
     * allow you to specify this.
     * @{
     */
 
    /// Set to read all the pixels in the file
    void setReadSize();
 
    /// Set to read only a subset of the pixels in the file
    /**
     *
     * \test Scenario: fitsFile calculating subimage sizes \ref tests_ioutils_fits_fitsFile_subimage_sizes "[test doc]"
     */
    void setReadSize( long x0,   ///< is the starting x-pixel to read
                      long y0,   ///< is the starting y-pixel to read
                      long xpix, ///< is the number of x-pixels to read
                      long ypix  ///< is the number of y-pixels to read
    );
 
    /// Set to read all frames from a cube.
    /**
     */
    void setCubeReadSize();
 
    /// Set the number of frames to read from a cube.
    /**
     *
     * \test Scenario: fitsFile calculating subimage sizes \ref tests_ioutils_fits_fitsFile_subimage_sizes "[test doc]"
     */
    void setCubeReadSize( long z0,     ///< is the starting frame to read
                          long zframes ///< is the number of frames to read
    );
 
    ///@}
 
}; // fitsFile
class fitsFile {…};
 
template <typename dataT>
void fitsFile<dataT>::construct()
{
}
void fitsFile<dataT>::construct() {…}
 
template <typename dataT>
fitsFile<dataT>::fitsFile()
{
    construct();
}
fitsFile<dataT>::fitsFile() {…}
 
template <typename dataT>
fitsFile<dataT>::fitsFile( const std::string &fname, bool doopen )
{
    construct();
    fileName( fname, doopen );
}
fitsFile<dataT>::fitsFile( const std::string &fname, bool doopen ) {…}
 
template <typename dataT>
fitsFile<dataT>::~fitsFile()
{
    if( m_isOpen )
        close();
 
    if( m_naxes )
        delete[] m_naxes;
}
fitsFile<dataT>::~fitsFile() {…}
 
template <typename dataT>
std::string fitsFile<dataT>::fileName()
{
    return m_fileName;
}
std::string fitsFile<dataT>::fileName() {…}
 
template <typename dataT>
int fitsFile<dataT>::fileName( const std::string &fname, bool doopen )
{
    if( m_isOpen )
    {
        close();
    }
 
    m_fileName = fname;
 
    if( doopen )
    {
        return open();
    }
 
    return 0;
}
int fitsFile<dataT>::fileName( const std::string &fname, bool doopen ) {…}
 
template <typename dataT>
int fitsFile<dataT>::naxis()
{
    return m_naxis;
}
int fitsFile<dataT>::naxis() {…}
 
template <typename dataT>
long fitsFile<dataT>::naxes( int dim )
{
    if( m_naxes == nullptr )
        return -1;
 
    if( dim >= m_naxis )
        return -1;
 
    return m_naxes[dim];
}
long fitsFile<dataT>::naxes( int dim ) {…}
 
template <typename dataT>
int fitsFile<dataT>::open()
{
    if( m_isOpen ) // no-op
    {
        return 0;
    }
 
    int fstatus = 0;
 
    fits_open_file( &m_fptr, m_fileName.c_str(), READONLY, &fstatus );
 
    if( fstatus )
    {
        std::string explan = "Error opening file";
        fitsErrText( explan, m_fileName, fstatus );
        //mxError( "fitsFile", MXE_FILEOERR, explan );
        mxThrowException(err::liberr, "fitsFile::open", explan);
 
        return -1;
    }
 
    fits_get_img_dim( m_fptr, &m_naxis, &fstatus );
    if( fstatus )
    {
        std::string explan = "Error getting number of axes in file";
        fitsErrText( explan, m_fileName, fstatus );
        //mxError( "fitsFile", MXE_FILERERR, explan );
        mxThrowException(err::liberr, "fitsFile::open", explan);
 
        return -1;
    }
 
    if( m_naxes )
        delete[] m_naxes;
    m_naxes = new long[m_naxis];
 
    fits_get_img_size( m_fptr, m_naxis, m_naxes, &fstatus );
    if( fstatus )
    {
        std::string explan = "Error getting dimensions in file";
        fitsErrText( explan, m_fileName, fstatus );
        //mxError( "fitsFile", MXE_FILERERR, explan );
        mxThrowException(err::liberr, "fitsFile::open", explan);
 
        return -1;
    }
 
    m_isOpen = true; // Only set this after opening is complete.
 
    return 0;
}
int fitsFile<dataT>::open() {…}
 
template <typename dataT>
int fitsFile<dataT>::open( const std::string &fname )
{
    return fileName( fname, true );
}
int fitsFile<dataT>::open( const std::string &fname ) {…}
 
template <typename dataT>
int fitsFile<dataT>::close()
{
    if( !m_isOpen )
    {
        return 0; // No error.
    }
 
    int fstatus = 0;
    fits_close_file( m_fptr, &fstatus );
 
    if( fstatus )
    {
        std::string explan = "Error closing file";
        fitsErrText( explan, m_fileName, fstatus );
        //mxError( "fitsFile", MXE_FILECERR, explan );
        mxThrowException(err::liberr, "fitsFile::close", explan);
 
        return -1;
    }
 
    m_isOpen = 0;
    fstatus = 0;
 
    if( m_naxes )
        delete[] m_naxes;
 
    m_naxes = nullptr;
 
    return 0;
}
int fitsFile<dataT>::close() {…}
 
template <typename dataT>
int fitsFile<dataT>::getDimensions()
{
    if( !m_isOpen or !m_naxes )
        return -1;
 
    return m_naxis;
}
int fitsFile<dataT>::getDimensions() {…}
 
template <typename dataT>
long fitsFile<dataT>::getSize()
{
    if( !m_isOpen or !m_naxes )
        return -1;
 
    long sz = 1;
 
    if( m_x0 > -1 && m_y0 > -1 && m_xpix > -1 && m_ypix > -1 && m_naxis == 2 )
    {
        return m_xpix * m_ypix;
    }
    else
    {
        for( int i = 0; i < m_naxis; ++i )
            sz *= m_naxes[i];
    }
 
    return sz;
}
long fitsFile<dataT>::getSize() {…}
 
template <typename dataT>
long fitsFile<dataT>::getSize( size_t axis )
{
    if( !m_isOpen or !m_naxes )
        return -1;
 
    if( m_x0 > -1 && m_y0 > -1 && m_xpix > -1 && m_ypix > -1 && m_naxis == 2 )
    {
        if( axis == 0 )
            return m_xpix;
        return m_ypix;
    }
    else
    {
        return m_naxes[axis];
    }
}
long fitsFile<dataT>::getSize( size_t axis ) {…}
 
template <typename dataT>
void fitsFile<dataT>::pixarrs( long **fpix, long **lpix, long **inc )
{
    *fpix = new long[m_naxis];
    *lpix = new long[m_naxis];
    *inc = new long[m_naxis];
 
    if( m_x0 > -1 && m_y0 > -1 && m_xpix > -1 && m_ypix > -1 && m_naxis == 2 )
    {
        ( *fpix )[0] = m_x0 + 1;
        ( *lpix )[0] = ( *fpix )[0] + m_xpix - 1;
        ( *fpix )[1] = m_y0 + 1;
        ( *lpix )[1] = ( *fpix )[1] + m_ypix - 1;
 
        ( *inc )[0] = 1;
        ( *inc )[1] = 1;
    }
    else
    {
        if( m_x0 < 0 && m_y0 < 0 && m_xpix < 0 && m_ypix < 0 && m_z0 < 0 && m_zframes < 0 )
        {
            for( int i = 0; i < m_naxis; i++ )
            {
                ( *fpix )[i] = 1;
                ( *lpix )[i] = m_naxes[i];
                ( *inc )[i] = 1;
            }
        }
        else
        {
            if( m_x0 > -1 && m_y0 > -1 && m_xpix > -1 && m_ypix > -1 )
            {
                ( *fpix )[0] = m_x0 + 1;
                ( *lpix )[0] = ( *fpix )[0] + m_xpix - 1;
                ( *fpix )[1] = m_y0 + 1;
                ( *lpix )[1] = ( *fpix )[1] + m_ypix - 1;
 
                ( *inc )[0] = 1;
                ( *inc )[1] = 1;
            }
            else
            {
                ( *fpix )[0] = 1;
                ( *lpix )[0] = m_naxes[0];
                ( *fpix )[1] = 1;
                ( *lpix )[1] = m_naxes[1];
 
                ( *inc )[0] = 1;
                ( *inc )[1] = 1;
            }
 
            if( m_z0 > -1 && m_zframes > -1 )
            {
                ( *fpix )[2] = m_z0 + 1;
                ( *lpix )[2] = ( *fpix )[2] + m_zframes - 1;
                ( *inc )[2] = 1;
            }
            else
            {
                ( *fpix )[2] = 1;
                ( *lpix )[2] = m_naxes[2];
                ( *inc )[2] = 1;
            }
        }
    }
}
void fitsFile<dataT>::pixarrs( long **fpix, long **lpix, long **inc ) {…}
 
/************************************************************/
/***                      Basic Arrays                    ***/
/************************************************************/
 
template <typename dataT>
int fitsFile<dataT>::read( dataT *data )
{
    int fstatus = 0;
 
    if( !m_isOpen )
    {
        if( open() < 0 )
            return -1;
    }
 
    long *fpix, *lpix, *inc;
    pixarrs( &fpix, &lpix, &inc );
 
    ///\todo test if there is a speed difference for full reads for fits_read_pix/subset
 
    // fits_read_pix(m_fptr, fitsType<dataT>(), fpix, nelements, (void *) &m_nulval,
    //(void *) data, &m_anynul, &fstatus);
    fits_read_subset(
        m_fptr, fitsType<dataT>(), fpix, lpix, inc, (void *)&m_nulval, (void *)data, &m_anynul, &fstatus );
 
    delete[] fpix;
    delete[] lpix;
    delete[] inc;
 
    if( fstatus && fstatus != 107 )
    {
        std::string explan = "Error reading data from file";
        fitsErrText( explan, m_fileName, fstatus );
        //mxError( "fitsFile", MXE_FILERERR, explan );
        mxThrowException(err::liberr, "fitsFile::read", explan);
 
        return -1;
    }
 
    return 0;
}
int fitsFile<dataT>::read( dataT *data ) {…}
 
template <typename dataT>
int fitsFile<dataT>::read( dataT *data, fitsHeader &head )
{
    if( read( data ) < 0 )
        return -1;
    if( readHeader( head ) < 0 )
        return -1;
 
    return 0;
}
int fitsFile<dataT>::read( dataT *data, fitsHeader &head ) {…}
 
template <typename dataT>
int fitsFile<dataT>::read( dataT *data, const std::string &fname )
{
    if( fileName( fname ) < 0 )
        return -1;
    if( read( data ) < 0 )
        return -1;
    return 0;
}
int fitsFile<dataT>::read( dataT *data, const std::string &fname ) {…}
 
template <typename dataT>
int fitsFile<dataT>::read( dataT *data, fitsHeader &head, const std::string &fname )
{
    if( fileName( fname ) < 0 )
        return -1;
    if( read( data ) < 0 )
        return -1;
    if( readHeader( head ) < 0 )
        return -1;
    return 0;
}
int fitsFile<dataT>::read( dataT *data, fitsHeader &head, const std::string &fname ) {…}
 
template <typename dataT>
int fitsFile<dataT>::read( dataT *im, const std::vector<std::string> &flist )
{
    if( flist.size() == 0 )
    {
        //mxError( "fitsFile", MXE_PARAMNOTSET,  );
        mxThrowException(err::paramnotset, "fitsFile::read", "Empty file list");
        return -1;
    }
 
    long sz0 = 0, sz1 = 0;
 
    for( int i = 0; i < flist.size(); ++i )
    {
        if( fileName( flist[i], 1 ) < 0 )
            return -1;
 
        if( read( im + i * sz0 * sz1 ) < 0 )
            return -1;
 
        sz0 = getSize( 0 );
        sz1 = getSize( 1 );
    }
 
    return 0;
}
int fitsFile<dataT>::read( dataT *im, const std::vector<std::string> &flist ) {…}
 
template <typename dataT>
int fitsFile<dataT>::read( dataT *im, std::vector<fitsHeader> &heads, const std::vector<std::string> &flist )
{
    if( flist.size() == 0 )
    {
        mxThrowException(err::paramnotset, "fitsFile::read", "Empty file list");
        return -1;
    }
 
    long sz0 = 0, sz1 = 0;
 
    for( size_t i = 0; i < flist.size(); ++i )
    {
        if( fileName( flist[i], 1 ) < 0 )
            return -1;
 
        if( read( im + i * sz0 * sz1 ) < 0 )
            return -1;
 
        if( readHeader( heads[i] ) < 0 )
            return -1;
 
        sz0 = getSize( 0 );
        sz1 = getSize( 1 );
    }
 
    return 0;
}
int fitsFile<dataT>::read( dataT *im, std::vector<fitsHeader> &heads, const std::vector<std::string> &flist ) {…}
 
/************************************************************/
/***                      Eigen Arrays                    ***/
/************************************************************/
 
template <typename arrT, bool isCube = improc::is_eigenCube<arrT>::value>
struct eigenArrResize
{
    // If it's a cube, always pass zsz
    void resize( arrT &arr, int xsz, int ysz, int zsz )
    {
        arr.resize( xsz, ysz, zsz );
    }
};
 
template <typename arrT>
struct eigenArrResize<arrT, false>
{
    // If it's not a cube, never pass zsz
    void resize( arrT &arr, int xsz, int ysz, int zsz )
    {
        static_cast<void>( zsz );
 
        arr.resize( xsz, ysz );
    }
};
 
template <typename dataT>
template <typename arrT>
int fitsFile<dataT>::read( arrT &im )
{
    ///\todo this can probably be made part of one read function (or call read(data *)) with a call to resize with
    ///SFINAE
    int fstatus = 0;
 
    if( !m_isOpen )
    {
        if( open() < 0 )
            return -1;
    }
 
    long *fpix, *lpix, *inc;
    pixarrs( &fpix, &lpix, &inc );
 
    eigenArrResize<arrT> arrresz;
    if( m_naxis > 2 )
    {
        arrresz.resize( im, lpix[0] - fpix[0] + 1, lpix[1] - fpix[1] + 1, lpix[2] - fpix[2] + 1 );
    }
    else if( m_naxis > 1 )
    {
        arrresz.resize( im, lpix[0] - fpix[0] + 1, lpix[1] - fpix[1] + 1, 1 );
    }
    else
    {
        arrresz.resize( im, lpix[0] - fpix[0] + 1, 1, 1 );
    }
 
    if( fits_read_subset( m_fptr,
                          fitsType<typename arrT::Scalar>(),
                          fpix,
                          lpix,
                          inc,
                          (void *)&m_nulval,
                          (void *)im.data(),
                          &m_anynul,
                          &fstatus ) < 0 )
        return -1;
 
    delete[] fpix;
    delete[] lpix;
    delete[] inc;
 
    if( fstatus && fstatus != 107 )
    {
        std::string explan = "Error reading data from file";
        fitsErrText( explan, m_fileName, fstatus );
        //mxError( "fitsFile", MXE_FILERERR, explan );
        mxThrowException(err::liberr, "fitsFile::read", explan);
        return -1;
    }
 
    return 0;
}
int fitsFile<dataT>::read( arrT &im ) {…}
 
template <typename dataT>
template <typename arrT>
int fitsFile<dataT>::read( arrT &data, fitsHeader &head )
{
    if( read( data ) < 0 )
        return -1;
    if( readHeader( head ) < 0 )
        return -1;
    return 0;
}
int fitsFile<dataT>::read( arrT &data, fitsHeader &head ) {…}
 
template <typename dataT>
template <typename arrT>
int fitsFile<dataT>::read( arrT &data, const std::string &fname )
{
    if( fileName( fname ) < 0 )
        return -1;
    if( read( data ) < 0 )
        return -1;
    return 0;
}
int fitsFile<dataT>::read( arrT &data, const std::string &fname ) {…}
 
template <typename dataT>
template <typename arrT>
int fitsFile<dataT>::read( arrT &data, fitsHeader &head, const std::string &fname )
{
    if( fileName( fname ) < 0 )
        return -1;
    if( read( data ) < 0 )
        return -1;
    if( readHeader( head ) < 0 )
        return -1;
 
    return 0;
}
int fitsFile<dataT>::read( arrT &data, fitsHeader &head, const std::string &fname ) {…}
 
template <typename dataT>
template <typename cubeT>
int fitsFile<dataT>::read( cubeT &cube, const std::vector<std::string> &flist, std::vector<fitsHeader> *heads )
{
    int fstatus = 0;
 
    if( flist.size() == 0 )
    {
        //mxError( "fitsFile", MXE_PARAMNOTSET, "Empty file list" );
        mxThrowException(err::paramnotset, "fitsFile::read", "Empty file list");
        return -1;
    }
 
    // Open the first file to get the dimensions.
    if( fileName( flist[0], 1 ) < 0 )
        return -1;
 
    long *fpix, *lpix, *inc;
    pixarrs( &fpix, &lpix, &inc );
 
    cube.resize( lpix[0] - fpix[0] + 1, lpix[1] - fpix[1] + 1, flist.size() );
 
    // Now read first image.
    fits_read_subset( m_fptr,
                      fitsType<typename cubeT::Scalar>(),
                      fpix,
                      lpix,
                      inc,
                      (void *)&m_nulval,
                      (void *)cube.image( 0 ).data(),
                      &m_anynul,
                      &fstatus );
 
    if( fstatus && fstatus != 107 )
    {
 
        std::string explan = "error reading data from file";
        fitsErrText( explan, m_fileName, fstatus );
 
        //mxError( "cfitsio", MXE_FILERERR, explan );
        mxThrowException(err::liberr, "fitsFile::read", explan);
 
        delete[] fpix;
        delete[] lpix;
        delete[] inc;
 
        return -1;
    }
 
    if( heads )
    {
        if( readHeader( ( *heads )[0] ) < 0 )
            return -1;
    }
 
    // Now read in the rest.
    for( int i = 1; i < flist.size(); ++i )
    {
        fileName( flist[i], 1 );
 
        // Now read image.
        fits_read_subset( m_fptr,
                          fitsType<typename cubeT::Scalar>(),
                          fpix,
                          lpix,
                          inc,
                          (void *)&m_nulval,
                          (void *)cube.image( i ).data(),
                          &m_anynul,
                          &fstatus );
 
        if( fstatus && fstatus != 107 )
        {
            std::string explan = "Error reading data from file";
            fitsErrText( explan, m_fileName, fstatus );
            //mxError( "cfitsio", MXE_FILERERR, explan );
            mxThrowException(err::liberr, "fitsFile::read", explan);
 
            delete[] fpix;
            delete[] lpix;
            delete[] inc;
 
            return -1;
        }
 
        if( heads )
        {
            if( readHeader( ( *heads )[i] ) < 0 )
                return -1;
        }
    }
 
    delete[] fpix;
    delete[] lpix;
    delete[] inc;
 
    return 0;
}
int fitsFile<dataT>::read( cubeT &cube, const std::vector<std::string> &flist, std::vector<fitsHeader> *heads ) {…}
 
template <typename dataT>
template <typename cubeT>
int fitsFile<dataT>::read( cubeT &cube, std::vector<fitsHeader> &heads, const std::vector<std::string> &flist )
{
    return read( cube, flist, &heads );
}
int fitsFile<dataT>::read( cubeT &cube, std::vector<fitsHeader> &heads, const std::vector<std::string> &flist ) {…}
 
template <typename dataT>
int fitsFile<dataT>::readHeader( fitsHeader &head )
{
    int fstatus = 0;
 
    char keyword[FLEN_KEYWORD];
    char value[FLEN_VALUE];
    char *comment;
 
    // The keys to look for if head is already populated
    std::list<fitsHeader::headerIterator> head_keys;
    std::list<fitsHeader::headerIterator>::iterator head_keys_it;
    //   int num_head_keys;
 
    bool head_keys_only = false;
    if( head.size() > 0 )
    {
        head_keys_only = true;
        fitsHeader::headerIterator headIt = head.begin();
        while( headIt != head.end() )
        {
            head_keys.push_back( headIt );
            ++headIt;
        }
        //      num_head_keys = head.size();
    }
 
    // If m_noComment is set, then we don't read in the comment
    if( m_noComment )
    {
        comment = 0;
    }
    else
    {
        comment = new char[FLEN_COMMENT];
    }
 
    int keysexist;
    int morekeys;
 
    if( !m_isOpen )
    {
        open();
    }
 
    // This gets the number of header keys to read
    fits_get_hdrspace( m_fptr, &keysexist, &morekeys, &fstatus );
 
    if( fstatus )
    {
        std::string explan = "Error reading header from file";
        fitsErrText( explan, m_fileName, fstatus );
        //mxError( "fitsFile", MXE_FILERERR, explan );
        mxThrowException(err::liberr, "fitsFile::readHeader", explan);
 
        if( comment )
            delete[] comment;
 
        return -1;
    }
 
    for( int i = 0; i < keysexist; i++ )
    {
        fits_read_keyn( m_fptr, i + 1, keyword, value, comment, &fstatus );
 
        if( fstatus )
        {
            std::string explan = "Error reading header from file";
            fitsErrText( explan, m_fileName, fstatus );
            //mxError( "fitsFile", MXE_FILERERR, explan );
            mxThrowException(err::liberr, "fitsFile::readHeader", explan);
 
            if( comment )
                delete[] comment;
 
            return -1;
        }
 
        if( !head_keys_only )
        {
            if( strcmp( keyword, "COMMENT" ) == 0 )
            {
                head.append<fitsCommentType>( keyword, fitsCommentType( value ), comment );
            }
            else if( strcmp( keyword, "HISTORY" ) == 0 )
            {
                head.append<fitsHistoryType>( keyword, fitsHistoryType( value ), comment );
            }
            else
            {
                // Otherwise we append it as an unknown type
                head.append( keyword, value, comment );
            }
        }
        else
        {
            head_keys_it = head_keys.begin();
            while( head_keys_it != head_keys.end() )
            {
                if( ( *( *head_keys_it ) ).keyword() == keyword )
                {
                    head[keyword].value( value );
                    if( comment )
                        head[keyword].comment( comment );
 
                    head_keys.erase( head_keys_it );
 
                    break;
                }
                ++head_keys_it;
            }
 
            // Quit if we're done.
            if( head_keys.empty() )
                break;
        }
    }
 
    if( comment )
        delete[] comment;
 
    return 0;
}
int fitsFile<dataT>::readHeader( fitsHeader &head ) {…}
 
template <typename dataT>
int fitsFile<dataT>::readHeader( fitsHeader &head, const std::string &fname )
{
    if( fileName( fname ) < 0 )
        return -1;
    if( readHeader( head ) < 0 )
        return -1;
}
int fitsFile<dataT>::readHeader( fitsHeader &head, const std::string &fname ) {…}
 
template <typename dataT>
int fitsFile<dataT>::readHeader( std::vector<fitsHeader> &heads, const std::vector<std::string> &flist )
{
    for( int i = 0; i < flist.size(); ++i )
    {
        fileName( flist[i], 1 );
 
        if( readHeader( heads[i] ) < 0 )
            return -1;
    }
 
    return 0;
}
int fitsFile<dataT>::readHeader( std::vector<fitsHeader> &heads, const std::vector<std::string> &flist ) {…}
 
template <typename dataT>
int fitsFile<dataT>::write( const dataT *im, int d1, int d2, int d3, fitsHeader *head )
{
    int fstatus = 0;
 
    if( m_isOpen )
        close();
    if( m_naxes )
        delete[] m_naxes;
 
    fstatus = 0;
    m_naxis = 1;
    if( d2 > 0 )
    {
        if( d3 > 1 )
            m_naxis = 3;
        else
            m_naxis = 2;
    }
 
    m_naxes = new long[m_naxis];
 
    m_naxes[0] = d1;
    if( m_naxis > 1 )
        m_naxes[1] = d2;
    if( m_naxis > 2 )
        m_naxes[2] = d3;
 
    std::string forceFileName = "!" + m_fileName;
 
    fits_create_file( &m_fptr, forceFileName.c_str(), &fstatus );
    if( fstatus )
    {
        std::string explan = "Error creating file";
        fitsErrText( explan, m_fileName, fstatus );
        //mxError( "fitsFile::write", MXE_FILEOERR, explan );
        mxThrowException(err::liberr, "fitsFile::write", explan);
 
        return -1;
    }
    m_isOpen = true;
 
    fits_create_img( m_fptr, fitsBITPIX<dataT>(), m_naxis, m_naxes, &fstatus );
    if( fstatus )
    {
        std::string explan = "Error creating image";
        fitsErrText( explan, m_fileName, fstatus );
        //mxError( "fitsFile::write", MXE_FILEWERR, explan );
        mxThrowException(err::liberr, "fitsFile::readHeader", explan);
 
        return -1;
    }
 
    long fpixel[3];
    fpixel[0] = 1;
    fpixel[1] = 1;
    fpixel[2] = 1;
 
    LONGLONG nelements = 1;
 
    for( int i = 0; i < m_naxis; ++i )
        nelements *= m_naxes[i];
 
    fits_write_pix( m_fptr, fitsType<dataT>(), fpixel, nelements, (void *)im, &fstatus );
    if( fstatus )
    {
        std::string explan = "Error writing data";
        fitsErrText( explan, m_fileName, fstatus );
        //mxError( "fitsFile::write", MXE_FILEWERR, explan );
        mxThrowException(err::liberr, "fitsFile::write", explan);
 
        return -1;
    }
 
    if( head != 0 )
    {
        fitsHeader::headerIterator it;
 
        for( it = head->begin(); it != head->end(); ++it )
        {
            int wrv = it->write( m_fptr );
            if( wrv != 0 )
            {
                std::string explan = "Error writing keyword";
                fitsErrText( explan, m_fileName, wrv );
                //mxError( "fitsFile::write", MXE_FILEWERR, explan );
                mxThrowException(err::liberr, "fitsFile::write", explan);
            }
        }
    }
 
    if( close() < 0 )
        return -1;
 
    return 0;
}
int fitsFile<dataT>::write( const dataT *im, int d1, int d2, int d3, fitsHeader *head ) {…}
 
template <typename dataT>
int fitsFile<dataT>::write( const dataT *im, int d1, int d2, int d3 )
{
    return write( im, d1, d2, d3, (fitsHeader *)0 );
}
int fitsFile<dataT>::write( const dataT *im, int d1, int d2, int d3 ) {…}
 
template <typename dataT>
int fitsFile<dataT>::write( const dataT *im, int d1, int d2, int d3, fitsHeader &head )
{
    return write( im, d1, d2, d3, &head );
}
int fitsFile<dataT>::write( const dataT *im, int d1, int d2, int d3, fitsHeader &head ) {…}
 
template <typename dataT>
int fitsFile<dataT>::write( const std::string &fname, const dataT *im, int d1, int d2, int d3 )
{
    if( fileName( fname, false ) < 0 )
        return -1;
    return write( im, d1, d2, d3, (fitsHeader *)0 );
}
int fitsFile<dataT>::write( const std::string &fname, const dataT *im, int d1, int d2, int d3 ) {…}
 
template <typename dataT>
int fitsFile<dataT>::write( const std::string &fname, const dataT *im, int d1, int d2, int d3, fitsHeader &head )
{
    if( fileName( fname, false ) < 0 )
        return -1;
    return write( im, d1, d2, d3, &head );
}
int fitsFile<dataT>::write( const std::string &fname, const dataT *im, int d1, int d2, int d3, fitsHeader &head ) {…}
 
template <typename dataT>
template <typename arrT>
int fitsFile<dataT>::write( const std::string &fname, const arrT &im )
{
    improc::eigenArrPlanes<arrT> planes;
 
    return write( fname, im.data(), im.rows(), im.cols(), planes( im ) );
}
int fitsFile<dataT>::write( const std::string &fname, const arrT &im ) {…}
 
template <typename dataT>
template <typename arrT>
int fitsFile<dataT>::write( const std::string &fname, const arrT &im, fitsHeader &head )
{
    improc::eigenArrPlanes<arrT> planes;
 
    return write( fname, im.data(), im.rows(), im.cols(), planes( im ), head );
}
int fitsFile<dataT>::write( const std::string &fname, const arrT &im, fitsHeader &head ) {…}
 
template <typename dataT>
void fitsFile<dataT>::setReadSize()
{
    m_x0 = -1;
    m_y0 = -1;
    m_xpix = -1;
    m_ypix = -1;
}
void fitsFile<dataT>::setReadSize() {…}
 
template <typename dataT>
void fitsFile<dataT>::setReadSize( long x0, long y0, long xpix, long ypix )
{
    m_x0 = x0;
    m_y0 = y0;
    m_xpix = xpix;
    m_ypix = ypix;
}
void fitsFile<dataT>::setReadSize( long x0, long y0, long xpix, long ypix ) {…}
 
template <typename dataT>
void fitsFile<dataT>::setCubeReadSize()
{
    m_z0 = -1;
    m_zframes = -1;
}
void fitsFile<dataT>::setCubeReadSize() {…}
 
template <typename dataT>
void fitsFile<dataT>::setCubeReadSize( long z0, long zframes )
{
    m_z0 = z0;
    m_zframes = zframes;
}
void fitsFile<dataT>::setCubeReadSize( long z0, long zframes ) {…}
 
/** \ingroup fits_processing_typedefs
 * @{
 */
 
/// A \ref fitsFile to work in signed characters
typedef fitsFile<char> fitsFilec;
 
/// A \ref fitsFile to work in unsigned characters
typedef fitsFile<unsigned char> fitsFileuc;
 
/// A \ref fitsFile to work in signed short integers
typedef fitsFile<short> fitsFiles;
 
/// A \ref fitsFile to work in unsigned short integers
typedef fitsFile<unsigned short> fitsFileus;
 
/// A \ref fitsFile to work in signed integers
typedef fitsFile<int> fitsFilei;
 
/// A \ref fitsFile to work in unsigned integers
typedef fitsFile<unsigned int> fitsFileui;
 
/// A \ref fitsFile to work in signed long integers
typedef fitsFile<long> fitsFilel;
 
/// A \ref fitsFile to work in single precision floats
typedef fitsFile<float> fitsFilef;
 
/// A \ref fitsFile to work in double precision
typedef fitsFile<double> fitsFiled;
 
///@}
 
} // namespace fits
} // namespace mx
 
#endif // ioutils_fits_fitsFile_hpp