/******************************************************************************
 * $Id: dted_ptstream.c 3b0bbf7a8a012d69a783ee1f9cfeb5c52b370021 2017-06-27 20:57:02Z Even Rouault $
 *
 * Project:  DTED Translator
 * Purpose:  DTED Point Stream Writer.
 * Author:   Frank Warmerdam, warmerdam@pobox.com
 *
 ******************************************************************************
 * Copyright (c) 2001, Frank Warmerdam
 *
 * 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 "dted_api.h"

CPL_CVSID("$Id: dted_ptstream.c 3b0bbf7a8a012d69a783ee1f9cfeb5c52b370021 2017-06-27 20:57:02Z Even Rouault $")

typedef struct {
    char     *pszFilename;
    DTEDInfo *psInfo;

    GInt16 **papanProfiles;

    int    nLLLong;
    int    nLLLat;
} DTEDCachedFile;

typedef struct {
    int nLevel;
    char *pszPath;

    double  dfPixelSize;

    int nOpenFiles;
    DTEDCachedFile *pasCF;

    int nLastFile;

    char *apszMetadata[DTEDMD_MAX+1];
} DTEDPtStream;

/************************************************************************/
/*                         DTEDCreatePtStream()                         */
/************************************************************************/

void *DTEDCreatePtStream( const char *pszPath, int nLevel )

{
    DTEDPtStream *psStream;
    int          i;
    VSIStatBuf   sStat;

/* -------------------------------------------------------------------- */
/*      Does the target directory already exist?  If not try to         */
/*      create it.                                                      */
/* -------------------------------------------------------------------- */
    if( CPLStat( pszPath, &sStat ) != 0 )
    {
        if( VSIMkdir( pszPath, 0755 ) != 0 )
        {
#ifndef AVOID_CPL
            CPLError( CE_Failure, CPLE_OpenFailed,
                      "Unable to find, or create directory `%s'.",
                      pszPath );
#endif
            return NULL;
        }
    }

/* -------------------------------------------------------------------- */
/*      Create the stream and initialize it.                            */
/* -------------------------------------------------------------------- */

    psStream = (DTEDPtStream *) CPLCalloc( sizeof(DTEDPtStream), 1 );
    psStream->nLevel = nLevel;
    psStream->pszPath = CPLStrdup( pszPath );
    psStream->nOpenFiles = 0;
    psStream->pasCF = NULL;
    psStream->nLastFile = -1;

    for( i = 0; i < DTEDMD_MAX+1; i++ )
        psStream->apszMetadata[i] = NULL;

    if( nLevel == 0 )
        psStream->dfPixelSize = 1.0 / 120.0;
    else if( nLevel == 1 )
        psStream->dfPixelSize = 1.0 / 1200.0;
    else /* if( nLevel == 2 ) */
        psStream->dfPixelSize = 1.0 / 3600.0;

    return (void *) psStream;
}

/************************************************************************/
/*                        DTEDPtStreamNewTile()                         */
/*                                                                      */
/*      Create a new DTED file file, add it to our list, and make it    */
/*      "current".                                                      */
/************************************************************************/

static int DTEDPtStreamNewTile( DTEDPtStream *psStream,
                                int nCrLong, int nCrLat )

{
    DTEDInfo        *psInfo;
    char            szFile[128];
    char            chNSHemi, chEWHemi;
    char            *pszFullFilename;
    const char      *pszError;

    /* work out filename */
    if( nCrLat < 0 )
        chNSHemi = 's';
    else
        chNSHemi = 'n';

    if( nCrLong < 0 )
        chEWHemi = 'w';
    else
        chEWHemi = 'e';

    snprintf( szFile, sizeof(szFile), "%c%03d%c%03d.dt%d",
             chEWHemi, ABS(nCrLong), chNSHemi, ABS(nCrLat),
             psStream->nLevel );

    pszFullFilename =
        CPLStrdup(CPLFormFilename( psStream->pszPath, szFile, NULL ));

    /* create the dted file */
    pszError = DTEDCreate( pszFullFilename, psStream->nLevel,
                           nCrLat, nCrLong );
    if( pszError != NULL )
    {
#ifndef AVOID_CPL
        CPLError( CE_Failure, CPLE_OpenFailed,
                  "Failed to create DTED file `%s'.\n%s",
                  pszFullFilename, pszError );
#endif
        return FALSE;
    }

    psInfo = DTEDOpen( pszFullFilename, "rb+", FALSE );

    if( psInfo == NULL )
    {
        CPLFree( pszFullFilename );
        return FALSE;
    }

    /* add cached file to stream */
    psStream->nOpenFiles++;
    psStream->pasCF =
        CPLRealloc(psStream->pasCF,
                   sizeof(DTEDCachedFile)*psStream->nOpenFiles);

    psStream->pasCF[psStream->nOpenFiles-1].psInfo = psInfo;
    psStream->pasCF[psStream->nOpenFiles-1].papanProfiles =
        CPLCalloc(sizeof(GInt16*),psInfo->nXSize);
    psStream->pasCF[psStream->nOpenFiles-1].pszFilename = pszFullFilename;
    psStream->pasCF[psStream->nOpenFiles-1].nLLLat = nCrLat;
    psStream->pasCF[psStream->nOpenFiles-1].nLLLong = nCrLong;

    psStream->nLastFile = psStream->nOpenFiles-1;

    return TRUE;
}

