利用GDAL实现影像的几何校正

一、概述


    遥感影像和地理坐标进行关联的方式一般有好几种，一种是直接给出了仿射变换系数，即6个参数，左上角地理坐标，纵横方向上的分辨率，以及旋转系数。在这种情况下，求出某一像素点的地理坐标非常容易，直接用公式可以求出，具体代码如下：
}
 
可
能有的读者额看到这个觉得非常熟悉，这个是GDAL里面实现的由像素点转换到地理坐标点的转换函数，这个函数的声明在gdal.h中，定义在
gdaltransformer.cpp文件中，注意，这个函数计算的是影像像元的左上角的地理坐标，也可以计算中心点坐标，具体看需求。
 
另外一种方式是记录影像控制点坐标的数据，一般控制点都记录着控制点自身的地理坐标，还有控制点在影像中的像素坐标点，还有控制点必须是在某一特定的地理坐标系下，在GDAL中控制点的结构体为：
} GDAL_GCP;
 
那读者会有疑问了，这个控制点的数据怎么获得呢，GDAL里面有函数接口可以获得，具体如下；
    const GDAL_GCP *pGCPList= poDatasetchangsha->GetGCPs();
 
我首先给出一个例子先有一个比较直观的感受，我手上有一个spot传感器的数据，它是带控制点信息的，控制点刚好是图像的四个角点。比如，用arcgis打开就已经是纠正显示了。没有实施纠正显示的情况和ARCGIS实现了纠正显示的情况如下图：


 
Arcgis显示的时候，边缘的一些像素由于在原始影像上找不到对应的采样点，所以就用白色填充。
 

二、几何校正方法


关于影像几何校正，一般可以采用下面的几种方法：
1、几何多项式校正模型
  这个模型比较简单，回避了遥感影像的成像过程，适用于覆盖面积不大喝地形比较平坦的区域。
2、有理函数数学模型
   可以根据地面控制点求算出仿射变换系数，可以调用GDAL的函数，其函数是：
                        double *padfGeoTransform, int bApproxOK );
这个函数的定义在gdal_misc.cpp文件中，改天可以去了解下。求出仿射变换系数后可以求出图像的大致范围和分辨率，这样就可以目标图像所占的行数和列数。
 
下面的代码是求出影像的envelope，大概的范围
            dbMaxy = max(max(max(dbY[0],dbY[1]),dbY[2]),dbY[3]);
 
但是上面的方法计算出来的范围有可能不是很准确，所以这里要寻找一种比较精确的方法，那就是求出每个边缘像素的边缘对应的地理坐标，然后逐个比较大小求出比较精确的范围。这种方法可以采用GDAL算法模块里面GDALSuggestedWarpOutput2函数。
                          double *padfExtent,int nOptions)；
这个函数里面的参数的意义解释可以参考GDAL官方文档，也可以直接看GDAL源代码，下面我就给出各个参数的意思。
                          int nReqOrder, int bReversed );函数生成。第四个参数是图像纠正后的仿射变换的6个参数，pnPixels是纠正后图像的列数，pnLines是纠正后图像的行数，padfExtent就是图像纠正后的地理范围，nOptions是相关选项的设置，默认情况下可以设为空。



三、重采样及其校正实现


到
此，我们已经完成了原始图像的像素坐标到校正后的地理空间坐标建立了关联关系，并且输出影像的大小、地理范围、分辨率、像素位深度都已经知道，这样我们就
可以创建一个空的图像了，但是校正后的图像的各个像素的值还没确定，接下来就是如何确定校正后图像中各个像素的值，也就是图像处理中的重采样技术，一般大
家听的比较多的重采样技术可能是最邻近方法、双线性内插重采样以及三次卷积内插重采样等方法。还有一些其他的重采样方法，比如平均值法，中值法等。
重采样技术我们可以直接调用GDAL的函数去实现，但是我这里决定自己写重采样的过程，以便今后做一些改进，这里我先用最邻近采样方法实现重采样。
 
整个过程的代码如下：


[cpp] view plaincopy

