/******************************************************************************

 * Project:  Selafin importer

 * Purpose:  Implementation of functions for reading records in Selafin files

 * Author:   François Hissel, francois.hissel@gmail.com

 *

 ******************************************************************************

 * Copyright (c) 2014,  François Hissel <francois.hissel@gmail.com>

 *

 * SPDX-License-Identifier: MIT

 ****************************************************************************/

#include "io_selafin.h"

#include "cpl_vsi.h"

#include "cpl_conv.h"

#include "cpl_port.h"

#include "cpl_error.h"

#include "cpl_quad_tree.h"

namespace Selafin

{

const char SELAFIN_ERROR_MESSAGE[] = "Error when reading Selafin file\n";

struct Point

{

    int nIndex;

    const Header *poHeader;

};

static void GetBoundsFunc(const void *hFeature, CPLRectObj *poBounds)

{

    const Point *poPoint = (const Point *)hFeature;

    poBounds->minx = poPoint->poHeader->paadfCoords[0][poPoint->nIndex];

    poBounds->maxx = poPoint->poHeader->paadfCoords[0][poPoint->nIndex];

    poBounds->miny = poPoint->poHeader->paadfCoords[1][poPoint->nIndex];

    poBounds->maxy = poPoint->poHeader->paadfCoords[1][poPoint->nIndex];

}

static int DumpFeatures(void *pElt, void * /* pUserData */)

{

    Point *poPoint = (Point *)pElt;

    delete poPoint;

    return TRUE;

}

/****************************************************************/

/*                         Header                               */

/****************************************************************/

Header::Header()

    : nHeaderSize(0), nStepSize(0), nMinxIndex(-1), nMaxxIndex(-1),

      nMinyIndex(-1), nMaxyIndex(-1), bTreeUpdateNeeded(true), nFileSize(0),

      fp(nullptr), pszFilename(nullptr), pszTitle(nullptr), nVar(0),

      papszVariables(nullptr), nPoints(0), nElements(0), nPointsPerElement(0),

      panConnectivity(nullptr), poTree(nullptr), panBorder(nullptr),