/************************************************************************/
/*                           DTEDWritePtLL()                            */
/************************************************************************/

static int DTEDWritePtLL( CPL_UNUSED DTEDPtStream *psStream,
                          DTEDCachedFile *psCF,
                          double dfLong,
                          double dfLat,
                          double dfElev )
{
/* -------------------------------------------------------------------- */
/*      Determine what profile this belongs in, and initialize the      */
/*      profile if it doesn't already exist.                            */
/* -------------------------------------------------------------------- */
    DTEDInfo *psInfo = psCF->psInfo;
    int iProfile, i, iRow;

    iProfile = (int) ((dfLong - psInfo->dfULCornerX) / psInfo->dfPixelSizeX);
    iProfile = MAX(0,MIN(psInfo->nXSize-1,iProfile));

    if( psCF->papanProfiles[iProfile] == NULL )
    {
        psCF->papanProfiles[iProfile] =
            CPLMalloc(sizeof(GInt16) * psInfo->nYSize);

        for( i = 0; i < psInfo->nYSize; i++ )
            psCF->papanProfiles[iProfile][i] = DTED_NODATA_VALUE;
    }

/* -------------------------------------------------------------------- */
/*      Establish where we fit in the profile.                          */
/* -------------------------------------------------------------------- */
    iRow = (int) ((psInfo->dfULCornerY-dfLat) / psInfo->dfPixelSizeY);
    iRow = MAX(0,MIN(psInfo->nYSize-1,iRow));

    psCF->papanProfiles[iProfile][iRow] = (GInt16) floor(dfElev+0.5);

    return TRUE;
}

/************************************************************************/
/*                            DTEDWritePt()                             */
/*                                                                      */
/*      Write a single point out, creating a new file if necessary      */
/*      to hold it.                                                     */
/************************************************************************/

int DTEDWritePt( void *hStream, double dfLong, double dfLat, double dfElev )

