1. 程式人生 > >cv2實現基於粒子濾波的目標跟蹤

cv2實現基於粒子濾波的目標跟蹤

目標跟蹤過程分為2部分,即目標特徵提取和目標跟蹤演算法。

      其中目標特徵提取又包括以下幾種:1. 各種色彩空間直方圖,利用色彩空間的直方圖分佈作為目標跟蹤的特徵,可以減少物體遠近距離的影響,因為其顏色分佈大致相同。2.輪廓特徵,提取目標的輪廓特徵,可以加快演算法的速度,且可以在目標有小部分影響的情況下同樣有效果。3. 紋理特徵,如果被跟蹤目標是有紋理的,則根據其紋理特徵來跟蹤效果會有所改善。

     目標跟蹤演算法目前大概分為以下4種:1. 基於meanshift演算法,即利用meanshift演算法可以快速找到領域目標最相似的地方,效果還不錯,但是其只能找到區域性最大值,且不能解決遮擋問題以及不能自適應跟蹤目標的形狀,方向等。其後面有學者對其做了改進,比如說camshift,就可以自適應物體的大小,方向,具有較好的跟蹤效果。2. Kalman濾波的思想,該思想是利用物體的運動模型來,即服從高斯模型,來對目標狀態進行預測,然後與觀察模型進行比較,根據2者之間的誤差來尋找運動目標的狀態,但是該演算法的精度不高,因為其高斯運動模型在現實生活中很多條件下並得不到滿足,並且該演算法對雜亂的背景也很敏感。3. 基於粒子濾波的思想,每次通過實驗可以重取樣粒子的分佈,根據該分佈對粒子進行擴散,然後通過擴散的結果來觀察目標的狀態,最後更新目標的狀態。該演算法最大的特點是跟蹤速度快,且能解決一部分遮擋問題,在實際應用過程中越來越多。4.基於目標建模的方法。該方法具有一定的針對性,需要提前知道所需跟蹤的目標是什麼,比如說車輛,人臉,行人等。由於已經知道了跟蹤目標,所以必須對目標進行建模,然後利用該模型來進行跟蹤。該方法的侷限性是必須提前知道所跟蹤的目標是什麼,因而其推廣性比較差。

 實現過程

  1)該階段要人工指定跟蹤目標,程式計算跟蹤目標的特徵,比如可以採用目標的顏色特徵。具體到Rob Hess的程式碼,開始時需要人工用滑鼠拖動出一個跟蹤區域,然後程式自動計算該區域色調(Hue)空間的直方圖,即為目標的特徵。直方圖可以用一個向量來表示,所以目標特徵就是一個N*1的向量V。

  • 2)搜尋階段-放狗
    好,我們已經掌握了目標的特徵,下面放出很多條狗,去搜索目標物件,這裡的狗就是粒子particle。狗有很多種放法。比如,a)均勻的放:即在整個影象平面均勻的撒粒子(uniform distribution);b)在上一幀得到的目標附近按照高斯分佈來放,可以理解成,靠近目標的地方多放,遠離目標的地方少放。Rob Hess的程式碼用的是後一種方法。狗放出去後,每條狗怎麼搜尋目標呢?就是按照初始化階段得到的目標特徵(色調直方圖,向量V)。每條狗計算它所處的位置處影象的顏色特徵,得到一個色調直方圖,向量Vi,計算該直方圖與目標直方圖的相似性。相似性有多種度量,最簡單的一種是計算sum(abs(Vi-V)).每條狗算出相似度後再做一次歸一化,使得所有的狗得到的相似度加起來等於1.
  • 3)決策階段
    我們放出去的一條條聰明的狗向我們發回報告,“一號狗處影象與目標的相似度是0.3”,“二號狗處影象與目標的相似度是0.02”,“三號狗處影象與目標的相似度是0.0003”,“N號狗處影象與目標的相似度是0.013”…那麼目標究竟最可能在哪裡呢?我們做次加權平均吧。設N號狗的影象畫素座標是(Xn,Yn),它報告的相似度是Wn,於是目標最可能的畫素座標X = sum(Xn*Wn),Y = sum(Yn*Wn).
  • 4)重取樣階段Resampling
    既然我們是在做目標跟蹤,一般說來,目標是跑來跑去亂動的。在新的一幀影象裡,目標可能在哪裡呢?還是讓我們放狗搜尋吧。但現在應該怎樣放狗呢?讓我們重溫下狗狗們的報告吧。“一號狗處影象與目標的相似度是0.3”,“二號狗處影象與目標的相似度是0.02”,“三號狗處影象與目標的相似度是0.0003”,“N號狗處影象與目標的相似度是0.013”…綜合所有狗的報告,一號狗處的相似度最高,三號狗處的相似度最低,於是我們要重新分佈警力,正所謂好鋼用在刀刃上,我們在相似度最高的狗那裡放更多條狗,在相似度最低的狗那裡少放狗,甚至把原來那條狗也撤回來。這就是Sampling Importance Resampling,根據重要性重取樣(更具重要性重新放狗)。
    (2)->(3)->(4)->(2)如是反覆迴圈,即完成了目標的動態跟蹤。

