【Python】OpenCVで特徴点の追跡【メダカの軌跡】

【Python】OpenCVで特徴点の追跡【メダカの軌跡】

タイトルのサムネイル写真は、メダカの動きをOpenCVでコーナー検出し、特徴点の追跡を行なって軌跡を描いたものです。この記事ではそのやり方を解説していきます。

Python上のOpenCVで calcOpticalFlowPyrLK 関数によって特徴点の追跡が可能です。つまりは、動画フレーム内で特徴点が移動した場所を追跡できるようになります。特徴点とは、コーナー検出などによって得られる画像内の物体の特徴的な座標になります。

はじめに

こちらのメダカビオトープの動画素材をOpenCVで加工していきます。

メダカビオトープの動画素材
メダカビオトープの動画素材

ここではPython3.xとOpenCVを使用しました。

$ pip list | grep opencv
opencv-contrib-python                             4.6.0.66
opencv-python                                     4.5.5.62

OpenCVがまだインストールされてない場合は、pipコマンドでインストールしましょう。numpyも使いますので合わせてインストールします。

$ pip install opencv-python
$ pip install opencv-contrib-python
$ pip install numpy

他にもPython x OpenCVに関する記事を書いてますので、合わせてご参考いただければと思います。

HSV色空間でメダカのみを抽出する

【Python】OpenCVで画像操作いろいろ(グレースケール・モノ・輪郭抽出・切り抜く・透過)で行なったような方法で、メダカのみをHSV色空間のフィルタリングで抽出してみます。

カラーピッカーでメダカのカラダのRGB色を取得して、HSV色空間へ変換、特定のHSV色のみを抽出した例です:

import cv2
import numpy as np
import itertools
import colorsys


def rgb2hsv(rgb):
    r,g,b = rgb
    h,s,v = colorsys.rgb_to_hsv(r, g, b)
    return int(h*180), int(s*255), int(v)

def hsv_maks(img, rgb_list=None):
    for i, hsv in enumerate([rgb2hsv(rgb) for rgb in rgb_list]):
        h,s,v = hsv
        lower =  np.array([h-3, s-4, v-4])
        upper =  np.array([h+3, s+4, v+4]) 
        src = cv2.inRange(img,lower,upper)

        if i == 0:
            mask_hsv = src
        else:
            mask_hsv = cv2.addWeighted(src1=src,alpha=1,src2=mask_hsv,beta=1,gamma=0)

    return mask_hsv

def main(inpath, outpath, rgb_list):
    # -------------------------------------------------
    # setup
    # -------------------------------------------------

    vidcap = cv2.VideoCapture(inpath)

    width = int(vidcap.get(cv2.CAP_PROP_FRAME_WIDTH))
    height = int(vidcap.get(cv2.CAP_PROP_FRAME_HEIGHT))
    fps = 30 #vidcap.get(cv2.CAP_PROP_FPS)

    fmt = cv2.VideoWriter_fourcc('m', 'p', '4', 'v')
    writer = cv2.VideoWriter(outpath, fmt, fps, (width, height)) # ライター作成

    counter = 0


    # -------------------------------------------------
    # loop
    # -------------------------------------------------    
    while True:

        grabbed_frame, frame = vidcap.read()
        if frame is None:
            break
        
        img_hsv = cv2.cvtColor(frame, cv2.COLOR_BGR2HSV)
        mask_hsv = hsv_maks(img_hsv, rgb_list=rgb_list)        
        img = cv2.bitwise_and(frame, frame, mask=mask_hsv)

        cv2.imshow('img', img)
        img = cv2.resize(img, (width, height))
        writer.write(img)

        key = cv2.waitKey(5)
        if key == 27: # Esc
            break

        counter += 1

    vidcap.release()
    writer.release()

if __name__ == '__main__':
    
    inpath = '../assets/biotope_short.mov' # 1280 x 720
    outpath = '../build/medaka_hsv.mp4'
    
    rgb_list = [
        [156, 134, 70], # 金色メダカ
        [144,146,80], # 金色メダカのお腹 
        [217, 186, 137], # 金色メダカ
        [122, 111, 48], # 金色メダカ
        [150,120,56], # 金色メダカ
        [133,106,42], # 金色メダカ
        [137,122,52], # 金色メダカ
        [144,145,80], # 金色メダカ
        [185,149,77], # 金色メダカ
        [130,117,69], # 金色メダカ        	
        	
        [89, 73, 47], # 黒メダカ
        [115, 105, 62], # 黒メダカ
        [79, 61, 29], # 黒メダカ
        [100, 97, 47], # 黒メダカ
        [120, 101, 46], # 黒メダカ
        [108, 110, 85], # 黒メダカのお腹
        [69,67,31], # 黒メダカ
        [87,77,16], # 黒メダカ
        [103,104,26], # 黒メダカ
        [80,75,43], # 黒メダカ
        [81,79,43], # 黒メダカ
        [98,87,54], # 黒メダカ	
        [147,156,156], # 黒メダカのお腹
        [133,135,139], # 黒メダカのお腹
        [140,133,135], # 黒メダカのお腹
        [53,46,12], # 黒メダカ	
        [78,69,18], # 黒メダカ	
        [77,85,42], # 黒メダカ	
        [90,82,40], # 黒メダカ	
    ]
    
    main(inpath, outpath, rgb_list)