{
    DTEDPtStream *psStream = (DTEDPtStream *) hStream;
    int          i;
    DTEDInfo     *psInfo;
    int          bOnBoundary = FALSE;

/* -------------------------------------------------------------------- */
/*      Determine if we are in a boundary region ... that is in the     */
/*      area of the edge "pixel" that is shared with adjacent           */
/*      tiles.                                                          */
/* -------------------------------------------------------------------- */
    if( (floor(dfLong - 0.5*psStream->dfPixelSize)
         != floor(dfLong + 0.5*psStream->dfPixelSize))
        || (floor(dfLat - 0.5*psStream->dfPixelSize)
            != floor(dfLat + 0.5*psStream->dfPixelSize)) )
    {
        bOnBoundary = TRUE;
        psStream->nLastFile = -1;
    }

/* ==================================================================== */
/*      Handle case where the tile is not on a boundary.  We only       */
/*      need one output tile.                                           */
/* ==================================================================== */
/* -------------------------------------------------------------------- */
/*      Is the last file used still applicable?                         */
/* -------------------------------------------------------------------- */
    if( !bOnBoundary )
    {
        if( psStream->nLastFile != -1 )
        {
            psInfo = psStream->pasCF[psStream->nLastFile].psInfo;

            if( dfLat > psInfo->dfULCornerY
                || dfLat < psInfo->dfULCornerY - 1.0 - psInfo->dfPixelSizeY
                || dfLong < psInfo->dfULCornerX
                || dfLong > psInfo->dfULCornerX + 1.0 + psInfo->dfPixelSizeX )
                psStream->nLastFile = -1;
        }

/* -------------------------------------------------------------------- */
/*      Search for the file to write to.                                */
/* -------------------------------------------------------------------- */
        for( i = 0; i < psStream->nOpenFiles && psStream->nLastFile == -1; i++ )
        {
            psInfo = psStream->pasCF[i].psInfo;

            if( !(dfLat > psInfo->dfULCornerY
                  || dfLat < psInfo->dfULCornerY - 1.0 - psInfo->dfPixelSizeY
                  || dfLong < psInfo->dfULCornerX
                  || dfLong > psInfo->dfULCornerX + 1.0 + psInfo->dfPixelSizeX) )
            {
                psStream->nLastFile = i;
            }
        }

/* -------------------------------------------------------------------- */
/*      If none found, create a new file.                               */
/* -------------------------------------------------------------------- */
        if( psStream->nLastFile == -1 )
        {
            int nCrLong, nCrLat;

            nCrLong = (int) floor(dfLong);
            nCrLat = (int) floor(dfLat);

            if( !DTEDPtStreamNewTile( psStream, nCrLong, nCrLat ) )
                return FALSE;
        }

/* -------------------------------------------------------------------- */
/*      Write data out to selected tile.                                */
/* -------------------------------------------------------------------- */
        return DTEDWritePtLL( psStream, psStream->pasCF + psStream->nLastFile,
                              dfLong, dfLat, dfElev );
    }

/* ==================================================================== */
/*      Handle case where we are on a boundary.  We may be writing      */
/*      the value to as many as four tiles.                             */
/* ==================================================================== */
    else
    {
        int nLatMin, nLatMax, nLongMin, nLongMax;
        int nCrLong, nCrLat;

        nLongMin = (int) floor( dfLong - 0.5*psStream->dfPixelSize );
        nLongMax = (int) floor( dfLong + 0.5*psStream->dfPixelSize );
        nLatMin = (int) floor( dfLat - 0.5*psStream->dfPixelSize );
        nLatMax = (int) floor( dfLat + 0.5*psStream->dfPixelSize );

        for( nCrLong = nLongMin; nCrLong <= nLongMax; nCrLong++ )
        {
            for( nCrLat = nLatMin; nCrLat <= nLatMax; nCrLat++ )
            {
                psStream->nLastFile = -1;

/* -------------------------------------------------------------------- */
/*      Find this tile in our existing list.                            */
/* -------------------------------------------------------------------- */
                for( i = 0; i < psStream->nOpenFiles; i++ )
                {
                    if( psStream->pasCF[i].nLLLong == nCrLong
                        && psStream->pasCF[i].nLLLat == nCrLat )
                    {
                        psStream->nLastFile = i;
                        break;
                    }
                }

/* -------------------------------------------------------------------- */
/*      Create the tile if not found.                                   */
/* -------------------------------------------------------------------- */
                if( psStream->nLastFile == -1 )
                {
                    if( !DTEDPtStreamNewTile( psStream, nCrLong, nCrLat ) )
                        return FALSE;
                }

/* -------------------------------------------------------------------- */
/*      Write to the tile.                                              */
/* -------------------------------------------------------------------- */
                if( !DTEDWritePtLL( psStream,
                                    psStream->pasCF + psStream->nLastFile,
                                    dfLong, dfLat, dfElev ) )
                    return FALSE;
            }
        }
    }

    return TRUE;
}

/************************************************************************/
/*                         DTEDClosePtStream()                          */
/************************************************************************/

void DTEDClosePtStream( void *hStream )

