iphoneで撮影した動画から心拍数を抽出してみた

#はじめに
近年カメラ映像から抽出した生体信号を利用して，様々な応用研究が行われています．
2020年には5Gが導入される予定ですし，遠隔医療やAIアシスタントでの利用など，今後の医療ヘルスケアの礎となる技術ではないでしょうか．
また，これらの生体信号は医療以外の分野でも利用価値があり，今後ますます注目されるはずです．

そこで今回，iPhoneで撮影した動画像から心拍数を取得してみました．
なおこの記事はあくまで僕の個人的な取り組みです，医学的に検証されたものではないのでご注意ください．

###参考文献
・小原一誠・他　「映像からの脈波情報抽出」　
　東北大学サイバーサイエンスセンター先端情報技術研究部 吉澤・杉田研究室
　http://www.sairct.idac.tohoku.ac.jp/wp-content/uploads/2016/07/160728_c2_yoshizawa_DLpre2.pdf

・北島利浩・他「生体信号を利用したプライバシー保護画像からの顔検出手法」
　電気学会論文誌 C(電子・情報・システム部門誌)
　IEEJ Transactions on Electronics, Information and Systems Vol.138 No.3 2018 pp.210–220

#解析アルゴリズム
１．検証データ取得（動画撮影）
２．動画をフレームに分割
３．画像からRGBのGの値のみを抽出
４．解析領域を指定し，領域内のGの時間変化を求める
５．パワースペクトルを算出し周波数でフィルタリング
６．ピーク周波数から脈拍数を算出

##1.検証データ取得
撮影機器：iphoneXR
撮影場所：室内
光源：蛍光灯，緑の照明（解析を容易にするため緑のライトも照射した）
撮影部位：顔
撮影方法：顔とカメラを固定
動画サイズ：1080p/60fpsの約50s

##2.動画像のフレーム分割

ライブラリを読み込んでおきます．

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import cv2
import seaborn as sns
import os
import shutil
from scipy import signal
import glob
from scipy import fftpack

フレーム分割する関数を定義します．

#フレーム分割
def frame_split(video_file='./image_wave.wmv', image_dir='./image_wave２/',
                   image_file='image_wave２-%s.png'):
    if os.path.exists(image_dir):
        shutil.rmtree(image_dir)

    if not os.path.exists(image_dir):
        os.makedirs(image_dir)

    i = 0
    
    cap = cv2.VideoCapture(video_file)
    while (cap.isOpened()):
        flag, frame = cap.read()  # Capture frame-by-frame
        if flag == False:  
            break
        cv2.imwrite(image_dir + image_file % str(i).zfill(6),
                    frame)  
        print('Save', image_dir + image_file % str(i).zfill(6))
        i += 1

    cap.release()

frame_splitを実行，今回のデータは50秒の動画なので，3000枚になります．

frame_split("/Users/IMG.mov")

##3.画像からRGBのGの値のみを抽出
###なぜGなのか？
これはこの分野ではもはや常識だと思いますが，一応簡単に書いておきます．
ヘモグロビンは緑の波長帯域に対して高い吸光度を持ちます．そのため，Gの値は血中ヘモグロビンの量に敏感に反応します．血中ヘモグロビンは心臓の動きよって増減するので，心拍数がわかります．
心拍数をとるなら緑を使うのが一般的です．心拍数を取得できるウェアブル端末ではだいたい緑が使われていると思います．（僕が確認した限りでは，fitbitとよくわからない中国製の端末では使われていました）

#画像ファイルを格納，時系列にソート
files = glob.glob("/Users/image_wave２/*.png")
files.sort()

#画像を読み込み，green成分を抽出
images = [cv2.imread(files[i]) for i in range(len(files))]
images_green = [pd.DataFrame(images[j][:,:,1]) for j in range(len(images))]

ここで注意が１つ．cv2.imreadで画像を読み込んだ場合，返り値のshapeは（高さ x 幅 x 色数）となりますが，色数はBGRの順番になっています．今回はGが必要なので[１]です．
images_greenには3000枚のG成分が格納されました．

##４．解析領域を指定し，領域内のGの時間変化を求める

解析領域は顔の頬の部分100×50pixelを指定しています．
メディアンフィルタで平滑化した後に平均値を算出します．

#parameter
in1 = 1100
in2 = 1200
co1 = 200
co2 = 250

#メディアンフィルタによる平滑化
images_green_median = [cv2.medianBlur(images_green[i].iloc[ in1 : in2 , co1 :co2 ].values,ksize=5)  for i in range(len(images_green))]

#平均値を算出
img_mean = [images_green_median[i].mean()  for i in range(len(images_green_median))]

プロットしてみると，

#サンプリング周波数
fs = 60

end_time = round(len(img_mean)/fs)
time = np.arange(0 , end_time , end_time/len(img_mean))
plt.plot(time, img_mean)
plt.xlabel("time(s)",size=15)

10sぐらいの周期性が確認できます，おそらく呼吸性の周期と思われます．

##５．パワースペクトルを算出し周波数でフィルタリング

###解析窓長について
脈拍数を40-150bpmとすると，その周波数は0.7-2.5Hz，周期は1.6-0.4sです．
よって0.7Hzを10周期分入れるとすれば、解析窓長は16s必要になります。

解析窓長を決める時に重要なのは周波数分解能と時間分解能はトレードオフであるということです．
今回のサンプリング周波数は60Hz，窓長は15sとしました．
よって周波数分解能は0.067Hz
なので心拍数の単位は0.067*60 = 4bpmです．

まとめると，15sに1回心拍数を算出します，算出される心拍数は4の倍数となります．50sの動画なので3回算出されます．

#パワースペクトルを求める関数を定義
def signal_fft (data,fs):
    n = len(data)-1
    y = fftpack.fft(data)/n
    y = y[0:round(n/2)]
    power = 2*(np.abs(y)**2)
    power = 10*np.log10(power)
    f = np.arange(0,fs/2,fs/n)
    return power , f

#15sの窓でパワースペクトル算出
time_interval = 15
loop_len = list(range(round(len(img_mean)/(fs*time_interval))))
fft_result = []
for n,i in enumerate(loop_len):
    if n ==0:
        img = img_mean[:time_interval*fs]
        fft_result.append(signal_fft(img,fs))
    else :
        img = img_mean[time_interval*fs*n : time_interval*fs*(n+1)]
        fft_result.append(signal_fft(img,fs)) 

#周波数フィルタリング（0.7-2.5Hzの値のみ抽出）
data = []
for j in range(len(fft_result)):
    fft_data = [fft_result[j][0][n] for n,i in enumerate(fft_result[j][1]) if 0.7< i < 2.5]
    data.append(fft_data)
frequ = [i for n,i in enumerate(fft_result[0][1]) if 0.7< i < 2.5]

##６．ピーク周波数から脈拍数を算出

data = np.array(data)
#ピークのindexを取得
peak_index  = [data[i].argmax() for i in range(len(data))]

#取得したindexから該当の周波数を取得し心拍数に変換
heart_rate = [round(frequ[i]*fs) for i in peak_index]

print("heart rate:{}".format(heart_rate))

heart rate:[76.0, 76.0, 72.0]

若干速い気がしますがこんなもんでしょうか．

#おわりに
今回は簡略化のため緑の光を照射していますが，太陽光には全ての波長が含まれているため，光を照射せずとも脈拍を取得できます．またヘモグロビンは近赤外光にも高い吸光特性を持っているので，緑以外でも心拍数を検知できるようです．