00001 /******************************************************************************
00002 * $Id: gdal_misc_cpp-source.html,v 1.13 2002/12/21 19:13:12 warmerda Exp $
00003 *
00004 * Project: GDAL Core
00005 * Purpose: Free standing functions for GDAL.
00006 * Author: Frank Warmerdam, warmerda@home.com
00007 *
00008 ******************************************************************************
00009 * Copyright (c) 1999, Frank Warmerdam
00010 *
00011 * Permission is hereby granted, free of charge, to any person obtaining a
00012 * copy of this software and associated documentation files (the "Software"),
00013 * to deal in the Software without restriction, including without limitation
00014 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
00015 * and/or sell copies of the Software, and to permit persons to whom the
00016 * Software is furnished to do so, subject to the following conditions:
00017 *
00018 * The above copyright notice and this permission notice shall be included
00019 * in all copies or substantial portions of the Software.
00020 *
00021 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
00022 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
00023 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
00024 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
00025 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
00026 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
00027 * DEALINGS IN THE SOFTWARE.
00028 ******************************************************************************
00029 *
00030 * $Log: gdal_misc_cpp-source.html,v $
00030 * Revision 1.13 2002/12/21 19:13:12 warmerda
00030 * updated
00030 *
00031 * Revision 1.40 2002/12/11 21:21:46 warmerda
00032 * fixed debug format problem
00033 *
00034 * Revision 1.39 2002/12/09 20:05:31 warmerda
00035 * fixed return flag from GDALReadTabFile
00036 *
00037 * Revision 1.38 2002/12/09 18:53:25 warmerda
00038 * GDALDecToDMS() now calls CPLDecToDMS()
00039 *
00040 * Revision 1.37 2002/12/05 17:55:30 warmerda
00041 * gdalreadtabfile should not be static
00042 *
00043 * Revision 1.36 2002/12/05 15:46:38 warmerda
00044 * added GDALReadTabFile()
00045 *
00046 * Revision 1.35 2002/07/09 20:33:12 warmerda
00047 * expand tabs
00048 *
00049 * Revision 1.34 2002/05/06 21:37:29 warmerda
00050 * added GDALGCPsToGeoTransform
00051 *
00052 * Revision 1.33 2002/04/25 16:18:41 warmerda
00053 * added extra checking
00054 *
00055 * Revision 1.32 2002/04/24 19:21:26 warmerda
00056 * Include <ctype.h> for toupper(), tolower().
00057 *
00058 * Revision 1.31 2002/04/19 12:22:05 dron
00059 * added GDALWriteWorldFile()
00060 *
00061 * Revision 1.30 2002/04/16 13:59:33 warmerda
00062 * added GDALVersionInfo
00063 *
00064 * Revision 1.29 2001/12/07 20:04:21 warmerda
00065 * fixed serious bug in random sampler
00066 *
00067 * Revision 1.28 2001/11/30 03:41:26 warmerda
00068 * Fixed bug with the block sampling rate being too low to satisfy large
00069 * sample count values. Fixed bug with tiled images including some uninitialized
00070 * or zero data in the sample set on partial edge tiles.
00071 *
00072 * Revision 1.27 2001/11/26 20:14:01 warmerda
00073 * added GDALProjDef stubs for old 'bridges'
00074 *
00075 * Revision 1.26 2001/11/19 16:03:16 warmerda
00076 * moved GDALDectoDMS here
00077 *
00078 * Revision 1.25 2001/08/15 15:05:44 warmerda
00079 * return magnitude for complex samples in random sampler
00080 *
00081 * Revision 1.24 2001/07/18 04:04:30 warmerda
00082 * added CPL_CVSID
00083 *
00084 * Revision 1.23 2001/05/01 18:09:25 warmerda
00085 * added GDALReadWorldFile()
00086 *
00087 * Revision 1.22 2000/12/04 20:45:14 warmerda
00088 * removed unused variable.
00089 *
00090 * Revision 1.21 2000/10/06 15:22:49 warmerda
00091 * added GDALDataTypeUnion
00092 *
00093 * Revision 1.20 2000/08/18 15:24:48 warmerda
00094 * added GDALTermProgress
00095 *
00096 * Revision 1.19 2000/08/09 16:25:42 warmerda
00097 * don't crash if block is null
00098 *
00099 * Revision 1.18 2000/07/11 14:35:43 warmerda
00100 * added documentation
00101 *
00102 * Revision 1.17 2000/07/05 17:53:33 warmerda
00103 * Removed unused code related to nXCheck.
00104 *
00105 * Revision 1.16 2000/06/27 17:21:26 warmerda
00106 * added GDALGetRasterSampleOverview
00107 *
00108 * Revision 1.15 2000/06/26 22:17:49 warmerda
00109 * added scaled progress support
00110 *
00111 * Revision 1.14 2000/06/05 17:24:05 warmerda
00112 * added real complex support
00113 *
00114 * Revision 1.13 2000/04/21 21:55:32 warmerda
00115 * made more robust if block read fails
00116 *
00117 * Revision 1.12 2000/04/17 20:59:40 warmerda
00118 * Removed printf.
00119 *
00120 * Revision 1.11 2000/04/17 20:59:14 warmerda
00121 * fixed sampling bug
00122 *
00123 * Revision 1.10 2000/03/31 13:41:45 warmerda
00124 * added gcp support functions
00125 *
00126 * Revision 1.9 2000/03/24 00:09:19 warmerda
00127 * added sort-of random sampling
00128 *
00129 * Revision 1.8 2000/03/23 16:53:21 warmerda
00130 * use overviews for approximate min/max
00131 *
00132 * Revision 1.7 2000/03/09 23:21:44 warmerda
00133 * added GDALDummyProgress
00134 *
00135 * Revision 1.6 2000/03/06 21:59:44 warmerda
00136 * added min/max calculate
00137 *
00138 * Revision 1.5 2000/03/06 02:20:15 warmerda
00139 * added getname functions for colour interpretations
00140 *
00141 * Revision 1.4 1999/07/23 19:35:47 warmerda
00142 * added GDALGetDataTypeName
00143 *
00144 * Revision 1.3 1999/05/17 02:00:45 vgough
00145 * made pure_virtual C linkage
00146 *
00147 * Revision 1.2 1999/05/16 19:32:13 warmerda
00148 * Added __pure_virtual.
00149 *
00150 * Revision 1.1 1998/12/06 02:50:16 warmerda
00151 * New
00152 *
00153 */
00154
00155 #include "gdal_priv.h"
00156 #include "cpl_string.h"
00157 #include <ctype.h>
00158
00159 CPL_CVSID("$Id: gdal_misc_cpp-source.html,v 1.13 2002/12/21 19:13:12 warmerda Exp $");
00160
00161 #ifdef HAVE_MITAB
00162 #include "ogr_spatialref.h"
00163 // from mitab component.
00164 OGRSpatialReference * MITABCoordSys2SpatialRef( const char * pszCoordSys );
00165 #endif
00166
00167 /************************************************************************/
00168 /* __pure_virtual() */
00169 /* */
00170 /* The following is a gross hack to remove the last remaining */
00171 /* dependency on the GNU C++ standard library. */
00172 /************************************************************************/
00173
00174 #ifdef __GNUC__
00175
00176 extern "C"
00177 void __pure_virtual()
00178
00179 {
00180 }
00181
00182 #endif
00183
00184 /************************************************************************/
00185 /* GDALDataTypeUnion() */
00186 /************************************************************************/
00187
00198 GDALDataType GDALDataTypeUnion( GDALDataType eType1, GDALDataType eType2 )
00199
00200 {
00201 int bFloating, bComplex, nBits, bSigned;
00202
00203 bComplex = GDALDataTypeIsComplex(eType1) | GDALDataTypeIsComplex(eType2);
00204
00205 switch( eType1 )
00206 {
00207 case GDT_Byte:
00208 nBits = 8;
00209 bSigned = FALSE;
00210 bFloating = FALSE;
00211 break;
00212
00213 case GDT_Int16:
00214 case GDT_CInt16:
00215 nBits = 16;
00216 bSigned = TRUE;
00217 bFloating = FALSE;
00218 break;
00219
00220 case GDT_UInt16:
00221 nBits = 16;
00222 bSigned = FALSE;
00223 bFloating = FALSE;
00224 break;
00225
00226 case GDT_Int32:
00227 case GDT_CInt32:
00228 nBits = 32;
00229 bSigned = TRUE;
00230 bFloating = FALSE;
00231 break;
00232
00233 case GDT_UInt32:
00234 nBits = 32;
00235 bSigned = FALSE;
00236 bFloating = FALSE;
00237 break;
00238
00239 case GDT_Float32:
00240 case GDT_CFloat32:
00241 nBits = 32;
00242 bSigned = TRUE;
00243 bFloating = TRUE;
00244 break;
00245
00246 case GDT_Float64:
00247 case GDT_CFloat64:
00248 nBits = 64;
00249 bSigned = TRUE;
00250 bFloating = TRUE;
00251 break;
00252
00253 default:
00254 CPLAssert( FALSE );
00255 return GDT_Unknown;
00256 }
00257
00258 switch( eType2 )
00259 {
00260 case GDT_Byte:
00261 break;
00262
00263 case GDT_Int16:
00264 nBits = MAX(nBits,16);
00265 bSigned = TRUE;
00266 break;
00267
00268 case GDT_UInt16:
00269 nBits = MAX(nBits,16);
00270 break;
00271
00272 case GDT_Int32:
00273 case GDT_CInt32:
00274 nBits = MAX(nBits,32);
00275 bSigned = TRUE;
00276 break;
00277
00278 case GDT_UInt32:
00279 nBits = MAX(nBits,32);
00280 break;
00281
00282 case GDT_Float32:
00283 case GDT_CFloat32:
00284 nBits = MAX(nBits,32);
00285 bSigned = TRUE;
00286 bFloating = TRUE;
00287 break;
00288
00289 case GDT_Float64:
00290 case GDT_CFloat64:
00291 nBits = MAX(nBits,64);
00292 bSigned = TRUE;
00293 bFloating = TRUE;
00294 break;
00295
00296 default:
00297 CPLAssert( FALSE );
00298 return GDT_Unknown;
00299 }
00300
00301 if( nBits == 8 )
00302 return GDT_Byte;
00303 else if( nBits == 16 && bComplex )
00304 return GDT_CInt16;
00305 else if( nBits == 16 && bSigned )
00306 return GDT_Int16;
00307 else if( nBits == 16 && !bSigned )
00308 return GDT_UInt16;
00309 else if( nBits == 32 && bFloating && bComplex )
00310 return GDT_CFloat32;
00311 else if( nBits == 32 && bFloating )
00312 return GDT_Float32;
00313 else if( nBits == 32 && bComplex )
00314 return GDT_CInt32;
00315 else if( nBits == 32 && bSigned )
00316 return GDT_Int32;
00317 else if( nBits == 32 && !bSigned )
00318 return GDT_UInt32;
00319 else if( nBits == 64 && bComplex )
00320 return GDT_CFloat64;
00321 else
00322 return GDT_Float64;
00323 }
00324
00325
00326 /************************************************************************/
00327 /* GDALGetDataTypeSize() */
00328 /************************************************************************/
00329 int GDALGetDataTypeSize( GDALDataType eDataType )
00330
00331 {
00332 switch( eDataType )
00333 {
00334 case GDT_Byte:
00335 return 8;
00336
00337 case GDT_UInt16:
00338 case GDT_Int16:
00339 return 16;
00340
00341 case GDT_UInt32:
00342 case GDT_Int32:
00343 case GDT_Float32:
00344 case GDT_CInt16:
00345 return 32;
00346
00347 case GDT_Float64:
00348 case GDT_CInt32:
00349 case GDT_CFloat32:
00350 return 64;
00351
00352 case GDT_CFloat64:
00353 return 128;
00354
00355 default:
00356 CPLAssert( FALSE );
00357 return 0;
00358 }
00359 }
00360
00361 /************************************************************************/
00362 /* GDALDataTypeIsComplex() */
00363 /************************************************************************/
00364
00365 int GDALDataTypeIsComplex( GDALDataType eDataType )
00366
00367 {
00368 switch( eDataType )
00369 {
00370 case GDT_CInt16:
00371 case GDT_CInt32:
00372 case GDT_CFloat32:
00373 case GDT_CFloat64:
00374 return TRUE;
00375
00376 default:
00377 return FALSE;
00378 }
00379 }
00380
00381 /************************************************************************/
00382 /* GDALGetDataTypeName() */
00383 /************************************************************************/
00384
00385 const char *GDALGetDataTypeName( GDALDataType eDataType )
00386
00387 {
00388 switch( eDataType )
00389 {
00390 case GDT_Byte:
00391 return "Byte";
00392
00393 case GDT_UInt16:
00394 return "UInt16";
00395
00396 case GDT_Int16:
00397 return "Int16";
00398
00399 case GDT_UInt32:
00400 return "UInt32";
00401
00402 case GDT_Int32:
00403 return "Int32";
00404
00405 case GDT_Float32:
00406 return "Float32";
00407
00408 case GDT_Float64:
00409 return "Float64";
00410
00411 case GDT_CInt16:
00412 return "CInt16";
00413
00414 case GDT_CInt32:
00415 return "CInt32";
00416
00417 case GDT_CFloat32:
00418 return "CFloat32";
00419
00420 case GDT_CFloat64:
00421 return "CFloat64";
00422
00423 default:
00424 return NULL;
00425 }
00426 }
00427
00428 /************************************************************************/
00429 /* GDALGetPaletteInterpretationName() */
00430 /************************************************************************/
00431
00432 const char *GDALGetPaletteInterpretationName( GDALPaletteInterp eInterp )
00433
00434 {
00435 switch( eInterp )
00436 {
00437 case GPI_Gray:
00438 return "Gray";
00439
00440 case GPI_RGB:
00441 return "RGB";
00442
00443 case GPI_CMYK:
00444 return "CMYK";
00445
00446 case GPI_HLS:
00447 return "HLS";
00448
00449 default:
00450 return "Unknown";
00451 }
00452 }
00453
00454 /************************************************************************/
00455 /* GDALGetColorInterpretationName() */
00456 /************************************************************************/
00457
00458 const char *GDALGetColorInterpretationName( GDALColorInterp eInterp )
00459
00460 {
00461 switch( eInterp )
00462 {
00463 case GCI_Undefined:
00464 return "Undefined";
00465
00466 case GCI_GrayIndex:
00467 return "Gray";
00468
00469 case GCI_PaletteIndex:
00470 return "Palette";
00471
00472 case GCI_RedBand:
00473 return "Red";
00474
00475 case GCI_GreenBand:
00476 return "Green";
00477
00478 case GCI_BlueBand:
00479 return "Blue";
00480
00481 case GCI_AlphaBand:
00482 return "Alpha";
00483
00484 case GCI_HueBand:
00485 return "Hue";
00486
00487 case GCI_SaturationBand:
00488 return "Saturation";
00489
00490 case GCI_LightnessBand:
00491 return "Lightness";
00492
00493 case GCI_CyanBand:
00494 return "Cyan";
00495
00496 case GCI_MagentaBand:
00497 return "Magenta";
00498
00499 case GCI_YellowBand:
00500 return "Yellow";
00501
00502 case GCI_BlackBand:
00503 return "Black";
00504
00505 default:
00506 return "Unknown";
00507 }
00508 }
00509
00510 /************************************************************************/
00511 /* GDALComputeRasterMinMax() */
00512 /************************************************************************/
00513
00532 void GDALComputeRasterMinMax( GDALRasterBandH hBand, int bApproxOK,
00533 double adfMinMax[2] )
00534
00535 {
00536 double dfMin=0.0, dfMax=0.0;
00537 GDALRasterBand *poBand;
00538
00539 /* -------------------------------------------------------------------- */
00540 /* Does the driver already know the min/max? */
00541 /* -------------------------------------------------------------------- */
00542 if( bApproxOK )
00543 {
00544 int bSuccessMin, bSuccessMax;
00545
00546 dfMin = GDALGetRasterMinimum( hBand, &bSuccessMin );
00547 dfMax = GDALGetRasterMaximum( hBand, &bSuccessMax );
00548
00549 if( bSuccessMin && bSuccessMax )
00550 {
00551 adfMinMax[0] = dfMin;
00552 adfMinMax[1] = dfMax;
00553 return;
00554 }
00555 }
00556
00557 /* -------------------------------------------------------------------- */
00558 /* If we have overview bands, use them for min/max. */
00559 /* -------------------------------------------------------------------- */
00560 if( bApproxOK )
00561 poBand = (GDALRasterBand *) GDALGetRasterSampleOverview( hBand, 2500 );
00562 else
00563 poBand = (GDALRasterBand *) hBand;
00564
00565 /* -------------------------------------------------------------------- */
00566 /* Figure out the ratio of blocks we will read to get an */
00567 /* approximate value. */
00568 /* -------------------------------------------------------------------- */
00569 int nBlockXSize, nBlockYSize;
00570 int nBlocksPerRow, nBlocksPerColumn;
00571 int nSampleRate;
00572 int bGotNoDataValue, bFirstValue = TRUE;
00573 double dfNoDataValue;
00574
00575 dfNoDataValue = poBand->GetNoDataValue( &bGotNoDataValue );
00576
00577 poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
00578 nBlocksPerRow = (poBand->GetXSize() + nBlockXSize - 1) / nBlockXSize;
00579 nBlocksPerColumn = (poBand->GetYSize() + nBlockYSize - 1) / nBlockYSize;
00580
00581 if( bApproxOK )
00582 nSampleRate =
00583 (int) MAX(1,sqrt((double) nBlocksPerRow * nBlocksPerColumn));
00584 else
00585 nSampleRate = 1;
00586
00587 for( int iSampleBlock = 0;
00588 iSampleBlock < nBlocksPerRow * nBlocksPerColumn;
00589 iSampleBlock += nSampleRate )
00590 {
00591 double dfValue = 0.0;
00592 int iXBlock, iYBlock, nXCheck, nYCheck;
00593 GDALRasterBlock *poBlock;
00594
00595 iYBlock = iSampleBlock / nBlocksPerRow;
00596 iXBlock = iSampleBlock - nBlocksPerRow * iYBlock;
00597
00598 poBlock = poBand->GetBlockRef( iXBlock, iYBlock );
00599 if( poBlock == NULL )
00600 continue;
00601
00602 if( (iXBlock+1) * nBlockXSize > poBand->GetXSize() )
00603 nXCheck = poBand->GetXSize() - iXBlock * nBlockXSize;
00604 else
00605 nXCheck = nBlockXSize;
00606
00607 if( (iYBlock+1) * nBlockYSize > poBand->GetYSize() )
00608 nYCheck = poBand->GetYSize() - iYBlock * nBlockYSize;
00609 else
00610 nYCheck = nBlockYSize;
00611
00612 /* this isn't the fastest way to do this, but is easier for now */
00613 for( int iY = 0; iY < nYCheck; iY++ )
00614 {
00615 for( int iX = 0; iX < nXCheck; iX++ )
00616 {
00617 int iOffset = iX + iY * nBlockXSize;
00618
00619 switch( poBlock->GetDataType() )
00620 {
00621 case GDT_Byte:
00622 dfValue = ((GByte *) poBlock->GetDataRef())[iOffset];
00623 break;
00624 case GDT_UInt16:
00625 dfValue = ((GUInt16 *) poBlock->GetDataRef())[iOffset];
00626 break;
00627 case GDT_Int16:
00628 dfValue = ((GInt16 *) poBlock->GetDataRef())[iOffset];
00629 break;
00630 case GDT_UInt32:
00631 dfValue = ((GUInt32 *) poBlock->GetDataRef())[iOffset];
00632 break;
00633 case GDT_Int32:
00634 dfValue = ((GInt32 *) poBlock->GetDataRef())[iOffset];
00635 break;
00636 case GDT_Float32:
00637 dfValue = ((float *) poBlock->GetDataRef())[iOffset];
00638 break;
00639 case GDT_Float64:
00640 dfValue = ((double *) poBlock->GetDataRef())[iOffset];
00641 break;
00642 case GDT_CInt16:
00643 dfValue = ((GInt16 *) poBlock->GetDataRef())[iOffset*2];
00644 break;
00645 case GDT_CInt32:
00646 dfValue = ((GInt32 *) poBlock->GetDataRef())[iOffset*2];
00647 break;
00648 case GDT_CFloat32:
00649 dfValue = ((float *) poBlock->GetDataRef())[iOffset*2];
00650 break;
00651 case GDT_CFloat64:
00652 dfValue = ((double *) poBlock->GetDataRef())[iOffset*2];
00653 break;
00654 default:
00655 CPLAssert( FALSE );
00656 }
00657
00658 if( bGotNoDataValue && dfValue == dfNoDataValue )
00659 continue;
00660
00661 if( bFirstValue )
00662 {
00663 dfMin = dfMax = dfValue;
00664 bFirstValue = FALSE;
00665 }
00666 else
00667 {
00668 dfMin = MIN(dfMin,dfValue);
00669 dfMax = MAX(dfMax,dfValue);
00670 }
00671 }
00672 }
00673 }
00674
00675 adfMinMax[0] = dfMin;
00676 adfMinMax[1] = dfMax;
00677 }
00678
00679 /************************************************************************/
00680 /* GDALDummyProgress() */
00681 /************************************************************************/
00682
00747 int GDALDummyProgress( double, const char *, void * )
00748
00749 {
00750 return TRUE;
00751 }
00752
00753 /************************************************************************/
00754 /* GDALScaledProgress() */
00755 /************************************************************************/
00756 typedef struct {
00757 GDALProgressFunc pfnProgress;
00758 void *pData;
00759 double dfMin;
00760 double dfMax;
00761 } GDALScaledProgressInfo;
00762
00763 int GDALScaledProgress( double dfComplete, const char *pszMessage,
00764 void *pData )
00765
00766 {
00767 GDALScaledProgressInfo *psInfo = (GDALScaledProgressInfo *) pData;
00768
00769 return psInfo->pfnProgress( dfComplete * (psInfo->dfMax - psInfo->dfMin)
00770 + psInfo->dfMin,
00771 pszMessage, psInfo->pData );
00772 }
00773
00774 /************************************************************************/
00775 /* GDALCreateScaledProgress() */
00776 /************************************************************************/
00777
00778 void *GDALCreateScaledProgress( double dfMin, double dfMax,
00779 GDALProgressFunc pfnProgress,
00780 void * pData )
00781
00782 {
00783 GDALScaledProgressInfo *psInfo;
00784
00785 psInfo = (GDALScaledProgressInfo *)
00786 CPLCalloc(sizeof(GDALScaledProgressInfo),1);
00787
00788 if( ABS(dfMin-dfMax) < 0.0000001 )
00789 dfMax = dfMin + 0.01;
00790
00791 psInfo->pData = pData;
00792 psInfo->pfnProgress = pfnProgress;
00793 psInfo->dfMin = dfMin;
00794 psInfo->dfMax = dfMax;
00795
00796 return (void *) psInfo;
00797 }
00798
00799 /************************************************************************/
00800 /* GDALDestroyScaledProgress() */
00801 /************************************************************************/
00802
00803 void GDALDestroyScaledProgress( void * pData )
00804
00805 {
00806 CPLFree( pData );
00807 }
00808
00809 /************************************************************************/
00810 /* GDALTermProgress() */
00811 /************************************************************************/
00812
00813 int GDALTermProgress( double dfComplete, const char *pszMessage, void * )
00814
00815 {
00816 static double dfLastComplete = -1.0;
00817
00818 if( dfLastComplete > dfComplete )
00819 {
00820 if( dfLastComplete > 1.0 )
00821 dfLastComplete = -1.0;
00822 else
00823 dfLastComplete = dfComplete;
00824 }
00825
00826 if( floor(dfLastComplete*10) != floor(dfComplete*10) )
00827 {
00828 int nPercent = (int) floor(dfComplete*100);
00829
00830 if( nPercent == 0 && pszMessage != NULL )
00831 fprintf( stdout, "%s:", pszMessage );
00832
00833 if( nPercent == 100 )
00834 fprintf( stdout, "%d - done.\n", (int) floor(dfComplete*100) );
00835 else
00836 {
00837 fprintf( stdout, "%d.", (int) floor(dfComplete*100) );
00838 fflush( stdout );
00839 }
00840 }
00841 else if( floor(dfLastComplete*30) != floor(dfComplete*30) )
00842 {
00843 fprintf( stdout, "." );
00844 fflush( stdout );
00845 }
00846
00847 dfLastComplete = dfComplete;
00848
00849 return TRUE;
00850 }
00851
00852 /************************************************************************/
00853 /* GDALGetRasterSampleOverview() */
00854 /************************************************************************/
00855
00871 GDALRasterBandH GDALGetRasterSampleOverview( GDALRasterBandH hBand,
00872 int nDesiredSamples )
00873
00874 {
00875 int nBestSamples;
00876 GDALRasterBandH hBestBand = hBand;
00877
00878 nBestSamples = GDALGetRasterBandXSize(hBand)
00879 * GDALGetRasterBandYSize(hBand);
00880
00881 for( int iOverview = 0;
00882 iOverview < GDALGetOverviewCount( hBand );
00883 iOverview++ )
00884 {
00885 GDALRasterBandH hOBand = GDALGetOverview( hBand, iOverview );
00886 int nOSamples;
00887
00888 nOSamples = GDALGetRasterBandXSize(hOBand)
00889 * GDALGetRasterBandYSize(hOBand);
00890
00891 if( nOSamples < nBestSamples && nOSamples > nDesiredSamples )
00892 {
00893 nBestSamples = nOSamples;
00894 hBestBand = hOBand;
00895 }
00896 }
00897
00898 return hBestBand;
00899 }
00900
00901 /************************************************************************/
00902 /* GDALGetRandomRasterSample() */
00903 /************************************************************************/
00904
00905 int GDALGetRandomRasterSample( GDALRasterBandH hBand, int nSamples,
00906 float *pafSampleBuf )
00907
00908 {
00909 GDALRasterBand *poBand;
00910
00911 poBand = (GDALRasterBand *) GDALGetRasterSampleOverview( hBand, nSamples );
00912
00913 /* -------------------------------------------------------------------- */
00914 /* Figure out the ratio of blocks we will read to get an */
00915 /* approximate value. */
00916 /* -------------------------------------------------------------------- */
00917 int nBlockXSize, nBlockYSize;
00918 int nBlocksPerRow, nBlocksPerColumn;
00919 int nSampleRate;
00920 int bGotNoDataValue;
00921 double dfNoDataValue;
00922 int nActualSamples = 0;
00923 int nBlockSampleRate;
00924 int nBlockPixels, nBlockCount;
00925
00926 dfNoDataValue = poBand->GetNoDataValue( &bGotNoDataValue );
00927
00928 poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
00929
00930 nBlocksPerRow = (poBand->GetXSize() + nBlockXSize - 1) / nBlockXSize;
00931 nBlocksPerColumn = (poBand->GetYSize() + nBlockYSize - 1) / nBlockYSize;
00932
00933 nBlockPixels = nBlockXSize * nBlockYSize;
00934 nBlockCount = nBlocksPerRow * nBlocksPerColumn;
00935
00936 if( nBlocksPerRow == 0 || nBlocksPerColumn == 0 || nBlockPixels == 0
00937 || nBlockCount == 0 )
00938 {
00939 CPLError( CE_Failure, CPLE_AppDefined,
00940 "GDALGetRandomSample(): returning because band"
00941 " appears degenerate." );
00942
00943 return FALSE;
00944 }
00945
00946 nSampleRate = (int) MAX(1,sqrt((double) nBlockCount)-2.0);
00947
00948 if( nSampleRate == nBlocksPerRow && nSampleRate > 1 )
00949 nSampleRate--;
00950
00951 while( nSampleRate > 1
00952 && ((nBlockCount-1) / nSampleRate + 1) * nBlockPixels < nSamples )
00953 nSampleRate--;
00954
00955 nBlockSampleRate =
00956 MAX(1,nBlockPixels / (nSamples / ((nBlockCount-1) / nSampleRate + 1)));
00957
00958 for( int iSampleBlock = 0;
00959 iSampleBlock < nBlockCount;
00960 iSampleBlock += nSampleRate )
00961 {
00962 double dfValue = 0.0, dfReal, dfImag;
00963 int iXBlock, iYBlock, iX, iY, iXValid, iYValid, iRemainder = 0;
00964 GDALRasterBlock *poBlock;
00965
00966 iYBlock = iSampleBlock / nBlocksPerRow;
00967 iXBlock = iSampleBlock - nBlocksPerRow * iYBlock;
00968
00969 poBlock = poBand->GetBlockRef( iXBlock, iYBlock );
00970 if( poBlock == NULL )
00971 continue;
00972
00973 if( iXBlock * nBlockXSize > poBand->GetXSize() )
00974 iXValid = poBand->GetXSize() - iXBlock * nBlockXSize;
00975 else
00976 iXValid = nBlockXSize;
00977
00978 if( iYBlock * nBlockYSize > poBand->GetYSize() )
00979 iYValid = poBand->GetYSize() - iYBlock * nBlockYSize;
00980 else
00981 iYValid = nBlockYSize;
00982
00983 for( iY = 0; iY < iYValid; iY++ )
00984 {
00985 for( iX = iRemainder; iX < iXValid; iX += nBlockSampleRate )
00986 {
00987 int iOffset;
00988
00989 iOffset = iX + iY * nBlockXSize;
00990 switch( poBlock->GetDataType() )
00991 {
00992 case GDT_Byte:
00993 dfValue = ((GByte *) poBlock->GetDataRef())[iOffset];
00994 break;
00995 case GDT_UInt16:
00996 dfValue = ((GUInt16 *) poBlock->GetDataRef())[iOffset];
00997 break;
00998 case GDT_Int16:
00999 dfValue = ((GInt16 *) poBlock->GetDataRef())[iOffset];
01000 break;
01001 case GDT_UInt32:
01002 dfValue = ((GUInt32 *) poBlock->GetDataRef())[iOffset];
01003 break;
01004 case GDT_Int32:
01005 dfValue = ((GInt32 *) poBlock->GetDataRef())[iOffset];
01006 break;
01007 case GDT_Float32:
01008 dfValue = ((float *) poBlock->GetDataRef())[iOffset];
01009 break;
01010 case GDT_Float64:
01011 dfValue = ((double *) poBlock->GetDataRef())[iOffset];
01012 break;
01013 case GDT_CInt16:
01014 dfReal = ((GInt16 *) poBlock->GetDataRef())[iOffset*2];
01015 dfImag = ((GInt16 *) poBlock->GetDataRef())[iOffset*2+1];
01016 dfValue = sqrt(dfReal*dfReal + dfImag*dfImag);
01017 break;
01018 case GDT_CInt32:
01019 dfReal = ((GInt32 *) poBlock->GetDataRef())[iOffset*2];
01020 dfImag = ((GInt32 *) poBlock->GetDataRef())[iOffset*2+1];
01021 dfValue = sqrt(dfReal*dfReal + dfImag*dfImag);
01022 break;
01023 case GDT_CFloat32:
01024 dfReal = ((float *) poBlock->GetDataRef())[iOffset*2];
01025 dfImag = ((float *) poBlock->GetDataRef())[iOffset*2+1];
01026 dfValue = sqrt(dfReal*dfReal + dfImag*dfImag);
01027 break;
01028 case GDT_CFloat64:
01029 dfReal = ((double *) poBlock->GetDataRef())[iOffset*2];
01030 dfImag = ((double *) poBlock->GetDataRef())[iOffset*2+1];
01031 dfValue = sqrt(dfReal*dfReal + dfImag*dfImag);
01032 break;
01033 default:
01034 CPLAssert( FALSE );
01035 }
01036
01037 if( bGotNoDataValue && dfValue == dfNoDataValue )
01038 continue;
01039
01040 if( nActualSamples < nSamples )
01041 pafSampleBuf[nActualSamples++] = dfValue;
01042 }
01043
01044 iRemainder = iX - iXValid;
01045 }
01046 }
01047
01048 return nActualSamples;
01049 }
01050
01051 /************************************************************************/
01052 /* GDALInitGCPs() */
01053 /************************************************************************/
01054
01055 void GDALInitGCPs( int nCount, GDAL_GCP * psGCP )
01056
01057 {
01058 for( int iGCP = 0; iGCP < nCount; iGCP++ )
01059 {
01060 memset( psGCP, 0, sizeof(GDAL_GCP) );
01061 psGCP->pszId = CPLStrdup("");
01062 psGCP->pszInfo = CPLStrdup("");
01063 psGCP++;
01064 }
01065 }
01066
01067 /************************************************************************/
01068 /* GDALDeinitGCPs() */
01069 /************************************************************************/
01070
01071 void GDALDeinitGCPs( int nCount, GDAL_GCP * psGCP )
01072
01073 {
01074 for( int iGCP = 0; iGCP < nCount; iGCP++ )
01075 {
01076 CPLFree( psGCP->pszId );
01077 CPLFree( psGCP->pszInfo );
01078 psGCP++;
01079 }
01080 }
01081
01082 /************************************************************************/
01083 /* GDALDuplicateGCPs() */
01084 /************************************************************************/
01085
01086 GDAL_GCP *GDALDuplicateGCPs( int nCount, const GDAL_GCP *pasGCPList )
01087
01088 {
01089 GDAL_GCP *pasReturn;
01090
01091 pasReturn = (GDAL_GCP *) CPLMalloc(sizeof(GDAL_GCP) * nCount);
01092 GDALInitGCPs( nCount, pasReturn );
01093
01094 for( int iGCP = 0; iGCP < nCount; iGCP++ )
01095 {
01096 CPLFree( pasReturn[iGCP].pszId );
01097 pasReturn[iGCP].pszId = CPLStrdup( pasGCPList[iGCP].pszId );
01098
01099 CPLFree( pasReturn[iGCP].pszInfo );
01100 pasReturn[iGCP].pszInfo = CPLStrdup( pasGCPList[iGCP].pszInfo );
01101
01102 pasReturn[iGCP].dfGCPPixel = pasGCPList[iGCP].dfGCPPixel;
01103 pasReturn[iGCP].dfGCPLine = pasGCPList[iGCP].dfGCPLine;
01104 pasReturn[iGCP].dfGCPX = pasGCPList[iGCP].dfGCPX;
01105 pasReturn[iGCP].dfGCPY = pasGCPList[iGCP].dfGCPY;
01106 pasReturn[iGCP].dfGCPZ = pasGCPList[iGCP].dfGCPZ;
01107 }
01108
01109 return pasReturn;
01110 }
01111
01112 /************************************************************************/
01113 /* GDALReadTabFile() */
01114 /* */
01115 /* Helper function for translator implementators wanting */
01116 /* support for MapInfo .tab-files. */
01117 /************************************************************************/
01118
01119 #define MAX_GCP 256
01120
01121 int GDALReadTabFile( const char * pszBaseFilename,
01122 double *padfGeoTransform, char **ppszWKT,
01123 int *pnGCPCount, GDAL_GCP **ppasGCPs )
01124
01125
01126 {
01127 const char *pszTAB;
01128 FILE *fpTAB;
01129 char **papszLines;
01130 char **papszTok=NULL;
01131 int bTypeRasterFound = FALSE;
01132 int bInsideTableDef = FALSE;
01133 int iLine, numLines=0;
01134 int nCoordinateCount = 0;
01135 GDAL_GCP asGCPs[MAX_GCP];
01136
01137
01138 /* -------------------------------------------------------------------- */
01139 /* Try lower case, then upper case. */
01140 /* -------------------------------------------------------------------- */
01141 pszTAB = CPLResetExtension( pszBaseFilename, "tab" );
01142
01143 fpTAB = VSIFOpen( pszTAB, "rt" );
01144
01145 #ifndef WIN32
01146 if( fpTAB == NULL )
01147 {
01148 pszTAB = CPLResetExtension( pszBaseFilename, "TAB" );
01149 fpTAB = VSIFOpen( pszTAB, "rt" );
01150 }
01151 #endif
01152
01153 if( fpTAB == NULL )
01154 return FALSE;
01155
01156 VSIFClose( fpTAB );
01157
01158 /* -------------------------------------------------------------------- */
01159 /* We found the file, now load and parse it. */
01160 /* -------------------------------------------------------------------- */
01161 papszLines = CSLLoad( pszTAB );
01162
01163 numLines = CSLCount(papszLines);
01164
01165 // Iterate all lines in the TAB-file
01166 for(iLine=0; iLine<numLines; iLine++)
01167 {
01168 CSLDestroy(papszTok);
01169 papszTok = CSLTokenizeStringComplex(papszLines[iLine], " \t(),;",
01170 TRUE, FALSE);
01171
01172 if (CSLCount(papszTok) < 2)
01173 continue;
01174
01175 // Did we find table definition
01176 if (EQUAL(papszTok[0], "Definition") && EQUAL(papszTok[1], "Table") )
01177 {
01178 bInsideTableDef = TRUE;
01179 }
01180 else if (bInsideTableDef && (EQUAL(papszTok[0], "Type")) )
01181 {
01182 // Only RASTER-type will be handled
01183 if (EQUAL(papszTok[1], "RASTER"))
01184 {
01185 bTypeRasterFound = TRUE;
01186 }
01187 else
01188 {
01189 CSLDestroy(papszTok);
01190 CSLDestroy(papszLines);
01191 return FALSE;
01192 }
01193 }
01194 else if (bTypeRasterFound && bInsideTableDef
01195 && CSLCount(papszTok) > 5
01196 && EQUAL(papszTok[4], "Label")
01197 && nCoordinateCount < MAX_GCP )
01198 {
01199 GDALInitGCPs( 1, asGCPs + nCoordinateCount );
01200
01201 asGCPs[nCoordinateCount].dfGCPPixel = atof(papszTok[2]);
01202 asGCPs[nCoordinateCount].dfGCPLine = atof(papszTok[3]);
01203 asGCPs[nCoordinateCount].dfGCPX = atof(papszTok[0]);
01204 asGCPs[nCoordinateCount].dfGCPY = atof(papszTok[1]);
01205 CPLFree( asGCPs[nCoordinateCount].pszId );
01206 asGCPs[nCoordinateCount].pszId = CPLStrdup(papszTok[5]);
01207
01208 nCoordinateCount++;
01209 }
01210 else if( bTypeRasterFound && bInsideTableDef
01211 && EQUAL(papszTok[0],"CoordSys")
01212 && ppszWKT != NULL )
01213 {
01214 #ifdef HAVE_MITAB
01215 OGRSpatialReference *poSRS = NULL;
01216
01217 poSRS = MITABCoordSys2SpatialRef( papszLines[iLine] );
01218 if( poSRS != NULL )
01219 {
01220 poSRS->exportToWkt( ppszWKT );
01221 delete poSRS;
01222 }
01223 #else
01224 CPLDebug( "GDAL", "GDALReadTabFile(): Found `%s',\n"
01225 "but GDALReadTabFile() not configured with MITAB callout.",
01226 papszLines[iLine] );
01227 #endif
01228 }
01229 }
01230
01231 /* -------------------------------------------------------------------- */
01232 /* Try to convert the GCPs into a geotransform definition, if */
01233 /* possible. Otherwise we will need to use them as GCPs. */
01234 /* -------------------------------------------------------------------- */
01235 if( !GDALGCPsToGeoTransform( nCoordinateCount, asGCPs, padfGeoTransform,
01236 FALSE ) )
01237 {
01238 CPLDebug( "GDAL",
01239 "GDALReadTabFile(%s) found file, wasn't able to derive a\n"
01240 "first order geotransform. Using points as GCPs.",
01241 pszTAB );
01242
01243 *ppasGCPs = (GDAL_GCP *)
01244 CPLCalloc(sizeof(GDAL_GCP),nCoordinateCount);
01245 memcpy( *ppasGCPs, asGCPs, sizeof(GDAL_GCP) * nCoordinateCount );
01246 *pnGCPCount = nCoordinateCount;
01247 }
01248 else
01249 {
01250 GDALDeinitGCPs( nCoordinateCount, asGCPs );
01251 }
01252
01253 CSLDestroy(papszTok);
01254 CSLDestroy(papszLines);
01255 return *pnGCPCount == 0;
01256 }
01257
01258 /************************************************************************/
01259 /* GDALReadWorldFile() */
01260 /* */
01261 /* Helper function for translator implementators wanting */
01262 /* support for ESRI world files. */
01263 /************************************************************************/
01264
01265 int GDALReadWorldFile( const char * pszBaseFilename, const char *pszExtension,
01266 double *padfGeoTransform )
01267
01268 {
01269 const char *pszTFW;
01270 char szExtUpper[32], szExtLower[32];
01271 int i;
01272 FILE *fpTFW;
01273 char **papszLines;
01274
01275 if( *pszExtension == '.' )
01276 pszExtension++;
01277
01278 /* -------------------------------------------------------------------- */
01279 /* Generate upper and lower case versions of the extension. */
01280 /* -------------------------------------------------------------------- */
01281 strncpy( szExtUpper, pszExtension, 32 );
01282 strncpy( szExtLower, pszExtension, 32 );
01283
01284 for( i = 0; szExtUpper[i] != '\0' && i < 32; i++ )
01285 {
01286 szExtUpper[i] = toupper(szExtUpper[i]);
01287 szExtLower[i] = tolower(szExtLower[i]);
01288 }
01289
01290 /* -------------------------------------------------------------------- */
01291 /* Try lower case, then upper case. */
01292 /* -------------------------------------------------------------------- */
01293 pszTFW = CPLResetExtension( pszBaseFilename, szExtLower );
01294
01295 fpTFW = VSIFOpen( pszTFW, "rt" );
01296
01297 #ifndef WIN32
01298 if( fpTFW == NULL )
01299 {
01300 pszTFW = CPLResetExtension( pszBaseFilename, szExtUpper );
01301 fpTFW = VSIFOpen( pszTFW, "rt" );
01302 }
01303 #endif
01304
01305 if( fpTFW == NULL )
01306 return FALSE;
01307
01308 VSIFClose( fpTFW );
01309
01310 /* -------------------------------------------------------------------- */
01311 /* We found the file, now load and parse it. */
01312 /* -------------------------------------------------------------------- */
01313 papszLines = CSLLoad( pszTFW );
01314 if( CSLCount(papszLines) >= 6
01315 && atof(papszLines[0]) != 0.0
01316 && atof(papszLines[3]) != 0.0 )
01317 {
01318 padfGeoTransform[0] = atof(papszLines[4]);
01319 padfGeoTransform[1] = atof(papszLines[0]);
01320 padfGeoTransform[2] = atof(papszLines[2]);
01321 padfGeoTransform[3] = atof(papszLines[5]);
01322 padfGeoTransform[4] = atof(papszLines[1]);
01323 padfGeoTransform[5] = atof(papszLines[3]);
01324
01325 // correct for center of pixel vs. top left of pixel
01326 padfGeoTransform[0] -= 0.5 * padfGeoTransform[1];
01327 padfGeoTransform[0] -= 0.5 * padfGeoTransform[2];
01328 padfGeoTransform[3] -= 0.5 * padfGeoTransform[4];
01329 padfGeoTransform[3] -= 0.5 * padfGeoTransform[5];
01330
01331 CSLDestroy(papszLines);
01332
01333 return TRUE;
01334 }
01335 else
01336 {
01337 CPLDebug( "GDAL",
01338 "GDALReadWorldFile(%s) found file, but it was corrupt.",
01339 pszTFW );
01340 CSLDestroy(papszLines);
01341 return FALSE;
01342 }
01343 }
01344
01345 /************************************************************************/
01346 /* GDALWriteWorldFile() */
01347 /* */
01348 /* Helper function for translator implementators wanting */
01349 /* support for ESRI world files. */
01350 /************************************************************************/
01351
01352 int GDALWriteWorldFile( const char * pszBaseFilename, const char *pszExtension,
01353 double *padfGeoTransform )
01354
01355 {
01356 const char *pszTFW;
01357 FILE *fpTFW;
01358
01359 pszTFW = CPLResetExtension( pszBaseFilename, pszExtension );
01360 fpTFW = VSIFOpen( pszTFW, "wt" );
01361 if( fpTFW == NULL )
01362 return FALSE;
01363
01364 /* -------------------------------------------------------------------- */
01365 /* We open the file, now fill it with the world data. */
01366 /* -------------------------------------------------------------------- */
01367 fprintf( fpTFW, "%.10f\n", padfGeoTransform[1] );
01368 fprintf( fpTFW, "%.10f\n", padfGeoTransform[4] );
01369 fprintf( fpTFW, "%.10f\n", padfGeoTransform[2] );
01370 fprintf( fpTFW, "%.10f\n", padfGeoTransform[5] );
01371 fprintf( fpTFW, "%.10f\n", padfGeoTransform[0]
01372 + 0.5 * padfGeoTransform[1]
01373 + 0.5 * padfGeoTransform[2] );
01374 fprintf( fpTFW, "%.10f\n", padfGeoTransform[3]
01375 + 0.5 * padfGeoTransform[4]
01376 + 0.5 * padfGeoTransform[5] );
01377
01378 VSIFClose( fpTFW );
01379 return TRUE;
01380 }
01381
01382 /************************************************************************/
01383 /* GDALVersionInfo() */
01384 /************************************************************************/
01385
01404 const char *GDALVersionInfo( const char *pszRequest )
01405
01406 {
01407 static char szResult[128];
01408
01409
01410 if( pszRequest == NULL || EQUAL(pszRequest,"VERSION_NUM") )
01411 sprintf( szResult, "%d", GDAL_VERSION_NUM );
01412 else if( EQUAL(pszRequest,"RELEASE_DATE") )
01413 sprintf( szResult, "%d", GDAL_RELEASE_DATE );
01414 else if( EQUAL(pszRequest,"RELEASE_NAME") )
01415 sprintf( szResult, "%s", GDAL_RELEASE_NAME );
01416 else // --version
01417 sprintf( szResult, "GDAL %s, released %d/%02d/%02d",
01418 GDAL_RELEASE_NAME,
01419 GDAL_RELEASE_DATE / 10000,
01420 (GDAL_RELEASE_DATE % 10000) / 100,
01421 GDAL_RELEASE_DATE % 100 );
01422
01423 return szResult;
01424 }
01425
01426 /************************************************************************/
01427 /* GDALDecToDMS() */
01428 /* */
01429 /* Translate a decimal degrees value to a DMS string with */
01430 /* hemisphere. */
01431 /************************************************************************/
01432
01433 const char *GDALDecToDMS( double dfAngle, const char * pszAxis,
01434 int nPrecision )
01435
01436 {
01437 return CPLDecToDMS( dfAngle, pszAxis, nPrecision );
01438 }
01439
01440 /************************************************************************/
01441 /* GDALGCPsToGeoTransform() */
01442 /************************************************************************/
01443
01460 int GDALGCPsToGeoTransform( int nGCPCount, const GDAL_GCP *pasGCPs,
01461 double *padfGeoTransform, int bApproxOK )
01462
01463 {
01464 int iAnchor=0, iPnt1, iPnt2, i;
01465 double adfDPixel[2], adfDLine[2], adfDX[2], adfDY[2];
01466
01467 /* -------------------------------------------------------------------- */
01468 /* Recognise a few special cases. */
01469 /* -------------------------------------------------------------------- */
01470 if( nGCPCount < 2 )
01471 return FALSE;
01472
01473 if( nGCPCount == 2 )
01474 {
01475 if( pasGCPs[1].dfGCPPixel == pasGCPs[0].dfGCPPixel
01476 || pasGCPs[1].dfGCPLine == pasGCPs[0].dfGCPLine )
01477 return FALSE;
01478
01479 padfGeoTransform[1] = (pasGCPs[1].dfGCPX - pasGCPs[0].dfGCPX)
01480 / (pasGCPs[1].dfGCPPixel - pasGCPs[0].dfGCPPixel);
01481 padfGeoTransform[2] = 0.0;
01482
01483 padfGeoTransform[4] = 0.0;
01484 padfGeoTransform[5] = (pasGCPs[1].dfGCPY - pasGCPs[0].dfGCPY)
01485 / (pasGCPs[1].dfGCPLine - pasGCPs[0].dfGCPLine);
01486
01487 padfGeoTransform[0] = pasGCPs[0].dfGCPX
01488 - pasGCPs[0].dfGCPPixel * padfGeoTransform[1]
01489 - pasGCPs[0].dfGCPLine * padfGeoTransform[2];
01490
01491 padfGeoTransform[3] = pasGCPs[0].dfGCPY
01492 - pasGCPs[0].dfGCPPixel * padfGeoTransform[4]
01493 - pasGCPs[0].dfGCPLine * padfGeoTransform[5];
01494
01495 return TRUE;
01496 }
01497
01498 /* -------------------------------------------------------------------- */
01499 /* We use the first point as our anchor. Select two other */
01500 /* points that don't have the same Pixel value to analyse. */
01501 /* -------------------------------------------------------------------- */
01502 iPnt1 = -1;
01503 iPnt2 = -1;
01504 for( i = 1; (iPnt1 == -1 || iPnt2 == -1 ) && i < nGCPCount; i++ )
01505 {
01506 double dfDPixel = pasGCPs[i].dfGCPPixel-pasGCPs[iAnchor].dfGCPPixel;
01507 double dfDLine = pasGCPs[i].dfGCPLine - pasGCPs[iAnchor].dfGCPLine;
01508 double dfDX = pasGCPs[i].dfGCPX - pasGCPs[iAnchor].dfGCPX;
01509 double dfDY = pasGCPs[i].dfGCPY - pasGCPs[iAnchor].dfGCPY;
01510
01511 if( iPnt1 == -1 && ABS(dfDPixel) > 0.001 )
01512 {
01513 iPnt1 = i;
01514 adfDPixel[0] = dfDPixel;
01515 adfDLine[0] = dfDLine;
01516 adfDX[0] = dfDX;
01517 adfDY[0] = dfDY;
01518 }
01519 else if( iPnt2 == -1 )
01520 {
01521 iPnt2 = i;
01522 adfDPixel[1] = dfDPixel;
01523 adfDLine[1] = dfDLine;
01524 adfDX[1] = dfDX;
01525 adfDY[1] = dfDY;
01526 }
01527 }
01528
01529 /* -------------------------------------------------------------------- */
01530 /* If necessary, scale one of the points to avoid divide by */
01531 /* zeros. */
01532 /* -------------------------------------------------------------------- */
01533 if( ABS((adfDLine[0] / adfDPixel[0] - adfDLine[1])) < 0.0001 )
01534 {
01535 adfDX[1] *= 2;
01536 adfDY[1] *= 2;
01537 adfDPixel[1] *= 2;
01538 adfDLine[1] *= 2;
01539 }
01540
01541 /* -------------------------------------------------------------------- */
01542 /* Compute X related coefficients. */
01543 /* -------------------------------------------------------------------- */
01544 padfGeoTransform[2] =
01545 (adfDX[1] - (adfDPixel[1] * adfDX[0]) / adfDPixel[0])
01546 / (adfDLine[1] - (adfDLine[0]*adfDPixel[1]) / adfDPixel[0]);
01547
01548 padfGeoTransform[1] = (adfDX[0] - adfDLine[0] * padfGeoTransform[2])
01549 / adfDPixel[0];
01550
01551 /* -------------------------------------------------------------------- */
01552 /* Compute Y related coefficients. */
01553 /* -------------------------------------------------------------------- */
01554 padfGeoTransform[5] =
01555 (adfDY[1] - (adfDPixel[1] * adfDY[0]) / adfDPixel[0])
01556 / (adfDLine[1] - (adfDLine[0]*adfDPixel[1]) / adfDPixel[0]);
01557
01558 padfGeoTransform[4] =
01559 (adfDY[0] - adfDLine[0] * padfGeoTransform[5]) / adfDPixel[0];
01560
01561 /* -------------------------------------------------------------------- */
01562 /* Compute top/left origin. */
01563 /* -------------------------------------------------------------------- */
01564
01565 padfGeoTransform[0] = pasGCPs[0].dfGCPX
01566 - pasGCPs[0].dfGCPPixel * padfGeoTransform[1]
01567 - pasGCPs[0].dfGCPLine * padfGeoTransform[2];
01568
01569 padfGeoTransform[3] = pasGCPs[0].dfGCPY
01570 - pasGCPs[0].dfGCPPixel * padfGeoTransform[4]
01571 - pasGCPs[0].dfGCPLine * padfGeoTransform[5];
01572
01573 /* -------------------------------------------------------------------- */
01574 /* Now check if any of the input points fit this poorly. */
01575 /* -------------------------------------------------------------------- */
01576 if( !bApproxOK )
01577 {
01578 double dfPixelSize = ABS(padfGeoTransform[1])
01579 + ABS(padfGeoTransform[2])
01580 + ABS(padfGeoTransform[4])
01581 + ABS(padfGeoTransform[5]);
01582
01583 for( i = 0; i < nGCPCount; i++ )
01584 {
01585 double dfErrorX, dfErrorY;
01586
01587 dfErrorX =
01588 (pasGCPs[i].dfGCPPixel * padfGeoTransform[1]
01589 + pasGCPs[i].dfGCPLine * padfGeoTransform[2]
01590 + padfGeoTransform[0])
01591 - pasGCPs[i].dfGCPX;
01592 dfErrorY =
01593 (pasGCPs[i].dfGCPPixel * padfGeoTransform[4]
01594 + pasGCPs[i].dfGCPLine * padfGeoTransform[5]
01595 + padfGeoTransform[3])
01596 - pasGCPs[i].dfGCPY;
01597
01598 if( ABS(dfErrorX) > 0.25 * dfPixelSize
01599 || ABS(dfErrorY) > 0.25 * dfPixelSize )
01600 return FALSE;
01601 }
01602 }
01603
01604 return TRUE;
01605 }
01606
01607 /************************************************************************/
01608 /* -------------------------------------------------------------------- */
01609 /* The following stubs are present to ensure that older GDAL */
01610 /* bridges don't fail with newer libraries. */
01611 /* -------------------------------------------------------------------- */
01612 /************************************************************************/
01613
01614 CPL_C_START
01615
01616 void *GDALCreateProjDef( const char * pszDef )
01617 {
01618 CPLDebug( "GDAL", "GDALCreateProjDef no longer supported." );
01619 return NULL;
01620 }
01621
01622 CPLErr GDALReprojectToLongLat( void *pDef, double *, double * )
01623 {
01624 CPLDebug( "GDAL", "GDALReprojectToLatLong no longer supported." );
01625 return CE_Failure;
01626 }
01627
01628 CPLErr GDALReprojectFromLongLat( void *pDef, double *, double * )
01629 {
01630 CPLDebug( "GDAL", "GDALReprojectFromLatLong no longer supported." );
01631 return CE_Failure;
01632 }
01633
01634 void GDALDestroyProjDef( void *pDef )
01635
01636 {
01637 CPLDebug( "GDAL", "GDALDestroyProjDef no longer supported." );
01638 }
01639
01640 CPL_C_END
1.2.3-20001105 written by Dimitri van Heesch,
© 1997-2000