{
    DTEDPtStream *psStream = (DTEDPtStream *) hStream;
    int           iFile, iMD;

/* -------------------------------------------------------------------- */
/*      Flush all DTED files.                                           */
/* -------------------------------------------------------------------- */
    for( iFile = 0; iFile < psStream->nOpenFiles; iFile++ )
    {
        int             iProfile;
        DTEDCachedFile  *psCF = psStream->pasCF + iFile;

        for( iProfile = 0; iProfile < psCF->psInfo->nXSize; iProfile++ )
        {
            if( psCF->papanProfiles[iProfile] != NULL )
            {
                DTEDWriteProfile( psCF->psInfo, iProfile,
                                  psCF->papanProfiles[iProfile] );
                CPLFree( psCF->papanProfiles[iProfile] );
            }
        }

        CPLFree( psCF->papanProfiles );

        for( iMD = 0; iMD <= DTEDMD_MAX; iMD++ )
        {
            if( psStream->apszMetadata[iMD] != NULL )
                DTEDSetMetadata( psCF->psInfo, (DTEDMetaDataCode)iMD,
                                 psStream->apszMetadata[iMD] );
        }

        DTEDClose( psCF->psInfo );
    }

/* -------------------------------------------------------------------- */
/*      Final cleanup.                                                  */
/* -------------------------------------------------------------------- */

    for( iMD = 0; iMD < DTEDMD_MAX+1; iMD++ )
        CPLFree( psStream->apszMetadata[iMD] );

    CPLFree( psStream->pasCF );
    CPLFree( psStream->pszPath );
    CPLFree( psStream );
}

/************************************************************************/
/*                           DTEDFillPixel()                            */
/************************************************************************/
static
void DTEDFillPixel( DTEDInfo *psInfo, GInt16 **papanProfiles,
                    GInt16 **papanDstProfiles, int iX, int iY,
                    int nPixelSearchDist, float *pafKernel )

{
    int nKernelWidth = 2 * nPixelSearchDist + 1;
    int nXMin, nXMax, nYMin, nYMax;
    double dfCoefSum = 0.0, dfValueSum = 0.0;
    int iXS, iYS;

    nXMin = MAX(0,iX - nPixelSearchDist);
    nXMax = MIN(psInfo->nXSize-1,iX + nPixelSearchDist);
    nYMin = MAX(0,iY - nPixelSearchDist);
    nYMax = MIN(psInfo->nYSize-1,iY + nPixelSearchDist);

    for( iXS = nXMin; iXS <= nXMax; iXS++ )
    {
        GInt16  *panThisProfile = papanProfiles[iXS];

        if( panThisProfile == NULL )
            continue;

        for( iYS = nYMin; iYS <= nYMax; iYS++ )
        {
            if( panThisProfile[iYS] != DTED_NODATA_VALUE )
            {
                int     iXK, iYK;
                float   fKernelCoef;

                iXK = iXS - iX + nPixelSearchDist;
                iYK = iYS - iY + nPixelSearchDist;

                fKernelCoef = pafKernel[iXK + iYK * nKernelWidth];
                dfCoefSum += fKernelCoef;
                dfValueSum += fKernelCoef * panThisProfile[iYS];
            }
        }
    }

    if( dfCoefSum == 0.0 )
        papanDstProfiles[iX][iY] = DTED_NODATA_VALUE;
    else
        papanDstProfiles[iX][iY] =
            (GInt16) floor(dfValueSum / dfCoefSum + 0.5);
}

/************************************************************************/
/*                          DTEDFillPtStream()                          */
/*                                                                      */
/*      Apply simple inverse distance interpolator to all no-data       */
/*      pixels based on available values within the indicated search    */
/*      distance (rectangular).                                         */
/************************************************************************/

void DTEDFillPtStream( void *hStream, int nPixelSearchDist )

