gdal几何校正源码 <含bug>

以下 代码输出结果并不正确。输出图像仅仅是将原始图像进行了重采样。欢迎高手指点。
另外，对于gdal到底是否可以作几何校正现在持怀疑态度。

#include "gdal_priv.h"
#include "gdalwarper.h"
#include "gdal.h"
#include "cpl_minixml.h"
#include "ogr_api.h"
#include "ogr_spatialref.h"

int main(int argc, char* argv[])
{
 GDALDatasetH  hSrcDS, hDstDS;

    // Open input and output files.

    GDALAllRegister();

    hSrcDS = GDALOpen( "E:\\spot_pan.tif", GA_ReadOnly );
 
    // Setup warp options.
   // GDALWarpOptions *psFirstWarpOptions = GDALCreateWarpOptions();
  

 // my transformer
 double* C = new double[30];
    double* R = new double[30];
    double* Xcontrol = new double[30];
    double* Ycontrol = new double[30];

  GDAL_GCP* pGCPs =new GDAL_GCP[30];
    FILE *pf_gcp;
 char *filename = "gcp_30_多项式正射.txt";

 pf_gcp = fopen(filename, "r");
   // assert(pf_gcp!=NULL);
  
 for(int i = 0; i < 30; i++)
 {
  fscanf(pf_gcp,"%lf %lf %lf %lf\n", &C[i], &R[i], &Xcontrol[i], &Ycontrol[i]);
     pGCPs[i].dfGCPLine = R[i];
  pGCPs[i].dfGCPPixel = C[i];
     pGCPs[i].dfGCPX = Xcontrol[i];
        pGCPs[i].dfGCPY = Ycontrol[i];
  pGCPs[i].dfGCPZ = 0.0;
       
  pGCPs[i].pszId = "1";
    //  _itoa(i,pGCPs[i].pszId,10);
  pGCPs[i].pszInfo = "";
 
 }
 fclose(pf_gcp);
  
   //* 控制点误差计算 Xcontrol 中返回的是其回归值
   void* yyy= GDALCreateGCPTransformer(30,pGCPs,3,0);
   double* z = new double[30];
   int *djkdjk = new int(1);
   int outresult = GDALGCPTransform(yyy,TRUE,30,Xcontrol,Ycontrol,z,djkdjk);
   i=30;
   for(i=0;i<30;i++)
    printf("%f %f\n",Xcontrol[i],Ycontrol[i]);
  // */

   void* pTransformerArg = GDALCreateGCPTransformer(30,pGCPs,3,0);
 
   
 // 求得预测的四致范围和分辨率
 double adfDstGeoTransform[6];
    int nPixels=0, nLines=0;
    CPLErr eErr;

    eErr = GDALSuggestedWarpOutput( hSrcDS,
                                    GDALGCPTransform, pTransformerArg,
                                    adfDstGeoTransform, &nPixels, &nLines );
    CPLAssert( eErr == CE_None );

    //创建输出影像
    GDALDriverH hDriver;
    GDALDataType eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1));;
 hDriver = GDALGetDriverByName( "GTiff" );
    CPLAssert( hDriver != NULL );
    hDstDS = GDALCreate( hDriver, "outRectified.tif", nPixels, nLines,
                         GDALGetRasterCount(hSrcDS), eDT, NULL );
    CPLAssert( hDstDS != NULL );

    // Write out the projection definition.
    OGRSpatialReference oSRS;
 char *pszSRS_WKT = NULL;
 
 oSRS.SetWellKnownGeogCS( "WGS84" ); //默认形式投影信息
 oSRS.SetUTM( 11, TRUE );
 oSRS.exportToWkt( &pszSRS_WKT );
    GDALSetProjection( hDstDS, pszSRS_WKT );
    GDALSetGeoTransform( hDstDS, adfDstGeoTransform );
 CPLFree( pszSRS_WKT );
  
   

 //确定输出影像的 warpoption的选项
 GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();

    psWarpOptions->hSrcDS = hSrcDS;
    psWarpOptions->nBandCount = 1;
    psWarpOptions->panSrcBands =
        (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );
    psWarpOptions->panSrcBands[0] = 1;
      
  psWarpOptions->hDstDS = hDstDS;
     psWarpOptions->panDstBands =
        (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );
     psWarpOptions->panDstBands[0] = 1;
     
   

  char* gcpid= new char[30];
 for( i = 0; i < 30; i++)
 {
 
     pGCPs[i].dfGCPLine = R[i];
  pGCPs[i].dfGCPPixel = C[i];
     pGCPs[i].dfGCPX = Xcontrol[i];
        pGCPs[i].dfGCPY = Ycontrol[i];
  pGCPs[i].dfGCPZ = 0.0;
        
  _itoa(i,gcpid,10);
  pGCPs[i].pszId = gcpid;
   
  pGCPs[i].pszInfo = "";

  gcpid++;
 
 }

  psWarpOptions->pTransformerArg = GDALCreateGCPTransformer(30,pGCPs,3,0);
   psWarpOptions->pfnTransformer = GDALGCPTransform;
     psWarpOptions->eResampleAlg = GRA_Bilinear ;
  psWarpOptions->pfnProgress = GDALTermProgress;

 
 // Initialize and execute the warp operation.
    GDALWarpOperation oOperation;

    eErr =  oOperation.Initialize( psWarpOptions );
 CPLAssert(eErr==CE_None );
     eErr = oOperation.ChunkAndWarpMulti( 0, 0, nPixels, nLines);
    CPLAssert(eErr==CE_None );
    GDALDestroyGCPTransformer( psWarpOptions->pTransformerArg );
    GDALDestroyWarpOptions( psWarpOptions );

    GDALClose( hDstDS );
    GDALClose( hSrcDS );

 
 return 0;
}

 

转载于:https://www.cnblogs.com/luoriyouhun/archive/2008/01/19/1045700.html