/******************************************************************************
 * $Id: mffdataset.cpp,v 1.27 2003/04/02 19:05:59 dron Exp $
 *
 * Project:  GView
 * Purpose:  Implementation of Atlantis MFF Support
 * Author:   Frank Warmerdam, warmerda@home.com
 *
 ******************************************************************************
 * Copyright (c) 2000, 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.
 ******************************************************************************
 *
 * $Log: mffdataset.cpp,v $
 * Revision 1.27  2003/04/02 19:05:59  dron
 * Large file support.
 *
 * Revision 1.26  2003/03/18 19:03:44  gwalter
 * Add flushing, fix coordinate interpretation.
 *
 * Revision 1.25  2003/03/13 18:34:02  gwalter
 * Fix data type determination in CreateCopy.
 *
 * Revision 1.24  2003/03/12 16:02:28  gwalter
 * Fix spheroid support.
 *
 * Revision 1.22  2003/03/03 20:10:06  gwalter
 * Updated MFF and HKV (MFF2) georeferencing support.
 *
 * Revision 1.21  2002/11/23 18:54:17  warmerda
 * added CREATIONDATATYPES metadata for drivers
 *
 * Revision 1.20  2002/09/04 06:50:37  warmerda
 * avoid static driver pointers
 *
 * Revision 1.19  2002/06/12 21:12:25  warmerda
 * update to metadata based driver info
 *
 * Revision 1.18  2001/11/11 23:51:00  warmerda
 * added required class keyword to friend declarations
 *
 * Revision 1.17  2001/07/18 04:51:57  warmerda
 * added CPL_CVSID
 *
 * Revision 1.16  2001/07/11 18:08:14  warmerda
 * Fixed pixel interpretation again.  TOP_LEFT is really 0.0,0.0 but
 * BOTTOM_RIGHT is the top left corner of the bottom right pixel!
 *
 * Revision 1.15  2001/07/11 18:00:57  warmerda
 * Fixed off-by-half-pixel issues with MFF GCPs.
 *
 * Revision 1.14  2000/10/11 16:06:51  warmerda
 * added read support for GCPn and NUM_GCPs keywords
 *
 * Revision 1.13  2000/09/25 21:20:13  warmerda
 * avoid initialization warnings
 *
 * Revision 1.12  2000/08/24 14:26:09  warmerda
 * improved handling of read-only file sets
 *
 * Revision 1.11  2000/08/15 19:28:26  warmerda
 * added help topic
 *
 * Revision 1.10  2000/07/26 21:05:21  warmerda
 * added support for reading tiled MFF files
 *
 * Revision 1.9  2000/07/12 19:21:34  warmerda
 * fixed cleanup of GCPs
 *
 * Revision 1.8  2000/06/05 17:24:06  warmerda
 * added real complex support
 *
 * Revision 1.7  2000/05/18 22:06:03  warmerda
 * added gcp and metadata support
 *
 * Revision 1.6  2000/05/15 14:18:27  warmerda
 * added COMPLEX_INTERPRETATION metadata
 *
 * Revision 1.5  2000/04/21 21:59:26  warmerda
 * added overview support
 *
 * Revision 1.4  2000/03/24 20:00:46  warmerda
 * Don't require IMAGE_FILE_FORMAT.
 *
 * Revision 1.3  2000/03/13 14:34:42  warmerda
 * avoid const problem on write
 *
 * Revision 1.2  2000/03/06 21:51:32  warmerda
 * fixed docs
 *
 * Revision 1.1  2000/03/06 19:23:20  warmerda
 * New
 *
 */

#include "rawdataset.h"
#include "cpl_string.h"
#include <ctype.h>
#include "ogr_spatialref.h"
#include "atlsci_spheroid.h"

CPL_CVSID("$Id: mffdataset.cpp,v 1.27 2003/04/02 19:05:59 dron Exp $");

CPL_C_START
void	GDALRegister_MFF(void);
CPL_C_END

typedef enum {
  MFFPRJ_NONE,
  MFFPRJ_LL,
  MFFPRJ_UTM,
  MFFPRJ_UNRECOGNIZED
} ;

static int         GetMFFProjectionType(const char * pszNewProjection);

/************************************************************************/
/* ==================================================================== */
/*				MFFDataset				*/
/* ==================================================================== */
/************************************************************************/

class MFFDataset : public RawDataset
{
    FILE	*fpImage;	// image data file.
    
    int         nGCPCount;
    GDAL_GCP    *pasGCPList;

    char *pszProjection;
    char *pszGCPProjection;
    double adfGeoTransform[6]; 

    void        ScanForGCPs();
    void        ScanForProjectionInfo();
 

  public:
    		MFFDataset();
    	        ~MFFDataset();
    
    char	**papszHdrLines;
    
    FILE        **pafpBandFiles;
    
    virtual int    GetGCPCount();
    virtual const char *GetGCPProjection();
    virtual const GDAL_GCP *GetGCPs();

    virtual const char *GetProjectionRef();
    virtual CPLErr GetGeoTransform( double * );

    static GDALDataset *Open( GDALOpenInfo * );
    static GDALDataset *Create( const char * pszFilename,
                                int nXSize, int nYSize, int nBands,
                                GDALDataType eType, char ** papszParmList );
    static GDALDataset *CreateCopy( const char * pszFilename, 
                                    GDALDataset *poSrcDS, 
                                    int bStrict, char ** papszOptions, 
                                    GDALProgressFunc pfnProgress, 
                                    void * pProgressData );

};

/************************************************************************/
/* ==================================================================== */
/*                            MFFTiledBand                              */
/* ==================================================================== */
/************************************************************************/

class MFFTiledBand : public GDALRasterBand
{
    friend class MFFDataset;

    FILE        *fpRaw;
    int         bNative;

  public:

                   MFFTiledBand( MFFDataset *, int, FILE *, int, int, 
                                 GDALDataType, int );
                   ~MFFTiledBand();

    virtual CPLErr IReadBlock( int, int, void * );
};


/************************************************************************/
/*                            MFFTiledBand()                            */
/************************************************************************/

MFFTiledBand::MFFTiledBand( MFFDataset *poDS, int nBand, FILE *fp, 
                            int nTileXSize, int nTileYSize, 
                            GDALDataType eDataType, int bNative )

{
    this->poDS = poDS;
    this->nBand = nBand;

    this->eDataType = eDataType; 

    this->bNative = bNative;

    this->nBlockXSize = nTileXSize;
    this->nBlockYSize = nTileYSize;

    this->fpRaw = fp;
}

/************************************************************************/
/*                           ~MFFTiledBand()                            */
/************************************************************************/

MFFTiledBand::~MFFTiledBand()

{
    VSIFCloseL( fpRaw );
}


/************************************************************************/
/*                             IReadBlock()                             */
/************************************************************************/

CPLErr MFFTiledBand::IReadBlock( int nBlockXOff, int nBlockYOff,
                                 void * pImage )