1. /**
2. * @brief 几何多项式校正
3. * @param pszSrcFile         输入文件路径
4. * @param pszDstFile         输出文件路径
5. * @param nGCPCount          GCP点对个数
6. * @param pasGCPList         GCP点对列表
7. * @param pszDstWKT          输出文件投影
8. * @param iOrder             多项式次数
9. * @param dResX              输出图像X方向分辨率
10. * @param dResY              输出图像Y方向分辨率
11. * @param eResampleMethod    重采样方式
12. * @param pszFormat          输出文件格式
13. * @return  返回值，目前只返回true或者false
14. */  
15. int ImageWarpByGCP(const char * pszSrcFile,  
16.                   const char * pszDstFile,  
17.                   int nGCPCount,  
18.                   const GDAL_GCP *pasGCPList,  
19.                   const char * pszDstWKT,  
20.                   int iOrder,  
21.                   double dResX,  
22.                   double dResY,  
23.                   GDALResampleAlg eResampleMethod,  
24.                   const char * pszFormat)  
25. {  
26.  
27.    GDALAllRegister();  
28.    CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");  
29.  
30.    // 打开原始图像  
31.    GDALDatasetH hSrcDS = GDALOpen(pszSrcFile, GA_ReadOnly);  
32.    if (NULL == hSrcDS)  
33.    {  
34.        return 0;  
35.    }  
36.  
37.    GDALDataType eDataType = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS, 1));  
38.    int nBandCount = GDALGetRasterCount(hSrcDS);  
39.  
40.    // 创建几何多项式坐标转换关系  
41.    void *hTransform = GDALCreateGCPTransformer(nGCPCount, pasGCPList, iOrder, FALSE );  
42.    if( NULL == hTransform )  
43.    {  
44.        GDALClose( hSrcDS );  
45.        return 0;  
46.    }  
47.  
48.    // 计算输出图像四至范围、大小、仿射变换六参数等信息  
49.    double adfGeoTransform[6] = {0};  
50.    double adfExtent[4] = {0};  
51.    int    nPixels = 0, nLines = 0;  
52.  
53.    if( GDALSuggestedWarpOutput2( hSrcDS, GDALGCPTransform, hTransform,  
54.        adfGeoTransform, &nPixels, &nLines, adfExtent, 0 ) != CE_None )  
55.    {  
56.        GDALClose( hSrcDS );  
57.        return 0;  
58.    }  
59.  
60.    // 下面开始根据用户指定的分辨率来反算输出图像的大小和六参数等信息  
61.    double dResXSize = dResX;  
62.    double dResYSize = dResY;  
63.  
64.    //如果为0，则默认为原始影像的分辨率  
65.    if(dResXSize == 0.0 && dResYSize == 0.0 )  
66.    {  
67.        double dbGeoTran[6] = {0};  
68.        GDALGCPsToGeoTransform(nGCPCount,pasGCPList,dbGeoTran,0);  
69.        dResXSize = fabs(dbGeoTran[1]);  
70.        dResYSize = fabs(dbGeoTran[5]);  
71.    }  
72.  
73.    // 如果用户指定了输出图像的分辨率  
74.    else if ( dResXSize != 0.0 || dResYSize != 0.0 )  
75.    {  
76.        if ( dResXSize == 0.0 ) dResXSize = adfGeoTransform[1];  
77.        if ( dResYSize == 0.0 ) dResYSize = adfGeoTransform[5];  
78.    }  
79.  
80.    if ( dResXSize < 0.0 ) dResXSize = -dResXSize;  
81.    if ( dResYSize > 0.0 ) dResYSize = -dResYSize;  
82.  
83.    // 计算输出图像的范围  
84.    double minX = adfGeoTransform[0];  
85.    double maxX = adfGeoTransform[0] + adfGeoTransform[1] * nPixels;  
86.    double maxY = adfGeoTransform[3];  
87.    double minY = adfGeoTransform[3] + adfGeoTransform[5] * nLines;  
88.  
89.    nPixels = ceil(( maxX - minX ) / dResXSize );  
90.    nLines  = ceil(( minY - maxY ) / dResYSize ) ;  
91.    adfGeoTransform[0] = minX;  
92.    adfGeoTransform[3] = maxY;  
93.    adfGeoTransform[1] = dResXSize;  
94.    adfGeoTransform[5] = dResYSize;  
95.  
96.    // 创建输出图像  
97.    GDALDriverH hDriver = GDALGetDriverByName( pszFormat );  
98.    if (NULL == hDriver)  
99.    {  
100.        return 0;  
101.    }  
102.    GDALDatasetH hDstDS = GDALCreate( hDriver, pszDstFile, nPixels, nLines, nBandCount, eDataType, NULL );  
103.    if (NULL == hDstDS)  
104.    {  
105.        return 0;  
106.    }  
107.    GDALSetProjection( hDstDS,  pszDstWKT);  
108.    GDALSetGeoTransform( hDstDS, adfGeoTransform );  
109.  
110.    //获得原始图像的行数和列数  
111.    int nXsize = GDALGetRasterXSize(hSrcDS);  
112.    int nYsize = GDALGetRasterYSize(hSrcDS);  
113.  
114.    //然后是图像重采样  
115.    int nFlag = 0;  
116.    float dfValue = 0;  
117.    CPLErr err = CE_Failure;  
118.    //最邻近采样  
119.    for (int nBandIndex = 0; nBandIndex < nBandCount; nBandIndex ++)  
120.    {  
121.        GDALRasterBandH hSrcBand = GDALGetRasterBand(hSrcDS,nBandIndex+1);  
122.        GDALRasterBandH hDstBand = GDALGetRasterBand(hDstDS,nBandIndex+1);  
123.        for (int nRow = 0; nRow < nLines; nRow ++)  
124.        {  
125.            for (int nCol = 0; nCol < nPixels; nCol ++)  
126.            {  
127.                double dbX = adfGeoTransform[0] + nCol * adfGeoTransform[1]  
128.                + nRow  * adfGeoTransform[2];  
129.                double dbY = adfGeoTransform[3] + nCol * adfGeoTransform[4]  
130.                + nRow  * adfGeoTransform[5];  
131.  
132.                //由输出的图像地理坐标系变换到原始的像素坐标系  
133.                GDALGCPTransform(hTransform,TRUE,1,&dbX,&dbY,NULL,&nFlag);  
134.                int nXCol = (int)(dbX + 0.5);  
135.                int nYRow = (int)(dbY + 0.5);  
136.  
137.                //超出范围的用0填充  
138.                if (nXCol < 0 || nXCol >= nXsize || nYRow < 0 || nYRow >= nYsize)  
139.                {  
140.                    dfValue = 0;  
141.                }  
142.  
143.                else  
144.                {  
145.                    err = GDALRasterIO(hSrcBand,GF_Read,nXCol,nYRow,1,1,&dfValue,1,1,eDataType,0,0);  
146.                      
147.                }  
148.                err = GDALRasterIO(hDstBand,GF_Write,nCol,nRow,1,1,&dfValue,1,1,eDataType,0,0);  
149.  
150.            }  
151.        }  
152.          
153.    }  
154.  
155.      
156.  
157.    if (hTransform != NULL)  
158.    {  
159.        GDALDestroyGCPTransformer(hTransform);  
160.        hTransform = NULL;  
161.    }  
162.  
163.    GDALClose(hSrcDS );  
164.    GDALClose(hDstDS );  
165.  
166.  
167.    return 1;  
168. }  

