gdalsimplewarp.cpp

/******************************************************************************
 *
 * Project:  Mapinfo Image Warper
 * Purpose:  Simple (source in memory) warp algorithm.
 * Author:   Frank Warmerdam, warmerdam@pobox.com
 *
 ******************************************************************************
 * Copyright (c) 2002, i3 - information integration and imaging, Fort Collin, CO
 *
 * 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 "cpl_port.h"
#include "gdal_alg.h"

#include <cstdlib>
#include <cstring>

#include <algorithm>

#include "cpl_conv.h"
#include "cpl_error.h"
#include "cpl_progress.h"
#include "cpl_string.h"
#include "cpl_vsi.h"
#include "gdal.h"
#include "gdal_priv.h"

CPL_CVSID("$Id$")

/************************************************************************/
/*                      GDALSimpleWarpRemapping()                       */
/*                                                                      */
/*      This function implements any raster remapping requested in      */
/*      the options list.  The remappings are applied to the source     */
/*      data before warping.  Two kinds are support ... REMAP           */
/*      commands which remap selected pixel values for any band and     */
/*      REMAP_MULTI which only remap pixels matching the input in       */
/*      all bands at once (i.e. to remap an RGB value to another).      */
/************************************************************************/

static void
GDALSimpleWarpRemapping( int nBandCount, GByte **papabySrcData,
                         int nSrcXSize, int nSrcYSize,
                         char **papszWarpOptions )

{

/* ==================================================================== */
/*      Process any and all single value REMAP commands.                */
/* ==================================================================== */
    char **papszRemaps = CSLFetchNameValueMultiple( papszWarpOptions,
                                                    "REMAP" );

    const int nRemaps = CSLCount(papszRemaps);
    for( int iRemap = 0; iRemap < nRemaps; iRemap++ )
    {

/* -------------------------------------------------------------------- */
/*      What are the pixel values to map from and to?                   */
/* -------------------------------------------------------------------- */
        char **papszTokens = CSLTokenizeString( papszRemaps[iRemap] );

        if( CSLCount(papszTokens) != 2 )
        {
            CPLError(CE_Warning, CPLE_AppDefined,
                     "Ill formed REMAP `%s' ignored in "
                     "GDALSimpleWarpRemapping()",
                     papszRemaps[iRemap] );
            CSLDestroy( papszTokens );
            continue;
        }

        const int nFromValue = atoi(papszTokens[0]);
        const int nToValue = atoi(papszTokens[1]);
        // TODO(schwehr): Why is it ok to narrow ints to byte without checking?
        const GByte byToValue = static_cast<GByte>(nToValue);

        CSLDestroy( papszTokens );

/* -------------------------------------------------------------------- */
/*      Pass over each band searching for matches.                      */
/* -------------------------------------------------------------------- */
        for( int iBand = 0; iBand < nBandCount; iBand++ )
        {
            GByte *pabyData = papabySrcData[iBand];
            int nPixelCount = nSrcXSize * nSrcYSize;

            while( nPixelCount != 0 )
            {
                if( *pabyData == nFromValue )
                    *pabyData = byToValue;

                pabyData++;
                nPixelCount--;
            }
        }
    }

    CSLDestroy( papszRemaps );

/* ==================================================================== */
/*      Process any and all REMAP_MULTI commands.                       */
/* ==================================================================== */
    papszRemaps = CSLFetchNameValueMultiple( papszWarpOptions,
                                             "REMAP_MULTI" );

    const int nRemapsMulti = CSLCount(papszRemaps);
    for( int iRemap = 0; iRemap < nRemapsMulti; iRemap++ )
    {
/* -------------------------------------------------------------------- */
/*      What are the pixel values to map from and to?                   */
/* -------------------------------------------------------------------- */
        char **papszTokens = CSLTokenizeString( papszRemaps[iRemap] );

        const int nTokens = CSLCount(papszTokens);
        if( nTokens % 2 == 1 ||
            nTokens == 0 ||
            nTokens > nBandCount * 2 )
        {
            CPLError(CE_Warning, CPLE_AppDefined,
                     "Ill formed REMAP_MULTI `%s' ignored in "
                     "GDALSimpleWarpRemapping()",
                     papszRemaps[iRemap]);
            CSLDestroy( papszTokens );
            continue;
        }

        const int nMapBandCount = nTokens / 2;

        int *panFromValue = static_cast<int *>(
            CPLMalloc(sizeof(int) * nMapBandCount ) );
        int *panToValue = static_cast<int *>(
            CPLMalloc(sizeof(int) * nMapBandCount ) );

        for( int iBand = 0; iBand < nMapBandCount; iBand++ )
        {
            panFromValue[iBand] = atoi(papszTokens[iBand]);
            panToValue[iBand] = atoi(papszTokens[iBand+nMapBandCount]);
        }

        CSLDestroy( papszTokens );

/* -------------------------------------------------------------------- */
/*      Search for matching values to replace.                          */
/* -------------------------------------------------------------------- */
        const int nPixelCount = nSrcXSize * nSrcYSize;

        for( int iPixel = 0; iPixel < nPixelCount; iPixel++ )
        {
            bool bMatch = true;

            // Always check band 0.
            for( int iBand = 0; bMatch && iBand < std::max(1, nMapBandCount);
                 iBand++ )
            {
                if( papabySrcData[iBand][iPixel] != panFromValue[iBand] )
                    bMatch = false;
            }

            if( !bMatch )
                continue;

            for( int iBand = 0; iBand < nMapBandCount; iBand++ )
                papabySrcData[iBand][iPixel] =
                    static_cast<GByte>( panToValue[iBand] );
        }

        CPLFree( panFromValue );
        CPLFree( panToValue );
    }

    CSLDestroy( papszRemaps );
}