{
    long    nOffset;
    int     nTilesPerRow;
    int     nWordSize, nBlockSize;

    nTilesPerRow = (nRasterXSize + nBlockXSize - 1) / nBlockXSize;
    nWordSize = GDALGetDataTypeSize( eDataType ) / 8;
    nBlockSize = nWordSize * nBlockXSize * nBlockYSize;

    nOffset = nBlockSize * (nBlockXOff + nBlockYOff*nTilesPerRow);

    if( VSIFSeekL( fpRaw, nOffset, SEEK_SET ) == -1 
        || VSIFReadL( pImage, 1, nBlockSize, fpRaw ) < 1 )
    {
        CPLError( CE_Failure, CPLE_FileIO, 
                  "Read of tile %d/%d failed with fseek or fread error.", 
                  nBlockXOff, nBlockYOff );
        return CE_Failure;
    }
    
    if( !bNative && nWordSize > 1 )
    {
        if( GDALDataTypeIsComplex( eDataType ) )
        {
            GDALSwapWords( pImage, nWordSize/2, nBlockXSize*nBlockYSize, 
                           nWordSize );
            GDALSwapWords( ((GByte *) pImage)+nWordSize/2, 
                           nWordSize/2, nBlockXSize*nBlockYSize, nWordSize );
        }
        else
            GDALSwapWords( pImage, nWordSize,
                           nBlockXSize * nBlockYSize, nWordSize );
    }
    
    return CE_None;
}

/************************************************************************/
/*                      MFF Spheroids                                   */
/************************************************************************/

class MFFSpheroidList : public SpheroidList
{

public:

  MFFSpheroidList();
  ~MFFSpheroidList();

};

MFFSpheroidList :: MFFSpheroidList()
{
  num_spheroids = 18;

  epsilonR = 0.1;
  epsilonI = 0.000001;   

  spheroids[0].SetValuesByRadii("SPHERE",6371007.0,6371007.0);
  spheroids[1].SetValuesByRadii("EVEREST",6377304.0,6356103.0);
  spheroids[2].SetValuesByRadii("BESSEL",6377397.0,6356082.0);
  spheroids[3].SetValuesByRadii("AIRY",6377563.0,6356300.0);
  spheroids[4].SetValuesByRadii("CLARKE_1858",6378294.0,6356621.0);
  spheroids[5].SetValuesByRadii("CLARKE_1866",6378206.4,6356583.8);
  spheroids[6].SetValuesByRadii("CLARKE_1880",6378249.0,6356517.0);
  spheroids[7].SetValuesByRadii("HAYFORD",6378388.0,6356915.0);
  spheroids[8].SetValuesByRadii("KRASOVSKI",6378245.0,6356863.0);
  spheroids[9].SetValuesByRadii("HOUGH",6378270.0,6356794.0);
  spheroids[10].SetValuesByRadii("FISHER_60",6378166.0,6356784.0);
  spheroids[11].SetValuesByRadii("KAULA",6378165.0,6356345.0);
  spheroids[12].SetValuesByRadii("IUGG_67",6378160.0,6356775.0);
  spheroids[13].SetValuesByRadii("FISHER_68",6378150.0,6356330.0);
  spheroids[14].SetValuesByRadii("WGS_72",6378135.0,6356751.0);
  spheroids[15].SetValuesByRadii("IUGG_75",6378140.0,6356755.0);
  spheroids[16].SetValuesByRadii("WGS_84",6378137.0,6356752.0);
  spheroids[17].SetValuesByRadii("HUGHES",6378273.0,6356889.4); 
}

MFFSpheroidList::~MFFSpheroidList()

{
}

/************************************************************************/
/* ==================================================================== */
/*				MFFDataset				*/
/* ==================================================================== */
/************************************************************************/

/************************************************************************/
/*                            MFFDataset()                             */
/************************************************************************/

MFFDataset::MFFDataset()
{
    papszHdrLines = NULL;
    pafpBandFiles = NULL;
    nGCPCount = 0;
    pasGCPList = NULL;

    pszProjection = CPLStrdup("");
    pszGCPProjection = CPLStrdup("");
    adfGeoTransform[0] = 0.0;
    adfGeoTransform[1] = 1.0;
    adfGeoTransform[2] = 0.0;
    adfGeoTransform[3] = 0.0;
    adfGeoTransform[4] = 0.0;
    adfGeoTransform[5] = 1.0;
}

/************************************************************************/
/*                            ~MFFDataset()                            */
/************************************************************************/

MFFDataset::~MFFDataset()

{
    FlushCache();
    CSLDestroy( papszHdrLines );
    if( pafpBandFiles != NULL )
    {
        for( int i = 0; i < GetRasterCount(); i++ )
        {
            if( pafpBandFiles[i] != NULL )
                VSIFCloseL( pafpBandFiles[i] );
        }
        CPLFree( pafpBandFiles );
    }

    if( nGCPCount > 0 )
    {
        GDALDeinitGCPs( nGCPCount, pasGCPList );
        CPLFree( pasGCPList );
    }
    CPLFree( pszProjection );
    CPLFree( pszGCPProjection );

}

/************************************************************************/
/*                            GetGCPCount()                             */
/************************************************************************/

int MFFDataset::GetGCPCount()

{
    return nGCPCount;
}

/************************************************************************/
/*                          GetGCPProjection()                          */
/************************************************************************/

const char *MFFDataset::GetGCPProjection()

{
    if( nGCPCount > 0 )
        return pszGCPProjection;
    else
        return "";
}

/************************************************************************/
/*                          GetProjectionRef()                          */
/************************************************************************/

const char *MFFDataset::GetProjectionRef()

{
   return ( pszProjection );
}

/************************************************************************/
/*                          GetGeoTransform()                           */
/************************************************************************/

CPLErr MFFDataset::GetGeoTransform( double * padfTransform )

{
    memcpy( padfTransform,  adfGeoTransform, sizeof(double) * 6 ); 
    return( CE_None );
}

/************************************************************************/
/*                               GetGCP()                               */
/************************************************************************/

const GDAL_GCP *MFFDataset::GetGCPs()

{
    return pasGCPList;
}

/************************************************************************/
/*                            ScanForGCPs()                             */
/************************************************************************/

void MFFDataset::ScanForGCPs()