{
    DTEDPtStream *psStream = (DTEDPtStream *) hStream;
    int           iFile, nKernelWidth;
    float         *pafKernel;
    int           iX, iY;

/* -------------------------------------------------------------------- */
/*      Setup inverse distance weighting kernel.                        */
/* -------------------------------------------------------------------- */
    nKernelWidth = 2 * nPixelSearchDist + 1;
    pafKernel = (float *) CPLMalloc(nKernelWidth*nKernelWidth*sizeof(float));

    for( iX = 0; iX < nKernelWidth; iX++ )
    {
        for( iY = 0; iY < nKernelWidth; iY++ )
        {
            pafKernel[iX + iY * nKernelWidth] = (float) (1.0 /
                sqrt( (nPixelSearchDist-iX) * (nPixelSearchDist-iX)
                      + (nPixelSearchDist-iY) * (nPixelSearchDist-iY) ));
        }
    }

/* ==================================================================== */
/*      Process each cached file.                                       */
/* ==================================================================== */
    for( iFile = 0; iFile < psStream->nOpenFiles; iFile++ )
    {
        DTEDInfo        *psInfo = psStream->pasCF[iFile].psInfo;
        GInt16          **papanProfiles = psStream->pasCF[iFile].papanProfiles;
        GInt16          **papanDstProfiles;

        papanDstProfiles = (GInt16 **)
            CPLCalloc(sizeof(GInt16*),psInfo->nXSize);

/* -------------------------------------------------------------------- */
/*      Setup output image.                                             */
/* -------------------------------------------------------------------- */
        for( iX = 0; iX < psInfo->nXSize; iX++ )
        {
            papanDstProfiles[iX] = (GInt16 *)
                CPLMalloc(sizeof(GInt16) * psInfo->nYSize);
        }

/* -------------------------------------------------------------------- */
/*      Interpolate all missing values, and copy over available values. */
/* -------------------------------------------------------------------- */
        for( iX = 0; iX < psInfo->nXSize; iX++ )
        {
            for( iY = 0; iY < psInfo->nYSize; iY++ )
            {
                if( papanProfiles[iX] == NULL
                    || papanProfiles[iX][iY] == DTED_NODATA_VALUE )
                {
                    DTEDFillPixel( psInfo, papanProfiles, papanDstProfiles,
                                   iX, iY, nPixelSearchDist, pafKernel );
                }
                else
                {
                    papanDstProfiles[iX][iY] = papanProfiles[iX][iY];
                }
            }
        }
/* -------------------------------------------------------------------- */
/*      Push new values back into cache.                                */
/* -------------------------------------------------------------------- */
        for( iX = 0; iX < psInfo->nXSize; iX++ )
        {
            CPLFree( papanProfiles[iX] );
            papanProfiles[iX] = papanDstProfiles[iX];
        }

        CPLFree( papanDstProfiles );
    }

    CPLFree( pafKernel );
}

/************************************************************************/
/*                      DTEDPtStreamSetMetadata()                       */
/************************************************************************/

void DTEDPtStreamSetMetadata( void *hStream, DTEDMetaDataCode eCode,
                              const char *pszValue )

{
    DTEDPtStream *psStream = (DTEDPtStream *) hStream;

    if( (int)eCode >= 0 && eCode < DTEDMD_MAX+1 )
    {
        CPLFree( psStream->apszMetadata[eCode] );
        psStream->apszMetadata[eCode] = CPLStrdup( pszValue );
    }
}

/************************************************************************/
/*                   DTEDPtStreamTrimEdgeOnlyTiles()                    */
/*                                                                      */
/*      Erase all tiles that only have boundary values set.             */
/************************************************************************/

void DTEDPtStreamTrimEdgeOnlyTiles( void *hStream )

{
    DTEDPtStream *psStream = (DTEDPtStream *) hStream;
    int iFile;

    for( iFile = psStream->nOpenFiles-1; iFile >= 0; iFile-- )
    {
        DTEDInfo        *psInfo = psStream->pasCF[iFile].psInfo;
        GInt16          **papanProfiles = psStream->pasCF[iFile].papanProfiles;
        int             iProfile, iPixel, bGotNonEdgeData = FALSE;

        for( iProfile = 1; iProfile < psInfo->nXSize-1; iProfile++ )
        {
            if( papanProfiles[iProfile] == NULL )
                continue;

            for( iPixel = 1; iPixel < psInfo->nYSize-1; iPixel++ )
            {
                if( papanProfiles[iProfile][iPixel] != DTED_NODATA_VALUE )
                {
                    bGotNonEdgeData = TRUE;
                    break;
                }
            }
        }

        if( bGotNonEdgeData )
            continue;

        /* Remove this tile */

        for( iProfile = 0; iProfile < psInfo->nXSize; iProfile++ )
        {
            if( papanProfiles[iProfile] != NULL )
                CPLFree( papanProfiles[iProfile] );
        }
        CPLFree( papanProfiles );

        DTEDClose( psInfo );

        VSIUnlink( psStream->pasCF[iFile].pszFilename );
        CPLFree( psStream->pasCF[iFile].pszFilename );

        memmove( psStream->pasCF + iFile,
                 psStream->pasCF + iFile + 1,
                 sizeof(DTEDCachedFile) * (psStream->nOpenFiles-iFile-1) );
        psStream->nOpenFiles--;
    }
}