【Datawhale】計算機視覺下 —— LBP特征描述算子

LBP特征描述算子

局部二值模型（Location Binary Pattern，LBP）是一種圖像紋理的描述算子，所以我們首先要知道什么是圖像的紋理特征，進而了解LBP算子的基本原理及其應用拓展。由於在原始的LBP提出后，研究人員還提出了各種改進方法，我們都將一一做介紹。最后使用opencv進行人臉識別。

圖像紋理特征

紋理是一種反映圖像中同質現象的視覺特征，它體現了物體表面的具有重復性和周期性變化的表面結構組織排列屬性。物體的紋理可以在其亮度和顏色上有所體現，而且幾乎所有圖像都包含了紋理信息，它可以表現這個物體的表面信息。有時候，它還能表達自身與外界環境的關系。

其中，紋理有三大特點，它們分別是：

1. 重復性：圖像可以看作是某種局部元素在全局區域的不斷重復出現。
2. 周期性：圖像中的元素並非隨機出現，而是按照一定的周期性重復出現。
3. 同質性：重復出現的元素在結構和尺寸上大致相同。

由上可見，紋理是某種局部序列性不斷重復、非隨機排列、在結構和尺寸上大致相同的統一體。紋理圖像示例如圖5-2所示，圖中第一行表示人工紋理，第二行是自然紋理。

紋理特征只是對物體表面特性進行描述，並不能反映物體的本質屬性，即圖像高層語義信息。

基本LBP

概念理解

LBP算子的基本思想是將中心像素的灰度值作為一個閾值，將其鄰域內的像素點灰度值進行比較，從而得到二進制編碼表示，來表示局部紋理特征。

鄰域的類型可分為四鄰域、D鄰域、八鄰域，四鄰域：該像素點的上下左右四個位置；

D領域：該像素點斜對角線上的四個相鄰位置。

八鄰域：四鄰域與D鄰域的並集。

基本的LBP算子考慮的是像素的八鄰域。

LBP表示方法有一個較為明顯的特點，它不容易收到圖像整體會讀線性變化的影響。也就是說，當圖像由於光線的影響使得整體灰度值發生線性均勻變化時，其LBP特征編碼是不變的。換句話說，它並不在意整體的灰度變化，而是關注像素之間的相對灰度改變。

例如，在某些情況下，陽光照射強度更低，導致拍攝圖像的整體亮度降低，但是實際上每個像素之間的差值仍然是固定的。那么在這種情況下，在圖片亮度對LBP特征編碼無影響。

主要步驟

使用一個\(3 \times 3\)的矩形形，處理待判斷像素點及其鄰域之間的關系。

1. 二值化過程

將像素點A的值與其八鄰域處的像素點逐一比較：

1. 如果A的像素值大於其鄰近點的像素值，則得到0。
2. 如果A的像素值小於其鄰近點的像素值，則得到1。

此過程以中心像素的灰度值為閾值，與鄰域的8個灰度值進行比較。

最后，根據順時針方向，將像素點A與周圍8個像素點比較所得到的0、1值連接起來，得到一個8位的二進制序列，然后將該二進制序列轉換為十進制數字作為點A的LBP值。

如下圖LBP原理示意圖所示，在左側\(3 \times 3\)的區域中，中心點的像素為76，並設置它為此次的閾值。然后現在我們對該中心點的8鄰域做進一步的處理。

1. 將中心點周圍的8個位置中灰度值大於76的像素點處理為1。例如，其鄰域中像素值為128、251、99、213的點，都被處理為1，填入對應的像素點位置上。
2. 將中心點周圍的8個位置中灰度值值小於76的像素點處理為0。例如，其鄰域中像素值為36、9、11、48的點，都被處理為0，填入對應的像素點位置上。
3. 最后得到的二值結果如右圖所示。

2. 中心點處理

完成上述的二值化過程之后，例如從像素點的正上方開始將得到的二值結果進行序列化，所以在上述例子中，二進制序列結果為“01011001”。

最后再將該二進制序列結果轉換為對應的十進制數“89”，作為當前中心點的像素值，如下圖所示。

3. 整合處理

上述過程僅僅是對一個像素點的處理過程，結果是一個LBP“編碼”。所以對於某幅圖像而言，需要進行逐行掃描完成每個像素點數值的更新。我們將采用分塊的形式進行編碼。其編碼過程如下：

假設此時有一幅100*100大小的圖像。

1. 在初始化時，將該圖划分為\(10 \times 10\) 個Block，其中每個Block的大小為\(10 \times 10\)。
2. 對每個Block的像素點提取其LBP特征，並建立一個計算某個“數字”出現的頻率統計直方圖。
3. 結束時，將生成\(10 \times 10\) 個統計直方圖，選擇性地對直方圖進行規范化處理。
4. 連接所有小塊的（規范化的）直方圖，整合后構成了整個窗口的特征向量，用來描述這幅圖片。
5. 得到特征向量之后，就可以使用各類算法對該圖像進行特定的處理了。

圓形鄰域的LBP算子

概念理解

基本LBP算子可以被進一步推廣到使用不同大小和形狀的鄰域。采用圓形的鄰域並結合雙線性插值運算使得我們可以獲得任意半徑和任意數目的鄰域像素點。該圓形鄰域可以用\(LBP^R_P\)表示，其中P表示圓形鄰域內參與運算的像素點個數，R表示鄰域的半徑。

計算方法

假設此時給出了一個半徑為2的8鄰域像素的圓形鄰域，圖中每個方格對應一個像素。

1. 處在方格中心的鄰域點（左、上、右、下4個黑點）：以該點所在方格的像素值作為它的值
2. 不在方格中心的鄰域點（斜45°方向的4個黑點）：線性插值法確定其值。

如下圖所示，若現在我們想計算左上角空心點的值，此時它並不在任何像素點內：

可以發現，在空心點的上方，距離\(2- \sqrt2\)處，有一個十叉點1，下方還有一個十叉點2。所以我們應首先分別計算這兩個十叉點1和2的水平插值。

其中點1的值根據與之處於同一行的\(I(i-2, j-2)\)以及\(I(i-2, j-2)\)的線形插值得到。

\[value(1) = I(i-2, j-2) + (2 - \sqrt2) \times (I(i-2, j-1) - I(i-2, j-2)) \]

同理計算出點2的值如下：

\[value(2) = I(i-1, j-2) + (2 - \sqrt2) \times (I(i-1, j-1) - I(i-1, j-2)) \]

再計算出點1和點2豎直的線性插值。

\[value = value(1) + (2 - \sqrt2) \times (value(2) - value(1)) \]

如此，計算出各鄰域像素點的灰度值之后，仍然使用基礎LBP算子的方法，得到二進制序列后轉換為十進制數等操作。

旋轉不變的LBP算子

根據上述介紹發現，像素點的二進制序列將由於圖像的旋轉而改變，因此它並不具有旋轉不變性。

那么，應該如何解決這樣的情況呢？

此后，Maenpaa等人又將 LBP 算子進行了擴展，提出了具有旋轉不變性的 LBP 算子，即不斷旋轉圓形鄰域得到一系列初始定義的 LBP 值，取其最小值作為該鄰域的 LBP 值。

如下圖所示，黑點代表0，白點代表1，假設初始的灰度值為255。經過7次旋轉，得到7個不同的二進制序列，並分別轉換為相應的十進制數，最后取到最小值“15”為最終中心點的LBP值，且序列為“00001111”。

統一化的LBP算子

理論基礎

由於LBP直方圖大多都是針對圖像中的各個分區分別計算的，對於一個普通大小的分塊區域，標准LBP算子得到的二進模式數目（LBP直方圖收集箱數目）較多，而實際的位於該分塊區域中的像素數目卻相對較少，這將會得到一個過於稀疏的直方圖，使得該直方圖失去統計的意義。

例如，在5×5鄰域內對20個像素點進行采樣，有\(2^20 = 1048576\)種二進制模式，數量過於龐大了。所以，我們想到了一種特殊的方法，來減少一些冗余的LBP模式。

研究者們提出了統一化模式（Uniform Patterns），再次改進了LBP算子的理論。

1. 跳變：

二進制序列中存在從1到0或者0到1的轉變，可以稱作是一次跳變。下面我們將舉例說明跳變次數的計算：

二進制序列	跳變次數
00000000 或 11111111	0次跳變
01000000 或 00010000	1次跳變
01010000	4次跳變

2. 統一化模式：

對於一個局部二進制模型而言，在將其二進制位串視為循環的情況下，如果其中包含的從0到1或者從1到0的轉變不多於2個，則稱為統一化模式。所以上例中的模式“01010000”就不屬於統一化模式。

3. 混合模式：

序列中包含的跳變為2次以上的，可以稱為混合模式。

統一化的意義

在隨后的LBP直方圖的計算過程中，只為統一化模式分配單獨的直方圖收集箱，而所有的非統一化模式都被放入一個公用收集箱中，使得LBP特征的數目大大減少。

一般來說，保留的統一化的模式往往是反映重要信息的那些模式，而那些非統一化模式中過多的轉變往往由噪聲引起，沒有良好的統計意義。

假設圖像分塊區域大小為\(18 \times 20\)，則像素的總數為360個。如果采用8鄰域像素的標准LBP算子，收集箱（特征）數目為\(2^8 = 256\)個，平均到每個收集箱的像素數目還不到2個。但是在統一化LBP算子的收集箱數目為59個（58個統一化模式收集箱加上1個非統一化模式收集箱），平均每個收集箱中將含有6個左右的像素點，因此更具有統計意義。

收集箱個數計算

將模式進行統一化后，實現了模式數量的降維，從之前的\(2^p\)轉換成為了\(p * (p - 1) + 2\)。筆者翻看了很多博客，並沒有詳細的說明最后的維度是怎么計算的，也就是我們收集箱的個數計算。那么這節我們就來討論一下。

接下來對8位二進制序列下的收集箱個數進行計算：

1. 0個轉變（2個）：

11111111，00000000。

**2. 1個轉變（7 x 2 = 14）： **01111111，00111111，00011111，00001111，00000111，00000011，00000001。

3. 2個轉變（42）：

• 6 x 2：01000000，00100000，00010000，00001000，00000100，00000010
• 5 x 2：01100000，00110000，00011000，00001100，00000110
• 4 x 2：01110000，00111000，00011100，00001110
• 3 x 2：01111000，00111100，00011110
• 2 x 2：01111100，00111110
• 1 x 2：01111110

值得注意的是，上面的“ \(K \times 2\) ”中的2表示有兩種情況，即反碼和原碼兩種，這樣便得到了58種統一化的編碼，比原來的256種減少了很多。

總結一下，經過上面的計算知道，這種統一化后的編碼個數可以用公式（4）表示。

\[T = P * (P - 1) + 2 \]

openCV實現人臉檢測

讀取圖片

這次仍然是把待讀取的照片放入img包中，圖片名為cv_3.jpeg。

import cv2
import matplotlib.pyplot as plt

filepath = "../img/cv_3.jpeg"

# 讀取圖片，路徑不能含有中文名，否則圖片讀取不出來
image = cv2.imread(filepath)
# 顯示圖片
plt.imshow('image', image)
plt.show()

這次的待檢測圖片：

加載級聯文件

我在這個地方是報錯了的，據說需要使用絕對路徑，但是隊里的其他人都不是像我這么玩的... 有點奇怪～

# 需要是Anaconda虛擬環境下的絕對路徑
cascade_path = "/Users/sonata/opt/anaconda3/share/opencv4/lbpcascades/lbpcascade_frontalface_improved.xml"  
# 下載到了本地使用

# 加載人臉級聯文件
faceCascade = cv2.CascadeClassifier(cascade_path)

人臉檢測

OpenCV給我們使用特征數據的方法：

def detectMultiScale(self, image, scaleFactor=None, minNeighbors=None, flags=None, minSize=None, maxSize=None)
params:
  1. scaleFactor: 指定每個圖像比例縮小多少圖像
  2. minNeighbors: 指定每個候選矩形必須保留多少個鄰居，值越大說明精度要求越高
  3. minSize：檢測到的最小矩形大小
  4. maxSize: 檢測到的最大矩形大小

所以我們使用此方法檢測圖片中的人臉

# 灰度轉換
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

# 人臉檢測
faces = faceCascade.detectMultiScale(gray, 1.1, 2, minSize=(100, 100))
print(faces)  # 識別的人臉信息

# 循環處理每一張臉
for x, y, w, h in faces:
    cv2.rectangle(img, pt1=(x, y), pt2=(x+w, y+h), color=[0, 0, 255], thickness=2)

img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
plt.imshow(img)
plt.show()

得到檢測結果如下：

處理臟數據

我們發現除了檢測到人臉數據，還有一些其他的臟數據，這個時候可以打印檢測出的人臉數據位置和大小結果如下：

[[ 306  357  177  177]
 [1264  410  116  116]
 [ 656  378  241  241]
 [1049  376  243  243]]

從大小中我們看到最大的兩個矩形，剛好是人臉數據，其余都是臟數據，那么繼續修改函數參數

faces = faceCascade.detectMultiScale(gray, 1.1, 2, minSize=(150, 150))
# 把最小矩形大小改成150

這樣我們就可以把臟數據給除去了：

C++實現人臉檢測

原始LBP算子