{
    int     nCorner;
    int	    NUM_GCPS = 0;
    
    if( CSLFetchNameValue(papszHdrLines, "NUM_GCPS") != NULL )
        NUM_GCPS = atoi(CSLFetchNameValue(papszHdrLines, "NUM_GCPS"));

    nGCPCount = 0;
    pasGCPList = (GDAL_GCP *) CPLCalloc(sizeof(GDAL_GCP),5+NUM_GCPS);

    for( nCorner = 0; nCorner < 5; nCorner++ )
    {
        const char * pszBase=NULL;
        double       dfRasterX=0.0, dfRasterY=0.0;
        char         szLatName[40], szLongName[40];

        if( nCorner == 0 )
        {
            dfRasterX = 0.5;
            dfRasterY = 0.5;
            pszBase = "TOP_LEFT_CORNER";
        }
        else if( nCorner == 1 )
        {
            dfRasterX = GetRasterXSize()-0.5;
            dfRasterY = 0.5;
            pszBase = "TOP_RIGHT_CORNER";
        }
        else if( nCorner == 2 )
        {
            dfRasterX = GetRasterXSize()-0.5;
            dfRasterY = GetRasterYSize()-0.5;
            pszBase = "BOTTOM_RIGHT_CORNER";
        }
        else if( nCorner == 3 )
        {
            dfRasterX = 0.5;
            dfRasterY = GetRasterYSize()-0.5;
            pszBase = "BOTTOM_LEFT_CORNER";
        }
        else if( nCorner == 4 )
        {
            dfRasterX = GetRasterXSize()/2.0;
            dfRasterY = GetRasterYSize()/2.0;
            pszBase = "CENTRE";
        }

        sprintf( szLatName, "%s_LATITUDE", pszBase );
        sprintf( szLongName, "%s_LONGITUDE", pszBase );
        
        if( CSLFetchNameValue(papszHdrLines, szLatName) != NULL
            && CSLFetchNameValue(papszHdrLines, szLongName) != NULL )
        {
            GDALInitGCPs( 1, pasGCPList + nGCPCount );
            
            CPLFree( pasGCPList[nGCPCount].pszId );

            pasGCPList[nGCPCount].pszId = CPLStrdup( pszBase );
                
            pasGCPList[nGCPCount].dfGCPX = 
                atof(CSLFetchNameValue(papszHdrLines, szLongName));
            pasGCPList[nGCPCount].dfGCPY = 
                atof(CSLFetchNameValue(papszHdrLines, szLatName));
            pasGCPList[nGCPCount].dfGCPZ = 0.0;

            pasGCPList[nGCPCount].dfGCPPixel = dfRasterX;
            pasGCPList[nGCPCount].dfGCPLine = dfRasterY;

            nGCPCount++;
        }
       
    }

/* -------------------------------------------------------------------- */
/*      Collect standalone GCPs.  They look like:                       */
/*                                                                      */
/*      GCPn = row, col, lat, long                                      */
/*      GCP1 = 1, 1, 45.0, -75.0                                        */
/* -------------------------------------------------------------------- */
    int i;

    for( i = 0; i < NUM_GCPS; i++ )
    {
        char	szName[25];
        char    **papszTokens;

        sprintf( szName, "GCP%d", i+1 );
        if( CSLFetchNameValue( papszHdrLines, szName ) == NULL )
            continue;

        papszTokens = CSLTokenizeStringComplex( 
            CSLFetchNameValue( papszHdrLines, szName ), 
            ",", FALSE, FALSE );
        if( CSLCount(papszTokens) == 4 )
        {
            GDALInitGCPs( 1, pasGCPList + nGCPCount );

            CPLFree( pasGCPList[nGCPCount].pszId );
            pasGCPList[nGCPCount].pszId = CPLStrdup( szName );

            pasGCPList[nGCPCount].dfGCPX = atof(papszTokens[3]);
            pasGCPList[nGCPCount].dfGCPY = atof(papszTokens[2]);
            pasGCPList[nGCPCount].dfGCPZ = 0.0;
            pasGCPList[nGCPCount].dfGCPPixel = atof(papszTokens[1])+0.5;
            pasGCPList[nGCPCount].dfGCPLine = atof(papszTokens[0])+0.5;

            nGCPCount++;
        }
    }
}

/************************************************************************/
/*                        ScanForProjectionInfo                         */
/************************************************************************/

void MFFDataset::ScanForProjectionInfo()
{
    const char *pszProjName, *pszOriginLong, *pszSpheroidName;
    const char *pszSpheroidEqRadius, *pszSpheroidPolarRadius;
    double eq_radius, polar_radius;
    OGRSpatialReference oProj;
    OGRSpatialReference oLL;
    MFFSpheroidList *mffEllipsoids;

    pszProjName = CSLFetchNameValue(papszHdrLines, 
                                    "PROJECTION_NAME");
    pszOriginLong = CSLFetchNameValue(papszHdrLines, 
                                      "PROJECTION_ORIGIN_LONGITUDE");
    pszSpheroidName = CSLFetchNameValue(papszHdrLines, 
                                      "SPHEROID_NAME");

    if (pszProjName == NULL)
    {
        CPLFree( pszProjection );
        CPLFree( pszGCPProjection );
        pszProjection=CPLStrdup("");
        pszGCPProjection=CPLStrdup("");
        return;
    }
    else if ((!EQUAL(pszProjName,"utm")) && (!EQUAL(pszProjName,"ll")))
    {
        CPLError(CE_Warning,CPLE_AppDefined,
                 "Warning- only utm and lat/long projections are currently supported.");
        CPLFree( pszProjection );
        CPLFree( pszGCPProjection );
        pszProjection=CPLStrdup("");
        pszGCPProjection=CPLStrdup("");
        return;
    }
    mffEllipsoids = new MFFSpheroidList;

    if( EQUAL(pszProjName,"utm") )
    {
        int nZone;

        if (pszOriginLong == NULL)
        {
            /* If origin not specified, assume 0.0 */
            CPLError(CE_Warning,CPLE_AppDefined,
                   "Warning- no projection origin longitude specified.  Assuming 0.0.");
            nZone = 31;
        }
        else
            nZone = 31 + (int) floor(atof(pszOriginLong)/6.0);


        if( pasGCPList[4].dfGCPY < 0 )
            oProj.SetUTM( nZone, 0 );
        else
            oProj.SetUTM( nZone, 1 );
     
        if (pszOriginLong != NULL)
            oProj.SetProjParm(SRS_PP_CENTRAL_MERIDIAN,atof(pszOriginLong));
        
    }

    if (pszOriginLong != NULL)
        oLL.SetProjParm(SRS_PP_LONGITUDE_OF_ORIGIN,atof(pszOriginLong));

    if ((pszSpheroidName == NULL))
    {
        CPLError(CE_Warning,CPLE_AppDefined,
            "Warning- unspecified ellipsoid.  Using wgs-84 parameters.\n");

        oProj.SetWellKnownGeogCS( "WGS84" );
        oLL.SetWellKnownGeogCS( "WGS84" );
    }
    else
    {
      if (mffEllipsoids->SpheroidInList(pszSpheroidName))
      { 
         oProj.SetGeogCS( "unknown","unknown",pszSpheroidName,
                         mffEllipsoids->GetSpheroidEqRadius(pszSpheroidName), 
                         mffEllipsoids->GetSpheroidInverseFlattening(pszSpheroidName)
                       );
         oLL.SetGeogCS( "unknown","unknown",pszSpheroidName,
                         mffEllipsoids->GetSpheroidEqRadius(pszSpheroidName), 
                         mffEllipsoids->GetSpheroidInverseFlattening(pszSpheroidName)
                      );
      }
      else if (EQUAL(pszSpheroidName,"USER_DEFINED"))
      {
          pszSpheroidEqRadius = CSLFetchNameValue(papszHdrLines,
                                                  "SPHEROID_EQUATORIAL_RADIUS");
          pszSpheroidPolarRadius = CSLFetchNameValue(papszHdrLines,
                                                  "SPHEROID_POLAR_RADIUS");
          if ((pszSpheroidEqRadius != NULL) && (pszSpheroidPolarRadius != NULL))
          {
            eq_radius = atof( pszSpheroidEqRadius );
            polar_radius = atof( pszSpheroidPolarRadius );
            oProj.SetGeogCS( "unknown","unknown","unknown",
                         eq_radius, eq_radius/(eq_radius - polar_radius));
            oLL.SetGeogCS( "unknown","unknown","unknown",
                         eq_radius, eq_radius/(eq_radius - polar_radius));          
          }
          else
          {
              CPLError(CE_Warning,CPLE_AppDefined,
                "Warning- radii not specified for user-defined ellipsoid. Using wgs-84 parameters. \n");
              oProj.SetWellKnownGeogCS( "WGS84" );
              oLL.SetWellKnownGeogCS( "WGS84" );
          }
      }
      else
      {
         CPLError(CE_Warning,CPLE_AppDefined,
            "Warning- unrecognized ellipsoid.  Using wgs-84 parameters.\n");
         oProj.SetWellKnownGeogCS( "WGS84" );
         oLL.SetWellKnownGeogCS( "WGS84" );
      }
    }  

    /* If a geotransform is sufficient to represent the GCP's (ie. each  */
    /* estimated gcp is within 0.25*pixel size of the actual value- this */
    /* is the test applied by GDALGCPsToGeoTransform), store the         */
    /* geotransform.                                                     */
    int transform_ok = FALSE;

    if (EQUAL(pszProjName,"LL"))
    {
        transform_ok = GDALGCPsToGeoTransform(nGCPCount,pasGCPList,adfGeoTransform,0);
    }
    else
    {
        OGRCoordinateTransformation *poTransform = NULL;
        double *dfPrjX, *dfPrjY; 
        int gcp_index;
        int    bSuccess = TRUE;

        dfPrjX = (double *) CPLMalloc(nGCPCount*sizeof(double));
        dfPrjY = (double *) CPLMalloc(nGCPCount*sizeof(double));


        poTransform = OGRCreateCoordinateTransformation( &oLL, &oProj );
        if( poTransform == NULL )
            bSuccess = FALSE;

        for(gcp_index=0;gcp_index<nGCPCount;gcp_index++)
        {
            dfPrjX[gcp_index] = pasGCPList[gcp_index].dfGCPX;
            dfPrjY[gcp_index] = pasGCPList[gcp_index].dfGCPY;

            if( bSuccess && !poTransform->Transform( 1, &(dfPrjX[gcp_index]), &(dfPrjY[gcp_index]) ) )
                bSuccess = FALSE;
 
        }

        if( bSuccess )
        {

            for(gcp_index=0;gcp_index<nGCPCount;gcp_index++)
            {
                pasGCPList[gcp_index].dfGCPX = dfPrjX[gcp_index];
                pasGCPList[gcp_index].dfGCPY = dfPrjY[gcp_index];

            }
            transform_ok = GDALGCPsToGeoTransform(nGCPCount,pasGCPList,adfGeoTransform,0);

        }
        CPLFree(dfPrjX);
        CPLFree(dfPrjY);

    }

    CPLFree( pszProjection );
    CPLFree( pszGCPProjection );
    pszProjection = NULL;
    pszGCPProjection = NULL;
    oProj.exportToWkt( &pszProjection );
    oProj.exportToWkt( &pszGCPProjection );

    if (transform_ok == FALSE)
    {
    /* transform is sufficient in some cases (slant range, standalone gcps) */
        adfGeoTransform[0] = 0.0;
        adfGeoTransform[1] = 1.0;
        adfGeoTransform[2] = 0.0;
        adfGeoTransform[3] = 0.0;
        adfGeoTransform[4] = 0.0;
        adfGeoTransform[5] = 1.0;
        CPLFree( pszProjection );
        pszProjection = CPLStrdup("");
    }

    delete mffEllipsoids;
  
}