/************************************************************************/
/*                        GDALSimpleImageWarp()                         */
/************************************************************************/

/**
 * Perform simple image warp.
 *
 * Copies an image from a source dataset to a destination dataset applying
 * an application defined transformation.   This algorithm is called simple
 * because it lacks many options such as resampling kernels (other than
 * nearest neighbour), support for data types other than 8bit, and the
 * ability to warp images without holding the entire source and destination
 * image in memory.
 *
 * The following option(s) may be passed in papszWarpOptions.
 * <ul>
 * <li> "INIT=v[,v...]": This option indicates that the output dataset should
 * be initialized to the indicated value in any area valid data is not written.
 * Distinct values may be listed for each band separated by columns.
 * </ul>
 *
 * @param hSrcDS the source image dataset.
 * @param hDstDS the destination image dataset.
 * @param nBandCount the number of bands to be warped.  If zero, all bands
 * will be processed.
 * @param panBandList the list of bands to translate.
 * @param pfnTransform the transformation function to call.  See
 * GDALTransformerFunc().
 * @param pTransformArg the callback handle to pass to pfnTransform.
 * @param pfnProgress the function used to report progress.  See
 * GDALProgressFunc().
 * @param pProgressArg the callback handle to pass to pfnProgress.
 * @param papszWarpOptions additional options controlling the warp.
 *
 * @return TRUE if the operation completes, or FALSE if an error occurs.
 */

int CPL_STDCALL
GDALSimpleImageWarp( GDALDatasetH hSrcDS, GDALDatasetH hDstDS,
                     int nBandCount, int *panBandList,
                     GDALTransformerFunc pfnTransform, void *pTransformArg,
                     GDALProgressFunc pfnProgress, void *pProgressArg,
                     char **papszWarpOptions )

{
    VALIDATE_POINTER1( hSrcDS, "GDALSimpleImageWarp", 0 );
    VALIDATE_POINTER1( hDstDS, "GDALSimpleImageWarp", 0 );

    bool bError = false;

/* -------------------------------------------------------------------- */
/*      If no bands provided assume we should process all bands.        */
/* -------------------------------------------------------------------- */
    if( nBandCount == 0 )
    {
        nBandCount = GDALGetRasterCount( hSrcDS );
        if( nBandCount == 0 )
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                     "No raster band in source dataset");
            return FALSE;
        }

        panBandList = static_cast<int *>(CPLCalloc(sizeof(int), nBandCount));

        for( int iBand = 0; iBand < nBandCount; iBand++ )
            panBandList[iBand] = iBand + 1;

        const int nResult =
            GDALSimpleImageWarp( hSrcDS, hDstDS, nBandCount, panBandList,
                                 pfnTransform, pTransformArg,
                                 pfnProgress, pProgressArg,
                                 papszWarpOptions );
        CPLFree( panBandList );
        return nResult;
    }

