/******************************************************************************
* $Id: idadataset.cpp,v 1.9 2005/09/14 01:12:00 fwarmerdam Exp $
*
* Project: IDA Raster Driver
* Purpose: Implemenents IDA driver/dataset/rasterband.
* Author: Frank Warmerdam, warmerdam@pobox.com
*
******************************************************************************
* Copyright (c) 2004, Frank Warmerdam <warmerdam@pobox.com>
*
* 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: idadataset.cpp,v $
* Revision 1.9 2005/09/14 01:12:00 fwarmerdam
* Fixed creation data types and help link.
*
* Revision 1.8 2005/08/23 16:21:25 fwarmerdam
* Fixed help link, dont barf in open if fp null.
*
* Revision 1.7 2005/05/05 13:55:42 fwarmerdam
* PAM Enable
*
* Revision 1.6 2005/02/08 20:22:39 fwarmerdam
* Added SetProjection support.
*
* Revision 1.5 2005/01/15 16:11:12 fwarmerdam
* lots of work on create support
*
* Revision 1.4 2005/01/04 22:14:33 fwarmerdam
* Nailed down open checking fairly closely.
*
* Revision 1.3 2004/12/26 21:25:36 fwarmerdam
* avoid spurious opens
*
* Revision 1.2 2004/12/26 17:37:46 fwarmerdam
* fix unix/dos format
*
* Revision 1.1 2004/12/26 16:17:00 fwarmerdam
* New
*
*/
#include "rawdataset.h"
#include "ogr_spatialref.h"
CPL_CVSID("$Id: idadataset.cpp,v 1.9 2005/09/14 01:12:00 fwarmerdam Exp $");
CPL_C_START
void GDALRegister_IDA(void);
CPL_C_END
// convert a Turbo Pascal real into a double
static double tp2c(GByte *r);
// convert a double into a Turbo Pascal real
static void c2tp(double n, GByte *r);
/************************************************************************/
/* ==================================================================== */
/* IDADataset */
/* ==================================================================== */
/************************************************************************/
class IDADataset : public RawDataset
{
friend class IDARasterBand;
int nImageType;
int nProjection;
char szTitle[81];
double dfLatCenter;
double dfLongCenter;
double dfXCenter;
double dfYCenter;
double dfDX;
double dfDY;
double dfParallel1;
double dfParallel2;
int nLower;
int nUpper;
int nMissing;
double dfM;
double dfB;
int nDecimals;
FILE *fpRaw;
char *pszProjection;
double adfGeoTransform[6];
void ProcessGeoref();
GByte abyHeader[512];
int bHeaderDirty;
public:
IDADataset();
~IDADataset();
virtual void FlushCache();
virtual const char *GetProjectionRef(void);
virtual CPLErr SetProjection( const char * );
virtual CPLErr GetGeoTransform( double * );
virtual CPLErr SetGeoTransform( double * );
static GDALDataset *Open( GDALOpenInfo * );
static GDALDataset *Create( const char * pszFilename,
int nXSize, int nYSize, int nBands,
GDALDataType eType,
char ** /* papszParmList */ );
};
/************************************************************************/
/* ==================================================================== */
/* IDARasterBand */
/* ==================================================================== */
/************************************************************************/
class IDARasterBand : public RawRasterBand
{
friend class IDADataset;
public:
IDARasterBand( IDADataset *poDSIn, FILE *fpRaw, int nXSize );
virtual ~IDARasterBand();
virtual double GetOffset( int *pbSuccess = NULL );
virtual CPLErr SetOffset( double dfNewValue );
virtual double GetScale( int *pbSuccess = NULL );
virtual CPLErr SetScale( double dfNewValue );
virtual double GetNoDataValue( int *pbSuccess = NULL );
};
/************************************************************************/
/* IDARasterBand */
/************************************************************************/
IDARasterBand::IDARasterBand( IDADataset *poDSIn,
FILE *fpRaw, int nXSize )
: RawRasterBand( poDSIn, 1, fpRaw, 512, 1, nXSize,
GDT_Byte, FALSE, FALSE )
{
}
/************************************************************************/
/* ~IDARasterBand() */
/************************************************************************/
IDARasterBand::~IDARasterBand()
{
}
/************************************************************************/
/* GetNoDataValue() */
/************************************************************************/
double IDARasterBand::GetNoDataValue( int *pbSuccess )
{
if( pbSuccess != NULL )
*pbSuccess = TRUE;
return ((IDADataset *) poDS)->nMissing;
}
/************************************************************************/
/* GetOffset() */
/************************************************************************/
double IDARasterBand::GetOffset( int *pbSuccess )
{
if( pbSuccess != NULL )
*pbSuccess = TRUE;
return ((IDADataset *) poDS)->dfB;
}
/************************************************************************/
/* SetOffset() */
/************************************************************************/
CPLErr IDARasterBand::SetOffset( double dfNewValue )
{
IDADataset *poIDS = (IDADataset *) poDS;
if( dfNewValue == poIDS->dfB )
return CE_None;
if( poIDS->nImageType != 200 )
{
CPLError( CE_Failure, CPLE_AppDefined,
"Setting explicit offset only support for image type 200.");
return CE_Failure;
}
poIDS->dfB = dfNewValue;
c2tp( dfNewValue, poIDS->abyHeader + 178 );
poIDS->bHeaderDirty = TRUE;
return CE_None;
}
/************************************************************************/
/* GetScale() */
/************************************************************************/
double IDARasterBand::GetScale( int *pbSuccess )
{
if( pbSuccess != NULL )
*pbSuccess = TRUE;
return ((IDADataset *) poDS)->dfM;
}
/************************************************************************/
/* SetScale() */
/************************************************************************/
CPLErr IDARasterBand::SetScale( double dfNewValue )
{
IDADataset *poIDS = (IDADataset *) poDS;
if( dfNewValue == poIDS->dfM )
return CE_None;
if( poIDS->nImageType != 200 )
{
CPLError( CE_Failure, CPLE_AppDefined,
"Setting explicit scale only support for image type 200.");
return CE_Failure;
}
poIDS->dfM = dfNewValue;
c2tp( dfNewValue, poIDS->abyHeader + 172 );
poIDS->bHeaderDirty = TRUE;
return CE_None;
}
/************************************************************************/
/* ==================================================================== */
/* IDADataset */
/* ==================================================================== */
/************************************************************************/
/************************************************************************/
/* IDADataset() */
/************************************************************************/
IDADataset::IDADataset()
{
fpRaw = NULL;
pszProjection = NULL;
adfGeoTransform[0] = 0.0;
adfGeoTransform[1] = 1.0;
adfGeoTransform[2] = 0.0;
adfGeoTransform[3] = 0.0;
adfGeoTransform[4] = 0.0;
adfGeoTransform[5] = 1.0;
bHeaderDirty = FALSE;
}
/************************************************************************/
/* ~IDADataset() */
/************************************************************************/
IDADataset::~IDADataset()
{
FlushCache();
if( fpRaw != NULL )
VSIFClose( fpRaw );
CPLFree( pszProjection );
}
/************************************************************************/
/* ProcessGeoref() */
/************************************************************************/
void IDADataset::ProcessGeoref()
{
OGRSpatialReference oSRS;
if( nProjection == 3 )
{
oSRS.SetWellKnownGeogCS( "WGS84" );
}
else if( nProjection == 4 )
{
oSRS.SetLCC( dfParallel1, dfParallel2,
dfLatCenter, dfLongCenter,
0.0, 0.0 );
oSRS.SetGeogCS( "Clarke 1866", "Clarke 1866", "Clarke 1866",
6378206.4, 293.97869821389662 );
}
else if( nProjection == 6 )
{
oSRS.SetLAEA( dfLatCenter, dfLongCenter, 0.0, 0.0 );
oSRS.SetGeogCS( "Sphere", "Sphere", "Sphere",
6370997.0, 0.0 );
}
else if( nProjection == 8 )
{
oSRS.SetACEA( dfParallel1, dfParallel2,
dfLatCenter, dfLongCenter,
0.0, 0.0 );
oSRS.SetGeogCS( "Clarke 1866", "Clarke 1866", "Clarke 1866",
6378206.4, 293.97869821389662 );
}
else if( nProjection == 9 )
{
oSRS.SetGH( dfLongCenter, 0.0, 0.0 );
oSRS.SetGeogCS( "Sphere", "Sphere", "Sphere",
6370997.0, 0.0 );
}
if( oSRS.GetRoot() != NULL )
{
CPLFree( pszProjection );
pszProjection = NULL;
oSRS.exportToWkt( &pszProjection );
}
adfGeoTransform[0] = 0 - dfDX * dfXCenter;
adfGeoTransform[1] = dfDX;
adfGeoTransform[2] = 0.0;
adfGeoTransform[3] = dfDY * dfYCenter;
adfGeoTransform[4] = 0.0;
adfGeoTransform[5] = -dfDY;
if( nProjection == 3 )
{
adfGeoTransform[0] += dfLongCenter;
adfGeoTransform[3] += dfLatCenter;
}
}
/************************************************************************/
/* FlushCache() */
/************************************************************************/
void IDADataset::FlushCache()
{
RawDataset::FlushCache();
if( bHeaderDirty )
{
VSIFSeek( fpRaw, 0, SEEK_SET );
VSIFWrite( abyHeader, 512, 1, fpRaw );
bHeaderDirty = FALSE;
}
}
/************************************************************************/
/* GetGeoTransform() */
/************************************************************************/
CPLErr IDADataset::GetGeoTransform( double *padfGeoTransform )
{
memcpy( padfGeoTransform, adfGeoTransform, sizeof(double) * 6 );
return CE_None;
}
/************************************************************************/
/* SetGeoTransform() */
/************************************************************************/
CPLErr IDADataset::SetGeoTransform( double *padfGeoTransform )
{
if( padfGeoTransform[2] != 0.0 || padfGeoTransform[4] != 0.0 )
return GDALPamDataset::SetGeoTransform( padfGeoTransform );
memcpy( adfGeoTransform, padfGeoTransform, sizeof(double) * 6 );
bHeaderDirty = TRUE;
dfDX = adfGeoTransform[1];
dfDY = -adfGeoTransform[5];
dfXCenter = -adfGeoTransform[0] / dfDX;
dfYCenter = adfGeoTransform[3] / dfDY;
c2tp( dfDX, abyHeader + 144 );
c2tp( dfDY, abyHeader + 150 );
c2tp( dfXCenter, abyHeader + 132 );
c2tp( dfYCenter, abyHeader + 138 );
return CE_None;
}
/************************************************************************/
/* GetProjectionRef() */
/************************************************************************/
const char *IDADataset::GetProjectionRef()
{
if( pszProjection )
return pszProjection;
else
return "";
}
/************************************************************************/
/* SetProjection() */
/************************************************************************/
CPLErr IDADataset::SetProjection( const char *pszWKTIn )
{
OGRSpatialReference oSRS;
oSRS.importFromWkt( (char **) &pszWKTIn );
if( !oSRS.IsGeographic() && !oSRS.IsProjected() )
GDALPamDataset::SetProjection( pszWKTIn );
/* -------------------------------------------------------------------- */
/* Clear projection parameters. */
/* -------------------------------------------------------------------- */
dfParallel1 = 0.0;
dfParallel2 = 0.0;
dfLatCenter = 0.0;
dfLongCenter = 0.0;
/* -------------------------------------------------------------------- */
/* Geographic. */
/* -------------------------------------------------------------------- */
if( oSRS.IsGeographic() )
{
// If no change, just return.
if( nProjection == 3 )
return CE_None;
nProjection = 3;
}
/* -------------------------------------------------------------------- */
/* Verify we don't have a false easting or northing as these */
/* will be ignored for the projections we do support. */
/* -------------------------------------------------------------------- */
if( oSRS.GetProjParm( SRS_PP_FALSE_EASTING ) != 0.0
|| oSRS.GetProjParm( SRS_PP_FALSE_NORTHING ) != 0.0 )
{
CPLError( CE_Failure, CPLE_AppDefined,
"Attempt to set a projection on an IDA file with a non-zero\n"
"false easting and/or northing. This is not supported." );
return CE_Failure;
}
/* -------------------------------------------------------------------- */
/* Lambert Conformal Conic. Note that we don't support false */
/* eastings or nothings. */
/* -------------------------------------------------------------------- */
const char *pszProjection = oSRS.GetAttrValue( "PROJECTION" );
if( pszProjection == NULL )
{
/* do nothing - presumably geographic */;
}
else if( EQUAL(pszProjection,SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP) )
{
nProjection = 4;
dfParallel1 = oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1,0.0);
dfParallel2 = oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2,0.0);
dfLatCenter = oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0);
dfLongCenter = oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0);
}
else if( EQUAL(pszProjection,SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA) )
{
nProjection = 6;
dfLatCenter = oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0);
dfLongCenter = oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0);
}
else if( EQUAL(pszProjection,SRS_PT_ALBERS_CONIC_EQUAL_AREA) )
{
nProjection = 8;
dfParallel1 = oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1,0.0);
dfParallel2 = oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2,0.0);
dfLatCenter = oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0);
dfLongCenter = oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0);
}
else if( EQUAL(pszProjection,SRS_PT_GOODE_HOMOLOSINE) )
{
nProjection = 9;
dfLongCenter = oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0);
}
else
{
return GDALPamDataset::SetProjection( pszWKTIn );
}
/* -------------------------------------------------------------------- */
/* Update header and mark it as dirty. */
/* -------------------------------------------------------------------- */
bHeaderDirty = TRUE;
abyHeader[23] = nProjection;
c2tp( dfLatCenter, abyHeader + 120 );
c2tp( dfLongCenter, abyHeader + 126 );
c2tp( dfParallel1, abyHeader + 156 );
c2tp( dfParallel2, abyHeader + 162 );
return CE_None;
}
/************************************************************************/
/* Open() */
/************************************************************************/
GDALDataset *IDADataset::Open( GDALOpenInfo * poOpenInfo )
{
/* -------------------------------------------------------------------- */
/* Is this an IDA file? */
/* -------------------------------------------------------------------- */
int nXSize, nYSize;
GIntBig nExpectedFileSize, nActualFileSize;
if( poOpenInfo->fp == NULL )
return NULL;
if( poOpenInfo->nHeaderBytes < 512 )
return NULL;
// projection legal?
if( poOpenInfo->pabyHeader[23] > 10 )
return NULL;
// imagetype legal?
if( (poOpenInfo->pabyHeader[22] > 14
&& poOpenInfo->pabyHeader[22] < 100)
|| (poOpenInfo->pabyHeader[22] > 114
&& poOpenInfo->pabyHeader[22] != 200 ) )
return NULL;
nXSize = poOpenInfo->pabyHeader[30] + poOpenInfo->pabyHeader[31] * 256;
nYSize = poOpenInfo->pabyHeader[32] + poOpenInfo->pabyHeader[33] * 256;
if( nXSize == 0 || nYSize == 0 )
return NULL;
// The file just be exactly the image size + header size in length.
nExpectedFileSize = nXSize * nYSize + 512;
VSIFSeek( poOpenInfo->fp, 0, SEEK_END );
nActualFileSize = VSIFTell( poOpenInfo->fp );
VSIRewind( poOpenInfo->fp );
if( nActualFileSize != nExpectedFileSize )
return NULL;
/* -------------------------------------------------------------------- */
/* Create the dataset. */
/* -------------------------------------------------------------------- */
IDADataset *poDS = new IDADataset();
memcpy( poDS->abyHeader, poOpenInfo->pabyHeader, 512 );
/* -------------------------------------------------------------------- */
/* Parse various values out of the header. */
/* -------------------------------------------------------------------- */
poDS->nImageType = poOpenInfo->pabyHeader[22];
poDS->nProjection = poOpenInfo->pabyHeader[23];
poDS->nRasterYSize = poOpenInfo->pabyHeader[30]
+ poOpenInfo->pabyHeader[31] * 256;
poDS->nRasterXSize = poOpenInfo->pabyHeader[32]
+ poOpenInfo->pabyHeader[33] * 256;
strncpy( poDS->szTitle, (const char *) poOpenInfo->pabyHeader+38, 80 );
poDS->szTitle[80] = '\0';
int nLastTitleChar = strlen(poDS->szTitle)-1;
while( nLastTitleChar > -1
&& (poDS->szTitle[nLastTitleChar] == 10
|| poDS->szTitle[nLastTitleChar] == 13
|| poDS->szTitle[nLastTitleChar] == ' ') )
poDS->szTitle[nLastTitleChar--] = '\0';
poDS->dfLatCenter = tp2c( poOpenInfo->pabyHeader + 120 );
poDS->dfLongCenter = tp2c( poOpenInfo->pabyHeader + 126 );
poDS->dfXCenter = tp2c( poOpenInfo->pabyHeader + 132 );
poDS->dfYCenter = tp2c( poOpenInfo->pabyHeader + 138 );
poDS->dfDX = tp2c( poOpenInfo->pabyHeader + 144 );
poDS->dfDY = tp2c( poOpenInfo->pabyHeader + 150 );
poDS->dfParallel1 = tp2c( poOpenInfo->pabyHeader + 156 );
poDS->dfParallel2 = tp2c( poOpenInfo->pabyHeader + 162 );
poDS->ProcessGeoref();
poDS->SetMetadataItem( "TITLE", poDS->szTitle );
/* -------------------------------------------------------------------- */
/* Handle various image types. */
/* -------------------------------------------------------------------- */
/*
GENERIC = 0
FEW S NDVI = 1
EROS NDVI = 6
ARTEMIS CUTOFF = 10
ARTEMIS RECODE = 11
ARTEMIS NDVI = 12
ARTEMIS FEWS = 13
ARTEMIS NEWNASA = 14
GENERIC DIFF = 100
FEW S NDVI DIFF = 101
EROS NDVI DIFF = 106
ARTEMIS CUTOFF DIFF = 110
ARTEMIS RECODE DIFF = 111
ARTEMIS NDVI DIFF = 112
ARTEMIS FEWS DIFF = 113
ARTEMIS NEWNASA DIFF = 114
CALCULATED =200
*/
poDS->nMissing = 0;
switch( poDS->nImageType )
{
case 1:
poDS->SetMetadataItem( "IMAGETYPE", "1, FEWS NDVI" );
poDS->dfM = 1/256.0;
poDS->dfB = -82/256.0;
break;
case 6:
poDS->SetMetadataItem( "IMAGETYPE", "6, EROS NDVI" );
poDS->dfM = 1/100.0;
poDS->dfB = -100/100.0;
break;
case 10:
poDS->SetMetadataItem( "IMAGETYPE", "10, ARTEMIS CUTOFF" );
poDS->dfM = 1.0;
poDS->dfB = 0.0;
poDS->nMissing = 254;
break;
case 11:
poDS->SetMetadataItem( "IMAGETYPE", "11, ARTEMIS RECODE" );
poDS->dfM = 4.0;
poDS->dfB = 0.0;
poDS->nMissing = 254;
break;
case 12: /* ANDVI */
poDS->SetMetadataItem( "IMAGETYPE", "12, ARTEMIS NDVI" );
poDS->dfM = 4/500.0;
poDS->dfB = -3/500.0 - 1.0;
poDS->nMissing = 254;
break;
case 13: /* AFEWS */
poDS->SetMetadataItem( "IMAGETYPE", "13, ARTEMIS FEWS" );
poDS->dfM = 1/256.0;
poDS->dfB = -82/256.0;
poDS->nMissing = 254;
break;
case 14: /* NEWNASA */
poDS->SetMetadataItem( "IMAGETYPE", "13, ARTEMIS NEWNASA" );
poDS->dfM = 0.75/250.0;
poDS->dfB = 0.0;
poDS->nMissing = 254;
break;
case 101: /* NDVI_DIFF (FEW S) */
poDS->dfM = 1/128.0;
poDS->dfB = -1.0;
poDS->nMissing = 0;
break;
case 106: /* EROS_DIFF */
poDS->dfM = 1/50.0;
poDS->dfB = -128/50.0;
poDS->nMissing = 0;
break;
case 110: /* CUTOFF_DIFF */
poDS->dfM = 2.0;
poDS->dfB = -128*2;
poDS->nMissing = 254;
break;
case 111: /* RECODE_DIFF */
poDS->dfM = 8;
poDS->dfB = -128*8;
poDS->nMissing = 254;
break;
case 112: /* ANDVI_DIFF */
poDS->dfM = 8/1000.0;
poDS->dfB = (-128*8)/1000.0;
poDS->nMissing = 254;
break;
case 113: /* AFEWS_DIFF */
poDS->dfM = 1/128.0;
poDS->dfB = -1;
poDS->nMissing = 254;
break;
case 114: /* NEWNASA_DIFF */
poDS->dfM = 0.75/125.0;
poDS->dfB = -128*poDS->dfM;
poDS->nMissing = 254;
break;
case 200:
/* we use the values from the header */
poDS->dfM = tp2c( poOpenInfo->pabyHeader + 171 );
poDS->dfB = tp2c( poOpenInfo->pabyHeader + 177 );
poDS->nMissing = poOpenInfo->pabyHeader[170];
break;
default:
poDS->dfM = 1.0;
poDS->dfB = 0.0;
break;
}
/* -------------------------------------------------------------------- */
/* Create the band. */
/* -------------------------------------------------------------------- */
if( poOpenInfo->eAccess == GA_ReadOnly )
{
poDS->fpRaw = poOpenInfo->fp;
poOpenInfo->fp = NULL;
}
else
{
poDS->fpRaw = VSIFOpen( poOpenInfo->pszFilename, "rb+" );
poDS->eAccess = GA_Update;
if( poDS->fpRaw == NULL )
{
CPLError( CE_Failure, CPLE_OpenFailed,
"Failed to open %s for write access.",
poOpenInfo->pszFilename );
return NULL;
}
}
poDS->SetBand( 1, new IDARasterBand( poDS, poDS->fpRaw,
poDS->nRasterXSize ) );
/* -------------------------------------------------------------------- */
/* Check for overviews. */
/* -------------------------------------------------------------------- */
poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
/* -------------------------------------------------------------------- */
/* Initialize any PAM information. */
/* -------------------------------------------------------------------- */
poDS->SetDescription( poOpenInfo->pszFilename );
poDS->TryLoadXML();
return( poDS );
}
/************************************************************************/
/* tp2c() */
/* */
/* convert a Turbo Pascal real into a double */
/************************************************************************/
static double tp2c(GByte *r)
{
double mant;
int sign, exp, i;
// handle 0 case
if (r[0] == 0)
return 0.0;
// extract sign: bit 7 of byte 5
sign = r[5] & 0x80 ? -1 : 1;
// extract mantissa from first bit of byte 1 to bit 7 of byte 5
mant = 0;
for (i = 1; i < 5; i++)
mant = (r[i] + mant) / 256;
mant = (mant + (r[5] & 0x7F)) / 128 + 1;
// extract exponent
exp = r[0] - 129;
// compute the damned number
return sign * ldexp(mant, exp);
}
/************************************************************************/
/* c2tp() */
/* */
/* convert a double into a Turbo Pascal real */
/************************************************************************/
static void c2tp(double x, GByte *r)
{
double mant, temp;
int negative, exp, i;
// handle 0 case
if (x == 0.0)
{
for (i = 0; i < 6; r[i++] = 0);
return;
}
// compute mantissa, sign and exponent
mant = frexp(x, &exp) * 2 - 1;
exp--;
negative = 0;
if (mant < 0)
{
mant = -mant;
negative = 1;
}
// stuff mantissa into Turbo Pascal real
mant = modf(mant * 128, &temp);
r[5] = (unsigned char) temp;
for (i = 4; i >= 1; i--)
{
mant = modf(mant * 256, &temp);
r[i] = (unsigned char) temp;
}
// add sign
if (negative)
r[5] |= 0x80;
// put exponent
r[0] = exp + 129;
}
/************************************************************************/
/* Create() */
/************************************************************************/
GDALDataset *IDADataset::Create( const char * pszFilename,
int nXSize, int nYSize, int nBands,
GDALDataType eType,
char ** /* papszParmList */ )
{
/* -------------------------------------------------------------------- */
/* Verify input options. */
/* -------------------------------------------------------------------- */
if( eType != GDT_Byte || nBands != 1 )
{
CPLError( CE_Failure, CPLE_AppDefined,
"Only 1 band, Byte datasets supported for IDA format." );
return NULL;
}
/* -------------------------------------------------------------------- */
/* Try to create the file. */
/* -------------------------------------------------------------------- */
FILE *fp;
fp = VSIFOpen( pszFilename, "wb" );
if( fp == NULL )
{
CPLError( CE_Failure, CPLE_OpenFailed,
"Attempt to create file `%s' failed.\n",
pszFilename );
return NULL;
}
/* -------------------------------------------------------------------- */
/* Prepare formatted header. */
/* -------------------------------------------------------------------- */
GByte abyHeader[512];
memset( abyHeader, 0, sizeof(abyHeader) );
abyHeader[22] = 200; /* image type - CALCULATED */
abyHeader[23] = 0; /* projection - NONE */
abyHeader[30] = nYSize % 256;
abyHeader[31] = nYSize / 256;
abyHeader[32] = nXSize % 256;
abyHeader[33] = nXSize / 256;
abyHeader[170] = 255; /* missing = 255 */
c2tp( 1.0, abyHeader + 171 ); /* slope = 1.0 */
c2tp( 0.0, abyHeader + 177 ); /* offset = 0 */
abyHeader[168] = 0; // lower limit
abyHeader[169] = 254; // upper limit
// pixel size = 1.0
c2tp( 1.0, abyHeader + 144 );
c2tp( 1.0, abyHeader + 150 );
if( VSIFWrite( abyHeader, 1, 512, fp ) != 512 )
{
CPLError( CE_Failure, CPLE_AppDefined,
"IO error writing %s.\n%s",
VSIStrerror( errno ) );
VSIFClose( fp );
return NULL;
}
/* -------------------------------------------------------------------- */
/* Now we need to extend the file to just the right number of */
/* bytes for the data we have to ensure it will open again */
/* properly. */
/* -------------------------------------------------------------------- */
if( VSIFSeek( fp, nXSize * nYSize - 1, SEEK_CUR ) != 0
|| VSIFWrite( abyHeader, 1, 1, fp ) != 1 )
{
CPLError( CE_Failure, CPLE_AppDefined,
"IO error writing %s.\n%s",
VSIStrerror( errno ) );
VSIFClose( fp );
return NULL;
}
VSIFClose( fp );
return (GDALDataset *) GDALOpen( pszFilename, GA_Update );
}
/************************************************************************/
/* GDALRegister_IDA() */
/************************************************************************/
void GDALRegister_IDA()
{
GDALDriver *poDriver;
if( GDALGetDriverByName( "IDA" ) == NULL )
{
poDriver = new GDALDriver();
poDriver->SetDescription( "IDA" );
poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
"Image Data and Analysis" );
poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
"frmt_various.html#IDA" );
poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES, "Byte" );
poDriver->pfnOpen = IDADataset::Open;
poDriver->pfnCreate = IDADataset::Create;
GetGDALDriverManager()->RegisterDriver( poDriver );
}
}
