特征匹配之：采用CPU指令集SSE快速计算特征点描述子之间的欧氏\汉明距离

摘要：最近在做倾斜摄影测量中的倾斜影像精确匹配，采用vlSift(做了相关实验，OpenCV中的Sift算法实现不如VLfeat中的Sift算法实现，具体原因还不明白，希望知道的博友告知)对影像提点，之后进行特征点匹配(省略匹配的详细过程)。由于需要匹配的特征点的数据量很大故需要采用SSE技术更快速的计算128维描述子之间的欧氏距离。经过试验，在采用SSE后倾斜影像精确匹配的速度提高了2倍(在经过3周不断去更换搜索匹配、策略和试验后，能得到这个结果很欣喜。好幸运，时间和精力没有白费！)下面分享SSE计算的代码，并对代码做简单说明。

1.先看main函数

#include "stdafx.h"
#include <iostream>
#include <emmintrin.h>
#include <xmmintrin.h>
#include <vector>

int main(int argc, _TCHAR* argv[])
{
	SSEComputerCharDescriptorArray();  //Char类型的描述子

	SSEComputerFloatDescriptorArray();  //float类型的描述子

	SSEComputerBinaryDescriptorArray();  //二进制类型的描述子

	return 0;
}

2.128维char类型描述子的SSE计算(可修改为64维)

