/******************************************************************************
 * $Id: gdalwarp.cpp,v 1.1 2003/08/19 21:32:24 jan Exp $
 *
 * Project:  Mapinfo Image Warper
 * Purpose:  Commandline program for doing a variety of image warps, including
 *           image reprojection.
 * 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.
 ******************************************************************************
 *
 * $Log: gdalwarp.cpp,v $
 * Revision 1.1  2003/08/19 21:32:24  jan
 * These files have been moved here from thuban/extensions/thuban/
 * See there in the Attic for the older history.
 *
 * Revision 1.7  2003/07/10 14:55:11  jonathan
 * (ProjectRasterFile): Make sure
 *         the function doesn't return NULL without first setting a Python
 *         Error.
 *
 * Revision 1.6  2003/06/26 17:00:59  jonathan
 * (get_gdal_version): New.
 *         Returns the version of gdal library being used as a string.
 *
 * Revision 1.5  2003/06/17 15:25:56  jonathan
 * (MFILENAME): Padding should be to twice sizeof(void*) because there are
 *         two digits for each hex byte.
 *
 * Revision 1.4  2003/06/12 17:33:25  jonathan
 * Removed debug printing as the image is rendered. Fixes RTbug #1937.
 *
 * Revision 1.3  2003/05/30 06:30:43  jonathan
 * Removed code that allowed gdalwarp
 *         to be compiled as a standalone app. Now the code can only be
 *         called from Python which simplifies the parameter passing.
 * (ProjectRasterFile): Handle Python arguments. Remove code that
 *         checks for a destination dataset. Add more clean up code.
 *
 * Revision 1.2  2003/05/21 17:25:16  jonathan
 * Improved the error handling
 *         and passed GDAL messages back up to the Python layer. Also
 *         tried to fix some memory leaks that were present in the original
 *         utility but didn't matter because the program aborted.
 *
 * Revision 1.1  2003/05/20 15:26:17  jonathan
 * New, but basically a direct copy
 *         of the gdalwarp utility provided in GDAL. Added function calls
 *         that can be accessed from python.
 *
 * Revision 1.7  2003/03/18 17:33:00  warmerda
 * Copy colortable if there is one on the source file.
 *
 * Revision 1.6  2002/12/17 18:23:19  warmerda
 * copy GCP projection if the main projection isn't set meaningfully
 *
 * Revision 1.5  2002/12/13 04:57:15  warmerda
 * implementing remaining commandline options
 *
 * Revision 1.4  2002/12/09 16:08:57  warmerda
 * added approximating transformer
 *
 * Revision 1.3  2002/12/07 22:58:09  warmerda
 * pass initialization warp option
 *
 * Revision 1.2  2002/12/07 17:09:13  warmerda
 * added -order flag
 *
 * Revision 1.1  2002/12/07 16:13:38  warmerda
 * New
 *
 */

#include "gdal_alg.h"
#include "cpl_string.h"
#include "ogr_srs_api.h"
#include "cpl_mfile.h"

CPL_C_START
void CPL_DLL GDALRegister_THUBANBMP(void);
CPL_C_END

                                                                                
/* MFILENAME takes the name of a variable to declare as being the
 * "filename", and the address of a MFILEReceiver.
 */
#define MFILENAME(name, ptr) \
char name[ 8 + 2 * sizeof( void* ) + 1]; \
{snprintf( name, sizeof(name), "\3\1\4MFILE%0*x", 2*sizeof(void*), (int)(ptr)); \
 memset( ptr, 0, sizeof( ptr ) );}


#include <Python.h>
#define PYTHON_ERR(x) \
    {const char *str = CPLGetLastErrorMsg(); \
        str != NULL ? PyErr_SetString(x, str) : PyErr_SetString(x, "");}

#define LEAVE_NOW(e) { err = e; goto getOut; }

CPL_CVSID("$Id: gdalwarp.cpp,v 1.1 2003/08/19 21:32:24 jan Exp $");

static GDALDatasetH 
GDALWarpCreateOutput( GDALDatasetH hSrcDS, const char *pszFilename, 
                      const char *pszFormat, const char *pszSourceSRS, 
                      const char *pszTargetSRS, int nOrder, 
                      char **papszCreateOptions );
static PyObject* get_gdal_version(PyObject *, PyObject * args);
static PyObject* ProjectRasterFile(PyObject *, PyObject * args);

static double	       dfMinX=0.0, dfMinY=0.0, dfMaxX=0.0, dfMaxY=0.0;
static double	       dfXRes=0.0, dfYRes=0.0;
static int             nForcePixels=0, nForceLines=0;

static PyMethodDef gdalwarp_methods[] = {
    { "ProjectRasterFile", ProjectRasterFile, METH_VARARGS },
    { "get_gdal_version", get_gdal_version, METH_VARARGS },
    {NULL, NULL}
};

