dted_ptstream.c

/******************************************************************************
 * $Id$
 *
 * 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"

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--;
    }
}