最近在做基于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 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.                        */
/* -------------------------------------------------------------------- */
//后面的代码省略

//此处是坐标正变换的时候新增的代码
	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;

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