改動GDAL庫支持RPC像方改正模型
近期在做基於RPC的像方改正模型。方便對數據進行測試，改動了GDAL庫中的RPC糾正模型，使之能夠支持RPC像方改正參數。
以下是RPC模型的公式，rn,cn為歸一化之後的圖像行列號坐標，PLH為歸一化後的經度緯度高程。
將上面的公式變形，使用偏移系數和縮放系數帶入，能夠得到圖像的行列號坐標與經緯度坐標之間的坐標轉換關系。整理後的公式例如以下所看到的。下標帶s的為縮放系數，下標為0的表示偏移系數，rc為圖像行列號，此處的PLH為地面經緯度坐標。P1～P4為有理函數的多項式系數。
使用像方改正模型的公式例如以下所看到的，Line和Sample為圖像的行列號，rc為通過RPC模型將地面點經緯度高程計算得到的行列號，deltaR和DeltaC為像方改正數。
deltaR和DeltaC像方改正數使用仿射變換模型，詳細公式例如以下，A0～A2為行方向的改正系數，B0～B2為列方向改正系數。無改正時這六個系數均為0.
將上面兩個公式合並之後，再將DeltaR和DeltaC移向等式坐標，合並同類項之後。就得到了終於的一個仿射變換系數。此時無改正時。六個系數變為0 1 0 0 0 1。該系數為以下是用的終於系數。
改動的代碼非常少，在GDAL源代碼中的alg目錄裏面的gdal_rpc.cpp中，詳細改動三處地方。
第一處。GDALRPCTransformInfo結構體，在結構體中添加兩個double [6]的數組。用於保存RPC像方改正系數。
改動後的代碼例如以下。最後兩個參數adfAffineTransform和adfReverseAffineTransform分別表示RPC像方改正系數及其逆變換系數。
typedef struct {
GDALTransformerInfo sTI;
GDALRPCInfo sRPC;
double adfPLToLatLongGeoTransform[6];
int bReversed;
double dfPixErrThreshold;
double dfHeightOffset;
double dfHeightScale;
char *pszDEMPath;
DEMResampleAlg eResampleAlg;
int bHasTriedOpeningDS;
GDALDataset *poDS;
OGRCoordinateTransformation *poCT;
double adfGeoTransform[6];
double adfReverseGeoTransform[6];
double adfAffineTransform[6]; //RPC adjustment affine transform
double adfReverseAffineTransform[6]; //RPC adjustment reverse affine transform
} GDALRPCTransformInfo;便於興許進行坐標轉換的時候使用。改動後的代碼例如以下，因為代碼太長，僅僅貼出我改動的部分代碼。以下代碼中The Affine transform parameters部分的代碼由我新加，主要是通過一個RPC_AFFINE來指定像方改正的六個系數，六個系數中間用空格隔開。然後將解析後的六個系數計算逆變換系數。默認系數是0 1 0 0 0 1，表示不進行像方改正。
//前面的代碼省略
/* -------------------------------------------------------------------- */
/* The DEM interpolation */
/* -------------------------------------------------------------------- */
const char *pszDEMInterpolation = CSLFetchNameValueDef( papszOptions, "RPC_DEMINTERPOLATION", "bilinear" );
if(EQUAL(pszDEMInterpolation, "near" ))
psTransform->eResampleAlg = DRA_NearestNeighbour;
else if(EQUAL(pszDEMInterpolation, "bilinear" ))
psTransform->eResampleAlg = DRA_Bilinear;
else if(EQUAL(pszDEMInterpolation, "cubic" ))
psTransform->eResampleAlg = DRA_Cubic;
else
psTransform->eResampleAlg = DRA_Bilinear;
/* -------------------------------------------------------------------- */
/* The Affine transform parameters */
/* -------------------------------------------------------------------- */
const char *pszRpcAffine = CSLFetchNameValueDef( papszOptions, "RPC_AFFINE", "0 1 0 0 0 1" );
if(pszRpcAffine != NULL) //解析RPC像方改正仿射變換參數
{
char** papszTokens = CSLTokenizeString2( pszRpcAffine, " ", 0 );
int nTokens = CSLCount(papszTokens);
if(nTokens == 6) //must be 6
{
for (int i=0; i<6; i++)
psTransform->adfAffineTransform[i] = atof(papszTokens[i]);
}
else
{
psTransform->adfAffineTransform[1] = 1;
psTransform->adfAffineTransform[5] = 1;
}
GDALInvGeoTransform( psTransform->adfAffineTransform, psTransform->adfReverseAffineTransform );
CSLDestroy(papszTokens);
}
/* -------------------------------------------------------------------- */
/* Establish a reference point for calcualating an affine */
/* geotransform approximate transformation. */
/* -------------------------------------------------------------------- */
//後面的代碼省略第三處改動的地方就是在進行坐標轉換的函數中。即函數GDALRPCTransform()中，這裏面主要改動兩部分。第一個部分就是坐標正變換的時候，第二個是坐標逆變換的時候，在坐標正變換的時候，也就是從經緯度坐標換算到原始圖像行列號坐標時，等使用RPC模型轉換完畢後。再使用仿射變換的逆變換系數進行改正。在坐標逆變換的時候，也就是從原始圖像行列號換算到經緯度坐標時，先使用仿射變換的正變換系數進行改正，然後將改正後的行列號帶入RPC模型進行轉換到經緯度。改動的代碼就兩行。可是原始的代碼太多。就貼出來我新增的部分。//此處是坐標正變換的時候新增的代碼
GDALApplyGeoTransform(psTransform->adfReverseAffineTransform, padfX[i], padfY[i], padfX + i, padfY + i );
panSuccess[i] = TRUE;
//此處是坐標逆變換的時候新增的代碼
GDALApplyGeoTransform(psTransform->adfAffineTransform, padfX[i], padfY[i], padfX + i, padfY + i );
double dfResultX, dfResultY;改動完之後，保存，然後又一次編譯GDAL庫就可以。之後我們就能夠使用gdalwarp.exe這個超牛的工具來進行校正了。詳細的命令就是在原來使用-rpc的命令基礎上，添加一個-to “RPC_AFFINE=0 1 0 0 0 1”就可以。當然這六個系數須要自己敲代碼使用控制點來進行反算，僅僅要三個控制點就可以，使用1個或兩個控制點僅僅能計算一個平移模型，即上面公式中的A0和B0。完整的命令行為：
gdalwarp.exe -rpc -to "RPC_AFFINE=-32.714672501057066 0.999199897235577 0.000158731686899 28.720843336473692 0.000589585516339 1.000068008511035" D:\rpctest\banda.tif D:\rpctest\banda_affine.tif
改動GDAL庫支持RPC像方改正模型