/* -------------------------------------------------------------------- */
/*      Post initial progress.                                          */
/* -------------------------------------------------------------------- */
    if( pfnProgress )
    {
        if( !pfnProgress( 0.0, "", pProgressArg ) )
            return FALSE;
    }

/* -------------------------------------------------------------------- */
/*      Load the source image band(s).                                  */
/* -------------------------------------------------------------------- */
    const int nSrcXSize = GDALGetRasterXSize(hSrcDS);
    const int nSrcYSize = GDALGetRasterYSize(hSrcDS);
    GByte **papabySrcData = static_cast<GByte **>(
        CPLCalloc(nBandCount, sizeof(GByte*)) );

    bool ok = true;
    for( int iBand = 0; iBand < nBandCount; iBand++ )
    {
        papabySrcData[iBand] = static_cast<GByte *>(
            VSI_MALLOC2_VERBOSE(nSrcXSize, nSrcYSize) );
        if( papabySrcData[iBand] == nullptr )
        {
            CPLError( CE_Failure, CPLE_OutOfMemory,
                      "GDALSimpleImageWarp out of memory." );
            ok = false;
            break;
        }

        if( GDALRasterIO(
                GDALGetRasterBand(hSrcDS,panBandList[iBand]), GF_Read,
                0, 0, nSrcXSize, nSrcYSize,
                papabySrcData[iBand], nSrcXSize, nSrcYSize, GDT_Byte,
                0, 0 ) != CE_None )
        {
            CPLError( CE_Failure, CPLE_FileIO,
                      "GDALSimpleImageWarp GDALRasterIO failure %s",
                      CPLGetLastErrorMsg() );
            ok = false;
            break;
        }
    }
    if( !ok )
    {
        for( int i=0; i <= nBandCount; i++ )
        {
            VSIFree(papabySrcData[i]);
        }
        CPLFree(papabySrcData);
        return FALSE;
    }

/* -------------------------------------------------------------------- */
/*      Check for remap request(s).                                     */
/* -------------------------------------------------------------------- */
    GDALSimpleWarpRemapping( nBandCount, papabySrcData, nSrcXSize, nSrcYSize,
                             papszWarpOptions );

/* -------------------------------------------------------------------- */
/*      Allocate scanline buffers for output image.                     */
/* -------------------------------------------------------------------- */
    const int nDstXSize = GDALGetRasterXSize( hDstDS );
    const int nDstYSize = GDALGetRasterYSize( hDstDS );
    GByte **papabyDstLine = static_cast<GByte **>(
        CPLCalloc(nBandCount, sizeof(GByte*)) );

    for( int iBand = 0; iBand < nBandCount; iBand++ )
        papabyDstLine[iBand] = static_cast<GByte *>(CPLMalloc( nDstXSize ));

/* -------------------------------------------------------------------- */
/*      Allocate x,y,z coordinate arrays for transformation ... one     */
/*      scanlines worth of positions.                                   */
/* -------------------------------------------------------------------- */
    double *padfX = static_cast<double *>(
        CPLMalloc(sizeof(double) * nDstXSize) );
    double *padfY = static_cast<double *>(
        CPLMalloc(sizeof(double) * nDstXSize) );
    double *padfZ = static_cast<double *>(
        CPLMalloc(sizeof(double) * nDstXSize) );
    int *pabSuccess = static_cast<int *>( CPLMalloc(sizeof(int) * nDstXSize) );