実行結果が次になります。

HSV色空間でメダカのみを抽出
HSV色空間でメダカのみを抽出

少しわかりづらいですが、メダカのカラダの色を捉えることができています。

hsv_maks() 関数内のlowerupperでHSV色空間の範囲を持たせています。範囲をゆるめればメダカのカラダの色がもう少しはっきりしますが、今度は葉っぱなど他の物体の色と被ってしまうので狭めに設定しました。

その代わり、メダカのカラダの色をカラーピッカーで取得してリスト化する作業が大変になります。

メダカのカラダのRGB色
メダカのカラダのRGB色

特徴点の検出(コーナー検出)

【Python】OpenCVでコーナーの検出【Harris/Shi-Tomasi】で行なったように、先ほどの映像から特徴点を検出します。

Shi-Tomasiコーナー検出のgoodFeaturesToTrack() 関数を使ってコーナー検出した例です:

import cv2
import numpy as np
import itertools
import colorsys


def rgb2hsv(rgb):
    r,g,b = rgb
    h,s,v = colorsys.rgb_to_hsv(r, g, b)
    return int(h*180), int(s*255), int(v)

def hsv_maks(img, rgb_list=None):
    for i, hsv in enumerate([rgb2hsv(rgb) for rgb in rgb_list]):
        h,s,v = hsv
        lower =  np.array([h-3, s-4, v-4])
        upper =  np.array([h+3, s+4, v+4]) 
        src = cv2.inRange(img,lower,upper)

        if i == 0:
            mask_hsv = src
        else:
            mask_hsv = cv2.addWeighted(src1=src,alpha=1,src2=mask_hsv,beta=1,gamma=0)

    return mask_hsv

def fill(size, color):
    w, h = size
    canvas =  np.zeros((h, w, 3), dtype="uint8")
    cv2.rectangle(canvas, (0,0), (w,h), color, -1)
    return canvas

def is_nearby(pt1, pt2, maxDistance, minDistance):
    r = (pt1[0] - pt2[0])**2 + (pt1[1] - pt2[1])**2
    if r < maxDistance**2 and r  > minDistance**2:
        return True
    return False

def clamp(n, smallest, largest):
    return max(smallest, min(n, largest))


def main(inpath, outpath, rgb_list):
    # -------------------------------------------------
    # setup
    # -------------------------------------------------

    vidcap = cv2.VideoCapture(inpath)

    width = int(vidcap.get(cv2.CAP_PROP_FRAME_WIDTH))
    height = int(vidcap.get(cv2.CAP_PROP_FRAME_HEIGHT))
    fps = 30 #vidcap.get(cv2.CAP_PROP_FPS)

    fmt = cv2.VideoWriter_fourcc('m', 'p', '4', 'v')
    writer = cv2.VideoWriter(outpath, fmt, fps, (width, height)) # ライター作成

    radius = 2

    counter = 0

    # -------------------------------------------------
    # loop
    # -------------------------------------------------    
    while True:

        grabbed_frame, frame = vidcap.read()
        if frame is None:
            break
        
        # HSVマスク作成
        img_hsv = cv2.cvtColor(frame, cv2.COLOR_BGR2HSV)
        mask_hsv = hsv_maks(img_hsv, rgb_list=rgb_list)        
        img_bgr = cv2.bitwise_and(frame, frame, mask=mask_hsv)

        # グレースケール変換
        gray = cv2.cvtColor(img_bgr, cv2.COLOR_BGR2GRAY)
        gray = cv2.dilate(gray,None,iterations = 1)

        # Shi-Tomasiコーナー検出
        corners = cv2.goodFeaturesToTrack(gray, maxCorners=180, qualityLevel=0.001, minDistance=12)
        corners = np.int0(corners)

        # キャンバスへ描画
        canvas = np.zeros_like(frame)
        
        for i in corners:
            pt = i.ravel()
            cv2.circle(canvas, pt, radius, (255,255,255), -1)

        for pair in itertools.combinations(corners, 2): # 組み合わせ
            d1, d2 = pair
            pt1 = d1.ravel()
            pt2 = d2.ravel()

            if is_nearby(pt1, pt2, 30, 0):
                cv2.line(canvas, pt1, pt2, (255,255,255))

#        canvas = cv2.addWeighted(src1=canvas,alpha=1, src2=frame,beta=0.5,gamma=0)


        # 書き出し
        cv2.imshow('Canvas', canvas)
        canvas = cv2.resize(canvas, (width, height))
        writer.write(canvas)

        key = cv2.waitKey(5)
        if key == 27: # Esc
            break

        counter += 1

    vidcap.release()
    writer.release()

if __name__ == '__main__':
    
    inpath = '../assets/biotope_short.mov' # 1280 x 720
    outpath = '../build/medaka_detected_corner.mp4'
    
    rgb_list = [
        [156, 134, 70], # 金色メダカ
        [144,146,80], # 金色メダカのお腹 
        [217, 186, 137], # 金色メダカ
        [122, 111, 48], # 金色メダカ
        [150,120,56], # 金色メダカ
        [133,106,42], # 金色メダカ
        [137,122,52], # 金色メダカ
        [144,145,80], # 金色メダカ
        [185,149,77], # 金色メダカ
        [130,117,69], # 金色メダカ        	
        	
        [89, 73, 47], # 黒メダカ
        [115, 105, 62], # 黒メダカ
        [79, 61, 29], # 黒メダカ
        [100, 97, 47], # 黒メダカ
        [120, 101, 46], # 黒メダカ
        [108, 110, 85], # 黒メダカのお腹
        [69,67,31], # 黒メダカ
        [87,77,16], # 黒メダカ
        [103,104,26], # 黒メダカ
        [80,75,43], # 黒メダカ
        [81,79,43], # 黒メダカ
        [98,87,54], # 黒メダカ	
        [147,156,156], # 黒メダカのお腹
        [133,135,139], # 黒メダカのお腹
        [140,133,135], # 黒メダカのお腹
        [53,46,12], # 黒メダカ	
        [78,69,18], # 黒メダカ	
        [77,85,42], # 黒メダカ	
        [90,82,40], # 黒メダカ	
    ]
    
    main(inpath, outpath, rgb_list)

実行結果が次になります。

メダカの特徴点の検出
メダカの特徴点の検出

Processingのジェネラティブアートっぽく、点と点を線で結んでみました。一定の距離以下の点と点を結ぶ方法は比較的カンタンでして、やり方は【Python】OpenCVで図形の描画からアニメーションまで【線・四角・丸・塗りつぶし】の記事で解説しています。興味ある方はぜひご覧ください。

特徴点を追跡する(モーション解析)

コーナー検出によって見つけられた特徴点を、動画フレームの前後で保持(追跡)する便利な機能がOpenCVに備えられてます。それが calcOpticalFlowPyrLK() 関数です。

calcOpticalFlowPyrLK()Lucas-Kanade法(英語論文) を使ったオプティカルフローと呼ばれるものです。入力には、以前と現在の二つのグレースケール画像と、以前の特徴点の情報が必要になります。結果として移動後の特徴点の配列を得られます。

次はメダカの特徴点(軌跡)を追跡した例です:

import cv2
import numpy as np
import colorsys


def rgb2hsv(rgb):
    r,g,b = rgb
    h,s,v = colorsys.rgb_to_hsv(r, g, b)
    return int(h*180), int(s*255), int(v)

def hsv_maks(img, rgb_list=None):
    for i, hsv in enumerate([rgb2hsv(rgb) for rgb in rgb_list]):
        h,s,v = hsv
        lower =  np.array([h-3, s-4, v-4])
        upper =  np.array([h+3, s+4, v+4]) 
        src = cv2.inRange(img,lower,upper)

        if i == 0:
            mask_hsv = src
        else:
            mask_hsv = cv2.addWeighted(src1=src,alpha=1,src2=mask_hsv,beta=1,gamma=0)

    return mask_hsv

def clamp(n, smallest, largest):
    return max(smallest, min(n, largest))