/************************************************************************/
/*                          get_gdal_version                            */
/************************************************************************/
static PyObject*
get_gdal_version(PyObject *, PyObject * args)
{
    PyObject *version = PyString_FromString(GDALVersionInfo("RELEASE_NAME"));

    Py_XINCREF(version);

    return version;
}

/************************************************************************/
/*                            initgdalwarp                              */
/************************************************************************/
#ifdef __cplusplus
extern "C"
#endif
void initgdalwarp(void)
{
    PyObject * shapelib = NULL;
    PyObject * c_api_func = NULL;
    PyObject * cobj = NULL;
                                                                                
    Py_InitModule("gdalwarp", gdalwarp_methods);
                                                                                
    Py_XDECREF(cobj);
    Py_XDECREF(c_api_func);
    Py_XDECREF(shapelib);
}

/************************************************************************/
/*                             SanitizeSRS                              */
/************************************************************************/

char *SanitizeSRS( const char *pszUserInput )

{
    OGRSpatialReferenceH hSRS;
    char *pszResult = NULL;

    hSRS = OSRNewSpatialReference( NULL );
    if( OSRSetFromUserInput( hSRS, pszUserInput ) == OGRERR_NONE )
        OSRExportToWkt( hSRS, &pszResult );

    OSRDestroySpatialReference( hSRS );

    return pszResult;
}