外部调用的实例


[cpp] view plaincopy

1. GDALDataset* poDataset = (GDALDataset*)GDALOpen("E\\248-279_5\\SCENE01\\IMAGERY.TIF",GA_ReadOnly);  
2.  
3.    char* pszWkt1 = (char*)poDataset->GetProjectionRef();  
4.    if (0 == strlen(pszWkt1))  
5.    {  
6.        pszWkt1 = (char*)poDataset->GetGCPProjection();  
7.    }  
8.  
9.    int nGCPCount = poDataset->GetGCPCount();  
10.    const GDAL_GCP *pGCPList = poDataset->GetGCPs();  
11.  
12.    GDAL_GCP *pGCPs = new GDAL_GCP[nGCPCount];  
13.    memcpy(pGCPs,pGCPList,sizeof(GDAL_GCP)*nGCPCount);  
14.  
15.    GDALClose(poDataset);  
16.  
17.    const char* pszInFile = "E:\\248-279_5\\SCENE01\\IMAGERY.TIF";  
18.    const char* pszOutFile = "E:\\248-279_5\\SCENE01\\IMAGERY-校正.TIF";  
19.    ImageWarpByGCP(pszInFile,pszOutFile,nGCPCount,pGCPs,pszWkt1,0,0,0,GRA_NearestNeighbour,"GTiff");  

 
对
于上述算法的函数，有些读者可能看出来了，算法的瓶颈就是在那个多重循环那，逐点计算是不是效率太低了，的确，这里可以用分块技术来优化一下，甚至有可能
的话，用GPU加速也是不错的选择，我这边只是向大家分享一下怎么去校正带控制点的影像，算是抛砖引玉吧。如果有什么不妥的地方也可以指出来，大家一起讨
论。
 

四、后续研究


1.  如何进行校正显示，不生成新到的校正影像。
2.      如何分块处理
3.      如何用硬件加速，GPU加速
 

参考文献

楼主，问一下您概述里面这一行const GDAL_GCP *pGCPList= poDatasetchangsha->GetGCPs();  poDatasetchangsha我总是报错ntelliSense:  "const GDAL_GCP *" 类型的值不能用于初始化 "const GDAL_GCP *" 类型的实体是怎么回事？头文件#include "stdafx.h"
#include "gdal_priv.h"
#include "cpl_string.h"