研究這個前前后后也有快兩三個月了，因為之前也一直在弄模板匹配方面的東西，所以偶爾還是有不少朋友咨詢或者問你有沒有研究過linemod這個演算法啊，那個效率啥的還不錯啊，有段時間一直不以為然，覺得我現在用的那個匹配因該很不錯的，沒必要深究了，后來呢，還是忍不住手癢，把論文打出來看了看，又找了點資料研究了下，結果沒想到一弄又是兩個月過去了，中間也折騰了很久，浪費了不少時間，總算還是有點識訓，稍微整理下做個交流，
至于linemod的具體數學原理，我也不需要詳談，畢竟論文和opencv的代碼就擺在那里， github上也有一些別人改進的版本，
我就覺得啊，linemod這個基于計算邊緣的模板匹配啊，他使用的是選中的特征點的梯度的角度方向作為特征，而不是梯度的值，而后計算模板和測驗不同位置角度的余弦的絕對值，這個都是常規的操作， 作者把這個角度線性量化為一些特定的值，這個本質上呢降低了演算法的精度，但是由于特征點較多，基本不會影響識別結果， 關鍵是這個量化啊，能夠帶來很多很多的好處，有些真的是意想不到，
論文里呢把360的角度量化為8個值，即以45度為間隔，分別用整數0、1、2、3、4、5、6、7表示，這樣呢不同的兩個角度之間的差異絕對值呢，只有0、1、2、3、4這5種可能，分別對應5個得分，比如模板的某個特征點的角度為210度，則量化值為4，目標中某個位置的角度值為52度，則量化值為1，這樣角度之間差異值為3，則對應的得分為1，
接著論文里說為了減少微小的變形引起的識別誤差，建議將量化后的值進行擴散，這個擴散也是設計的非常有技巧，很有意思，充分利用了或運算的優異特性，
后續說為了減少計算量呢，可以提前計算出8個回應圖，這樣匹配計算時就可以直接查表，而無需實時計算，
再后續還有一個線性化記憶體，算了，我已經沒看那個了，到前面這一步就已經打止了，因為我已經開始編程了，
第一步呢，我就是在考慮演算法的優化問題，我看了下opencv的代碼，寫的很好，又很不好，讓你讀的很難受，但是寫的確實穩健，考慮到了很多不同的硬體配置，這也許就是大工程的特性吧，
角度量化的問題和代碼方面我不想提，也有很多可優化的地方，大家可自行考慮，
在談到提速之前，我說一個重點，那就是所謂的梯度擴散、計算回應圖都是在查找模板時進行的， 對不同的圖都要有重新計算，而不是離線玩，所以這里的每個耗時，都和檢測速度有關，
第一：那個梯度擴散，CV的代碼有下面這一大堆：
/****************************************************************************************\
* Response maps *
\****************************************************************************************
static void orUnaligned8u(const uchar * src, const int src_stride,
uchar * dst, const int dst_stride,
const int width, const int height)
{
#if CV_SSE2
volatile bool haveSSE2 = checkHardwareSupport(CPU_SSE2);
#if CV_SSE3
volatile bool haveSSE3 = checkHardwareSupport(CPU_SSE3);
#endif
bool src_aligned = reinterpret_cast<unsigned long long>(src) % 16 == 0;
#endif
for (int r = 0; r < height; ++r)
{
int c = 0;
#if CV_SSE2
// Use aligned loads if possible
if (haveSSE2 && src_aligned)
{
for ( ; c < width - 15; c += 16)
{
const __m128i* src_ptr = reinterpret_cast<const __m128i*>(src + c);
__m128i* dst_ptr = reinterpret_cast<__m128i*>(dst + c);
*dst_ptr = _mm_or_si128(*dst_ptr, *src_ptr);
}
}
#if CV_SSE3
// Use LDDQU for fast unaligned load
else if (haveSSE3)
{
for ( ; c < width - 15; c += 16)
{
__m128i val = _mm_lddqu_si128(reinterpret_cast<const __m128i*>(src + c));
__m128i* dst_ptr = reinterpret_cast<__m128i*>(dst + c);
*dst_ptr = _mm_or_si128(*dst_ptr, val);
}
}
#endif
// Fall back to MOVDQU
else if (haveSSE2)
{
for ( ; c < width - 15; c += 16)
{
__m128i val = _mm_loadu_si128(reinterpret_cast<const __m128i*>(src + c));
__m128i* dst_ptr = reinterpret_cast<__m128i*>(dst + c);
*dst_ptr = _mm_or_si128(*dst_ptr, val);
}
}
#endif
for ( ; c < width; ++c)
dst[c] |= src[c];
// Advance to next row
src += src_stride;
dst += dst_stride;
}
}
/**
* \brief Spread binary labels in a quantized image.
*
* Implements section 2.3 "Spreading the Orientations."
*
* \param[in] src The source 8-bit quantized image.
* \param[out] dst Destination 8-bit spread image.
* \param T Sampling step. Spread labels T/2 pixels in each direction.
*/
static void spread(const Mat& src, Mat& dst, int T)
{
// Allocate and zero-initialize spread (OR'ed) image
dst = Mat::zeros(src.size(), CV_8U);
// Fill in spread gradient image (section 2.3)
for (int r = 0; r < T; ++r)
{
int height = src.rows - r;
for (int c = 0; c < T; ++c)
{
orUnaligned8u(&src.at<unsigned char>(r, c), static_cast<int>(src.step1()), dst.ptr(),
static_cast<int>(dst.step1()), src.cols - c, height);
}
}
}
我翻譯成我容易接受的代碼，并且剔除一些對硬體環境的判斷的陳述句，如下所示：
void IM_Spread_Ori(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int T)
{
memset(Dest, 0, Height * Stride);
for (int J = 0; J < T; J++)
{
int H = Height - J;
for (int I = 0; I < T; I++)
{
int W = Width - I;
int BlockSize = 16, Block = W / BlockSize;
unsigned char *SrcP = Src + J * Stride + I;
unsigned char *DestP = Dest;
for (int Y = 0; Y < H; Y++)
{
for (int X = 0; X < Block * BlockSize; X += BlockSize)
{
__m128i SrcV = _mm_loadu_si128((__m128i*)(SrcP + X));
__m128i DstV = _mm_loadu_si128((__m128i*)(DestP + X));
_mm_storeu_si128((__m128i *)(DestP + X), _mm_or_si128(SrcV, DstV));
}
for (int X = Block * BlockSize; X < W; X++)
{
DestP[X] |= SrcP[X];
}
SrcP += Stride;
DestP += Stride;
}
}
}
}
要說啊，這個代碼本身來說有是個比較高效的代碼了，但是，我一想到論文中的T=8，那就意味著差不多是8*8=64次全圖這樣的資料or操作，哪怕就算or操作再快， 這個也不太可能快過 64次memcpy的，特別當一個圖比較大的時候，這個就有點明顯了，我測驗了下，對于一個3000*3000的灰度圖（工業上遇到這么大的圖應該不算離譜吧），初步測驗了下，居然需要大概200ms的時間，對于模板匹配這種需要高頻操作的需求來說，單獨這一步的耗時還是大了點，
有么有在不改變效果的情況下，進一步提高這個演算法的方法呢，其實是有的，我們知道or操作時不分前后順序的，即多個資料or可以隨便誰和誰先操作，因此，我們可以安排T行之間先or，然后再對結果記性T列之間or操作，這樣則只需要2*T次or計算，而且有一些額外的好處就是避免很多cache miss，這是隱藏的速度提升，
改寫后的代碼如下所示：
void IM_Spread(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int T)
{
// 利用了或運算的行和列分離特性
memset(Dest, 0, Height * Stride);
int BlockSize = 16, Block = Width / BlockSize;
for (int Y = 0; Y < Height; Y++)
{
unsigned char *LinePS = Src + Y * Stride;
unsigned char *LinePD = Dest + Y * Stride;
for (int J = 0; J < ((Y + T > Height) ? Height - Y : T); J++)
{
// 高度方向依次向下進行T次或運算
for (int X = 0; X < Block * BlockSize; X += BlockSize)
{
__m128i SrcV = _mm_loadu_si128((__m128i*)(LinePS + X));
__m128i DstV = _mm_loadu_si128((__m128i*)(LinePD + X));
_mm_storeu_si128((__m128i *)(LinePD + X), _mm_or_si128(SrcV, DstV));
}
for (int X = Block * BlockSize; X < Width; X++)
{
LinePD[X] |= LinePS[X];
}
LinePS += Stride; // 源資料向下移動，目標資料不動
}
}
BlockSize = 16, Block = (Width - T) / BlockSize;
for (int Y = 0; Y < Height; Y++)
{
unsigned char *LinePD = Dest + Y * Stride;
for (int X = 0; X < Block * BlockSize; X += BlockSize)
{
__m128i Value = _mm_setzero_si128();
for (int I = 0; I < T; I++)
{
__m128i SrcV = _mm_loadu_si128((__m128i*)(LinePD + X + I));
Value = _mm_or_si128(Value, SrcV);
}
// 這個讀取和寫入是沒有重疊的，所以可以不利用中間記憶體
_mm_storeu_si128((__m128i *)(LinePD + X), Value);
}
for (int X = Block * BlockSize; X < Width; X++)
{
int Value = https://www.cnblogs.com/Imageshop/p/0;
for (// 此處洗掉部分代碼，供讀者自行補充 )
{
Value = Value | LinePD[X + I];
}
LinePD[X] = Value;
}
}
}
同樣的測驗圖，同樣的T=8，速度一下子提升到了45ms左右，有近5倍的速度提升，
為什么我會分享這段代碼呢，因為后面我發現他根本沒什么卵用，
第二：那個計算回應圖的代碼，也可以繼續優化，
// Auto-generated by create_similarity_lut.py
CV_DECL_ALIGNED(16) static const unsigned char SIMILARITY_LUT[256] = {0, 4, 3, 4, 2, 4, 3, 4, 1, 4, 3, 4, 2, 4, 3, 4, 0, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 0, 3, 4, 4, 3, 3, 4, 4, 2, 3, 4, 4, 3, 3, 4, 4, 0, 1, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 0, 2, 3, 3, 4, 4, 4, 4, 3, 3, 3, 3, 4, 4, 4, 4, 0, 2, 1, 2, 0, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 0, 3, 2, 3, 1, 3, 2, 3, 0, 3, 2, 3, 1, 3, 2, 3, 0, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 0, 4, 3, 4, 2, 4, 3, 4, 1, 4, 3, 4, 2, 4, 3, 4, 0, 1, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 0, 3, 4, 4, 3, 3, 4, 4, 2, 3, 4, 4, 3, 3, 4, 4, 0, 2, 1, 2, 0, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 0, 2, 3, 3, 4, 4, 4, 4, 3, 3, 3, 3, 4, 4, 4, 4, 0, 3, 2, 3, 1, 3, 2, 3, 0, 3, 2, 3, 1, 3, 2, 3, 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4};
/**
* \brief Precompute response maps for a spread quantized image.
*
* Implements section 2.4 "Precomputing Response Maps."
*
* \param[in] src The source 8-bit spread quantized image.
* \param[out] response_maps Vector of 8 response maps, one for each bit label.
*/
static void computeResponseMaps(const Mat& src, std::vector<Mat>& response_maps)
{
CV_Assert((src.rows * src.cols) % 16 == 0);
// Allocate response maps
response_maps.resize(8);
for (int i = 0; i < 8; ++i)
response_maps[i].create(src.size(), CV_8U);
Mat lsb4(src.size(), CV_8U);
Mat msb4(src.size(), CV_8U);
for (int r = 0; r < src.rows; ++r)
{
const uchar* src_r = src.ptr(r);
uchar* lsb4_r = lsb4.ptr(r);
uchar* msb4_r = msb4.ptr(r);
for (int c = 0; c < src.cols; ++c)
{
// Least significant 4 bits of spread image pixel
lsb4_r[c] = src_r[c] & 15;
// Most significant 4 bits, right-shifted to be in [0, 16)
msb4_r[c] = (src_r[c] & 240) >> 4;
}
}
#if CV_SSSE3
volatile bool haveSSSE3 = checkHardwareSupport(CV_CPU_SSSE3);
if (haveSSSE3)
{
const __m128i* lut = reinterpret_cast<const __m128i*>(SIMILARITY_LUT);
for (int ori = 0; ori < 8; ++ori)
{
__m128i* map_data = https://www.cnblogs.com/Imageshop/p/response_maps[ori].ptr<__m128i>();
__m128i* lsb4_data = https://www.cnblogs.com/Imageshop/p/lsb4.ptr<__m128i>();
__m128i* msb4_data = https://www.cnblogs.com/Imageshop/p/msb4.ptr<__m128i>();
// Precompute the 2D response map S_i (section 2.4)
for (int i = 0; i < (src.rows * src.cols) / 16; ++i)
{
// Using SSE shuffle for table lookup on 4 orientations at a time
// The most/least significant 4 bits are used as the LUT index
__m128i res1 = _mm_shuffle_epi8(lut[2*ori + 0], lsb4_data[i]);
__m128i res2 = _mm_shuffle_epi8(lut[2*ori + 1], msb4_data[i]);
// Combine the results into a single similarity score
map_data[i] = _mm_max_epu8(res1, res2);
}
}
}
else
#endif
{
// For each of the 8 quantized orientations...
for (int ori = 0; ori < 8; ++ori)
{
uchar* map_data = https://www.cnblogs.com/Imageshop/p/response_maps[ori].ptr();
uchar* lsb4_data = https://www.cnblogs.com/Imageshop/p/lsb4.ptr();
uchar* msb4_data = https://www.cnblogs.com/Imageshop/p/msb4.ptr();
const uchar* lut_low = SIMILARITY_LUT + 32*ori;
const uchar* lut_hi = lut_low + 16;
for (int i = 0; i < src.rows * src.cols; ++i)
{
map_data[i] = std::max(lut_low[ lsb4_data[i] ], lut_hi[ msb4_data[i] ]);
}
}
}
}
看上去又是一大堆代碼，簡化后如下所示：
void IM_ComputeResponseMaps_Slow(unsigned char *Src, unsigned char **ResponseMaps, int Width, int Height)
{
static const unsigned char SIMILARITY_LUT[256] = { 0, 4, 3, 4, 2, 4, 3, 4, 1, 4, 3, 4, 2, 4, 3, 4, 0, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 0, 3, 4, 4, 3, 3, 4, 4, 2, 3, 4, 4, 3, 3, 4, 4, 0, 1, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 0, 2, 3, 3, 4, 4, 4, 4, 3, 3, 3, 3, 4, 4, 4, 4, 0, 2, 1, 2, 0, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 0, 3, 2, 3, 1, 3, 2, 3, 0, 3, 2, 3, 1, 3, 2, 3, 0, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 0, 4, 3, 4, 2, 4, 3, 4, 1, 4, 3, 4, 2, 4, 3, 4, 0, 1, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 0, 3, 4, 4, 3, 3, 4, 4, 2, 3, 4, 4, 3, 3, 4, 4, 0, 2, 1, 2, 0, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 0, 2, 3, 3, 4, 4, 4, 4, 3, 3, 3, 3, 4, 4, 4, 4, 0, 3, 2, 3, 1, 3, 2, 3, 0, 3, 2, 3, 1, 3, 2, 3, 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4 };
unsigned char *lsb4 = (unsigned char *)malloc(Width * Height);
unsigned char *msb4 = (unsigned char *)malloc(Width * Height);
for (int Y = 0; Y < Height * Width; Y++)
{
lsb4[Y] = Src[Y] & 15;
msb4[Y] = (Src[Y] & 240) >> 4;
}
int BlockSize = 16, Block = (Width * Height) / BlockSize;
for (int Z = 0; Z < 8; Z++)
{
for (int Y = 0; Y < Block * BlockSize; Y += BlockSize)
{
__m128i Res1 = _mm_shuffle_epi8(_mm_loadu_si128((__m128i *)(SIMILARITY_LUT + 32 * Z + 0)), _mm_loadu_si128((__m128i *)(lsb4 + Y)));
__m128i Res2 = _mm_shuffle_epi8(_mm_loadu_si128((__m128i *)(SIMILARITY_LUT + 32 * Z + 16)), _mm_loadu_si128((__m128i *)(msb4 + Y)));
_mm_storeu_si128((__m128i *)(ResponseMaps[Z] + Y), _mm_max_epu8(Res1, Res2));
}
for (int Y = Block * BlockSize; Y < Width * Height; Y++)
{
ResponseMaps[Z][Y] = IM_Max(SIMILARITY_LUT[lsb4[Y] + 32 * Z], SIMILARITY_LUT[msb4[Y] + 32 * Z + 16]);
}
}
free(lsb4);
free(msb4);
}
同樣3000*3000的測驗圖，這個函式平均耗時80ms，也算是非常快的，
那有沒有改進空間呢，其實是有的，下面這兩句明顯也是可以用SIMD指令優化的嘛，&操作直接由對應的_mm_loadu_si128指令，至于byte型別資料的移位，確實沒有直接指令可以使用，但是自己寫個又有什么難度呢.
lsb4[Y] = Src[Y] & 15;
msb4[Y] = (Src[Y] & 240) >> 4;
比如這樣就可以：
// 無符號位元組資料右移四位
inline __m128i _mm_srli4_epu8(__m128i v)
{
v = _mm_srli_epi16(v, 4);
v = _mm_and_si128(v, _mm_set1_epi8(0x0f));
return v;
}
好了，其他的代碼就不需要我寫了吧，優化后這個速度能夠提高到50ms，
我從這個代碼里最大的識訓不是其他的，就是_mm_shuffle_epi8這個陳述句，利用這個很巧妙的實作了一個查找表的程序，其實我想起來了，在我的博文【演算法隨記七】巧用SIMD指令實作急速的位元組流按位反轉演算法， 里就已經使用了這個技巧，他能輕松的實作少于16個元素的位元組型別的查找表，而且效率比普通的C語言查表方式不知道高了多少倍，后面文章我們還會說道這個指令的一個更為優異的特性，
說到這里，大家也許會認為我會繼續談下后續的線性記憶體方面的優化，但是可惜了，后面的我就沒有看了，因為我覺得到了這一步，我就已經有了我自己的路可以走了，不需要后續的那個東西，那個可能對我來說還是累贅，于是我耐著性子，在我以前的大框架的基礎上，修改區域函式，終于能跑出了初步的效果， 比如下圖，我們取T=8時，得到的匹配結果如下：
整個的基本都錯位了，
而當我們的T取為3時，結果就比較好了，
但是仔細觀察，可以明顯的發現目標還是會有一到2個像素的偏移，如果我們T取值為1時，結果如下：
整體的準確度明顯有所提高，為了證明這個結果，在T=1時，我還測驗了很多其他影像，結果都表面是完全的，比如一下幾幅影像：
這是兩幅很具有代表性的測驗影像，一個是高重疊，第二個是強邊緣，可以看到效果很準確，
T=1意味這什么，即不需要進行擴散，所以我們由此沒有理由不懷疑論文里結論的準確性，為此，我嘗試分析深層次的原因，
我的測驗方法是：從一副影像中剪切一塊小的影像，然后直接按照上述代碼計算這個小塊影像在原圖中的各個不同位置的得分，注意到我們得分通常是按照計算后的總分值除以4*特征點的數量，
具體的測驗代碼詳見附件工程： 擴散結果驗證（生成的RAW影像保存D盤根目錄下）， 這里對算了做一些簡單的簡化，但不影響實質，
以下面兩幅圖為例：
小圖在大圖中的準確坐標位置為【 104，76】，
為了顯示方便，我們把計算得到的得分值，量化到0和32768之間，然后保存為RAW檔案，這樣就可以用PS打開查看其像素結果了，當T=1時，可得到如下視覺結果（下左圖），明顯圖中有個最亮的點，我們放大后查看其結果(下右圖），
可以看到準確的定位到了104,76這個坐標，
當我們選擇T=3時，同樣的結果如下面所示：
我們可以看到，坐標104,76此處依舊是最大值，但是由于擴散的作用，使得周邊也出現了同樣的最大值，由于PS只能同時有四個取樣點，所以實際上還不止，如果把T擴大到8，那么將會有更多的同樣的最大值，這就導致了一個問題，程式無法確定這些最大值那個才是真正的準確的位置了，而又必須確定一個，否則無法得到最終的結果，因此，在上述的T=8的匹配中，由于程式無原則的去了一個最值，然后逐層像金字塔上層傳遞，每層傳遞都有可能出現錯誤，所以導致最終誤差越來越大，
那么結論來了，要想準確的匹配，根本就不需要擴散程序，
好希望我的結論是錯誤的啊，本來還想用這擴散來解決模型的建立慢的問題，以及實作可同時檢查帶有縮放和旋轉的匹配呢，可惜，都是泡影，
雖然如此，但是這個演算法還是有很好的價值的，下一篇文章將講述基于T=1時改演算法的進一步擴展和優化，以及如何實作更高效率的演算法效果，先分享一個測驗工具了：16角度高速模板匹配，
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/445899.html
標籤：其他
上一篇：開篇_Blog