利用GDAL實現影像的幾何校正
一、概述
遙感影像和地理坐標進行關聯的方式一般有好幾種，一種是直接給出了仿射變換系數，即6個參數，左上角地理坐標，縱橫方向上的分辨率，以及旋轉系數。在這樣的情況下，求出某一像素點的地理坐標非常easy，直接用公式能夠求出，詳細代碼例如以下：
void CPL_STDCALL GDALApplyGeoTransform(double *padfGeoTransform,
double dfPixel, double dfLine,
double *pdfGeoX, double *pdfGeoY
{
*pdfGeoX = padfGeoTransform[0]+ dfPixel * padfGeoTransform[1]
+ dfLine * padfGeoTransform[2];
*pdfGeoY = padfGeoTransform[3]+ dfPixel * padfGeoTransform[4]
+ dfLine * padfGeoTransform[5];
}
可能有的讀者額看到這個覺得非常熟悉。這個是GDAL裏面實現的由像素點轉換到地理坐標點的轉換函數，這個函數的聲明在gdal.h中，定義在gdaltransformer.cpp文件裏，註意。這個函數計算的是影像像元的左上角的地理坐標，也能夠計算中心點坐標，詳細看需求。
第二種方式是記錄影像控制點坐標的數據，一般控制點都記錄著控制點自身的地理坐標。還有控制點在影像中的像素坐標點。還有控制點必須是在某一特定的地理坐標系下，在GDAL中控制點的結構體為：
/** Ground Control Point */
typedef struct
{
/** Unique identifier, often numeric */
char *pszId;
/** Informational message or "" */
char *pszInfo;
/** Pixel (x) location of GCP on raster */
double
/** Line (y) location of GCP on raster */
double dfGCPLine;
/** X position of GCP in georeferenced space*/
double dfGCPX;
/** Y position of GCP in georeferenced space*/
double dfGCPY;
/** Elevation of GCP, or zero if not known */
double dfGCPZ;
} GDAL_GCP;
那讀者會有疑問了，這個控制點的數據怎麽獲得呢。GDAL裏面有函數接口能夠獲得，詳細例如以下；
char* pszWkt1 = (char*)poDatasetchangsha->GetProjectionRef();
if (0 == strlen(pszWkt1))
{
pszWkt1= (char*)poDatasetchangsha->GetGCPProjection();
}
int nGCPCount = poDatasetchangsha->GetGCPCount();
const GDAL_GCP *pGCPList= poDatasetchangsha->GetGCPs();
我首先給出一個樣例先有一個比較直觀的感受，我手上有一個spot傳感器的數據，它是帶控制點信息的，控制點剛好是圖像的四個角點。比方。用arcgis打開就已經是糾正顯示了。
沒有實施糾正顯示的情況和ARCGIS實現了糾正顯示的情況例如以下圖：
Arcgis顯示的時候，邊緣的一些像素因為在原始影像上找不到相應的採樣點，所以就用白色填充。
二、幾何校正方法
關於影像幾何校正，一般能夠採用以下的幾種方法：
1、幾何多項式校正模型
這個模型比較簡單，回避了遙感影像的成像過程，適用於覆蓋面積不大喝地形比較平坦的區域。
2、有理函數數學模型
建立像素和地面位置相應關系的簡單數學模型。理論上能夠達到跟嚴格衛星傳感器模型相當的定位精度。其長處在於多項式包括高程信息,能夠提高校正精度。
這個模型自帶了高程信息，一般須要影像提供RPC文件，或者用戶自己選擇地面控制點。
3、局部區域校正模型
這個模型的基本思想是利用控制點建立不規則三角網。然後分區域利用幾何多項式校正，可是這樣的模型須要非常多的控制點，對於地形起伏非常大的區域須要的控制點很多其它，往往實施難度比較大。
4、衛星傳感器模型
這樣的模型是依據衛星遙感成像的幾何關系，利用攝影測量學中成像瞬間的地面點、透視中心以及相應的像素點三點共線建立起來的一種模型，可是此模型須要DEM數據，長處是精度比較高。
因為本文主要關註幾何多項式模型，所以僅僅利用GDAL中幾何多項式的模型來校正。
能夠依據地面控制點求算出仿射變換系數，能夠調用GDAL的函數。其函數是：
int CPL_DLL CPL_STDCALL
GDALGCPsToGeoTransform(int nGCPCount,const GDAL_GCP*pasGCPs,
double *padfGeoTransform, int bApproxOK );
這個函數的定義在gdal_misc.cpp文件裏，改天能夠去了解下。
求出仿射變換系數後能夠求出圖像的大致範圍和分辨率，這樣就能夠目標圖像所占的行數和列數。
以下的代碼是求出影像的envelope，大概的範圍
double dbX[4];
double dbY[4];
dbX[0] = dblGeoTransform[0]+ 0*dblGeoTransform[1] + 0*dblGeoTransform[2];
dbY[0] = dblGeoTransform[3]+ 0*dblGeoTransform[4] + 0*dblGeoTransform[5];
//右上角坐標
dbX[1] = dblGeoTransform[0]+ (lWidth/*-0.5*/)*dblGeoTransform[1] + 0*dblGeoTransform[2];
dbY[1] = dblGeoTransform[3]+ (lWidth/*-0.5*/)*dblGeoTransform[4] + 0*dblGeoTransform[5];
//右下角點坐標
dbX[2] = dblGeoTransform[0]+ (lWidth/*-0.5*/)*dblGeoTransform[1] + (lHeight/*-0.5*/)*dblGeoTransform[2];
dbY[2] = dblGeoTransform[3]+ (lWidth/*-0.5*/)*dblGeoTransform[4] + (lHeight/*-0.5*/)*dblGeoTransform[5];
//左下角坐標
dbX[3] = dblGeoTransform[0]+ 0*dblGeoTransform[1] + (lHeight/*-0.5*/)*dblGeoTransform[2];
dbY[3] = dblGeoTransform[3]+ 0*dblGeoTransform[4] + (lHeight/*-0.5*/)*dblGeoTransform[5];
double dbMinx = 0;
double dbMaxx = 0;
double dbMiny = 0;
double dbMaxy = 0;
dbMinx = min(min(min(dbX[0],dbX[1]),dbX[2]),dbX[3]);
dbMaxx = max(max(max(dbX[0],dbX[1]),dbX[2]),dbX[3]);
dbMiny = min(min(min(dbY[0],dbY[1]),dbY[2]),dbY[3]);
dbMaxy = max(max(max(dbY[0],dbY[1]),dbY[2]),dbY[3]);
可是上面的方法計算出來的範圍有可能不是非常準確，所以這裏要尋找一種比較精確的方法。那就是求出每一個邊緣像素的邊緣相應的地理坐標，然後逐個比較大小求出比較精確的範圍。這樣的方法能夠採用GDAL算法模塊裏面GDALSuggestedWarpOutput2函數。
CPLErr CPL_STDCALL
GDALSuggestedWarpOutput2(GDALDatasetH hSrcDS,
GDALTransformerFunc pfnTransformer,
void *pTransformArg,
double *padfGeoTransformOut,
int *pnPixels, int *pnLines,
double *padfExtent,int nOptions)。
這個函數裏面的參數的意義解釋能夠參考GDAL官方文檔，也能夠直接看GDAL源代碼，以下我就給出各個參數的意思。
第一個參數是數據源，第二個參數是轉換的回調函數，比方是GCP控制點轉換的函數GDALGCPTransform。或者是RPC校正的函數等等。第三個參數是代表轉換的關系，這裏的轉換關系是控制點幾何校正的轉換關系，此轉換關系能夠由void CPL_DLL *
GDALCreateGCPTransformer(int nGCPCount,const GDAL_GCP*pasGCPList,
int nReqOrder, int bReversed );函數生成。第四個參數是圖像糾正後的仿射變換的6個參數，pnPixels是糾正後圖像的列數，pnLines是糾正後圖像的行數。padfExtent就是圖像糾正後的地理範圍。nOptions是相關選項的設置，默認情況下能夠設為空。
三、重採樣及其校正實現
到此，我們已經完畢了原始圖像的像素坐標到校正後的地理空間坐標建立了關聯關系，而且輸出影像的大小、地理範圍、分辨率、像素位深度都已經知道。這樣我們就能夠創建一個空的圖像了。可是校正後的圖像的各個像素的值還沒確定。接下來就是怎樣確定校正後圖像中各個像素的值，也就是圖像處理中的重採樣技術，一般大家聽的比較多的重採樣技術可能是最鄰近方法、雙線性內插重採樣以及三次卷積內插重採樣等方法。另一些其它的重採樣方法。比方平均值法，中值法等。
重採樣技術我們能夠直接調用GDAL的函數去實現。可是我這裏決定自己寫重採樣的過程，以便今後做一些改進，這裏我先用最鄰近採樣方法實現重採樣。
整個過程的代碼例如以下：
/**
* @brief 幾何多項式校正
* @param pszSrcFile 輸入文件路徑
* @param pszDstFile 輸出文件路徑
* @param nGCPCount GCP點對個數
* @param pasGCPList GCP點對列表
* @param pszDstWKT 輸出文件投影
* @param iOrder 多項式次數
* @param dResX 輸出圖像X方向分辨率
* @param dResY 輸出圖像Y方向分辨率
* @param eResampleMethod 重採樣方式
* @param pszFormat 輸出文件格式
* @return 返回值，眼下僅僅返回true或者false
*/
int ImageWarpByGCP(const char * pszSrcFile,
const char * pszDstFile,
int nGCPCount,
const GDAL_GCP *pasGCPList,
const char * pszDstWKT,
int iOrder,
double dResX,
double dResY,
GDALResampleAlg eResampleMethod,
const char * pszFormat)
{
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
// 打開原始圖像
GDALDatasetH hSrcDS = GDALOpen(pszSrcFile, GA_ReadOnly);
if (NULL == hSrcDS)
{
return 0;
}
GDALDataType eDataType = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS, 1));
int nBandCount = GDALGetRasterCount(hSrcDS);
// 創建幾何多項式坐標轉換關系
void *hTransform = GDALCreateGCPTransformer(nGCPCount, pasGCPList, iOrder, FALSE );
if( NULL == hTransform )
{
GDALClose( hSrcDS );
return 0;
}
// 計算輸出圖像四至範圍、大小、仿射變換六參數等信息
double adfGeoTransform[6] = {0};
double adfExtent[4] = {0};
int nPixels = 0, nLines = 0;
if( GDALSuggestedWarpOutput2( hSrcDS, GDALGCPTransform, hTransform,
adfGeoTransform, &nPixels, &nLines, adfExtent, 0 ) != CE_None )
{
GDALClose( hSrcDS );
return 0;
}
// 以下開始依據用戶指定的分辨率來反算輸出圖像的大小和六參數等信息
double dResXSize = dResX;
double dResYSize = dResY;
//假設為0，則默覺得原始影像的分辨率
if(dResXSize == 0.0 && dResYSize == 0.0 )
{
double dbGeoTran[6] = {0};
GDALGCPsToGeoTransform(nGCPCount,pasGCPList,dbGeoTran,0);
dResXSize = fabs(dbGeoTran[1]);
dResYSize = fabs(dbGeoTran[5]);
}
// 假設用戶指定了輸出圖像的分辨率
else if ( dResXSize != 0.0 || dResYSize != 0.0 )
{
if ( dResXSize == 0.0 ) dResXSize = adfGeoTransform[1];
if ( dResYSize == 0.0 ) dResYSize = adfGeoTransform[5];
}
if ( dResXSize < 0.0 ) dResXSize = -dResXSize;
if ( dResYSize > 0.0 ) dResYSize = -dResYSize;
// 計算輸出圖像的範圍
double minX = adfGeoTransform[0];
double maxX = adfGeoTransform[0] + adfGeoTransform[1] * nPixels;
double maxY = adfGeoTransform[3];
double minY = adfGeoTransform[3] + adfGeoTransform[5] * nLines;
nPixels = ceil(( maxX - minX ) / dResXSize );
nLines = ceil(( minY - maxY ) / dResYSize ) ;
adfGeoTransform[0] = minX;
adfGeoTransform[3] = maxY;
adfGeoTransform[1] = dResXSize;
adfGeoTransform[5] = dResYSize;
// 創建輸出圖像
GDALDriverH hDriver = GDALGetDriverByName( pszFormat );
if (NULL == hDriver)
{
return 0;
}
GDALDatasetH hDstDS = GDALCreate( hDriver, pszDstFile, nPixels, nLines, nBandCount, eDataType, NULL );
if (NULL == hDstDS)
{
return 0;
}
GDALSetProjection( hDstDS, pszDstWKT);
GDALSetGeoTransform( hDstDS, adfGeoTransform );
//獲得原始圖像的行數和列數
int nXsize = GDALGetRasterXSize(hSrcDS);
int nYsize = GDALGetRasterYSize(hSrcDS);
//然後是圖像重採樣
int nFlag = 0;
float dfValue = 0;
CPLErr err = CE_Failure;
//最鄰近採樣
for (int nBandIndex = 0; nBandIndex < nBandCount; nBandIndex ++)
{
GDALRasterBandH hSrcBand = GDALGetRasterBand(hSrcDS,nBandIndex+1);
GDALRasterBandH hDstBand = GDALGetRasterBand(hDstDS,nBandIndex+1);
for (int nRow = 0; nRow < nLines; nRow ++)
{
for (int nCol = 0; nCol < nPixels; nCol ++)
{
double dbX = adfGeoTransform[0] + nCol * adfGeoTransform[1]
+ nRow * adfGeoTransform[2];
double dbY = adfGeoTransform[3] + nCol * adfGeoTransform[4]
+ nRow * adfGeoTransform[5];
//由輸出的圖像地理坐標系變換到原始的像素坐標系
GDALGCPTransform(hTransform,TRUE,1,&dbX,&dbY,NULL,&nFlag);
int nXCol = (int)(dbX + 0.5);
int nYRow = (int)(dbY + 0.5);
//超出範圍的用0填充
if (nXCol < 0 || nXCol >= nXsize || nYRow < 0 || nYRow >= nYsize)
{
dfValue = 0;
}
else
{
err = GDALRasterIO(hSrcBand,GF_Read,nXCol,nYRow,1,1,&dfValue,1,1,eDataType,0,0);
}
err = GDALRasterIO(hDstBand,GF_Write,nCol,nRow,1,1,&dfValue,1,1,eDataType,0,0);
}
}
}
if (hTransform != NULL)
{
GDALDestroyGCPTransformer(hTransform);
hTransform = NULL;
}
GDALClose(hSrcDS );
GDALClose(hDstDS );
return 1;
}外部調用的實例
GDALDataset* poDataset = (GDALDataset*)GDALOpen("E\\248-279_5\\SCENE01\\IMAGERY.TIF",GA_ReadOnly);
char* pszWkt1 = (char*)poDataset->GetProjectionRef();
if (0 == strlen(pszWkt1))
{
pszWkt1 = (char*)poDataset->GetGCPProjection();
}
int nGCPCount = poDataset->GetGCPCount();
const GDAL_GCP *pGCPList = poDataset->GetGCPs();
GDAL_GCP *pGCPs = new GDAL_GCP[nGCPCount];
memcpy(pGCPs,pGCPList,sizeof(GDAL_GCP)*nGCPCount);
GDALClose(poDataset);
const char* pszInFile = "E:\\248-279_5\\SCENE01\\IMAGERY.TIF";
const char* pszOutFile = "E:\\248-279_5\\SCENE01\\IMAGERY-校正.TIF";
ImageWarpByGCP(pszInFile,pszOutFile,nGCPCount,pGCPs,pszWkt1,0,0,0,GRA_NearestNeighbour,"GTiff");
對於上述算法的函數。有些讀者可能看出來了，算法的瓶頸就是在那個多重循環那，逐點計算是不是效率太低了。的確，這裏能夠用分塊技術來優化一下。甚至有可能的話。用GPU加速也是不錯的選擇，我這邊僅僅是向大家分享一下怎麽去校正帶控制點的影像，算是拋磚引玉吧。假設有什麽不妥的地方也能夠指出來，大家一起討論。
四、興許研究
1. 怎樣進行校正顯示，不生成新到的校正影像。
2. 怎樣分塊處理
3. 怎樣用硬件加速，GPU加速
參考文獻
張永生,劉軍. 高分辨率遙感衛星立體影像RPC模型定位的算法及其優化測繪project, 2004,
劉軍,王冬紅,毛國苗. 基於RPC模型的IKONOS衛星影像高精度立體定位. 測繪通報, 2004
劉軍,張永生,範永弘. 基於通用成像模型———有理函數模型的攝影測量定位方法. 測繪通報, 2003
鞏丹超,鄧雪清,張雲彬. 新型遙感衛星傳感器幾何模型—有理函數模型. 海洋測繪, 2003
GDAL官方文檔以及源代碼
阮秋琦、數字圖像處理（第三版）
利用GDAL實現影像的幾何校正