/* -------------------------------------------------------------------- */
/*      Establish the value we will use to initialize the bands.  We    */
/*      default to -1 indicating the initial value should be read       */
/*      and preserved from the source file, but allow this to be        */
/*      overridden by passed                                            */
/*      option(s).                                                      */
/* -------------------------------------------------------------------- */
    int * const panBandInit =
        static_cast<int *>( CPLCalloc(sizeof(int), nBandCount) );
    if( CSLFetchNameValue( papszWarpOptions, "INIT" ) )
    {
        char **papszTokens =
            CSLTokenizeStringComplex( CSLFetchNameValue( papszWarpOptions,
                                                         "INIT" ),
                                      " ,", FALSE, FALSE );

        const int nTokenCount = CSLCount(papszTokens);

        for( int iBand = 0; iBand < nBandCount; iBand++ )
        {
            if( nTokenCount == 0 )
                panBandInit[iBand] = 0;
            else
                panBandInit[iBand] =
                    atoi(papszTokens[std::min(iBand, nTokenCount- 1)]);
        }

        CSLDestroy(papszTokens);
    }

/* -------------------------------------------------------------------- */
/*      Loop over all the scanlines in the output image.                */
/* -------------------------------------------------------------------- */
    for( int iDstY = 0; iDstY < nDstYSize; iDstY++ )
    {
        // Clear output buffer to "transparent" value.  Should not we
        // really be reading from the destination file to support overlay?
        for( int iBand = 0; iBand < nBandCount; iBand++ )
        {
            if( panBandInit[iBand] == -1 )
            {
                if( GDALRasterIO(
                        GDALGetRasterBand(hDstDS,iBand + 1), GF_Read,
                        0, iDstY, nDstXSize, 1,
                        papabyDstLine[iBand], nDstXSize, 1, GDT_Byte,
                        0, 0 ) != CE_None )
                {
                    bError = TRUE;
                    break;
                }
            }
            else
            {
                memset( papabyDstLine[iBand], panBandInit[iBand], nDstXSize );
            }
        }

        // Set point to transform.
        for( int iDstX = 0; iDstX < nDstXSize; iDstX++ )
        {
            padfX[iDstX] = iDstX + 0.5;
            padfY[iDstX] = iDstY + 0.5;
            padfZ[iDstX] = 0.0;
        }

        // Transform the points from destination pixel/line coordinates
        // to source pixel/line coordinates.
        pfnTransform( pTransformArg, TRUE, nDstXSize,
                      padfX, padfY, padfZ, pabSuccess );

        // Loop over the output scanline.
        for( int iDstX = 0; iDstX < nDstXSize; iDstX++ )
        {
            if( !pabSuccess[iDstX] )
                continue;

            // We test against the value before casting to avoid the
            // problem of asymmetric truncation effects around zero.  That is
            // -0.5 will be 0 when cast to an int.
            if( padfX[iDstX] < 0.0 || padfY[iDstX] < 0.0 )
                continue;

            const int iSrcX = static_cast<int>( padfX[iDstX] );
            const int iSrcY = static_cast<int>( padfY[iDstX] );

            if( iSrcX >= nSrcXSize || iSrcY >= nSrcYSize )
                continue;

            const int iSrcOffset = iSrcX + iSrcY * nSrcXSize;

            for( int iBand = 0; iBand < nBandCount; iBand++ )
                papabyDstLine[iBand][iDstX] = papabySrcData[iBand][iSrcOffset];
        }

        // Write scanline to disk.
        for( int iBand = 0; iBand < nBandCount; iBand++ )
        {
            if( GDALRasterIO(
                    GDALGetRasterBand(hDstDS,iBand+1), GF_Write,
                    0, iDstY, nDstXSize, 1,
                    papabyDstLine[iBand], nDstXSize,
                    1, GDT_Byte, 0, 0 ) != CE_None )
            {
                bError = TRUE;
                break;
            }
        }

        if( pfnProgress != nullptr )
        {
            if( !pfnProgress( (iDstY + 1) / static_cast<double>(nDstYSize),
                              "", pProgressArg ) )
            {
                CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
                bError = TRUE;
                break;
            }
        }
    }

/* -------------------------------------------------------------------- */
/*      Cleanup working buffers.                                        */
/* -------------------------------------------------------------------- */
    for( int iBand = 0; iBand < nBandCount; iBand++ )
    {
        CPLFree( papabyDstLine[iBand] );
        CPLFree( papabySrcData[iBand] );
    }

    CPLFree( panBandInit );
    CPLFree( papabyDstLine );
    CPLFree( papabySrcData );
    CPLFree( padfX );
    CPLFree( padfY );
    CPLFree( padfZ );
    CPLFree( pabSuccess );

    return !bError;
}