/*****************************************************************************************
 *                                                                                       *
 * OpenSpace                                                                             *
 *                                                                                       *
 * Copyright (c) 2014-2019                                                               *
 *                                                                                       *
 * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
 * software and associated documentation files (the "Software"), to deal in the Software *
 * without restriction, including without limitation the rights to use, copy, modify,    *
 * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
 * permit persons to whom the Software is furnished to do so, subject to the following   *
 * conditions:                                                                           *
 *                                                                                       *
 * The above copyright notice and this permission notice shall be included in all copies *
 * or substantial portions of the Software.                                              *
 *                                                                                       *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
 * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
 * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
 * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
 * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
 * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
 ****************************************************************************************/

#include <modules/fitsfilereader/include/fitsfilereader.h>

#include <openspace/util/distanceconversion.h>
#include <ghoul/fmt.h>
#include <ghoul/logging/logmanager.h>
#include <ghoul/misc/dictionary.h>
#include <CCfits>
#include <fstream>

using namespace CCfits;

namespace {
    constexpr const char* _loggerCat = "FitsFileReader";
} // namespace

namespace openspace {

FitsFileReader::FitsFileReader(bool verboseMode) {
    _verboseMode = verboseMode;
    FITS::setVerboseMode(_verboseMode);
}

FitsFileReader::~FitsFileReader() {
    if (_infile) {
        _infile->destroy();
        _infile = nullptr;
    }
}

bool FitsFileReader::isPrimaryHDU() {
    return _infile->extension().empty();
}

template <typename T>
std::shared_ptr<ImageData<T>> FitsFileReader::readImage(const std::string& path) {
    try {
        _infile = std::make_unique<FITS>(path, Read, true);
        // Primary HDU Object
        if (isPrimaryHDU()) {
            return readImageInternal<T>(_infile->pHDU());
        }
        // Extension HDU Object
        return readImageInternal<T>(_infile->currentExtension());
    } catch (const FitsException& e){
        LERROR("Could not read FITS image from table. " + e.message() );
    }

    return nullptr;
}

template <typename T>
std::shared_ptr<std::unordered_map<std::string, T>> FitsFileReader::readHeader(
                                                       std::vector<std::string>& keywords)
{
    try {
        HDU& image = isPrimaryHDU() ?
            static_cast<HDU&>(_infile->pHDU()) :
            static_cast<HDU&>(_infile->currentExtension());

        std::vector<T> values;
        image.readKeys(keywords, values);

        if (values.size() != keywords.size()) {
            LERROR("Number of keywords does not match number of values");
        }

        std::unordered_map<std::string, T> result;
        std::transform(
            keywords.begin(),
            keywords.end(),
            values.begin(),
            std::inserter(
                result,
                result.end()
            ),
            [](std::string key, T value) { return std::make_pair(key, value); }
        );
        return std::make_shared<std::unordered_map<std::string, T>>(std::move(result));
    } catch (const FitsException& e) {
        LERROR("Could not read FITS header. " + e.message() );
    }
    return nullptr;
}

template <typename T>
std::shared_ptr<T> FitsFileReader::readHeaderValue(const std::string key) {
    try {
        HDU& image = isPrimaryHDU() ?
            static_cast<HDU&>(_infile->pHDU()) :
            static_cast<HDU&>(_infile->currentExtension());

        T value;
        image.readKey(key, value);
        return std::make_unique<T>(value);
    } catch (FitsException& e) {
        LERROR("Could not read FITS key. " + e.message() );
    }
    return nullptr;
}

template<typename T>
std::shared_ptr<TableData<T>> FitsFileReader::readTable(std::string& path,
                                              const std::vector<std::string>& columnNames,
                                                                             int startRow,
                                                                               int endRow,
                                                                               int hduIdx,
                                                                             bool readAll)
{
    // We need to lock reading when using multithreads because CCfits can't handle
    // multiple I/O drivers.
    std::lock_guard g(_mutex);

    try {
        _infile = std::make_unique<FITS>(path, Read, readAll);

        // Make sure FITS file is not a Primary HDU Object (aka an image).
        if (!isPrimaryHDU()) {
            ExtHDU& table = _infile->extension(hduIdx);
            int numCols = static_cast<int>(columnNames.size());
            int numRowsInTable = static_cast<int>(table.rows());
            std::unordered_map<string, std::vector<T>> contents;
            //LINFO("Read file: " + _infile->name());

            int firstRow = std::max(startRow, 1);

            if (endRow < firstRow) {
                endRow = numRowsInTable;
            }

            for (int i = 0; i < numCols; ++i) {
                std::vector<T> columnData;
                //LINFO("Read column: " + columnNames[i]);
                table.column(columnNames[i]).read(columnData, firstRow, endRow);
                contents[columnNames[i]] = columnData;
            }

            // Create TableData object of table contents.
            TableData<T> loadedTable = {
                std::move(contents),
                static_cast<int>(table.rows()),
                table.getRowsize(),
                table.name()
            };

            return std::make_shared<TableData<T>>(loadedTable);
        }
    }
    catch (FitsException& e) {
        LERROR(fmt::format(
            "Could not read FITS table from file '{}'. Make sure it's not an image file.",
            e.message()
        ));
    }
    return nullptr;
}

std::vector<float> FitsFileReader::readFitsFile(std::string filePath, int& nValuesPerStar,
                                                int firstRow, int lastRow,
                                               std::vector<std::string> filterColumnNames,
                                                                           int multiplier)
{
    std::vector<float> fullData;
    srand(1234567890);
    if (firstRow <= 0) {
        firstRow = 1;
    }

    // Define what columns to read.
    std::vector<std::string> allColumnNames = {
        "Position_X",
        "Position_Y",
        "Position_Z",
        "Velocity_X",
        "Velocity_Y",
        "Velocity_Z",
        "Gaia_Parallax",
        "Gaia_G_Mag",
        "Tycho_B_Mag",
        "Tycho_V_Mag",
        "Gaia_Parallax_Err",
        "Gaia_Proper_Motion_RA",
        "Gaia_Proper_Motion_RA_Err",
        "Gaia_Proper_Motion_Dec",
        "Gaia_Proper_Motion_Dec_Err",
        "Tycho_B_Mag_Err",
        "Tycho_V_Mag_Err"
    };

    // Append additional filter parameters to default rendering parameters.
    allColumnNames.insert(
        allColumnNames.end(),
        filterColumnNames.begin(),
        filterColumnNames.end()
    );

    std::string allNames = "Columns to read: \n";
    for (const std::string& colName : allColumnNames) {
        allNames += colName + "\n";
    }
    LINFO(allNames);

    // Read columns from FITS file. If rows aren't specified then full table will be read.
    std::shared_ptr<TableData<float>> table = readTable<float>(
        filePath,
        allColumnNames,
        firstRow,
        lastRow
    );

    if (!table) {
        throw ghoul::RuntimeError(fmt::format("Failed to open Fits file '{}'", filePath));
    }

    int nStars = table->readRows - firstRow + 1;

    int nNullArr = 0;
    int nColumnsRead = static_cast<int>(allColumnNames.size());
    int defaultCols = 17; // Number of columns that are copied by predefined code.
    if (nColumnsRead != defaultCols) {
        LINFO("Additional columns will be read! Consider add column in code for "
            "significant speedup!");
    }
    // Declare how many values to save per star
    nValuesPerStar = nColumnsRead + 1; // +1 for B-V color value.

    // Copy columns to local variables.
    std::unordered_map<std::string, std::vector<float>>& tableContent = table->contents;

    // Default render parameters!
    std::vector<float> posXcol = std::move(tableContent[allColumnNames[0]]);
    std::vector<float> posYcol = std::move(tableContent[allColumnNames[1]]);
    std::vector<float> posZcol = std::move(tableContent[allColumnNames[2]]);
    std::vector<float> velXcol = std::move(tableContent[allColumnNames[3]]);
    std::vector<float> velYcol = std::move(tableContent[allColumnNames[4]]);
    std::vector<float> velZcol = std::move(tableContent[allColumnNames[5]]);
    std::vector<float> parallax = std::move(tableContent[allColumnNames[6]]);
    std::vector<float> magCol = std::move(tableContent[allColumnNames[7]]);
    std::vector<float> tycho_b = std::move(tableContent[allColumnNames[8]]);
    std::vector<float> tycho_v = std::move(tableContent[allColumnNames[9]]);

    // Default filter parameters
    // Additional filter parameters are handled as well but slows down reading
    std::vector<float> parallax_err = std::move(tableContent[allColumnNames[10]]);
    std::vector<float> pr_mot_ra = std::move(tableContent[allColumnNames[11]]);
    std::vector<float> pr_mot_ra_err = std::move(tableContent[allColumnNames[12]]);
    std::vector<float> pr_mot_dec = std::move(tableContent[allColumnNames[13]]);
    std::vector<float> pr_mot_dec_err = std::move(tableContent[allColumnNames[14]]);
    std::vector<float> tycho_b_err = std::move(tableContent[allColumnNames[15]]);
    std::vector<float> tycho_v_err = std::move(tableContent[allColumnNames[16]]);

    // Construct data array. OBS: ORDERING IS IMPORTANT! This is where slicing happens.
    for (int i = 0; i < nStars * multiplier; ++i) {
        std::vector<float> values(nValuesPerStar);
        size_t idx = 0;

        // Default order for rendering:
        // Position [X, Y, Z]
        // Absolute Magnitude
        // B-V Color
        // Velocity [X, Y, Z]

        // Store positions.
        values[idx++] = posXcol[i % nStars];
        values[idx++] = posYcol[i % nStars];
        values[idx++] = posZcol[i % nStars];

        // Return early if star doesn't have a measured position.
        if (values[0] == -999 && values[1] == -999 && values[2] == -999) {
            nNullArr++;
            continue;
        }

        // Store color values.
        values[idx++] = magCol[i % nStars] == -999 ? 20.f : magCol[i % nStars];
        values[idx++] = tycho_b[i % nStars] - tycho_v[i % nStars];

        // Store velocity. Convert it to m/s with help by parallax.
        values[idx++] = convertMasPerYearToMeterPerSecond(
            velXcol[i % nStars],
            parallax[i % nStars]
        );
        values[idx++] = convertMasPerYearToMeterPerSecond(
            velYcol[i % nStars],
            parallax[i % nStars]
        );
        values[idx++] = convertMasPerYearToMeterPerSecond(
            velZcol[i % nStars],
            parallax[i % nStars]
        );

        // Store additional parameters to filter by.
        values[idx++] = parallax[i % nStars];
        values[idx++] = parallax_err[i % nStars];
        values[idx++] = pr_mot_ra[i % nStars];
        values[idx++] = pr_mot_ra_err[i % nStars];
        values[idx++] = pr_mot_dec[i % nStars];
        values[idx++] = pr_mot_dec_err[i % nStars];
        values[idx++] = tycho_b[i % nStars];
        values[idx++] = tycho_b_err[i % nStars];
        values[idx++] = tycho_v[i % nStars];
        values[idx++] = tycho_v_err[i % nStars];

        // Read extra columns, if any. This will slow down the sorting tremendously!
        for (size_t col = defaultCols; col < nColumnsRead; ++col) {
            std::vector<float> vecData = std::move(tableContent[allColumnNames[col]]);
            values[idx++] = vecData[i];
        }

        for (int j = 0; j < nValuesPerStar; ++j) {
            // The astronomers in Vienna use -999 as default value. Change it to 0.
            if (values[j] == -999) {
                values[j] = 0.f;
            }
            else if (multiplier > 1) {
                values[j] *= static_cast<float>(rand()) / static_cast<float>(RAND_MAX);
            }
        }

        fullData.insert(fullData.end(), values.begin(), values.end());
    }

    // Define what columns to read.
    /*auto allColumnNames = std::vector<std::string>({
        "ra",
        "dec",
        "parallax",
        "pmra",
        "pmdec",
        "phot_g_mean_mag",
        "phot_bp_mean_mag",
        "phot_rp_mean_mag",
        "radial_velocity",
        });
    // Append additional filter parameters to default rendering parameters.
    allColumnNames.insert(allColumnNames.end(), filterColumnNames.begin(),
        filterColumnNames.end());

    std::string allNames = "Columns to read: \n";
    for (auto colName : allColumnNames) {
        allNames += colName + "\n";
    }
    LINFO(allNames);

    // Read columns from FITS file. If rows aren't specified then full table will be read.
    auto table = readTable<float>(filePath, allColumnNames, firstRow, lastRow);

    if (!table) {
        throw ghoul::RuntimeError(fmt::format("Failed to open Fits file '{}'", filePath));
    }

    int nStars = table->readRows - firstRow + 1;

    int nNullArr = 0;
    size_t nColumnsRead = allColumnNames.size();
    size_t defaultCols = 9; // Number of columns that are copied by predefined code.
    if (nColumnsRead != defaultCols) {
        LINFO("Additional columns will be read! Consider add column in code for "
            "significant speedup!");
    }
    // Declare how many values to save per star
    nValuesPerStar = 8;

    // Copy columns to local variables.
    std::unordered_map<std::string, std::vector<float>>& tableContent = table->contents;

    std::vector<float> ra = std::move(tableContent[allColumnNames[0]]);
    std::vector<float> dec = std::move(tableContent[allColumnNames[1]]);
    std::vector<float> parallax = std::move(tableContent[allColumnNames[2]]);
    std::vector<float> pmra = std::move(tableContent[allColumnNames[3]]);
    std::vector<float> pmdec = std::move(tableContent[allColumnNames[4]]);
    std::vector<float> meanMagG = std::move(tableContent[allColumnNames[5]]);
    std::vector<float> meanMagBp = std::move(tableContent[allColumnNames[6]]);
    std::vector<float> meanMagRp = std::move(tableContent[allColumnNames[7]]);
    std::vector<float> radial_vel = std::move(tableContent[allColumnNames[8]]);

    // Construct data array. OBS: ORDERING IS IMPORTANT! This is where slicing happens.
    for (int i = 0; i < nStars; ++i) {
        std::vector<float> values(nValuesPerStar);
        size_t idx = 0;

        // Default order for rendering:
        // Position [X, Y, Z]
        // Mean G-band Magnitude
        // Bp-Rp Color
        // Velocity [X, Y, Z]

        // Return early if star doesn't have a measured position.
        if (std::isnan(ra[i]) || std::isnan(dec[i])) {
            nNullArr++;
            continue;
        }

        // Store positions. Set to a default distance if parallax doesn't exist.
        float radiusInKiloParsec = 9.0;
        if (!std::isnan(parallax[i])) {
            // Parallax is in milliArcseconds -> distance in kiloParsecs
            // https://gea.esac.esa.int/archive/documentation/GDR2/Gaia_archive/
            // chap_datamodel/sec_dm_main_tables/ssec_dm_gaia_source.html
            radiusInKiloParsec = 1.0 / parallax[i];
        }
        // Convert to Galactic Coordinates from Galactic Lon & Lat.
        // https://gea.esac.esa.int/archive/documentation/GDR2/Data_processing/
        // chap_cu3ast/sec_cu3ast_intro/ssec_cu3ast_intro_tansforms.html#SSS1
        //values[idx++] = radiusInKiloParsec * cos(glm::radians(b_latitude[i])) *
        //cos(glm::radians(l_longitude[i])); // Pos X
        //values[idx++] = radiusInKiloParsec * cos(glm::radians(b_latitude[i])) *
        //sin(glm::radians(l_longitude[i])); // Pos Y
        //values[idx++] = radiusInKiloParsec * sin(glm::radians(b_latitude[i])); // Pos Z


        // Convert ICRS Equatorial Ra and Dec to Galactic latitude and longitude.
        glm::mat3 aPrimG = glm::mat3(
            // Col 0
            glm::vec3(-0.0548755604162154, 0.4941094278755837, -0.8676661490190047),
            // Col 1
            glm::vec3(-0.8734370902348850, -0.4448296299600112, -0.1980763734312015),
            // Col 2
            glm::vec3(-0.4838350155487132, 0.7469822444972189, 0.4559837761750669)
        );
        glm::vec3 rICRS = glm::vec3(
            cos(glm::radians(ra[i])) * cos(glm::radians(dec[i])),
            sin(glm::radians(ra[i])) * cos(glm::radians(dec[i])),
            sin(glm::radians(dec[i]))
        );
        glm::vec3 rGal = aPrimG * rICRS;
        values[idx++] = radiusInKiloParsec * rGal.x; // Pos X
        values[idx++] = radiusInKiloParsec * rGal.y; // Pos Y
        values[idx++] = radiusInKiloParsec * rGal.z; // Pos Z

        // Store magnitude render value. (Set default to high mag = low brightness)
        values[idx++] = std::isnan(meanMagG[i]) ? 20.f : meanMagG[i]; // Mean G-band Mag

        // Store color render value. (Default value is bluish stars)
        values[idx++] = std::isnan(meanMagBp[i]) && std::isnan(meanMagRp[i]) ? 0.f
            : meanMagBp[i] - meanMagRp[i]; // Bp-Rp Color

        // Store velocity.
        if (std::isnan(pmra[i])) pmra[i] = 0.f;
        if (std::isnan(pmdec[i])) pmdec[i] = 0.f;

        // Convert Proper Motion from ICRS [Ra,Dec] to Galactic Tanget Vector [l,b].
        glm::vec3 uICRS = glm::vec3(
            -sin(glm::radians(ra[i])) * pmra[i] -
            cos(glm::radians(ra[i])) * sin(glm::radians(dec[i])) * pmdec[i],
            cos(glm::radians(ra[i])) * pmra[i] -
            sin(glm::radians(ra[i])) * sin(glm::radians(dec[i])) * pmdec[i],
            cos(glm::radians(dec[i]))  * pmdec[i]
        );
        glm::vec3 pmVecGal = aPrimG * uICRS;

        // Convert to Tangential vector [m/s] from Proper Motion vector [mas/yr]
        float tanVelX = 1000.0 * 4.74 * radiusInKiloParsec * pmVecGal.x;
        float tanVelY = 1000.0 * 4.74 * radiusInKiloParsec * pmVecGal.y;
        float tanVelZ = 1000.0 * 4.74 * radiusInKiloParsec * pmVecGal.z;

        // Calculate True Space Velocity [m/s] if we have the radial velocity
        if (!std::isnan(radial_vel[i])) {
            // Calculate Radial Velocity in the direction of the star.
            // radial_vel is given in [km/s] -> convert to [m/s].
            float radVelX = 1000.0 * radial_vel[i] * rGal.x;
            float radVelY = 1000.0 * radial_vel[i] * rGal.y;
            float radVelZ = 1000.0 * radial_vel[i] * rGal.z;

            // Use Pythagoras theorem for the final Space Velocity [m/s].
            values[idx++] = sqrt(pow(radVelX, 2) + pow(tanVelX, 2)); // Vel X [U]
            values[idx++] = sqrt(pow(radVelY, 2) + pow(tanVelY, 2)); // Vel Y [V]
            values[idx++] = sqrt(pow(radVelZ, 2) + pow(tanVelZ, 2)); // Vel Z [W]
        }
        // Otherwise use the vector [m/s] we got from proper motion.
        else {
            radial_vel[i] = 0.f;
            values[idx++] = tanVelX; // Vel X [U]
            values[idx++] = tanVelY; // Vel Y [V]
            values[idx++] = tanVelZ; // Vel Z [W]
        }

        fullData.insert(fullData.end(), values.begin(), values.end());
    }*/

    LINFO(fmt::format("{} out of {} read stars were null arrays", nNullArr, nStars));
    LINFO(fmt::format("Multiplier: {}", multiplier));

    return fullData;
}

std::vector<float> FitsFileReader::readSpeckFile(const std::string& filePath,
                                                 int& nRenderValues)
{
    std::vector<float> fullData;

    std::ifstream fileStream(filePath);

    if (!fileStream.good()) {
        LERROR(fmt::format("Failed to open Speck file '{}'", filePath));
        return fullData;
    }

    int nValuesPerStar = 0;
    int nNullArr = 0;
    size_t nStars = 0;

    // The beginning of the speck file has a header that either contains comments
    // (signaled by a preceding '#') or information about the structure of the file
    // (signaled by the keywords 'datavar', 'texturevar', 'texture' and 'maxcomment')
    std::string line;
    while (true) {
        std::streampos position = fileStream.tellg();
        std::getline(fileStream, line);

        if (line[0] == '#' || line.empty()) {
            continue;
        }

        if (line.substr(0, 7) != "datavar" && line.substr(0, 10) != "texturevar" &&
            line.substr(0, 7) != "texture" && line.substr(0, 10) != "maxcomment")
        {
            // We read a line that doesn't belong to the header, so we have to jump back
            // before the beginning of the current line.
            fileStream.seekg(position);
            break;
        }

        if (line.substr(0, 7) == "datavar") {
            // datavar lines are structured as follows:
            // datavar # description
            // where # is the index of the data variable; so if we repeatedly overwrite
            // the 'nValues' variable with the latest index, we will end up with the total
            // number of values (+3 since X Y Z are not counted in the Speck file index)
            std::stringstream str(line);

            std::string dummy;
            str >> dummy;
            str >> nValuesPerStar;
            nValuesPerStar += 1; // We want the number, but the index is 0 based
        }
    }

    nValuesPerStar += 3; // X Y Z are not counted in the Speck file indices

    // Order in DR1 file:       DR2 - GaiaGroupMembers:
    // 0 BVcolor                0 color
    // 1 lum                    1 lum
    // 2 Vabsmag                2 absmag
    // 3 Vappmag                3 Gmag
    // 4 distly                 4 distpc
    // 5 distpcPctErr           5 plx
    // 6 U                      6 ra
    // 7 V                      7 dec
    // 8 W                      8 RadVel
    // 9 speed                  9 Teff
    // 10 sptypeindex           10 vx
    // 11 lumclassindex         11 vy
    // 12 catsource             12 vz
    // 13 texture               13 speed
    //                          14 texture

    do {
        std::vector<float> readValues(nValuesPerStar);
        nStars++;

        std::getline(fileStream, line);
        std::stringstream str(line);

        // Read values.
        for (int i = 0; i < nValuesPerStar; ++i) {
            str >> readValues[i];
        }

        // Check if star is a nullArray.
        bool nullArray = true;
        for (float f : readValues) {
            if (f != 0.0) {
                nullArray = false;
                break;
            }
        }

        // Insert to data if we found some values.
        if (!nullArray) {
            // Re-order data here because Octree expects the data in correct order when
            // read.
            // Default order for rendering:
            // Position [X, Y, Z]
            // Absolute Magnitude
            // B-V Color
            // Velocity [X, Y, Z]

            nRenderValues = 8;
            std::vector<float> renderValues(nRenderValues);

            // Gaia DR1 data from AMNH measures positions in Parsec, but
            // RenderableGaiaStars expects kiloParsec (because fits file from Vienna had
            // in kPc).
            // Thus we need to convert positions twice atm.
            renderValues[0] = readValues[0] / 1000.0; // PosX
            renderValues[1] = readValues[1] / 1000.0; // PosY
            renderValues[2] = readValues[2] / 1000.0; // PosZ
            renderValues[3] = readValues[6]; // AbsMag
            renderValues[4] = readValues[3]; // color
            renderValues[5] = readValues[13] * readValues[16]; // Vel X
            renderValues[6] = readValues[14] * readValues[16]; // Vel Y
            renderValues[7] = readValues[15] * readValues[16]; // Vel Z

            fullData.insert(fullData.end(), renderValues.begin(), renderValues.end());
        }
        else {
            nNullArr++;
        }

    } while (!fileStream.eof());

    LINFO(fmt::format("{} out of {} read stars were null arrays", nNullArr, nStars));

    return fullData;
}

// This is pretty annoying, the read method is not derived from the HDU class
// in CCfits - need to explicitly cast to the sub classes to access read
template<typename T>
const std::shared_ptr<ImageData<T>> FitsFileReader::readImageInternal(ExtHDU& image) {
   try {
        std::valarray<T> contents;
        image.read(contents);
        ImageData<T> im = { std::move(contents), image.axis(0), image.axis(1) };
        return std::make_shared<ImageData<T>>(im);
    } catch (const FitsException& e){
        LERROR("Could not read FITS image EXTHDU. " + e.message() );
    }
    return nullptr;
}

template<typename T>
const std::shared_ptr<ImageData<T>> FitsFileReader::readImageInternal(PHDU& image) {
    try {
        std::valarray<T> contents;
        image.read(contents);
        ImageData<T> im = { std::move(contents), image.axis(0), image.axis(1) };
        return std::make_shared<ImageData<T>>(im);
    } catch (const FitsException& e){
        LERROR("Could not read FITS image PHDU. " + e.message() );
    }
    return nullptr;
}

} // namespace openspace