Kullback–Leibler divergence(KL情報量)から画像の類似度を推定してみた(ブレ検出付き)


こんにちは！！[Aster](https://twitter.com/Beehibiki)です！！！

2022年9月23日-25日にかけて[IwakenLab](https://iwakenlab.jp/)で熱海に温泉開発合宿を行いました！！
まず、支援者としてたくさんの方々から金銭面中心に支援をいただきました。本当にありがとうございました！！！

開発合宿に参加したメンバーで「開発合宿ブログリレー」として、10月発信していくことにしました。(アドカレのような企画の10月版です。10月1日~10月31日までブログリレーしていきます。

https://iwaken.notion.site/XR-a2faaa2ebc7c4b5381d48e1d2aebd11e


## 開発合宿で開発したもの(Auto pic in video)
今回2泊3日でKullback–Leibler divergence(KL情報量)から画像の類似度を推定を行い、自動で撮影するシステムを開発いたしました。
なぜKL情報量から類似度の推定を行ったかというと、パターン認識と機械学習(PRML)の本を読んでおり、ただ単純にKL情報量を使ってみたかったからです。

このシステムはまず人が理想な人の配置を決めそのボックスを配置します。
次にYOLO v5を用いて人の検出を行います。
人が決めたボックスとYOLO v5で検出したボックスの中心座標と幅と高さを求めます。
中心座標を平均、幅と高さを分散と仮定し、グラフを算出します。
算出したグラフをKullback–Leibler divergence(KL情報量)を用いてグラフの類似度を求めます。
類似度が高いものだけを検出し、動画から写真を保存する。

という流れで動画から理想な配置の写真を保存しています。

![IwakenLab開発合宿成果発表会 - Google スライド — Mozilla Firefox 2022_10_01 21_17_56.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1340659/f1099ac6-790c-1e03-e87a-41842f076fd5.png)


全体のコードは最後に記載しています。
2泊3日でコードを書いたので、殴り書きになってしまいました。そのため、もっと良い書き方がありましたら、コメントくれると嬉しいです。

## Kullback–Leibler divergence(KL情報量)とは
Kullback-Leibler divergence ( KLダイバージェンス、KL情報量 )は、２つの確率分布がどの程度似ているかを表す尺度です。
定義は以下のようになります。

```math
KL(p||q) = \int_{-\infty}^{\infty}p(x)\ln \frac{p(x)}{q(x)}dx

```

```math
KL(p||p) = \int_{-\infty}^{\infty}p(x)\ln \frac{p(x)}{p(x)}dx
         = \int_{-\infty}^{\infty}p(x)\ln(1)dx
         = 0
```
上記の式より、同じ確率分布では0となることがわかります。
そして、常に0を含む正の値となり、確率分布が似ていない程、大きな値となることもわかります。
グラフで見ると

![high.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1340659/f99746e6-15c2-0ce9-d165-6e7b863d2e76.png)
グラフの類似度が高い
![low.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1340659/d483ec80-8667-2192-2eb3-f21b8fc79212.png)
グラフの類似度が低い

このグラフの比較を用いることで画像の類似度を比較しています。

コードは以下のようになります。(類似度が高いほうの例)

```python:graph.py
import matplotlib.pyplot as plt
import numpy as np

#上のgraphの例

# サンプル数
N=10000

# サンプルをN個生成
x = np.random.normal(size=N)
x2 = np.random.normal(size=N)

# 生成した点のヒストグラムをプロット
plt.figure(figsize=(10, 5))
hist, _ = np.histogram(x)
plt.plot(hist, label="normal")
hist, _ = np.histogram(x2)
plt.plot(hist, label="normal2")

plt.grid()
#--------------------------
# Kullback-Leibler divergence
def KLD(a, b, bins=10, epsilon=.00001):

    a_hist, _ = np.histogram(a, bins=bins) 
    b_hist, _ = np.histogram(b, bins=bins)
    

    a_hist = (a_hist+epsilon)/np.sum(a_hist)
    b_hist = (b_hist+epsilon)/np.sum(b_hist)
    
    return np.sum([ai * np.log(ai / bi) for ai, bi in zip(a_hist, b_hist)])

print(KLD(x, x2))
```

## ラプラシアンフィルタを用いたブレ検出
今回動画から写真を保存する中で一つ問題が発生しました。
それがブレた写真も保存してしまう現象です。
これは、1フレームごと類似度を推定してしまっているからです。
そこでこれを解決できないかと考え、ラプラシアンフィルタならブレているか検出できるのではないかと考えました。

KL情報量で類似度が高い画像をラプラシアンフィルタにかけぶれているか検出するだけです。


 ![IwakenLab開発合宿成果発表会 - Google スライド — Mozilla Firefox 2022_10_01 21_18_04.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1340659/55391eb0-ce15-1646-fa78-75d4ff80d101.png)

## ラプラシアンフィルタとは
ラプラシアンフィルタ（Laplacian Filter）は、二次微分を利用して画像から輪郭を抽出するフィルタです。

![無題のプレゼンテーション - Google スライド — Mozilla Firefox 2022_10_15 22_37_33.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1340659/72a500b6-6cd3-b9ac-2d57-dbef8b7c3cef.png)


上の図の赤枠の中のように、全ての画素値が同じ領域では、出力される画素の値は0になります。
対象画素と周囲画素に差がある場合には、出力される画素の値の絶対値が大きくなります。
結果として、ラプラシアンフィルタを用いると画像のエッジ抽出が可能となります。
この特性を生かし、輪郭がぼけている写真を取り除くことができます。

コードは以下のようになります。

```python:laplacian.py

n = 0

dir_path = "save_image"
dir_train = "/train"
dir_test = "/test"
basename = "/evolution_frame"
base_path = os.path.join(dir_path, basename)

gray = cv2.cvtColor(img_val, cv2.COLOR_BGR2GRAY)
laplacian = cv2.Laplacian(gray,cv2.CV_64F)

if laplacian.var() > 60:
    print("Not Blurry")
    cv2.imwrite( dir_path+ dir_train + basename +str(n) + '.jpg', imgs)
    cv2.imwrite(dir_path+ dir_test + basename +str(n) + '.jpg', img_val)
    n += 1
else:
    print("Blurry")
```

# 全体コード

KL情報量を用いた画像の類似度の比較方法です。
```python:KL_evaluation.py
from cmath import log
import math
import sys
from tokenize import group
import cv2
from pandas import value_counts
import torch
import numpy as np


def calkLD_simgle(objects):
    
    mid_object = []
    box = []

    if len(objects) == 1:
        #理想の配置のボックスを決定(今回は真ん中)
        '''
        webcam

        x = 640
        y = 420
        w = 400
        h = 600
        '''

        '''
        iPhone13 mini
        '''
        x = 960
        y = 540
        w = 800
        h = 1000

        box.append(x)
        box.append(y)
        box.append(w)
        box.append(h)
        mid_object.append(box)
        
        #objectsはYoloで検出した人のボックス
        #mid_objectは理想の配置

        kl = kl_calculation(objects,mid_object)
        return kl

#KL情報量による類似度比較
def kl_calculation(object1,object2):
    mean1, vcm1 = transfer(object1)
    mean2, vcm2 = transfer(object2)

    value1 = det_k(vcm1)
    value2 = det_k(vcm2)
    ratio_value = value2/value1
    log_value = math.log(ratio_value)
    inverse = np.linalg.inv(vcm2)

    multi = np.dot(inverse,vcm1)

    trace_value = trace(multi)

    delt = []
    delt.append(mean1[0] - mean2[0])
    delt.append(mean1[1] - mean2[1])
    smt = mulVecMul_k(delt,inverse)

    third_value = smt[0] * delt[0] + smt[1] * delt[1]

    kl = 0.5 * log_value + 0.5 * trace_value + 0.5  * third_value -1
  
    return kl

def trace(a):
    result = a[0][0] + a[1][1]
    return result

def det_k(vcm):
    return vcm[0][0] * vcm[1][1] - vcm[0][1] * vcm[1][0]

def mulVecMul_k(ve,mu):

    re = []
    value0 = ve[0] * mu[0][0] + ve[1] * mu[1][0]
    value1 = ve[0] * mu[0][1] + ve[1] * mu[1][1]

    re.append(value0)
    re.append(value1)
    return re

# 平均、分散共分散行列を出力
def transfer(objects):
    mean = []
    vcm = []
    mean.append(objects[0][0])
    mean.append(objects[0][1])

    sigmax = objects[0][2] / 2.58 / 2.0
    sigmay = objects[0][3] / 2.58 / 2.0

    s_sigx = sigmax * sigmax
    s_sigy = sigmay * sigmay

    first = []
    second = []
    first.append(s_sigx)
    first.append(0)
    second.append(0)
    second.append(s_sigy)

    vcm.append(first)
    vcm.append(second)

    return mean,vcm
```


yolov5での人の検出とラプラシアンフィルタでのブレ検出のコードです。
```python:yolov5_detect.py
import math
import sys
from tokenize import group
import cv2
import torch
from kl_evaluation import calkLD_multi
import os

model = torch.hub.load("ultralytics/yolov5", "yolov5s", pretrained=True)
model.conf = 0.5 #--- 検出の下限値（<1）。設定しなければすべて検出
model.classes = [0] #--- 0:person クラスだけ検出する。設定しなければすべて検出

camera = cv2.VideoCapture(0)

'''
EpocCamを使う場合の設定
camera = cv2.VideoCapture(1)
'''

framewidth = camera.get(3)
frameheight = camera.get(4)

framemiddlewidth = framewidth/2
framemiddleheight = frameheight/2

dir_path = "save_image"
dir_train = "/train"
dir_test = "/test"

basename = "/evolution_frame"

base_path = os.path.join(dir_path, basename)

n = 0

while True:

    ret, imgs = camera.read()

    results = model(imgs)
    bboxs = []

    cc2 = (128,0,0)

    img_val = imgs.copy()

    cv2.rectangle(
        imgs,
        (int(960 - 800/2), int(540 - 1000/2)),
        (int(960 + 800/2), int(540 + 1000/2)),
        color=cc2,
        thickness=2,
        )

    for *box, conf, cls in results.xyxy[0]:

        bbox = []
        s = model.names[int(cls)]+":"+'{:.1f}'.format(float(conf)*100)

        cc = (255,255,0)
        
        minx,miny,maxx,maxy = box[0],box[1],box[2],box[3]

        w = maxx - minx
        h = maxy - miny
        x = minx + w/2
        y = miny + h/2

        bbox.append(x)
        bbox.append(y)
        bbox.append(w)
        bbox.append(h)
        bboxs.append(bbox)

        center_x_sum = 0
        center_y_sum = 0
        grounpcenter_x = 0
        grounpcenter_y = 0

        for i in range(len(bboxs)):
            center_x = (bboxs[i][0] + bboxs[i][2]) /2
            center_y = (bboxs[i][1] + bboxs[i][3]) /2

            center_x_sum = center_x_sum + center_x
            center_y_sum = center_y_sum + center_y
        
        grounpcenter_x = center_x_sum/len(bboxs)
        grounpcenter_y = center_y_sum/len(bboxs)

        xMiddle = grounpcenter_x
        yMiddle = grounpcenter_y
        
        xDifference = xMiddle - framemiddlewidth
        yDifference = yMiddle - framemiddleheight

        xDifferenceAbsolute = abs(xMiddle - framemiddlewidth)
        yDifferenceAbsolute = abs(yMiddle - framemiddleheight)

        cv2.rectangle(
            imgs,
            (int(box[0]), int(box[1])),
            (int(box[2]), int(box[3])),
            color=cc,
            thickness=2,
            )

        cv2.rectangle(imgs, (int(box[0]), int(box[1])-20), (int(box[0])+len(s)*10, int(box[1])), cc, -1)
        cv2.putText(imgs, s, (int(box[0]), int(box[1])-5), cv2.FONT_HERSHEY_PLAIN, 1, cc2, 1, cv2.LINE_AA)

    if len(bboxs) != 0 :
        #KL情報量から類似度比較
        calk_val =  calkLD_multi(bboxs)
        if calk_val != None: 
            test = calk_val.to('cpu').detach().numpy().copy()
            if test <= 0.1:
                print(test)
                
                #ラプラシアンフィルタによるブレ検出
                gray = cv2.cvtColor(img_val, cv2.COLOR_BGR2GRAY)

                laplacian = cv2.Laplacian(gray,cv2.CV_64F)

                if laplacian.var() > 60:

                    print("Not Blurry")
                    cv2.imwrite( dir_path+ dir_train + basename +str(n) + '.jpg', imgs)
                    cv2.imwrite(dir_path+ dir_test + basename +str(n) + '.jpg', img_val)
                    n += 1
                else:
                    print("Blurry")



    cv2.imshow('color',imgs)


    if cv2.waitKey(1) & 0xFF == ord('q'):
        break
```

# まとめ
今回は1人のみしか配置することができませんでしたが、今後は複数人でも対応できるように開発を進めていきます。

このようにKL情報量を用いることで理想な配置の写真をとることができ、ラプラシアンフィルタを用いてぶれていない写真のみを抽出することが可能となります。
これで、動画からこのフレームの写真を取り除きたいと思ったら自動で抽出してくれるという便利なシステムができます！！