示例程式碼:

import cv2
import numpy as np


# 權重計算函式
def likelihood(x, y, func, image, w=30, h=30):
    x1 = np.int32(max(0, x - w / 2))
    y1 = np.int32(max(0, y - h / 2))
    x2 = np.int32(min(image.shape[1], x + w / 2))
    y2 = np.int32(min(image.shape[0], y + h / 2))
    # 擷取圖片區域
    region = image[y1: y2, x1: x2]
    # 統計符合色值的畫素點個數
    count = region[func(region)].size
    return (float(count) / image.size) if count > 0 else 0.0001


# 粒子(particle)的初始化函式
def init_particle(func, image):
    mask = image.copy()
    mask[func(mask) == False] = 0
    # 查詢物體輪廓
    _, contours, _= cv2.findContours(mask, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
    if len(contours) <=0:
        return None
    # 面積最大輪廓
    max_contour = max(contours, key=cv2.contourArea)
    # 獲取外部矩形邊界(返回值:返回四個值,分別是x, y, w, h;x, y是矩陣左上方的座標,w,h是矩陣的寬和高
    max_rect = np.array(cv2.boundingRect(max_contour))
    max_rect = max_rect[: 2] + max_rect[2:] / 2
    weight = likelihood(max_rect[0], max_rect[1], func, image)
    particles = np.ndarray((500, 3), dtype=np.float32)
    particles[:] = [max_rect[0], max_rect[1], weight]
    return particles


# 根據粒子的權重的後驗概率分佈重新取樣
def resample(particles):
    tmp_particles = particles.copy()
    # 元素相加
    weights = particles[:, 2].cumsum()
    last_weight = weights[weights.shape[0] - 1]
    for i in range(particles.shape[0]):
        weight = np.random.rand() * last_weight
        # 只返回第一次出現的關係表示式為真的索引
        particles[i] = tmp_particles[(weights > weight).argmax()]
        particles[i][2] = 1.0


# 預測
def predict(particles, variance=13.0):
    # np.random.randn:以給定的形狀建立一個數組,陣列元素符合標準正態分佈N(0,1)
    # 13位半徑,向外隨機擴散
    particles[:, 0] += np.random.randn((particles.shape[0])) * variance
    particles[:, 1] += np.random.randn((particles.shape[0])) * variance


# 權重處理
def weight(particles, func, image):
    for i in range(particles.shape[0]):
        particles[i][2] = likelihood(particles[i][0], particles[i][1], func, image)
        # 權重相加
        sum_weight = particles[:, 2].sum()
        particles[:, 2] *= (particles.shape[0] / sum_weight)


# 測定座標
def measure(particles):
    x = (particles[:, 0] * particles[:, 2]).sum()
    y = (particles[:, 0] * particles[:, 2]).sum()
    weight = particles[:, 2].sum()
    return x / weight, y / weight


# 粒子(particle)座標獲取
particle_filter_cur_frame = 0
def particle_filter(particles, func, image, max_frame=10):
    global particle_filter_cur_frame

    if image[func(image)].size <= 0:
        if particle_filter_cur_frame >= max_frame:
            return None, -1, -1
    else:
        particle_filter_cur_frame = 0
        if particles is None:
            particles = init_particle(func, image)

    if particles is None:
        return None, -1, -1

    resample(particles)
    predict(particles)
    weight(particles, func, image)

    x, y = measure(particles)
    return particles, x, y


if __name__ == '__main__':
    def is_color(region):
        # 色調範圍限制
        return (region >= 11) & (region < 34)

    # 開啟攝像頭
    cap = cv2.VideoCapture(0)
    particles = None

    while cv2.waitKey(30) < 0:
        # 讀取影象
        _, frame = cap.read()
        # 影象轉化為HSV格式
        frame_hsv = cv2.cvtColor(frame, cv2.COLOR_BGR2HSV)
        frame_h = frame_hsv[:, :, 0]
        _, frame_s = cv2.threshold(frame_hsv[:, :, 1], 0, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)
        _, frame_v = cv2.threshold(frame_hsv[:, :, 2], 0, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)

        # s或v為0的畫素點,h賦值為0
        frame_h[(frame_s == 0) | (frame_v == 0)] = 0

        # 獲取粒子點,中心座標
        particles, x, y = particle_filter(particles, is_color, frame_h)

        if particles is not None:
            # 防止出邊界
            valid_particles = np.int32(particles[(particles[:, 0] >= 0) & (particles[:, 0] < frame.shape[1])
                                       & (particles[:, 1] >= 0) & (particles[:, 1] < frame.shape[0])])

            # 修改粒子點的顏色
            for i in range(valid_particles.shape[0]):
                frame[valid_particles[i][1], valid_particles[i][0]] = [255, 0, 0]
            p = np.array([x, y], dtype=np.int32)
            # 根據中心點畫框
            cv2.rectangle(frame, tuple(p - 15), tuple(p + 15), (0, 0, 255), thickness=2)
        cv2.imshow('green', frame)

    cap.release()
    cv2.destroyAllWindows()

測試結果: