【工程應用五】 opencv中linemod模板匹配演算法諸多疑惑和自我解讀。
- 2022 年 3 月 22 日
- 筆記
- ，邊緣特徵, [03] 影像演算法優化, linemod, 模板匹配
研究這個前前後後也有快兩三個月了，因為之前也一直在弄模板匹配方面的東西，所以偶爾還是有不少朋友諮詢或者問你有沒有研究過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 = 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 = response_maps[ori].ptr<__m128i>();
__m128i* lsb4_data = lsb4.ptr<__m128i>();
__m128i* msb4_data = 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 = response_maps[ori].ptr<uchar>();
uchar* lsb4_data = lsb4.ptr<uchar>();
uchar* msb4_data = msb4.ptr<uchar>();
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角度高速模板匹配。
