近期在做基于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;第二处。在函数GDALCreateRPCTransformer()中，主要是将參数papszOptions中的像方改正系数进行解析，然后给结构体中新加的两个数组赋值。便于兴许进行坐标转换的时候使用。改动后的代码例如以下，因为代码太长，仅仅贴出我改动的部分代码。以下代码中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: pctestanda.tif D: pctestanda_affine.tif