#define N 128
void SSEComputerCharDescriptorArray()
{
	   // 模拟输入的unsigned char*
	    unsigned char charArray_1[N], charArray_2[N];
	    for (int i = 0; i < N; i++)
	    {
		  charArray_1[i] = (unsigned char)i;
		  charArray_2[i] = (unsigned char)(i + 254);
           }

		unsigned char *offset_1, *offset_2;
		offset_1 = charArray_1;   //拿到描述子数组指针
		offset_2 = charArray_2 ;
	
        // 开始SSE处理
		__m128i indexJ3 = _mm_set_epi32(0, 0, 0, 0);

		register __m128i nullValue8 = _mm_set_epi8(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
		register __m128i nullValue16 = _mm_set_epi16(0, 0, 0, 0, 0, 0, 0, 0);

		register __m128i charSets_1, charSets_2;
		register __m128i src_1, src_2;  
		register __m128i subResult, mulResult, mulResult_1, mulResult_2;
		for (int i = 0; i < N / 8; i++, offset_1+= 8, offset_2 += 8)
		{
		  //每次load16个char，但后边解包时只用8个，故+8
		  // 每次8个char同时计算
		  charSets_1 = _mm_loadu_si128((const __m128i*) (offset_1));   
		  charSets_2 = _mm_loadu_si128((const __m128i*) (offset_2));

		  src_1 = _mm_unpacklo_epi8(charSets_1, nullValue8);
		  src_2 = _mm_unpacklo_epi8(charSets_2, nullValue8);

		  subResult = _mm_sub_epi16(src_2, src_1);    // 相减，产生16位的结果
		  mulResult = _mm_mullo_epi16(subResult, subResult);   // 求平方，产生32位的结果

		  /*为防止在后边累加求和时发生溢出，将平方结果由16位扩展到32位, 即变成每次4个数求和*/
		  mulResult_1 = _mm_unpacklo_epi16(mulResult, nullValue16);  
		  mulResult_2 = _mm_unpackhi_epi16(mulResult , nullValue16);

		  indexJ3 = _mm_add_epi32(indexJ3, mulResult_1);  // 累加计算
		  indexJ3 = _mm_add_epi32(indexJ3, mulResult_2);
		}  
		const __int32 *q = (const __int32 *)&indexJ3;
		int dest = *(q + 0) + *(q + 1) + *(q + 2) + *(q + 3);
		//输出结果
		std::cout << "dest = " << dest << std::endl;
}

具体，SSE函数不明白可以查Intel的官方说明，URL:
SSE函数用法的官方说明

3.128维float类型描述子的SSE计算(可修改为64维)

也直接上代码吧：

void SSEComputerFloatDescriptorArray()
{
	// 模拟输入的 float类型的描述子
	    float floatArray_1[N], floatArray_2[N];
	    for (int i = 0; i < N; i++)
	    {
		  floatArray_1[i] = i;
		  floatArray_2[i] = i + 2.0;
       }

     	// 开始SSE处理
		float *offset_1, *offset_2;
		offset_1 = floatArray_1;   // 获得描述子的指针
		offset_2 = floatArray_2;
		
		__m128 squareSum = _mm_setzero_ps(); // 累加器初始化为0

		register __m128 subResult, squareResult;   // 相减和平方结果的临时变量
		register __m128 floatSets_1, floatSets_2; //每4个float数一组，同时计算

		for (int i = 0; i < N / 4; i++, offset_1 += 4, offset_2 += 4)
		{
			floatSets_1 = _mm_loadu_ps(offset_1);  //将每4个float数导入成一个__m128类型
			floatSets_2 = _mm_loadu_ps(offset_2);
			
			subResult = _mm_sub_ps( floatSets_1, floatSets_2); // 4个float同时求差
			squareResult = _mm_mul_ps (subResult, subResult); // 4个float同时平方，此处可能会数据溢出
			squareSum = _mm_add_ps(squareSum, squareResult);  // 累加计算
		}  
		// 将两个描述子之间的欧氏距离的平方！，，，，是平方，输出
		const float *q = (const float *)&squareSum;
		double  dist = *(q + 0) + *(q + 1) + *(q + 2) + *(q + 3);
		std::cout << "dest = " << dist << std::endl;
}

另外，用SSE时load数据比set数据快，在计算char类型描述子的距离时，我刚开始是set数据，最后的结果是SSE计算非但没快，反而更慢了!!!

3.128维二进制类型描述子的SSE计算(可修改为64维)

此处计算的是汉明距离，在surf特征点匹配中用到多。上代码：

void SSEComputerBinaryDescriptorArray()
{
	/*计算两个描述子之间的汉明距离
	   汉明距离：两个描述子的二进制对应位不相同的个数，如 (十进制数52)0 0 1 1  0 1 0 0 与 (十进制数38)0 0 1 0  0 1 1 0 的汉明距离为2*/

	// 汉明距离查表数组     
	static const unsigned char popCountTable[] =
	{
		0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
		1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
		1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
		2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
		1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
		2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
		2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
		3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, 4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8
	};    // 这段你可以完整复制，我也是复制自OpenCV，hahahaha..............

       // 相同的套路，模拟两个128维度的描述子
	unsigned char charArray_1[N], charArray_2[N];
	for (int i = 0; i < N; i++)
	{
		charArray_1[i] = (unsigned char)0;
		charArray_2[i] = (unsigned char)2;
	}

	unsigned char *offset_1, *offset_2;
	offset_1 = charArray_1;
	offset_2 = charArray_2;
	register __m128i charSets_1, charSets_2;
	register __m128i tempDist;

	int dist = 0;
	for (int i = 0; i < N / 16; i++, offset_1 += 16, offset_2 += 16)
	{
		charSets_1 = _mm_loadu_si128((const __m128i*)offset_1);
		charSets_2 = _mm_loadu_si128((const __m128i*)offset_2);

		tempDist = _mm_xor_si128(charSets_1, charSets_2);
		const __int8 *q = (const __int8 *)&tempDist;
		// 查表求距离
		dist += popCountTable[*q]  + popCountTable[*(q + 1)]  + popCountTable[*(q + 2)]  + popCountTable[*(q + 3)] + popCountTable[*(q + 4)]
			        + popCountTable[*(q + 5)]  + popCountTable[*(q + 6)]  + popCountTable[*(q + 7)] + popCountTable[*(q + 8)]
			       + popCountTable[*(q + 9)] + popCountTable[*(q + 10)]  + popCountTable[*(q + 11)] + popCountTable[*(q + 12)]
			       + popCountTable[*(q + 13)]  + popCountTable[*(q + 14)]  + popCountTable[*(q + 15)];
	}
	// 输出两个描述子之间的汉明距离
	std::cout << "dest = " << dist << std::endl;
}

至于如何将SIFT/Surf的float类型的描述子转为char类型(char类型运算代价低，有助于提高特征匹配的速度)，我们下篇见。

分享到此结束，希望对你有所帮助！