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

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

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

はじめに
HSV色空間でメダカのみを抽出する
特徴点の検出（コーナー検出）
特徴点を追跡する（モーション解析）
1. calcOpticalFlowPyrLK

はじめに

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

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

shell

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

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

shell

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

他にもPython x OpenCVに関する記事を書いてます。

【Python】OpenCVで画像操作いろいろ（グレースケール・モノ・輪郭抽出・切り抜く・透過）

2022-07-21

【Python】OpenCVでコーナーの検出【Harris/Shi-Tomasi】

2022-07-26

【Python】OpenCVで図形の描画からアニメーションまで【線・四角・丸・塗りつぶし】

2022-07-27

Pythonで始めるOpenCV入門

Kindleで見る楽天市場で見る

Pythonで始めるOpenCV 4プログラミング

Amazonで見る楽天市場で見る

OpenCVではじめようディープラーニングによる画像認識

Amazonで見る Kindleで見る楽天市場で見る

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

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

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

medaka_hsv.py

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_maks() 関数内のlower と upperでHSV色空間の範囲を持たせてます。範囲をゆるめればメダカのカラダの色がもう少しはっきりしますが、今度は葉っぱなど他の物体の色と被ってしまうので狭めに設定しました。その代わり、メダカのカラダの色をカラーピッカーで取得してリスト化する作業が大変になります。

特徴点の検出（コーナー検出）

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

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

medaka_detected_corner.py

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

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

medaka_optical_flow.py

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)

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

▼ ぜひこちらの映像もご覧ください。

PythonとOpenCVを使ったメダカで動かすジェネラティブアート〜その2 - YouTube

calcOpticalFlowPyrLK

calcOpticalFlowPyrLK(prevImg, nextImg, prevPts, nextPts, status=..., err=..., winSize=..., maxLevel=..., criteria=..., flags: int = ..., minEigThreshold=...)

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