static PyObject*
ProjectRasterFile(PyObject *, PyObject * args)
{
    GDALDatasetH	hSrcDS, hDstDS;
    const char         *pszFormat = "GTiff";
    char               *pszTargetSRS = NULL;
    char               *pszSourceSRS = NULL;
    const char         *pszSrcFilename = NULL;
    int                 bCreateOutput = FALSE, nOrder = 0;
    void               *hTransformArg, *hGenImgProjArg=NULL, *hApproxArg=NULL;
    char               **papszWarpOptions = NULL;
    double             dfErrorThreshold = 0.125;
    GDALTransformerFunc pfnTransformer = NULL;
    char                **papszCreateOptions = NULL;
    int                 err = 0;

    PyObject * returnFileData = NULL;

    PyObject * filename;
    PyObject * srcImageArgs;
    PyObject * dstImageArgs;
    PyObject * extents;
    PyObject * resolution;
    PyObject * imageRes;

    MFILEReceiver receiver;
    MFILENAME(pszDstFilename, &receiver);

    //dummycall();

    if (!PyArg_ParseTuple(args, "OOOOOO", &filename, &srcImageArgs, 
                &dstImageArgs, 
                &extents, 
                &resolution, &imageRes))
    {
        return NULL;
    }


    dfXRes=0.0; dfYRes=0.0;

    pszSrcFilename = PyString_AsString( filename );

    pszSourceSRS = SanitizeSRS( PyString_AsString( srcImageArgs ) );
    pszTargetSRS = SanitizeSRS( PyString_AsString( dstImageArgs ) );

    dfMinX = PyFloat_AsDouble( PyTuple_GetItem( extents, 0 ) );
    dfMinY = PyFloat_AsDouble( PyTuple_GetItem( extents, 1 ) );
    dfMaxX = PyFloat_AsDouble( PyTuple_GetItem( extents, 2 ) );
    dfMaxY = PyFloat_AsDouble( PyTuple_GetItem( extents, 3 ) );

    nForcePixels = ( int )PyInt_AsLong( PyTuple_GetItem( imageRes, 0 ) );
    nForceLines  = ( int )PyInt_AsLong( PyTuple_GetItem( imageRes, 1 ) );

    pszFormat = "THUBANBMP";

    bCreateOutput = TRUE;

    GDALAllRegister();
    GDALRegister_THUBANBMP();

/* -------------------------------------------------------------------- */
/*      Open source dataset.                                            */
/* -------------------------------------------------------------------- */

    CPLPushErrorHandler( CPLQuietErrorHandler );

    hSrcDS = GDALOpen( pszSrcFilename, GA_ReadOnly );
    
    if( hSrcDS == NULL )
    {
        PYTHON_ERR( PyExc_IOError );
        LEAVE_NOW( CPLGetLastErrorNo() );
    }

    if( pszSourceSRS == NULL )
    {
        if( GDALGetProjectionRef( hSrcDS ) != NULL 
            && strlen(GDALGetProjectionRef( hSrcDS )) > 0 )
            pszSourceSRS = CPLStrdup(GDALGetProjectionRef( hSrcDS ));

        else if( GDALGetGCPProjection( hSrcDS ) != NULL
                 && strlen(GDALGetGCPProjection(hSrcDS)) > 0 
                 && GDALGetGCPCount( hSrcDS ) > 1 )
            pszSourceSRS = CPLStrdup(GDALGetGCPProjection( hSrcDS ));
        else
            pszSourceSRS = CPLStrdup("");
    }

    if( pszTargetSRS == NULL )
        pszTargetSRS = CPLStrdup(pszSourceSRS);

    hDstDS = GDALWarpCreateOutput( hSrcDS, pszDstFilename, pszFormat, 
                                   pszSourceSRS, pszTargetSRS, nOrder,
                                   papszCreateOptions );
    papszWarpOptions = CSLSetNameValue( papszWarpOptions, "INIT", "0" );
    CSLDestroy( papszCreateOptions );
    papszCreateOptions = NULL;

    if( hDstDS == NULL )
    {
        PYTHON_ERR( PyExc_IOError );
        LEAVE_NOW( CPLE_FileIO );
    }

/* -------------------------------------------------------------------- */
/*      Create a transformation object from the source to               */
/*      destination coordinate system.                                  */
/* -------------------------------------------------------------------- */
    hTransformArg = hGenImgProjArg = 
        GDALCreateGenImgProjTransformer( hSrcDS, pszSourceSRS, 
                                         hDstDS, pszTargetSRS, 
                                         TRUE, 1000.0, nOrder );

    if( hTransformArg == NULL )
    {
        PYTHON_ERR( PyExc_ValueError );
        LEAVE_NOW( CPLE_IllegalArg );
    }

    pfnTransformer = GDALGenImgProjTransform;

/* -------------------------------------------------------------------- */
/*      Warp the transformer with a linear approximator unless the      */
/*      acceptable error is zero.                                       */
/* -------------------------------------------------------------------- */
    if( dfErrorThreshold != 0.0 )
    {
        hTransformArg = hApproxArg = 
            GDALCreateApproxTransformer( GDALGenImgProjTransform, 
                                         hGenImgProjArg, dfErrorThreshold );
        pfnTransformer = GDALApproxTransform;
    }

/* -------------------------------------------------------------------- */
/*      Now actually invoke the warper to do the work.                  */
/* -------------------------------------------------------------------- */
    GDALSimpleImageWarp( hSrcDS, hDstDS, 0, NULL, 
                         pfnTransformer, hTransformArg,
                         GDALDummyProgress, NULL, papszWarpOptions );

    CSLDestroy( papszWarpOptions );
    papszWarpOptions = NULL;

    if( hApproxArg != NULL )
        GDALDestroyApproxTransformer( hApproxArg );

    if( hGenImgProjArg != NULL )
        GDALDestroyGenImgProjTransformer( hGenImgProjArg );

/* -------------------------------------------------------------------- */
/*      Cleanup.                                                        */
/* -------------------------------------------------------------------- */
    CPLErrorReset();
    err = 0;

getOut:

    GDALClose( hDstDS );
    GDALClose( hSrcDS );

    if ( receiver.data != NULL )
    {
        returnFileData = PyString_FromStringAndSize( ( char * )receiver.data, 
                                                 receiver.len );
#if 0
        FILE *fp;
        fp = fopen("dummy.bmp", "w+b");
        fwrite(receiver->data, receiver->len, 1, fp);
        fclose(fp);
        free( receiver->data );
#endif
    }

    if ( papszWarpOptions ) CSLDestroy( papszWarpOptions );
    if ( pszTargetSRS != NULL ) CPLFree(pszTargetSRS);
    if ( pszSourceSRS != NULL ) CPLFree(pszSourceSRS);

    GDALDumpOpenDatasets( stderr );
    GDALDestroyDriverManager();
    CPLPopErrorHandler();
    
    if ( !err && CPLGetLastErrorNo() )
    {
        PYTHON_ERR( PyExc_StandardError );
        err = CPLGetLastErrorNo();
    }

    if ( err ) 
    {
        CPLError( CE_Failure, err, "" );
        return NULL;
    }
    else
    {
        Py_XINCREF( returnFileData );
        return returnFileData;
    }
    
}

/************************************************************************/
/*                        GDALWarpCreateOutput()                        */
/*                                                                      */
/*      Create the output file based on various commandline options,    */
/*      and the input file.                                             */
/************************************************************************/

static GDALDatasetH 
GDALWarpCreateOutput( GDALDatasetH hSrcDS, const char *pszFilename, 
                      const char *pszFormat, const char *pszSourceSRS, 
                      const char *pszTargetSRS, int nOrder,
                      char **papszCreateOptions )

