Surf算子匹配

/**
* @brief 计算获得金字塔影像的surf算子
* @param strImg 影像路径
* @param surfPnts surf算子容器
* @param iMaxSearchSize 金字塔影像宽高的最大阈值
* @param dblThres surf阈值
* @param Surf匹配使用的影像金字塔层号
* SURF算法计算过程缓慢，读取金字塔影像计算，读取金字塔的层数根据
* iMaxSearchSize确定
*/
bool OnExtarctPyramidSurfPnt(string strImg, vector<Ipoint> &surfPnts,
							 int iMaxSearchSize, double dblThres, 
							 int &iLevel)

OpenSURF基于OpenCV开发，使用该开源代码，需要同时下载OpenCV和OpenSURF

（1）影像建立金字塔

遥感影像通常尺寸较大，直接利用OpenSURF提取Surf算子运算速度较慢（利用OpenSurf的源代码直接做大图像的匹配根本是不可能的，速度太慢），因此需要基于金字塔影像进行提取

使用GDAL创建金字塔

/**
 * @brief 构建影像金字塔
 */
bool COMRasterDataset::BuildPyramid()
{
	long lWidth, lHeight;
	if (!GetSize(lWidth, lHeight))
	{
		return false;
	}

	int iLevelCount = 0;
	long lMinSize = lWidth<lHeight?lWidth:lHeight;
	while (lMinSize > iMaxLevelPyramidSize)
	{
		lMinSize /= 2;
		iLevelCount++;
	}
	int *iOverviewLists = new int[iLevelCount];
	for (int i=0; i<iLevelCount; i++)
	{
		iOverviewLists[i] = pow(2.0, i+1);
	}
	if (mpDataset->BuildOverviews("GAUSS", iLevelCount, iOverviewLists,
		0, NULL, GDALDummyProgress, NULL) != CE_None)
	{
		delete []iOverviewLists;
		return false;
	}

	delete []iOverviewLists;

	return true;
}

（2）设置匹配数据的尺寸的最大值iMaxSearchSize（如1024），确定读取金字塔的第几层。（在GDAL中pBand->GetOverview(iLevel),iLevel=0获取的并不是原始数据层）

iLevel = -1;
long lWidth, lHeight;
if (!pRasterDataset->GetSize(lWidth, lHeight))
{
	return false;
}
long lMinSize = lWidth<lHeight?lWidth:lHeight;
while (lMinSize>iMaxSearchSize)
{
	iLevel++;
	lMinSize /= 2;
	lWidth /= 2;
	lHeight /= 2;
}

（3）读取金字塔层的影像数据，并进行正规化（将数据规则化到0-1之间）
对于多波段影像的彩色数据，最好能够按照一定的比例合并为一个波段进行计算，例如pPanBand = 0.2*pBand1 + 0.5*pBand2 + 0.3*pBand3;
读取金字塔层的数据：

GDALRasterBand *pBand = 
		(GDALRasterBand*)mpDataset->GetRasterBand(iBand+1);
if (!pBand)
{
	return false;
}
GDALRasterBand *pOverviewBand = 
	pBand->GetOverview(iLevel);
if (!pOverviewBand)
{
	return false;
}

CPLErr err = pOverviewBand->RasterIO(GF_Read, iStartCol, iStartRow,
	iWidth, iHeight, pBuffer, iWidth, iHeight,
	pOverviewBand->GetRasterDataType(), 0, 0);
if (err != CE_None)
{
	return false;
}

将数据进行归一化，以unsigned char型为例：
for (int i=0; i<iHeight; i++)
{
	for (int j=0; j<iWidth; j++)
	{
		pDest[i*iWidth+j] = float(pSrc[i*iWidth+j]/255.0);
	}
}

（4）OpenSURF中计算用到的是积分影像，对获取的正规化数据进行积分计算

m,n位置的像素值等于所有横坐标小于等于m，纵坐标小于等n的值

if (!pSrc)
{
	return false;
}

float fTemp = 0;
for (int i=0; i<iWidth; i++)
{
	fTemp += pSrc[i];
	pIntergral[i] = fTemp;
}

for (int j=1; j<iHeight; j++)
{
	fTemp = 0;
	for (int i=0; i<iWidth; i++)
	{
		fTemp += pSrc[j*iWidth+i];
		pIntergral[j*iWidth+i] = pIntergral[(j-1)*iWidth+i] + fTemp;
	}
}
return true;

（5）改写OpenSURF的源代码

OpenSURF中的参数都是OpenCV的影像参数，如果想使用可以有两种方法：第一，需要将数据构造成为OpenCV的影像；第二，改造OpenSURF的代码，传递数据指针

//! Library function builds vector of described interest points
inline void surfDetDes(IplImage *img,  /* image to find Ipoints in */
                       std::vector<Ipoint> &ipts, /* reference to vector of Ipoints */
                       bool upright = false, /* run in rotation invariant mode? */
                       int octaves = OCTAVES, /* number of octaves to calculate */
                       int intervals = INTERVALS, /* number of intervals per octave */
                       int init_sample = INIT_SAMPLE, /* initial sampling step */
                       float thres = THRES /* blob response threshold */)

修改第一个参数IplImage *img为float *pImageBuffer

 

（6）在对原始影像以及参考影像进行SURF算子提取后，对两张影像的算子点进行匹配

//! Populate IpPairVec with matched ipts 
void getMatches(IpVec &ipts1, IpVec &ipts2, IpPairVec &matches)

（7）剔除粗差点，计算两张影像的仿射变换关系

放大显示：

因为是基于金字塔影像进行匹配，因此匹配点的坐标会存在一定的误差，计算得到的仿射变换关系也仅仅是两张相片的粗略位置关系。