def main(inpath, outpath, rgb_list):
    # -------------------------------------------------
    # setup
    # -------------------------------------------------

    vidcap = cv2.VideoCapture(inpath)

    width = int(vidcap.get(cv2.CAP_PROP_FRAME_WIDTH))
    height = int(vidcap.get(cv2.CAP_PROP_FRAME_HEIGHT))
    fps = 30 #vidcap.get(cv2.CAP_PROP_FPS)

    fmt = cv2.VideoWriter_fourcc('m', 'p', '4', 'v')
    writer = cv2.VideoWriter(outpath, fmt, fps, (width, height)) # ライター作成

    radius = 2

    counter = 0
    alpha = 0.0
    beta = 1.0
    
    
    # Parameters for lucas kanade optical flow
    lk_params = dict( winSize  = (15,15),
                    maxLevel = 2,
                    criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, 10, 0.03))

    grabbed_frame, frame = vidcap.read()
    old_masks = []

    # -------------------------------------------------
    # loop
    # -------------------------------------------------    
    while True:

        grabbed_frame, frame = vidcap.read()
        if frame is None:
            break
        mask = np.zeros_like(frame)
        
        img_hsv = cv2.cvtColor(frame, cv2.COLOR_BGR2HSV)
        mask_hsv = hsv_maks(img_hsv, rgb_list=rgb_list)        
        img_bgr = cv2.bitwise_and(frame, frame, mask=mask_hsv)
        gray = cv2.cvtColor(img_bgr, cv2.COLOR_BGR2GRAY)                
        corners = cv2.goodFeaturesToTrack(gray, maxCorners=180, qualityLevel=0.001, minDistance=12)

        if counter > 0:

            new_pts, st, err = cv2.calcOpticalFlowPyrLK(old_gray, gray, old_pts, None, **lk_params)
            good_new = new_pts[st==1]
            good_old = old_pts[st==1]

            for new,old in zip(good_new,good_old):
                pt0 = np.int0(new.ravel())
                pt1 = np.int0(old.ravel())
                cv2.line(mask, pt0, pt1, (255,128,0), 2)
                cv2.circle(mask, pt1, radius, (255,0,255), -1)

            old_masks.insert(0, mask)
            if len(old_masks) > 15:
                old_masks = old_masks[:15]
            a = 1
            for m in old_masks:
                mask = cv2.addWeighted(src1=m,alpha=a, src2=mask,beta=1,gamma=0)
                a -= 0.05

            alpha = clamp(alpha, 0, 1)
            beta = clamp(beta, 0, 1)

            cv2.imshow('img', mask)
            mask = cv2.resize(mask, (width, height))
            writer.write(mask)


        key = cv2.waitKey(5)
        if key == 27: # Esc
            break

        counter += 1
        old_pts = corners.copy()
        old_gray = gray

    vidcap.release()
    writer.release()

if __name__ == '__main__':
    
    inpath = '../assets/biotope_short.mov' # 1280 x 720
    outpath = '../build/medaka_flow_short.mp4'
    
    rgb_list = [
        [156, 134, 70], # 金色メダカ
        [144,146,80], # 金色メダカのお腹 
        [217, 186, 137], # 金色メダカ
        [122, 111, 48], # 金色メダカ
        [150,120,56], # 金色メダカ
        [133,106,42], # 金色メダカ
        [137,122,52], # 金色メダカ
        [144,145,80], # 金色メダカ
        [185,149,77], # 金色メダカ
        [130,117,69], # 金色メダカ        	
        	
        [89, 73, 47], # 黒メダカ
        [115, 105, 62], # 黒メダカ
        [79, 61, 29], # 黒メダカ
        [100, 97, 47], # 黒メダカ
        [120, 101, 46], # 黒メダカ
        [108, 110, 85], # 黒メダカのお腹
        [69,67,31], # 黒メダカ
        [87,77,16], # 黒メダカ
        [103,104,26], # 黒メダカ
        [80,75,43], # 黒メダカ
        [81,79,43], # 黒メダカ
        [98,87,54], # 黒メダカ	
        [147,156,156], # 黒メダカのお腹
        [133,135,139], # 黒メダカのお腹
        [140,133,135], # 黒メダカのお腹
        [53,46,12], # 黒メダカ	
        [78,69,18], # 黒メダカ	
        [77,85,42], # 黒メダカ	
        [90,82,40], # 黒メダカ	
    ]
    
#    print(rgb2hsv([156, 134, 70]))
    
    main(inpath, outpath, rgb_list)

実行結果が次になります。

メダカの特徴点の検出
メダカの特徴点の検出

ぜひ続きの映像をご覧ください。

calcOpticalFlowPyrLK

calcOpticalFlowPyrLK(prevImg, nextImg, prevPts, nextPts, status=..., err=..., winSize=..., maxLevel=..., criteria=..., flags: int = ..., minEigThreshold=...)
パラメータ 意味
prevImg 1番目の入力画像
nextImg 2番目の入力画像
prevPts 1番目の入力画像の特徴点
winSize ピラミッドレベルにおける探索窓のサイズ
maxLevel 画像ピラミッドの最大レベル数
criteria 反復探索アルゴリズムの停止基準
derivLambda 画像の空間微分の相対的な重み

(大きな動きの場合はピラミッドのスケールをアップする)

記事に関するご質問などがあれば、
Twitter または お問い合わせ までご連絡ください。
Python学習にオススメの本をご紹介!
Pandasでデータサイエンスはじめよう!
関連記事