{
    GDALDriverH hDriver;
    GDALDatasetH hDstDS;
    void *hTransformArg;
    double adfDstGeoTransform[6];
    int nPixels=0, nLines=0;
    GDALColorTableH hCT;

/* -------------------------------------------------------------------- */
/*      Find the output driver.                                         */
/* -------------------------------------------------------------------- */
    hDriver = GDALGetDriverByName( pszFormat );
    if( hDriver == NULL 
        || GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) == NULL )
    {
        int	iDr;
        
        /*
        printf( "Output driver `%s' not recognised or does not support\n", 
                pszFormat );
        printf( "direct output file creation.  The following format drivers are configured\n"
                "and support direct output:\n" );
        */

        for( iDr = 0; iDr < GDALGetDriverCount(); iDr++ )
        {
            GDALDriverH hDriver = GDALGetDriver(iDr);

            if( GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL) != NULL )
            {
                /*
                printf( "  %s: %s\n",
                        GDALGetDriverShortName( hDriver  ),
                        GDALGetDriverLongName( hDriver ) );
                */
            }
        }
        //printf( "\n" );
        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Create a transformation object from the source to               */
/*      destination coordinate system.                                  */
/* -------------------------------------------------------------------- */
    hTransformArg = 
        GDALCreateGenImgProjTransformer( hSrcDS, pszSourceSRS, 
                                         NULL, pszTargetSRS, 
                                         TRUE, 1000.0, nOrder );

    if( hTransformArg == NULL )
        return NULL;

/* -------------------------------------------------------------------- */
/*      Get approximate output definition.                              */
/* -------------------------------------------------------------------- */
    if( GDALSuggestedWarpOutput( hSrcDS, 
                                 GDALGenImgProjTransform, hTransformArg, 
                                 adfDstGeoTransform, &nPixels, &nLines )
        != CE_None )
        return NULL;

    GDALDestroyGenImgProjTransformer( hTransformArg );

/* -------------------------------------------------------------------- */
/*      Did the user override some parameters?                          */
/* -------------------------------------------------------------------- */
    if( dfXRes != 0.0 && dfYRes != 0.0 )
    {
        CPLAssert( nPixels == 0 && nLines == 0 );
        if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
        {
            dfMinX = adfDstGeoTransform[0];
            dfMaxX = adfDstGeoTransform[0] + adfDstGeoTransform[1] * nPixels;
            dfMaxY = adfDstGeoTransform[3];
            dfMinY = adfDstGeoTransform[3] + adfDstGeoTransform[5] * nLines;
        }

        nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
        nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);
        adfDstGeoTransform[0] = dfMinX;
        adfDstGeoTransform[3] = dfMaxY;
        adfDstGeoTransform[1] = dfXRes;
        adfDstGeoTransform[5] = -dfYRes;
    }

    else if( nForcePixels != 0 && nForceLines != 0 )
    {
        if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
        {
            dfMinX = adfDstGeoTransform[0];
            dfMaxX = adfDstGeoTransform[0] + adfDstGeoTransform[1] * nPixels;
            dfMaxY = adfDstGeoTransform[3];
            dfMinY = adfDstGeoTransform[3] + adfDstGeoTransform[5] * nLines;
        }

        dfXRes = (dfMaxX - dfMinX) / nForcePixels;
        dfYRes = (dfMaxY - dfMinY) / nForceLines;

        adfDstGeoTransform[0] = dfMinX;
        adfDstGeoTransform[3] = dfMaxY;
        adfDstGeoTransform[1] = dfXRes;
        adfDstGeoTransform[5] = -dfYRes;

        nPixels = nForcePixels;
        nLines = nForceLines;
    }

    else if( dfMinX != 0.0 || dfMinY != 0.0 || dfMaxX != 0.0 || dfMaxY != 0.0 )
    {
        dfXRes = adfDstGeoTransform[1];
        dfYRes = fabs(adfDstGeoTransform[5]);

        nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
        nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);

        adfDstGeoTransform[0] = dfMinX;
        adfDstGeoTransform[3] = dfMaxY;
    }

/* -------------------------------------------------------------------- */
/*      Create the output file.                                         */
/* -------------------------------------------------------------------- */
    //printf( "Creating output file is that %dP x %dL.\n", nPixels, nLines );

    hDstDS = GDALCreate( hDriver, pszFilename, nPixels, nLines, 
                         GDALGetRasterCount(hSrcDS),
                         GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1)),
                         papszCreateOptions );

    if( hDstDS == NULL )
        return NULL;

/* -------------------------------------------------------------------- */
/*      Write out the projection definition.                            */
/* -------------------------------------------------------------------- */
    GDALSetProjection( hDstDS, pszTargetSRS );
    GDALSetGeoTransform( hDstDS, adfDstGeoTransform );

/* -------------------------------------------------------------------- */
/*      Copy the color table, if required.                              */
/* -------------------------------------------------------------------- */
    hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS,1) );
    if( hCT != NULL )
        GDALSetRasterColorTable( GDALGetRasterBand(hDstDS,1), hCT );

    return hDstDS;
}
    