/************************************************************************/
/*                                Open()                                */
/************************************************************************/

GDALDataset *MFFDataset::Open( GDALOpenInfo * poOpenInfo )

{
    int		i, bNative = TRUE;
    char        **papszHdrLines;

/* -------------------------------------------------------------------- */
/*      We assume the user is pointing to the header file.              */
/* -------------------------------------------------------------------- */
    if( poOpenInfo->nHeaderBytes < 17 || poOpenInfo->fp == NULL )
        return NULL;

    if( !EQUAL(CPLGetExtension(poOpenInfo->pszFilename),"hdr") )
        return NULL;

/* -------------------------------------------------------------------- */
/*      Load the .hdr file, and compress white space out around the     */
/*      equal sign.                                                     */
/* -------------------------------------------------------------------- */
    papszHdrLines = CSLLoad( poOpenInfo->pszFilename );
    if( papszHdrLines == NULL )
        return NULL;

    for( i = 0; papszHdrLines[i] != NULL; i++ )
    {
        int       bAfterEqual = FALSE;
        int       iSrc, iDst;
        char     *pszLine = papszHdrLines[i];

        for( iSrc=0, iDst=0; pszLine[iSrc] != '\0'; iSrc++ )
        {
            if( bAfterEqual || pszLine[iSrc] != ' ' )
            {
                pszLine[iDst++] = pszLine[iSrc];
            }

            if( iDst > 0 && pszLine[iDst-1] == '=' )
                bAfterEqual = FALSE;
        }
        pszLine[iDst] = '\0';
    }

/* -------------------------------------------------------------------- */
/*      Verify it is an MFF file.                                       */
/* -------------------------------------------------------------------- */
    if( CSLFetchNameValue( papszHdrLines, "IMAGE_FILE_FORMAT" ) != NULL
        && !EQUAL(CSLFetchNameValue(papszHdrLines,"IMAGE_FILE_FORMAT"),"MFF") )
    {
        CSLDestroy( papszHdrLines );
        return NULL;
    }

    if( (CSLFetchNameValue( papszHdrLines, "IMAGE_LINES" ) == NULL 
         || CSLFetchNameValue(papszHdrLines,"LINE_SAMPLES") == NULL)
        && (CSLFetchNameValue( papszHdrLines, "no_rows" ) == NULL 
            || CSLFetchNameValue(papszHdrLines,"no_columns") == NULL) )
    {
        CSLDestroy( papszHdrLines );
        return NULL;
    }
    
/* -------------------------------------------------------------------- */
/*      Create a corresponding GDALDataset.                             */
/* -------------------------------------------------------------------- */
    MFFDataset 	*poDS;

    poDS = new MFFDataset();

    poDS->papszHdrLines = papszHdrLines;
    
/* -------------------------------------------------------------------- */
/*      Set some dataset wide information.                              */
/* -------------------------------------------------------------------- */
    if( CSLFetchNameValue(papszHdrLines,"no_rows") != NULL
        && CSLFetchNameValue(papszHdrLines,"no_columns") != NULL )
    {
        poDS->RasterInitialize( 
            atoi(CSLFetchNameValue(papszHdrLines,"no_columns")),
            atoi(CSLFetchNameValue(papszHdrLines,"no_rows")) );
    }
    else
    {
        poDS->RasterInitialize( 
            atoi(CSLFetchNameValue(papszHdrLines,"LINE_SAMPLES")),
            atoi(CSLFetchNameValue(papszHdrLines,"IMAGE_LINES")) );
    }

    if( CSLFetchNameValue( papszHdrLines, "BYTE_ORDER" ) != NULL )
    {
#ifdef CPL_MSB
        bNative = EQUAL(CSLFetchNameValue(papszHdrLines,"BYTE_ORDER"),"MSB");
#else
        bNative = EQUAL(CSLFetchNameValue(papszHdrLines,"BYTE_ORDER"),"LSB");
#endif
    }

/* -------------------------------------------------------------------- */
/*      Get some information specific to APP tiled files.               */
/* -------------------------------------------------------------------- */
    int bTiled, nTileXSize=0, nTileYSize=0;
    const char *pszRefinedType = NULL;

    pszRefinedType = CSLFetchNameValue(papszHdrLines, "type" );

    bTiled = CSLFetchNameValue(papszHdrLines,"no_rows") != NULL;
    if( bTiled )
    {
        if( CSLFetchNameValue(papszHdrLines,"tile_size_rows") )
            nTileYSize = 
                atoi(CSLFetchNameValue(papszHdrLines,"tile_size_rows"));
        if( CSLFetchNameValue(papszHdrLines,"tile_size_columns") )
            nTileXSize = 
                atoi(CSLFetchNameValue(papszHdrLines,"tile_size_columns"));
    }

/* -------------------------------------------------------------------- */
/*      Read the directory to find matching band files.                 */
/* -------------------------------------------------------------------- */
    char       **papszDirFiles;
    char       *pszTargetBase, *pszTargetPath;
    int        nRawBand, nSkipped=0;

    pszTargetPath = CPLStrdup(CPLGetPath(poOpenInfo->pszFilename));
    pszTargetBase = CPLStrdup(CPLGetBasename( poOpenInfo->pszFilename ));
    papszDirFiles = CPLReadDir( CPLGetPath( poOpenInfo->pszFilename ) );
    if( papszDirFiles == NULL )
        return NULL;

    for( nRawBand = 0; TRUE; nRawBand++ )
    {
        const char  *pszExtension;
        int          nBand;
        GDALDataType eDataType;

        /* Find the next raw band file. */
        for( i = 0; papszDirFiles[i] != NULL; i++ )
        {
            if( !EQUAL(CPLGetBasename(papszDirFiles[i]),pszTargetBase) )
                continue;

            pszExtension = CPLGetExtension(papszDirFiles[i]);
            if( isdigit(pszExtension[1])
                && atoi(pszExtension+1) == nRawBand 
                && strchr("bBcCiIjJrRxXzZ",pszExtension[0]) != NULL )
                break;
        }

        if( papszDirFiles[i] == NULL  )
            break;

        /* open the file for required level of access */
        FILE     *fpRaw;
        const char *pszRawFilename = CPLFormFilename(pszTargetPath, 
                                                     papszDirFiles[i], NULL );

        if( poOpenInfo->eAccess == GA_Update )
            fpRaw = VSIFOpenL( pszRawFilename, "rb+" );
        else
            fpRaw = VSIFOpenL( pszRawFilename, "rb" );
        
        if( fpRaw == NULL )
        {
            CPLError( CE_Warning, CPLE_OpenFailed, 
                      "Unable to open %s ... skipping.\n", 
                      pszRawFilename );
            nSkipped++;
            continue;
        }

        pszExtension = CPLGetExtension(papszDirFiles[i]);
        if( pszRefinedType != NULL )
        {
            if( EQUAL(pszRefinedType,"C*4") )
                eDataType = GDT_CFloat32;
            else if( EQUAL(pszRefinedType,"C*8") )
                eDataType = GDT_CFloat64;
            else if( EQUAL(pszRefinedType,"R*4") )
                eDataType = GDT_Float32;
            else if( EQUAL(pszRefinedType,"R*8") )
                eDataType = GDT_Float64;
            else if( EQUAL(pszRefinedType,"I*1") )
                eDataType = GDT_Byte;
            else if( EQUAL(pszRefinedType,"I*2") )
                eDataType = GDT_Int16;
            else if( EQUAL(pszRefinedType,"I*4") )
                eDataType = GDT_Int32;
            else if( EQUAL(pszRefinedType,"U*2") )
                eDataType = GDT_UInt16;
            else if( EQUAL(pszRefinedType,"U*4") )
                eDataType = GDT_UInt32;
            else if( EQUAL(pszRefinedType,"J*1") )
                continue; /* we don't support 1 byte complex */
            else if( EQUAL(pszRefinedType,"J*2") )
                eDataType = GDT_CInt16;
            else if( EQUAL(pszRefinedType,"K*4") )
                eDataType = GDT_CInt32;
            else
                continue;
        }
        else if( EQUALN(pszExtension,"b",1) )
        {
            eDataType = GDT_Byte;
        }
        else if( EQUALN(pszExtension,"i",1) )
        {
            eDataType = GDT_UInt16;
        }
        else if( EQUALN(pszExtension,"j",1) )
        {
            eDataType = GDT_CInt16;
        }
        else if( EQUALN(pszExtension,"r",1) )
        {
            eDataType = GDT_Float32;
        }
        else if( EQUALN(pszExtension,"x",1) )
        {
            eDataType = GDT_CFloat32;
        }
        else
            continue;

        nBand = poDS->GetRasterCount() + 1;

        int nPixelOffset = GDALGetDataTypeSize(eDataType)/8;
        GDALRasterBand *poBand = NULL;
        
        if( bTiled )
        {
            poBand = 
                new MFFTiledBand( poDS, nBand, fpRaw, nTileXSize, nTileYSize,
                                  eDataType, bNative );
        }
        else
        {
            poBand = 
                new RawRasterBand( poDS, nBand, fpRaw, 0, nPixelOffset,
                                   nPixelOffset * poDS->GetRasterXSize(),
                                   eDataType, bNative, TRUE );
        }

        poDS->SetBand( nBand, poBand );
    }

    CPLFree(pszTargetPath);
    CPLFree(pszTargetBase);
    CSLDestroy(papszDirFiles);

/* -------------------------------------------------------------------- */
/*      Check if we have bands.                                         */
/* -------------------------------------------------------------------- */
    if( poDS->GetRasterCount() == 0 )
    {
        if( nSkipped > 0 && poOpenInfo->eAccess )
        {
            CPLError( CE_Failure, CPLE_OpenFailed, 
                      "Failed to open %d files that were apparently bands.\n"
                      "Perhaps this dataset is readonly?\n", 
                      nSkipped );
            delete poDS;
            return NULL;
        }
        else
        {
            CPLError( CE_Failure, CPLE_OpenFailed, 
                      "MFF header file read successfully, but no bands\n"
                      "were successfully found and opened." );
            delete poDS;
            return NULL;
        }
    }
    
/* -------------------------------------------------------------------- */
/*      Check for overviews.                                            */
/* -------------------------------------------------------------------- */
    poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );

/* -------------------------------------------------------------------- */
/*      Set all information from the .hdr that isn't well know to be    */
/*      metadata.                                                       */
/* -------------------------------------------------------------------- */
    for( i = 0; papszHdrLines[i] != NULL; i++ )
    {
        const char *pszValue;
        char       *pszName;

        pszValue = CPLParseNameValue(papszHdrLines[i], &pszName);
        if( pszName == NULL || pszValue == NULL )
            continue;

        if( !EQUAL(pszName,"END") 
            && !EQUAL(pszName,"FILE_TYPE") 
            && !EQUAL(pszName,"BYTE_ORDER") 
            && !EQUAL(pszName,"no_columns") 
            && !EQUAL(pszName,"no_rows") 
            && !EQUAL(pszName,"type") 
            && !EQUAL(pszName,"tile_size_rows") 
            && !EQUAL(pszName,"tile_size_columns") 
            && !EQUAL(pszName,"IMAGE_FILE_FORMAT") 
            && !EQUAL(pszName,"IMAGE_LINES") 
            && !EQUAL(pszName,"LINE_SAMPLES") )
        {
            poDS->SetMetadataItem( pszName, pszValue );
        }

        CPLFree( pszName );
    }

/* -------------------------------------------------------------------- */
/*      Any GCPs in header file?                                        */
/* -------------------------------------------------------------------- */
    poDS->ScanForGCPs();
    poDS->ScanForProjectionInfo();

    return( poDS );
}