// 原始LBP算子只是計算8鄰域內的局部二值模式
template <typename _Tp> static
void olbp_(InputArray _src, OutputArray _dst) {
    // get matrices
    Mat src = _src.getMat();
    // allocate memory for result
    _dst.create(src.rows-2, src.cols-2, CV_8UC1);
    Mat dst = _dst.getMat();
    // zero the result matrix
    dst.setTo(0);
    // calculate patterns
    for(int i=1;i<src.rows-1;i++) 
    {
        for(int j=1;j<src.cols-1;j++)
        {
            _Tp center = src.at<_Tp>(i,j);
            unsigned char code = 0;
            code |= (src.at<_Tp>(i-1,j-1) >= center) << 7;
            code |= (src.at<_Tp>(i-1,j) >= center) << 6;
            code |= (src.at<_Tp>(i-1,j+1) >= center) << 5;
            code |= (src.at<_Tp>(i,j+1) >= center) << 4;
            code |= (src.at<_Tp>(i+1,j+1) >= center) << 3;
            code |= (src.at<_Tp>(i+1,j) >= center) << 2;
            code |= (src.at<_Tp>(i+1,j-1) >= center) << 1;
            code |= (src.at<_Tp>(i,j-1) >= center) << 0;
            dst.at<unsigned char>(i-1,j-1) = code;
        }
    }
}
// 外部接口，根據不同的數據類型調用模板函數
void cv::olbp(InputArray src, OutputArray dst) {
    switch (src.getMat().type()) {
    case CV_8SC1:   olbp_<char>(src,dst); break;
    case CV_8UC1:   olbp_<unsigned char>(src,dst); break;
    case CV_16SC1:  olbp_<short>(src,dst); break;
    case CV_16UC1:  olbp_<unsigned short>(src,dst); break;
    case CV_32SC1:  olbp_<int>(src,dst); break;
    case CV_32FC1:  olbp_<float>(src,dst); break;
    case CV_64FC1:  olbp_<double>(src,dst); break;
    default: 
        string error_msg = format("Using Original Local Binary Patterns for feature extraction only works                                     
        on single-channel images (given %d). Please pass the image data as a grayscale image!", type);
        CV_Error(CV_StsNotImplemented, error_msg);
        break;
    }
}
Mat cv::olbp(InputArray src) {
    Mat dst;
    olbp(src, dst);
    return dst;
}

圓形LBP算子

// src為輸入圖像，dst為輸出圖像，radius為半徑，neighbor為計算當前點LBP所需的鄰域像素點數，也就是樣本點個數
template <typename _Tp> static // 模板函數，根據不同的原始數據類型得到不同的結果
inline void elbp_(InputArray _src, OutputArray _dst, int radius, int neighbors) 
{
    //get matrices
    Mat src = _src.getMat();
    // allocate memory for result因此不用在外部給_dst分配內存空間，輸出數據類型都是int
    _dst.create(src.rows-2*radius, src.cols-2*radius, CV_32SC1);
    Mat dst = _dst.getMat();
    // zero
    dst.setTo(0);
    for(int n=0; n<neighbors; n++) 
    {
        // sample points 獲取當前采樣點
        float x = static_cast<float>(-radius) * sin(2.0*CV_PI*n/static_cast<float>(neighbors));
        float y = static_cast<float>(radius) * cos(2.0*CV_PI*n/static_cast<float>(neighbors));
        // relative indices 下取整和上取整
        int fx = static_cast<int>(floor(x)); // 向下取整
        int fy = static_cast<int>(floor(y));
        int cx = static_cast<int>(ceil(x));  // 向上取整
        int cy = static_cast<int>(ceil(y));
        // fractional part 小數部分
        float tx = x - fx;
        float ty = y - fy;
        // set interpolation weights 設置四個點的插值權重
        float w1 = (1 - tx) * (1 - ty);
        float w2 =      tx  * (1 - ty);
        float w3 = (1 - tx) *      ty;
        float w4 =      tx  *      ty;
        // iterate through your data 循環處理圖像數據
        for(int i=radius; i < src.rows-radius;i++) 
        {
            for(int j=radius;j < src.cols-radius;j++) 
            {
                // calculate interpolated value 計算插值，t表示四個點的權重和
                float t = w1*src.at<_Tp>(i+fy,j+fx) + 
w2*src.at<_Tp>(i+fy,j+cx) + 
w3*src.at<_Tp>(i+cy,j+fx) + 
w4*src.at<_Tp>(i+cy,j+cx);
                // floating point precision, so check some machine-dependent epsilon
                // std::numeric_limits<float>::epsilon()=1.192092896e-07F
                // 當t>=src(i,j)的時候取1，並進行相應的移位
                dst.at<int>(i-radius,j-radius) += ((t > src.at<_Tp>(i,j)) || 
                            (std::abs(t-src.at<_Tp>(i,j)) < std::numeric_limits<float>::epsilon())) << n;
            }
        }
    }
}

// 外部接口，根據不同的數據類型調用模板函數
static void elbp(InputArray src, OutputArray dst, int radius, int neighbors)
{
    int type = src.type();
    switch (type) {
    case CV_8SC1:   elbp_<char>(src,dst, radius, neighbors); break;
    case CV_8UC1:   elbp_<unsigned char>(src, dst, radius, neighbors); break;
    case CV_16SC1:  elbp_<short>(src,dst, radius, neighbors); break;
    case CV_16UC1:  elbp_<unsigned short>(src,dst, radius, neighbors); break;
    case CV_32SC1:  elbp_<int>(src,dst, radius, neighbors); break;
    case CV_32FC1:  elbp_<float>(src,dst, radius, neighbors); break;
    case CV_64FC1:  elbp_<double>(src,dst, radius, neighbors); break;
    default:
        string error_msg = format("Using Circle Local Binary Patterns for feature extraction only works                                    
         on single-channel images (given %d). Please pass the image data as a grayscale image!", type);
        CV_Error(CV_StsNotImplemented, error_msg);
        break;
    }
}
Mat cv::elbp(InputArray src, int radius, int neighbors) {
    Mat dst;
    elbp(src, dst, radius, neighbors);
    return dst;
}

總結

在這一節的學習中，比較困擾我的是LBP算子它對於人臉檢測過程中的用處，以及優化后的LBP算子是如何進行特征降維的，它的計算過程是什么。不過總算都是搞清楚了，也花了很多時間去理解，希望依舊能評個優秀～

參考資料

《數字圖像處理與python實現》

https://blog.csdn.net/saw009/article/details/80105871

https://blog.csdn.net/zouxy09/article/details/7929531

https://blog.csdn.net/dujian996099665/article/details/8886576