きっかけ
近年、スマートウォッチやスマートリングの普及で誰でも、自分の健康状態をモニタリングすることができるようになった。
だが、非接触のセンサーでそれっぽい値がとれないだろうか？ とふと考えた。
ちなみに非接触センサーで心拍数を測定することを
rPPG（remote Photoplethysmography）と略して言うことが多いらしい。
調査
先行研究として、やられているものがあればそのアルゴリズムを使ってみようと考えた。
上記では緑成分に着目して、時間変化量をみてやることにしている。
実際にApplewatchを持っている方はわかるだろうが、心拍を測定しているときは時計の下部から緑の光が出ている。
POS法について
論文等の調査をしてみるとPOS(Plane-Orthogonal-to-Skin)法というものが出てくる。
自分は知らなかったが以下論文で発表された方法
Algorithmic Principles of Remote PPG
簡単に説明するとRGB成分から動的に心拍成分の変動を抽出するもの。
RGB成分→X,Y成分の合成信号→比率に応じて線形結合→rPPG という感じ
数式を使って説明すると
- RGB信号の正規化
顔領域の平均RGB値を正規化：
$$
\bar{C}(t) = \frac{C(t)}{\mathrm{mean}(C(t))} =
\begin{bmatrix}
\bar{R}(t) \
\bar{G}(t) \
\bar{B}(t)
\end{bmatrix}
$$
2. 2成分信号の生成
$$
X_s(t) = \bar{G}(t) - \bar{B}(t)
$$
$$
Y_s(t) = -2 \cdot \bar{R}(t) + \bar{G}(t) + \bar{B}(t)
$$
3. 時間的標準偏差の比率の計算
$$
\alpha = \frac{\sigma(X_s)}{\sigma(Y_s)}
$$
4. rPPG信号の合成
$$
s(t) = X_s(t) + \alpha \cdot Y_s(t)
$$
pythonで書くと
import numpy as np
def pos_method(rgb_signals):
"""
POS法によるrPPG信号抽出
Parameters:
rgb_signals (np.ndarray): shape=(T, 3) のRGB信号 (R, G, B)
各行が1フレームの平均RGB
Returns:
s (np.ndarray): shape=(T,) のrPPG信号
"""
# Step 1: 正規化（各チャンネルを平均で割る）
C = rgb_signals
mean_color = np.mean(C, axis=0)
C_bar = C / mean_color # shape: (T, 3)
# Step 2: Xs, Ys の計算
Xs = C_bar[:, 1] - C_bar[:, 2] # G - B
Ys = -2 * C_bar[:, 0] + C_bar[:, 1] + C_bar[:, 2] # -2R + G + B
# Step 3: ゼロ平均化（中心化）
Xs = Xs - np.mean(Xs)
Ys = Ys - np.mean(Ys)
# Step 4: α の計算（標準偏差の比率）
std_X = np.std(Xs)
std_Y = np.std(Ys)
alpha = std_X / std_Y if std_Y != 0 else 0
# Step 5: 合成信号
s = Xs + alpha * Ys
return s
PCのカメラで動かしてみる
originalの論文では1.6秒ごとにPOSを計算しているので、それに合わせる。
meadiapipeなどはpipでinstallしておく
import cv2
import numpy as np
import mediapipe as mp
from scipy.signal import find_peaks
import time
def pos_method(rgb_window):
C = rgb_window / np.mean(rgb_window, axis=0)
Xs = C[:, 1] - C[:, 2]
Ys = -2 * C[:, 0] + C[:, 1] + C[:, 2]
Xs -= np.mean(Xs)
Ys -= np.mean(Ys)
std_Y = np.std(Ys)
alpha = np.std(Xs) / std_Y if std_Y != 0 else 0
s = Xs + alpha * Ys
return s
def estimate_bpm(rppg_signal, fps):
peaks, _ = find_peaks(rppg_signal, distance=fps * 0.5) # 0.5秒以上の間隔
peak_intervals = np.diff(peaks) / fps
if len(peak_intervals) == 0:
return None
avg_period = np.mean(peak_intervals)
bpm = 60 / avg_period
return bpm
# 初期設定
buffer_size = 48 # 1.6秒 @30fps
fps = 30 # 固定（または cap.get(cv2.CAP_PROP_FPS) で取得）
rgb_buffer = []
bpm_text = "--"
last_bpm_time = time.time()
# Mediapipe 初期化
mp_face = mp.solutions.face_mesh
face_mesh = mp_face.FaceMesh(static_image_mode=False, max_num_faces=1, refine_landmarks=True)
cap = cv2.VideoCapture(0)
while True:
ret, frame = cap.read()
if not ret:
break
rgb = cv2.cvtColor(frame, cv2.COLOR_BGR2RGB)
results = face_mesh.process(rgb)
h, w, _ = frame.shape
if results.multi_face_landmarks:
landmarks = results.multi_face_landmarks[0].landmark
cx = int(landmarks[1].x * w)
cy = int(landmarks[1].y * h)
size = 30
roi = frame[max(0, cy - size):cy + size, max(0, cx - size):cx + size]
if roi.size > 0:
avg_rgb = np.mean(roi.reshape(-1, 3), axis=0)
rgb_buffer.append(avg_rgb)
# FIFO バッファ更新
if len(rgb_buffer) > buffer_size:
rgb_buffer.pop(0)
if len(rgb_buffer) == buffer_size:
rppg = pos_method(np.array(rgb_buffer))
if time.time() - last_bpm_time > 1.0:
bpm = estimate_bpm(rppg, fps)
if bpm is not None and 40 < bpm < 180:
bpm_text = f"{int(bpm)} BPM"
last_bpm_time = time.time()
# ROI枠の表示
cv2.rectangle(frame, (cx - size, cy - size), (cx + size, cy + size), (0, 255, 0), 2)
# BPM表示
cv2.putText(frame, f"Heart Rate: {bpm_text}", (30, 50),
cv2.FONT_HERSHEY_SIMPLEX, 1, (0, 0, 255), 2)
cv2.imshow("POS-based Heart Rate Estimation (48 frames)", frame)
if cv2.waitKey(1) & 0xFF == ord('q'):
break
cap.release()
cv2.destroyAllWindows()
精度はあまりよくない感じでした・・・