      panStartDate(nullptr), nSteps(0), nEpsg(0)

{

    paadfCoords[0] = nullptr;

    paadfCoords[1] = nullptr;

    for (size_t i = 0; i < 7; ++i)

        anUnused[i] = 0;

    adfOrigin[0] = 0.0;

    adfOrigin[1] = 0.0;

}

Header::~Header()

{

    CPLFree(pszFilename);

    CPLFree(pszTitle);

    if (papszVariables != nullptr)

    {

        for (int i = 0; i < nVar; ++i)

            CPLFree(papszVariables[i]);

        CPLFree(papszVariables);

    }

    CPLFree(panConnectivity);

    CPLFree(panBorder);

    if (poTree != nullptr)

    {

        CPLQuadTreeForeach(poTree, DumpFeatures, nullptr);

        CPLQuadTreeDestroy(poTree);

    }

    CPLFree(panStartDate);

    for (size_t i = 0; i < 2; ++i)

        CPLFree(paadfCoords[i]);

    if (fp != nullptr)

        VSIFCloseL(fp);

}

void Header::setUpdated()

{

    nHeaderSize = 88 + 16 + nVar * 40 + 12 * 4 +

                  ((panStartDate == nullptr) ? 0 : 32) + 24 +

                  (nElements * nPointsPerElement + 2) * 4 + (nPoints + 2) * 12;

    nStepSize = 12 + nVar * (nPoints + 2) * 4;

}

int Header::getPosition(int nStep, int nFeature, int nAttribute) const

{

    int a = (nFeature != -1 || nAttribute != -1)

                ? (12 + nAttribute * (nPoints + 2) * 4 + 4 + nFeature * 4)

                : 0;

    int b = nStep * nStepSize;

    return nHeaderSize + b + a;

}

CPLRectObj *Header::getBoundingBox() const

{

    CPLRectObj *poBox = new CPLRectObj;

    poBox->minx = paadfCoords[0][nMinxIndex];

    poBox->maxx = paadfCoords[0][nMaxxIndex];

    poBox->miny = paadfCoords[1][nMinyIndex];

    poBox->maxy = paadfCoords[1][nMaxyIndex];

    return poBox;

}

void Header::updateBoundingBox()

{

    if (nPoints > 0)

    {

        nMinxIndex = 0;

        for (int i = 1; i < nPoints; ++i)

            if (paadfCoords[0][i] < paadfCoords[0][nMinxIndex])

                nMinxIndex = i;

        nMaxxIndex = 0;

        for (int i = 1; i < nPoints; ++i)

            if (paadfCoords[0][i] > paadfCoords[0][nMaxxIndex])

                nMaxxIndex = i;

        nMinyIndex = 0;

        for (int i = 1; i < nPoints; ++i)

            if (paadfCoords[1][i] < paadfCoords[1][nMinyIndex])

                nMinyIndex = i;

        nMaxyIndex = 0;

        for (int i = 1; i < nPoints; ++i)

            if (paadfCoords[1][i] > paadfCoords[1][nMaxyIndex])

                nMaxyIndex = i;

    }

}

int Header::getClosestPoint(const double &dfx, const double &dfy,

                            const double &dfMax)

{

    // If there is no quad-tree of the points, build it now

    if (bTreeUpdateNeeded)

    {

        if (poTree != nullptr)

        {

            CPLQuadTreeForeach(poTree, DumpFeatures, nullptr);

            CPLQuadTreeDestroy(poTree);

        }

    }

    if (bTreeUpdateNeeded || poTree == nullptr)

    {

        bTreeUpdateNeeded = false;

        CPLRectObj *poBB = getBoundingBox();

        poTree = CPLQuadTreeCreate(poBB, GetBoundsFunc);

        delete poBB;

        CPLQuadTreeSetBucketCapacity(poTree, 2);

        for (int i = 0; i < nPoints; ++i)

        {

            Point *poPoint = new Point;

            poPoint->poHeader = this;

            poPoint->nIndex = i;

            CPLQuadTreeInsert(poTree, (void *)poPoint);

        }

    }

    // Now we can look for the nearest neighbour using this tree

    int nIndex = -1;

    CPLRectObj poObj;

    poObj.minx = dfx - dfMax;

    poObj.maxx = dfx + dfMax;

    poObj.miny = dfy - dfMax;

    poObj.maxy = dfy + dfMax;

    int nFeatureCount = 0;

    void **phResults = CPLQuadTreeSearch(poTree, &poObj, &nFeatureCount);

    if (nFeatureCount <= 0)

        return -1;

    double dfMin = dfMax * dfMax;

    for (int i = 0; i < nFeatureCount; ++i)

    {

        Point *poPoint = (Point *)(phResults[i]);

        double dfa = dfx - poPoint->poHeader->paadfCoords[0][poPoint->nIndex];

        dfa *= dfa;

        if (dfa >= dfMin)

            continue;

        const double dfb =

            dfy - poPoint->poHeader->paadfCoords[1][poPoint->nIndex];

        const double dfc = dfa + dfb * dfb;

        if (dfc < dfMin)

        {

            dfMin = dfc;

            nIndex = poPoint->nIndex;

        }

    }

    CPLFree(phResults);

    return nIndex;

}

void Header::addPoint(const double &dfx, const double &dfy)

{

    // We add the point to all the tables

    nPoints++;

    for (size_t i = 0; i < 2; ++i)

        paadfCoords[i] =

            (double *)CPLRealloc(paadfCoords[i], sizeof(double) * nPoints);

    paadfCoords[0][nPoints - 1] = dfx;

    paadfCoords[1][nPoints - 1] = dfy;

    panBorder = (int *)CPLRealloc(panBorder, sizeof(int) * nPoints);

    panBorder[nPoints - 1] = 0;

    // We update the bounding box

    if (nMinxIndex == -1 || dfx < paadfCoords[0][nMinxIndex])

        nMinxIndex = nPoints - 1;

    if (nMaxxIndex == -1 || dfx > paadfCoords[0][nMaxxIndex])

        nMaxxIndex = nPoints - 1;

    if (nMinyIndex == -1 || dfy < paadfCoords[1][nMinyIndex])

        nMinyIndex = nPoints - 1;

    if (nMaxyIndex == -1 || dfy > paadfCoords[1][nMaxyIndex])

        nMaxyIndex = nPoints - 1;

    // We update some parameters of the header

    bTreeUpdateNeeded = true;

    setUpdated();

}

void Header::removePoint(int nIndex)

{

    // We remove the point from all the tables

    nPoints--;

    for (size_t i = 0; i < 2; ++i)

    {

        for (int j = nIndex; j < nPoints; ++j)

            paadfCoords[i][j] = paadfCoords[i][j + 1];

        paadfCoords[i] =

            (double *)CPLRealloc(paadfCoords[i], sizeof(double) * nPoints);

    }

    for (int j = nIndex; j < nPoints; ++j)

        panBorder[j] = panBorder[j + 1];

    panBorder = (int *)CPLRealloc(panBorder, sizeof(int) * nPoints);

    // We must also remove all the elements referencing the deleted feature,

    // otherwise the file will not be consistent any inter

    int nOldElements = nElements;

    for (int i = 0; i < nElements; ++i)

    {

        bool bReferencing = false;

        int *panTemp = panConnectivity + i * nPointsPerElement;

        for (int j = 0; j < nPointsPerElement; ++j)

            bReferencing |= (panTemp[j] == nIndex + 1);

        if (bReferencing)

        {

            nElements--;

            for (int j = i; j < nElements; ++j)

                for (int k = 0; k < nPointsPerElement; ++k)

                    panConnectivity[j * nPointsPerElement + k] =

                        panConnectivity[(j + 1) * nPointsPerElement + k];

            --i;

        }

    }

    if (nOldElements != nElements)

        panConnectivity = (int *)CPLRealloc(

            panConnectivity, sizeof(int) * nElements * nPointsPerElement);

    // Now we update the bounding box if needed

    if (nPoints == 0)

    {

        nMinxIndex = -1;

        nMaxxIndex = -1;

        nMinyIndex = -1;

        nMaxyIndex = -1;

    }

    else

    {

        if (nIndex == nMinxIndex)

        {

            nMinxIndex = 0;

            for (int i = 1; i < nPoints; ++i)

                if (paadfCoords[0][i] < paadfCoords[0][nMinxIndex])

                    nMinxIndex = i;

        }

        if (nIndex == nMaxxIndex)

        {

            nMaxxIndex = 0;

            for (int i = 1; i < nPoints; ++i)

                if (paadfCoords[0][i] > paadfCoords[0][nMaxxIndex])

                    nMaxxIndex = i;

        }

        if (nIndex == nMinyIndex)

        {

            nMinyIndex = 0;

            for (int i = 1; i < nPoints; ++i)

                if (paadfCoords[1][i] < paadfCoords[1][nMinyIndex])

                    nMinyIndex = i;

        }

        if (nIndex == nMaxyIndex)

        {

            nMaxyIndex = 0;

            for (int i = 1; i < nPoints; ++i)

                if (paadfCoords[1][i] > paadfCoords[1][nMaxyIndex])

                    nMaxyIndex = i;

        }

    }

    // We update some parameters of the header

    bTreeUpdateNeeded = true;

    setUpdated();

}

#ifdef notdef

/****************************************************************/

/*                         TimeStep                             */

/****************************************************************/

TimeStep::TimeStep(int nRecordsP, int nFieldsP) : nFields(nFieldsP)

{

    papadfData = (double **)VSI_MALLOC2_VERBOSE(sizeof(double *), nFieldsP);

    for (int i = 0; i < nFieldsP; ++i)

        papadfData[i] =

            (double *)VSI_MALLOC2_VERBOSE(sizeof(double), nRecordsP);

}

TimeStep::~TimeStep()

{

    for (int i = 0; i < nFields; ++i)

        CPLFree(papadfData[i]);

    CPLFree(papadfData);

}

/****************************************************************/

/*                       TimeStepList                           */

/****************************************************************/

TimeStepList::~TimeStepList()

{

    TimeStepList *poFirst = this;

    while (poFirst != 0)

    {

        TimeStepList *poTmp = poFirst->poNext;

        delete poFirst->poStep;

        delete poFirst;

        poFirst = poTmp;

    }

}

#endif

/****************************************************************/

/*                     General functions                        */

/****************************************************************/

int read_integer(VSILFILE *fp, int &nData, bool bDiscard)

{

    unsigned char anb[4];

    if (VSIFReadL(anb, 1, 4, fp) < 4)

    {

        CPLError(CE_Failure, CPLE_FileIO, "%s", SELAFIN_ERROR_MESSAGE);

        return 0;

    };

    if (!bDiscard)

    {

        memcpy(&nData, anb, 4);

        CPL_MSBPTR32(&nData);

    }

    return 1;

}

int write_integer(VSILFILE *fp, int nData)

{

    unsigned char anb[4];

    CPL_MSBPTR32(&nData);

    memcpy(anb, &nData, 4);

    if (VSIFWriteL(anb, 1, 4, fp) < 4)

    {

        CPLError(CE_Failure, CPLE_FileIO, "%s", SELAFIN_ERROR_MESSAGE);

        return 0;

    };

    return 1;

}

int read_string(VSILFILE *fp, char *&pszData, vsi_l_offset nFileSize,

                bool bDiscard)

{

    int nLength = 0;

    read_integer(fp, nLength);

    if (nLength <= 0 || nLength == INT_MAX ||

        static_cast<unsigned>(nLength) > nFileSize)

    {

        CPLError(CE_Failure, CPLE_FileIO, "%s", SELAFIN_ERROR_MESSAGE);

        return 0;

    }

    if (bDiscard)

    {

        if (VSIFSeekL(fp, nLength + 4, SEEK_CUR) != 0)

        {

            CPLError(CE_Failure, CPLE_FileIO, "%s", SELAFIN_ERROR_MESSAGE);

            return 0;

        }

    }

    else

    {

        pszData = (char *)VSI_MALLOC_VERBOSE(nLength + 1);

        if (pszData == nullptr)

        {

            return 0;

        }

        if ((int)VSIFReadL(pszData, 1, nLength, fp) < (int)nLength)

        {

            CPLError(CE_Failure, CPLE_FileIO, "%s", SELAFIN_ERROR_MESSAGE);

            VSIFree(pszData);

            pszData = nullptr;

            return 0;

        }

        pszData[nLength] = 0;

        if (VSIFSeekL(fp, 4, SEEK_CUR) != 0)

        {

            CPLError(CE_Failure, CPLE_FileIO, "%s", SELAFIN_ERROR_MESSAGE);

            VSIFree(pszData);

            pszData = nullptr;

            return 0;

        }

    }

    return nLength;

}

int write_string(VSILFILE *fp, char *pszData, size_t nLength)

{

    if (nLength == 0)

        nLength = strlen(pszData);

    if (write_integer(fp, static_cast<int>(nLength)) == 0)

        return 0;

    if (VSIFWriteL(pszData, 1, nLength, fp) < nLength)

    {

        CPLError(CE_Failure, CPLE_FileIO, "%s", SELAFIN_ERROR_MESSAGE);

        return 0;

    }

    if (write_integer(fp, static_cast<int>(nLength)) == 0)

        return 0;

    return 1;

}

int read_intarray(VSILFILE *fp, int *&panData, vsi_l_offset nFileSize,

                  bool bDiscard)

{

    int nLength = 0;

    read_integer(fp, nLength);

    panData = nullptr;

    if (nLength < 0 || static_cast<unsigned>(nLength) / 4 > nFileSize)

    {

        CPLError(CE_Failure, CPLE_FileIO, "%s", SELAFIN_ERROR_MESSAGE);

        return -1;

    }

    if (bDiscard)

    {

        if (VSIFSeekL(fp, nLength + 4, SEEK_CUR) != 0)

        {

            CPLError(CE_Failure, CPLE_FileIO, "%s", SELAFIN_ERROR_MESSAGE);

            return -1;

        }

    }

    else

    {

        if (nLength == 0)

            panData = nullptr;

        else

        {

            panData = (int *)VSI_MALLOC2_VERBOSE(nLength / 4, sizeof(int));

            if (panData == nullptr)

                return -1;

        }

        for (int i = 0; i < nLength / 4; ++i)

            if (read_integer(fp, panData[i]) == 0)

            {

                CPLFree(panData);

                panData = nullptr;

                CPLError(CE_Failure, CPLE_FileIO, "%s", SELAFIN_ERROR_MESSAGE);

                return -1;

            }

        if (VSIFSeekL(fp, 4, SEEK_CUR) != 0)

        {

            CPLFree(panData);

            panData = nullptr;

            CPLError(CE_Failure, CPLE_FileIO, "%s", SELAFIN_ERROR_MESSAGE);

            return -1;

        }

    }

    return nLength / 4;

}

int write_intarray(VSILFILE *fp, int *panData, size_t nLength)

{

    if (write_integer(fp, static_cast<int>(nLength * 4)) == 0)

        return 0;

    for (size_t i = 0; i < nLength; ++i)

    {

        if (write_integer(fp, panData[i]) == 0)

        {

            CPLError(CE_Failure, CPLE_FileIO, "%s", SELAFIN_ERROR_MESSAGE);

            return 0;

        }

    }

    if (write_integer(fp, static_cast<int>(nLength * 4)) == 0)

        return 0;

    return 1;

}

int read_float(VSILFILE *fp, double &dfData, bool bDiscard)

{

    float dfVal = 0.0;

    if (VSIFReadL(&dfVal, 1, 4, fp) < 4)

    {

        CPLError(CE_Failure, CPLE_FileIO, "%s", SELAFIN_ERROR_MESSAGE);

        return 0;

    };

    if (!bDiscard)

    {

        CPL_MSBPTR32(&dfVal);

        dfData = dfVal;

    }

    return 1;

}

int write_float(VSILFILE *fp, double dfData)

{

    float dfVal = (float)dfData;

    CPL_MSBPTR32(&dfVal);

    if (VSIFWriteL(&dfVal, 1, 4, fp) < 4)

    {

        CPLError(CE_Failure, CPLE_FileIO, "%s", SELAFIN_ERROR_MESSAGE);

        return 0;

    };

    return 1;

}

int read_floatarray(VSILFILE *fp, double **papadfData, vsi_l_offset nFileSize,

                    bool bDiscard)

{

    int nLength = 0;

    read_integer(fp, nLength);

    if (nLength < 0 || static_cast<unsigned>(nLength) / 4 > nFileSize)

    {

        CPLError(CE_Failure, CPLE_FileIO, "%s", SELAFIN_ERROR_MESSAGE);

        return -1;

    }

    if (bDiscard)

    {

        if (VSIFSeekL(fp, nLength + 4, SEEK_CUR) != 0)

        {

            CPLError(CE_Failure, CPLE_FileIO, "%s", SELAFIN_ERROR_MESSAGE);

            return -1;

        }

    }

    else

    {

        if (nLength == 0)

            *papadfData = nullptr;

        else

        {

            *papadfData =

                (double *)VSI_MALLOC2_VERBOSE(sizeof(double), nLength / 4);

            if (*papadfData == nullptr)

                return -1;

        }

        for (int i = 0; i < nLength / 4; ++i)

            if (read_float(fp, (*papadfData)[i]) == 0)

            {

                CPLFree(*papadfData);

                *papadfData = nullptr;

                CPLError(CE_Failure, CPLE_FileIO, "%s", SELAFIN_ERROR_MESSAGE);

                return -1;

            }

        if (VSIFSeekL(fp, 4, SEEK_CUR) != 0)

        {

            CPLFree(*papadfData);

            *papadfData = nullptr;

            CPLError(CE_Failure, CPLE_FileIO, "%s", SELAFIN_ERROR_MESSAGE);

            return -1;

        }

    }

    return nLength / 4;

}

int write_floatarray(VSILFILE *fp, double *papadfData, size_t nLength)

{

    if (write_integer(fp, static_cast<int>(nLength * 4)) == 0)

        return 0;

    for (size_t i = 0; i < nLength; ++i)

    {

        if (write_float(fp, papadfData[i]) == 0)

        {

            CPLError(CE_Failure, CPLE_FileIO, "%s", SELAFIN_ERROR_MESSAGE);

            return 0;

        }

    }

    if (write_integer(fp, static_cast<int>(nLength * 4)) == 0)

        return 0;

    return 1;

}

void Header::UpdateFileSize()

{

    VSIFSeekL(fp, 0, SEEK_END);

    nFileSize = VSIFTellL(fp);

    VSIRewindL(fp);

}

Header *read_header(VSILFILE *fp, const char *pszFilename)

{

    // Save the filename

    Header *poHeader = new Header();

    poHeader->fp = fp;

    poHeader->UpdateFileSize();

    poHeader->pszFilename = CPLStrdup(pszFilename);

    int *panTemp = nullptr;

    // Read the title

    int nLength = read_string(fp, poHeader->pszTitle, poHeader->nFileSize);

    if (nLength == 0)

    {

        delete poHeader;

        return nullptr;

    }

    // Read the array of 2 integers, with the number of variables at the first

    // position

    nLength = read_intarray(fp, panTemp, poHeader->nFileSize);

    if (nLength != 2)

    {

        delete poHeader;

        CPLFree(panTemp);

        return nullptr;

    }

    poHeader->nVar = panTemp[0];

    poHeader->anUnused[0] = panTemp[1];

    CPLFree(panTemp);

    if (poHeader->nVar < 0)

    {

        poHeader->nVar = 0;

        delete poHeader;

        return nullptr;

    }

    if (poHeader->nVar > 1000000 && poHeader->nFileSize / sizeof(int) <

                                        static_cast<unsigned>(poHeader->nVar))

    {

        poHeader->nVar = 0;

        delete poHeader;

        return nullptr;

    }

    // For each variable, read its name as a string of 32 characters

    poHeader->papszVariables =

        (char **)VSI_MALLOC2_VERBOSE(sizeof(char *), poHeader->nVar);

    if (poHeader->nVar > 0 && poHeader->papszVariables == nullptr)

    {

        poHeader->nVar = 0;

        delete poHeader;

        return nullptr;

    }

    for (int i = 0; i < poHeader->nVar; ++i)

    {

        nLength =

            read_string(fp, poHeader->papszVariables[i], poHeader->nFileSize);

        if (nLength == 0)

        {

            poHeader->nVar = i;

            delete poHeader;

            return nullptr;

        }

        // We eliminate quotes in the names of the variables because SQL

        // requests don't seem to appreciate them

        char *pszc = poHeader->papszVariables[i];

        while (*pszc != 0)

        {

            if (*pszc == '\'')

                *pszc = ' ';

            pszc++;

        }

    }

    // Read an array of 10 integers

    nLength = read_intarray(fp, panTemp, poHeader->nFileSize);

    if (nLength < 10)

    {

        delete poHeader;

        CPLFree(panTemp);

        return nullptr;

    }

    poHeader->anUnused[1] = panTemp[0];

    poHeader->nEpsg = panTemp[1];

    poHeader->adfOrigin[0] = panTemp[2];

    poHeader->adfOrigin[1] = panTemp[3];

    for (size_t i = 4; i < 9; ++i)

        poHeader->anUnused[i - 2] = panTemp[i];

    // If the last integer was 1, read an array of 6 integers with the starting

    // date

    if (panTemp[9] == 1)

    {

        nLength =

            read_intarray(fp, poHeader->panStartDate, poHeader->nFileSize);

        if (nLength < 6)

        {

            delete poHeader;

            CPLFree(panTemp);

            return nullptr;

        }

    }

    CPLFree(panTemp);

    // Read an array of 4 integers with the number of elements, points and

    // points per element

    nLength = read_intarray(fp, panTemp, poHeader->nFileSize);

    if (nLength < 4)

    {

        delete poHeader;

        CPLFree(panTemp);

        return nullptr;

    }

    poHeader->nElements = panTemp[0];

    poHeader->nPoints = panTemp[1];

    poHeader->nPointsPerElement = panTemp[2];

    if (poHeader->nElements < 0 || poHeader->nPoints < 0 ||

        poHeader->nPointsPerElement < 0 || panTemp[3] != 1)

    {

        delete poHeader;

        CPLFree(panTemp);

        return nullptr;

    }

    CPLFree(panTemp);

    // Read the connectivity table as an array of nPointsPerElement*nElements

    // integers, and check if all point numbers are valid

    nLength = read_intarray(fp, poHeader->panConnectivity, poHeader->nFileSize);

    if (poHeader->nElements != 0 &&

        nLength / poHeader->nElements != poHeader->nPointsPerElement)

    {

        delete poHeader;

        return nullptr;

    }

    for (int i = 0; i < poHeader->nElements * poHeader->nPointsPerElement; ++i)

    {

        if (poHeader->panConnectivity[i] <= 0 ||

            poHeader->panConnectivity[i] > poHeader->nPoints)

        {

            delete poHeader;

            return nullptr;

        }

    }

    // Read the array of nPoints integers with the border points

    nLength = read_intarray(fp, poHeader->panBorder, poHeader->nFileSize);

    if (nLength != poHeader->nPoints)

    {

        delete poHeader;

        return nullptr;

    }

    // Read two arrays of nPoints floats with the coordinates of each point

    for (size_t i = 0; i < 2; ++i)

    {

        read_floatarray(fp, poHeader->paadfCoords + i, poHeader->nFileSize);

        if (nLength < poHeader->nPoints)

        {

            delete poHeader;

            return nullptr;

        }

        if (poHeader->nPoints != 0 && poHeader->paadfCoords[i] == nullptr)

        {

            delete poHeader;

            return nullptr;

        }

        for (int j = 0; j < poHeader->nPoints; ++j)

            poHeader->paadfCoords[i][j] += poHeader->adfOrigin[i];

    }

    // Update the boundinx box

    poHeader->updateBoundingBox();

    // Update the size of the header and calculate the number of time steps

    poHeader->setUpdated();

    int nPos = poHeader->getPosition(0);

    if (static_cast<vsi_l_offset>(nPos) > poHeader->nFileSize)

    {

        delete poHeader;

        return nullptr;

    }

    vsi_l_offset nStepsBig =

        poHeader->nVar != 0

            ? (poHeader->nFileSize - nPos) / (poHeader->getPosition(1) - nPos)

            : 0;

    if (nStepsBig > INT_MAX)

        poHeader->nSteps = INT_MAX;

    else

        poHeader->nSteps = static_cast<int>(nStepsBig);

    return poHeader;

}

int write_header(VSILFILE *fp, Header *poHeader)

{

    VSIRewindL(fp);

    if (write_string(fp, poHeader->pszTitle, 80) == 0)

        return 0;

    int anTemp[10] = {0};

    anTemp[0] = poHeader->nVar;

    anTemp[1] = poHeader->anUnused[0];

    if (write_intarray(fp, anTemp, 2) == 0)

        return 0;

    for (int i = 0; i < poHeader->nVar; ++i)

        if (write_string(fp, poHeader->papszVariables[i], 32) == 0)

            return 0;

    anTemp[0] = poHeader->anUnused[1];

    anTemp[1] = poHeader->nEpsg;

    anTemp[2] = (int)poHeader->adfOrigin[0];

    anTemp[3] = (int)poHeader->adfOrigin[1];

    for (size_t i = 4; i < 9; ++i)

        anTemp[i] = poHeader->anUnused[i - 2];

    anTemp[9] = (poHeader->panStartDate != nullptr) ? 1 : 0;

    if (write_intarray(fp, anTemp, 10) == 0)

        return 0;

    if (poHeader->panStartDate != nullptr &&

        write_intarray(fp, poHeader->panStartDate, 6) == 0)

        return 0;

    anTemp[0] = poHeader->nElements;

    anTemp[1] = poHeader->nPoints;

    anTemp[2] = poHeader->nPointsPerElement;

    anTemp[3] = 1;

    if (write_intarray(fp, anTemp, 4) == 0)

        return 0;

    if (write_intarray(fp, poHeader->panConnectivity,

                       static_cast<size_t>(poHeader->nElements) *

                           poHeader->nPointsPerElement) == 0)

        return 0;

    if (write_intarray(fp, poHeader->panBorder, poHeader->nPoints) == 0)

        return 0;

    double *dfVals =

        (double *)VSI_MALLOC2_VERBOSE(sizeof(double), poHeader->nPoints);

    if (poHeader->nPoints > 0 && dfVals == nullptr)

        return 0;

    for (size_t i = 0; i < 2; ++i)

    {

        for (int j = 0; j < poHeader->nPoints; ++j)

            dfVals[j] = poHeader->paadfCoords[i][j] - poHeader->adfOrigin[i];

        if (write_floatarray(fp, dfVals, poHeader->nPoints) == 0)

        {

            CPLFree(dfVals);

            return 0;

        }

    }

    CPLFree(dfVals);

    return 1;

}

#ifdef notdef

int read_step(VSILFILE *fp, const Header *poHeader, TimeStep *&poStep)

{

    poStep = new TimeStep(poHeader->nPoints, poHeader->nVar);

    int nLength = 0;

    if (read_integer(fp, nLength) == 0 || nLength != 1)

    {

        delete poStep;

        return 0;

    }

    if (read_float(fp, poStep->dfDate) == 0)

    {

        delete poStep;

        return 0;

    }

    if (read_integer(fp, nLength) == 0 || nLength != 1)

    {

        delete poStep;

        return 0;

    }

    for (int i = 0; i < poHeader->nVar; ++i)

    {

        nLength = read_floatarray(fp, &(poStep->papadfData[i]));

        if (nLength != poHeader->nPoints)

        {

            delete poStep;

            return 0;

        }

    }

    return 1;

}

int write_step(VSILFILE *fp, const Header *poHeader, const TimeStep *poStep)

{

    if (write_integer(fp, 1) == 0)

        return 0;

    if (write_float(fp, poStep->dfDate) == 0)

        return 0;

    if (write_integer(fp, 1) == 0)

        return 0;

    for (int i = 0; i < poHeader->nVar; ++i)

    {

        if (write_floatarray(fp, poStep->papadfData[i]) == 0)

            return 0;

    }

    return 1;

}

int read_steps(VSILFILE *fp, const Header *poHeader, TimeStepList *&poSteps)

{

    poSteps = 0;

    TimeStepList *poCur, *poNew;

    for (int i = 0; i < poHeader->nSteps; ++i)

    {

        poNew = new TimeStepList(0, 0);

        if (read_step(fp, poHeader, poNew->poStep) == 0)

        {

            delete poSteps;

            return 0;

        }

        if (poSteps == 0)

            poSteps = poNew;

        else

            poCur->poNext = poNew;

        poCur = poNew;

    }

    return 1;

}

int write_steps(VSILFILE *fp, const Header *poHeader,

                const TimeStepList *poSteps)

{

    const TimeStepList *poCur = poSteps;

    while (poCur != 0)

    {

        if (write_step(fp, poHeader, poCur->poStep) == 0)

            return 0;

        poCur = poCur->poNext;

    }

    return 1;

}

#endif

}  // namespace Selafin