Main Page   Class Hierarchy   Compound List   File List   Compound Members   File Members   Related Pages  

gdal_misc.cpp

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

Generated at Sat Dec 21 14:01:58 2002 for GDAL by doxygen1.2.3-20001105 written by Dimitri van Heesch, © 1997-2000