int GetMFFProjectionType(const char *pszNewProjection)
{
    OGRSpatialReference *oSRS;
    char *modifiableProjection = NULL;

    if( !EQUALN(pszNewProjection,"GEOGCS",6)
       && !EQUALN(pszNewProjection,"PROJCS",6)
       && !EQUAL(pszNewProjection,"") )
      {
          return MFFPRJ_UNRECOGNIZED;       
      }
      else if (EQUAL(pszNewProjection,""))
      { 
          return MFFPRJ_NONE;  
      }
      else
      {
          /* importFromWkt updates the pointer, so don't use pszNewProjection directly */
             modifiableProjection=CPLStrdup(pszNewProjection);

             oSRS = new OGRSpatialReference;
             oSRS->importFromWkt(&modifiableProjection);

             if ((oSRS->GetAttrValue("PROJECTION") != NULL) && 
                 (EQUAL(oSRS->GetAttrValue("PROJECTION"),SRS_PT_TRANSVERSE_MERCATOR)))
             {
               return MFFPRJ_UTM;
             }
             else if ((oSRS->GetAttrValue("PROJECTION") == NULL) && (oSRS->IsGeographic()))
             {
                  return MFFPRJ_LL;
             }
             else
             {
                  return MFFPRJ_UNRECOGNIZED;
             }
      }
}

/************************************************************************/
/*                               Create()                               */
/************************************************************************/

GDALDataset *MFFDataset::Create( const char * pszFilenameIn,
                                 int nXSize, int nYSize, int nBands,
                                 GDALDataType eType,
                                 char ** papszParmList )

{
/* -------------------------------------------------------------------- */
/*      Verify input options.                                           */
/* -------------------------------------------------------------------- */
    if( eType != GDT_Byte && eType != GDT_Float32 && eType != GDT_UInt16 
        && eType != GDT_CInt16 && eType != GDT_CFloat32 )
    {
        CPLError( CE_Failure, CPLE_AppDefined,
              "Attempt to create MFF file with currently unsupported\n"
              "data type (%s).\n",
              GDALGetDataTypeName(eType) );

        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Establish the base filename (path+filename, less extension).    */
/* -------------------------------------------------------------------- */
    char	*pszBaseFilename;
    int         i;

    pszBaseFilename = (char *) CPLMalloc(strlen(pszFilenameIn)+5);
    strcpy( pszBaseFilename, pszFilenameIn );
    
    for( i = strlen(pszBaseFilename)-1; i > 0; i-- )
    {
        if( pszBaseFilename[i] == '.' )
        {
            pszBaseFilename[i] = '\0';
            break;
        }

        if( pszBaseFilename[i] == '/' || pszBaseFilename[i] == '\\' )
            break;
    }
    
/* -------------------------------------------------------------------- */
/*      Create the header file.                                         */
/* -------------------------------------------------------------------- */
    FILE       *fp;
    const char *pszFilename;

    pszFilename = CPLFormFilename( NULL, pszBaseFilename, "hdr" );

    fp = VSIFOpen( pszFilename, "wt" );
    if( fp == NULL )
    {
        CPLError( CE_Failure, CPLE_OpenFailed, 
                  "Couldn't create %s.\n", pszFilename );
        return NULL;
    }

    fprintf( fp, "IMAGE_FILE_FORMAT = MFF\n" );
    fprintf( fp, "FILE_TYPE = IMAGE\n" );
    fprintf( fp, "IMAGE_LINES = %d\n", nYSize );
    fprintf( fp, "LINE_SAMPLES = %d\n", nXSize );
#ifdef CPL_MSB     
    fprintf( fp, "BYTE_ORDER = MSB\n" );
#else
    fprintf( fp, "BYTE_ORDER = LSB\n" );
#endif

    if (CSLFetchNameValue(papszParmList,"NO_END") == NULL)
        fprintf( fp, "END\n" );
    
    VSIFClose( fp );
   
/* -------------------------------------------------------------------- */
/*      Create the data files, but don't bother writing any data to them.*/
/* -------------------------------------------------------------------- */
    for( int iBand = 0; iBand < nBands; iBand++ )
    {
        char       szExtension[4];

        if( eType == GDT_Byte )
            sprintf( szExtension, "b%02d", iBand );
        else if( eType == GDT_UInt16 )
            sprintf( szExtension, "i%02d", iBand );
        else if( eType == GDT_Float32 )
            sprintf( szExtension, "r%02d", iBand );
        else if( eType == GDT_CInt16 )
            sprintf( szExtension, "j%02d", iBand );
        else if( eType == GDT_CFloat32 )
            sprintf( szExtension, "x%02d", iBand );

        pszFilename = CPLFormFilename( NULL, pszBaseFilename, szExtension );
        fp = VSIFOpen( pszFilename, "wb" );
        if( fp == NULL )
        {
            CPLError( CE_Failure, CPLE_OpenFailed, 
                      "Couldn't create %s.\n", pszFilename );
            return NULL;
        }

        VSIFWrite( (void *) "", 1, 1, fp );
        VSIFClose( fp );
    }

/* -------------------------------------------------------------------- */
/*      Open the dataset normally.                                      */
/* -------------------------------------------------------------------- */
    GDALDataset *poDS;

    strcat( pszBaseFilename, ".hdr" );
    poDS = (GDALDataset *) GDALOpen( pszBaseFilename, GA_Update );
    CPLFree( pszBaseFilename );
    
    return poDS;
}

/************************************************************************/
/*                             CreateCopy()                             */
/************************************************************************/

GDALDataset *
MFFDataset::CreateCopy( const char * pszFilename, GDALDataset *poSrcDS, 
                        int bStrict, char ** papszOptions, 
                        GDALProgressFunc pfnProgress, void * pProgressData )

{
    MFFDataset	*poDS;
    GDALDataType eType = poSrcDS->GetRasterBand(1)->GetRasterDataType();
    int          iBand;
    char **newpapszOptions=NULL;

    if( !pfnProgress( 0.0, NULL, pProgressData ) )
        return NULL;

    /* check that other bands match type- sets type */
    /* to unknown if they differ.                  */
    for( iBand = 1; iBand < poSrcDS->GetRasterCount(); iBand++ )
     {
         GDALRasterBand *poBand = poSrcDS->GetRasterBand( iBand+1 );
         eType = GDALDataTypeUnion( eType, poBand->GetRasterDataType() );
     }

    newpapszOptions=CSLDuplicate(papszOptions);
    newpapszOptions=CSLSetNameValue(newpapszOptions,"NO_END","TRUE");

    poDS = (MFFDataset *) Create( pszFilename, 
                                  poSrcDS->GetRasterXSize(), 
                                  poSrcDS->GetRasterYSize(), 
                                  poSrcDS->GetRasterCount(), 
                                  eType, newpapszOptions );
    
    CSLDestroy(newpapszOptions);
   

    /* Check that Create worked- return Null if it didn't */
    if (poDS == NULL)
        return NULL;


/* -------------------------------------------------------------------- */
/*      Copy the image data.                                            */
/* -------------------------------------------------------------------- */
    int         nXSize = poDS->GetRasterXSize();
    int         nYSize = poDS->GetRasterYSize();
    int  	nBlockXSize, nBlockYSize, nBlockTotal, nBlocksDone;

    poDS->GetRasterBand(1)->GetBlockSize( &nBlockXSize, &nBlockYSize );

    nBlockTotal = ((nXSize + nBlockXSize - 1) / nBlockXSize)
        * ((nYSize + nBlockYSize - 1) / nBlockYSize)
        * poSrcDS->GetRasterCount();

    nBlocksDone = 0;
    for( iBand = 0; iBand < poSrcDS->GetRasterCount(); iBand++ )
    {
        GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand( iBand+1 );
        GDALRasterBand *poDstBand = poDS->GetRasterBand( iBand+1 );
        int	       iYOffset, iXOffset;
        void           *pData;
        CPLErr  eErr;


        pData = CPLMalloc(nBlockXSize * nBlockYSize
                          * GDALGetDataTypeSize(eType) / 8);

        for( iYOffset = 0; iYOffset < nYSize; iYOffset += nBlockYSize )
        {
            for( iXOffset = 0; iXOffset < nXSize; iXOffset += nBlockXSize )
            {
                int	nTBXSize, nTBYSize;

                if( !pfnProgress( (nBlocksDone++) / (float) nBlockTotal,
                                  NULL, pProgressData ) )
                {
                    CPLError( CE_Failure, CPLE_UserInterrupt, 
                              "User terminated" );
                    delete poDS;

                    GDALDriver *poMFFDriver = 
                        (GDALDriver *) GDALGetDriverByName( "MFF" );
                    poMFFDriver->Delete( pszFilename );
                    return NULL;
                }

                nTBXSize = MIN(nBlockXSize,nXSize-iXOffset);
                nTBYSize = MIN(nBlockYSize,nYSize-iYOffset);

                eErr = poSrcBand->RasterIO( GF_Read, 
                                            iXOffset, iYOffset, 
                                            nTBXSize, nTBYSize,
                                            pData, nTBXSize, nTBYSize,
                                            eType, 0, 0 );
                if( eErr != CE_None )
                {
                    return NULL;
                }
            
                eErr = poDstBand->RasterIO( GF_Write, 
                                            iXOffset, iYOffset, 
                                            nTBXSize, nTBYSize,
                                            pData, nTBXSize, nTBYSize,
                                            eType, 0, 0 );

                if( eErr != CE_None )
                {
                    return NULL;
                }
            }
        }

        CPLFree( pData );
    }

/* -------------------------------------------------------------------- */
/*      Copy georeferencing information, if enough is available.        */
/* -------------------------------------------------------------------- */


/* -------------------------------------------------------------------- */
/*      Establish the base filename (path+filename, less extension).    */
/* -------------------------------------------------------------------- */
    char	*pszBaseFilename;
    int         i;
    FILE       *fp;
    const char *pszFilenameGEO;

    pszBaseFilename = (char *) CPLMalloc(strlen(pszFilename)+5);
    strcpy( pszBaseFilename, pszFilename );
    
    for( i = strlen(pszBaseFilename)-1; i > 0; i-- )
    {
        if( pszBaseFilename[i] == '.' )
        {
            pszBaseFilename[i] = '\0';
            break;
        }

        if( pszBaseFilename[i] == '/' || pszBaseFilename[i] == '\\' )
            break;
    }

    pszFilenameGEO = CPLFormFilename( NULL, pszBaseFilename, "hdr" );

    fp = VSIFOpen( pszFilenameGEO, "at" );
    if( fp == NULL )
    {
        CPLError( CE_Failure, CPLE_OpenFailed, 
                  "Couldn't open %s for appending.\n", pszFilenameGEO );
        return NULL;
    }  
    

    /* MFF requires corner and center gcps */
    double	*padfTiepoints;
    int         src_prj;
    int         georef_created = FALSE;

    padfTiepoints = (double *) CPLMalloc(2*sizeof(double)*5);

    src_prj = GetMFFProjectionType(poSrcDS->GetProjectionRef());

    if ((src_prj != MFFPRJ_NONE) && (src_prj != MFFPRJ_UNRECOGNIZED))
    {
      double    *tempGeoTransform = NULL; 

      tempGeoTransform = (double *) CPLMalloc(6*sizeof(double));

      if (( poSrcDS->GetGeoTransform( tempGeoTransform ) == CE_None)
          && (tempGeoTransform[0] != 0.0 || tempGeoTransform[1] != 1.0
          || tempGeoTransform[2] != 0.0 || tempGeoTransform[3] != 0.0
              || tempGeoTransform[4] != 0.0 || ABS(tempGeoTransform[5]) != 1.0 ))
      {
          OGRSpatialReference oUTMorLL;
          OGRSpatialReference oLL;
          OGRCoordinateTransformation *poTransform = NULL;          
          char *srcProjection=NULL;
          char *newGCPProjection=NULL;

          padfTiepoints[0]=tempGeoTransform[0] + tempGeoTransform[1]*0.5 +\
                           tempGeoTransform[2]*0.5;

          padfTiepoints[1]=tempGeoTransform[3] + tempGeoTransform[4]*0.5 +\
                           tempGeoTransform[5]*0.5;

          padfTiepoints[2]=tempGeoTransform[0] + tempGeoTransform[2]*0.5 +\
                           tempGeoTransform[1]*(poSrcDS->GetRasterXSize()-0.5);

          padfTiepoints[3]=tempGeoTransform[3] + tempGeoTransform[5]*0.5 +\
                           tempGeoTransform[4]*(poSrcDS->GetRasterXSize()-0.5);

          padfTiepoints[4]=tempGeoTransform[0] + tempGeoTransform[1]*0.5 +\
                           tempGeoTransform[2]*(poSrcDS->GetRasterYSize()-0.5);

          padfTiepoints[5]=tempGeoTransform[3] + tempGeoTransform[4]*0.5 +\
                           tempGeoTransform[5]*(poSrcDS->GetRasterYSize()-0.5);

          padfTiepoints[6]=tempGeoTransform[0] +\
                           tempGeoTransform[1]*(poSrcDS->GetRasterXSize()-0.5) +\
                           tempGeoTransform[2]*(poSrcDS->GetRasterYSize()-0.5);

          padfTiepoints[7]=tempGeoTransform[3]+\
                           tempGeoTransform[4]*(poSrcDS->GetRasterXSize()-0.5)+\
                           tempGeoTransform[5]*(poSrcDS->GetRasterYSize()-0.5);

          padfTiepoints[8]=tempGeoTransform[0]+\
                           tempGeoTransform[1]*(poSrcDS->GetRasterXSize())/2.0+\
                           tempGeoTransform[2]*(poSrcDS->GetRasterYSize())/2.0;

          padfTiepoints[9]=tempGeoTransform[3]+\
                           tempGeoTransform[4]*(poSrcDS->GetRasterXSize())/2.0+\
                           tempGeoTransform[5]*(poSrcDS->GetRasterYSize())/2.0;

          srcProjection = CPLStrdup(poSrcDS->GetProjectionRef());
          oUTMorLL.importFromWkt(&srcProjection);
          (oUTMorLL.GetAttrNode("GEOGCS"))->exportToWkt(&newGCPProjection);
          oLL.importFromWkt(&newGCPProjection);
          if EQUALN(poSrcDS->GetProjectionRef(),"PROJCS",6)
          {
            // projected coordinate system- need to translate gcps */
            int bSuccess=TRUE;
            int index;

            poTransform = OGRCreateCoordinateTransformation( &oUTMorLL, &oLL );
            if( poTransform == NULL )
                bSuccess = FALSE;

            for (index=0;index<5;index++)
            {
                if( !bSuccess || !poTransform->Transform( 1, &(padfTiepoints[index*2]), &(padfTiepoints[index*2+1]) ) )
                  bSuccess = FALSE;
            }
            if (bSuccess == TRUE)
               georef_created = TRUE;
          }
          else
          {
            georef_created = TRUE;
          }
      }
      CPLFree(tempGeoTransform);
    } 
  
    if (georef_created == TRUE)
    {
          char *szValue;
      
          szValue = (char *) CPLMalloc(255);
    /* -------------------------------------------------------------------- */
    /*      top left                                                        */
    /* -------------------------------------------------------------------- */
          sprintf( szValue, "TOP_LEFT_CORNER_LATITUDE = %.10f\n", padfTiepoints[1] );
          fprintf( fp, szValue );
          sprintf( szValue, "TOP_LEFT_CORNER_LONGITUDE = %.10f\n", padfTiepoints[0] );
          fprintf( fp, szValue );
    /* -------------------------------------------------------------------- */
    /*      top_right                                                       */
    /* -------------------------------------------------------------------- */
          sprintf( szValue, "TOP_RIGHT_CORNER_LATITUDE = %.10f\n", padfTiepoints[3] );
          fprintf( fp, szValue );
          sprintf( szValue, "TOP_RIGHT_CORNER_LONGITUDE = %.10f\n", padfTiepoints[2] );
          fprintf( fp, szValue );
    /* -------------------------------------------------------------------- */
    /*      bottom_left                                                     */
    /* -------------------------------------------------------------------- */
          sprintf( szValue, "BOTTOM_LEFT_CORNER_LATITUDE = %.10f\n", padfTiepoints[5] );
          fprintf( fp, szValue );
          sprintf( szValue, "BOTTOM_LEFT_CORNER_LONGITUDE = %.10f\n", padfTiepoints[4] );
          fprintf( fp, szValue );
    /* -------------------------------------------------------------------- */
    /*      bottom_right                                                    */
    /* -------------------------------------------------------------------- */
          sprintf( szValue, "BOTTOM_RIGHT_CORNER_LATITUDE = %.10f\n", padfTiepoints[7] );
          fprintf( fp, szValue );
          sprintf( szValue, "BOTTOM_RIGHT_CORNER_LONGITUDE = %.10f\n", padfTiepoints[6] );
          fprintf( fp, szValue );
    /* -------------------------------------------------------------------- */
    /*      Center                                                          */
    /* -------------------------------------------------------------------- */
          sprintf( szValue, "CENTRE_LATITUDE = %.10f\n", padfTiepoints[9] );
          fprintf( fp, szValue );
          sprintf( szValue, "CENTRE_LONGITUDE = %.10f\n", padfTiepoints[8] );
          fprintf( fp, szValue );


          CPLFree(szValue);
    /* ------------------------------------------------------------------- */
    /*     Ellipsoid/projection                                            */
    /* --------------------------------------------------------------------*/

          
          OGRSpatialReference *oSRS;
          MFFSpheroidList *mffEllipsoids;
          double eq_radius, inv_flattening;
          OGRErr ogrerrorEq=OGRERR_NONE;
          OGRErr ogrerrorInvf=OGRERR_NONE;
          OGRErr ogrerrorOl=OGRERR_NONE;
          char *modifiableProjection = NULL;
          char *pszNewProjection;
          char *spheroid_name = NULL;

          pszNewProjection = CPLStrdup( poSrcDS->GetProjectionRef() );
 
          if( !EQUALN(pszNewProjection,"GEOGCS",6)
           && !EQUALN(pszNewProjection,"PROJCS",6)
           && !EQUAL(pszNewProjection,"") )
          {
            CPLError( CE_Warning, CPLE_AppDefined,
                    "Only OGC WKT Projections supported for writing to MFF.\n"
                    "%s not supported.",
                      pszNewProjection );
          }
          else if (!EQUAL(pszNewProjection,""))
          {
           

          /* importFromWkt updates the pointer, so don't use pszNewProjection directly */
             modifiableProjection=CPLStrdup(pszNewProjection);

             oSRS = new OGRSpatialReference;
             oSRS->importFromWkt(&modifiableProjection);

             if ((oSRS->GetAttrValue("PROJECTION") != NULL) && 
                 (EQUAL(oSRS->GetAttrValue("PROJECTION"),SRS_PT_TRANSVERSE_MERCATOR)))
             {
                 char *ol_txt;
    
                 ol_txt=(char *) CPLMalloc(255);
                 fprintf(fp,"PROJECTION_NAME = UTM\n");
                 sprintf(ol_txt,"PROJECTION_ORIGIN_LONGITUDE = %f\n",
                         oSRS->GetProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0,&ogrerrorOl));
                 fprintf(fp,ol_txt);
                 CPLFree(ol_txt);
             }
             else if ((oSRS->GetAttrValue("PROJECTION") == NULL) && (oSRS->IsGeographic()))
             {
                  fprintf(fp,"PROJECTION_NAME = LL\n");
             }
             else
             {
                  CPLError( CE_Warning, CPLE_AppDefined,
                  "Unrecognized projection- no georeferencing information transferred.");
                  fprintf(fp,"PROJECTION_NAME = LL\n");
             }
             eq_radius = oSRS->GetSemiMajor(&ogrerrorEq);
             inv_flattening = oSRS->GetInvFlattening(&ogrerrorInvf);
             if ((ogrerrorEq == OGRERR_NONE) && (ogrerrorInvf == OGRERR_NONE)) 
             {
                 char *ol_txt;
    
                 ol_txt=(char *) CPLMalloc(255);
                 mffEllipsoids = new MFFSpheroidList;
                 spheroid_name = mffEllipsoids->GetSpheroidNameByEqRadiusAndInvFlattening(eq_radius,inv_flattening);
                 if (spheroid_name != NULL)
                 {
                     sprintf(ol_txt,"SPHEROID_NAME = %s\n",spheroid_name );
                     fprintf(fp,ol_txt);
                 } 
                 else
                 {
                     sprintf(ol_txt,
       "SPHEROID_NAME = USER_DEFINED\nSPHEROID_EQUATORIAL_RADIUS = %.10f\nSPHEROID_POLAR_RADIUS = %.10f\n",
                     eq_radius,eq_radius*(1-1.0/inv_flattening) );
                     fprintf(fp,ol_txt);
                 }
                 delete mffEllipsoids;
                 CPLFree(ol_txt);
                 CPLFree(spheroid_name);
              }
          } 
          CPLFree(pszNewProjection);         
    } 
      
    CPLFree( padfTiepoints );
    fprintf( fp, "END\n" );
    VSIFClose( fp );
   
    /* End of georeferencing stuff */

    /* Make sure image data gets flushed */
    for( iBand = 0; iBand < poDS->GetRasterCount(); iBand++ )
    {
        RawRasterBand *poDstBand =  (RawRasterBand *) poDS->GetRasterBand( iBand+1 );
        poDstBand->FlushCache();
    }


    if( !pfnProgress( 1.0, NULL, pProgressData ) )
    {
        CPLError( CE_Failure, CPLE_UserInterrupt, 
                  "User terminated" );
        delete poDS;

        GDALDriver *poMFFDriver = 
            (GDALDriver *) GDALGetDriverByName( "MFF" );
        poMFFDriver->Delete( pszFilename );
        return NULL;
    }

    return poDS;
}


/************************************************************************/
/*                         GDALRegister_MFF()                          */
/************************************************************************/

void GDALRegister_MFF()

{
    GDALDriver	*poDriver;

    if( GDALGetDriverByName( "MFF" ) == NULL )
    {
        poDriver = new GDALDriver();
        
        poDriver->SetDescription( "MFF" );
        poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, 
                                   "Atlantis MFF Raster" );
        poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, 
                                   "frmt_various.html#MFF" );
        poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "hdr" );
        poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES, 
                                   "Byte UInt16 Float32 CInt16 CFloat32" );

        poDriver->pfnOpen = MFFDataset::Open;
        poDriver->pfnCreate = MFFDataset::Create;
        poDriver->pfnCreateCopy = MFFDataset::CreateCopy;

        GetGDALDriverManager()->RegisterDriver( poDriver );
    }
}

