Last active
July 3, 2024 14:09
-
-
Save kongmunist/ba659019a483117a846dc2101e27f13d to your computer and use it in GitHub Desktop.
Webcam Photoplethysmyography
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# Must be run in console to work properly | |
import numpy as np | |
import cv2 | |
import time | |
from scipy import signal | |
import matplotlib.pyplot as plt | |
import matplotlib | |
import threading | |
import scipy.signal as sig | |
face_cascade = cv2.CascadeClassifier('haarcascade_frontalface_default.xml') | |
eye_cascade = cv2.CascadeClassifier('haarcascade_eye.xml') | |
## Filtering function | |
def applyFF(data, sampleFreq): | |
if sampleFreq > 3: | |
sos = sig.iirdesign([.66, 3.0], [.5, 4.0], 1.0, 40.0, fs=sampleFreq, output='sos') | |
return sig.sosfiltfilt(sos, data) | |
else: | |
return data | |
## Calculate the eye landmarks | |
def getFaceXYHWAndEyeXYHW(im): | |
gray = cv2.cvtColor(im, cv2.COLOR_BGR2GRAY) | |
# Only take one face, the first | |
faces = face_cascade.detectMultiScale(gray, 1.3, 5) | |
if len(faces) < 1: | |
return None | |
(x, y, w, h) = faces[0] | |
# Crop out faces and detect eyes in that window. | |
roi_gray = gray[y:y + h, x:x + w] | |
eyes, numDetects = eye_cascade.detectMultiScale2(roi_gray, minNeighbors=10) | |
if len(numDetects) < 2: | |
return None | |
# Change eye coords to be in image coordinates instead of face coordinates | |
eyes[0][0] += x | |
eyes[1][0] += x | |
eyes[0][1] += y | |
eyes[1][1] += y | |
return [faces[0], eyes[0], eyes[1]] | |
def getFace(im): | |
gray = cv2.cvtColor(im, cv2.COLOR_BGR2GRAY) | |
# Only take one face, the first | |
faces = face_cascade.detectMultiScale(gray, 1.3, 5) | |
if len(faces) < 1: | |
return None | |
return faces[0] | |
## Set up data vars for the face/intensity tracking | |
dataLen = 120 | |
camTimes = [0]*dataLen | |
intensities = [] | |
x = list(range(len(intensities))) | |
fig,axes = plt.subplots(2, 1) | |
# Add waveform of intensity | |
ax = axes[0] | |
line1, = ax.plot(x, intensities, 'r-') # Returns a tuple of line objects, thus the comma | |
ax.set_xlim(0,dataLen) | |
ax.set_ylim(0, 255) | |
# Add FFT graph | |
ax2 = axes[1] | |
line2, = ax2.plot(x,intensities, 'b-') | |
plt.show() | |
def updatefig(self, *args): | |
if len(intensities) < dataLen//2: | |
return [line1, ax, line2, ax2] | |
fig.tight_layout() | |
fig.canvas.draw() # This one isn't needed for some reason, but I left it in anyway | |
fig.canvas.flush_events() | |
fs = 1 / (sum(camTimes) / dataLen) | |
tmpIntens = sig.detrend(applyFF(intensities, fs)) | |
line1.set_data(list(range(len(intensities))), intensities) | |
ax.set_ylim(min(intensities)-2, max(intensities)+2) | |
freqs, pows = signal.welch(tmpIntens, fs = fs, nperseg=256) | |
line2.set_data(freqs, pows) | |
ax2.set_ylim(0, max(pows)) | |
ax2.set_xlim(freqs.min(), freqs.max()) | |
print("output BPM: ", round(freqs[np.argmax(pows)]*60,2), fs) | |
return [line1, ax, line2, ax2] | |
ani = matplotlib.animation.FuncAnimation(fig, updatefig, | |
interval=1, blit=True, | |
repeat_delay=1) | |
def getHeadboxFromHead(face): | |
return face[0] + face[2] // 4, face[1] + face[3] // 2, face[0] + 3 * face[ | |
2] // 4, face[1] + face[3] // 2 + 50 | |
# Set up webcam | |
cap = cv2.VideoCapture(0) | |
cap.set(cv2.CAP_PROP_AUTO_EXPOSURE, 0) | |
def readIntensity(intensities, curFrame, cropBoxBounds): | |
now = 0 | |
eyeleft = 0 | |
headTop = 0 | |
eyeright = 0 | |
eyeTop = 0 | |
while True: | |
ret, frame = cap.read() | |
scaleFactor = 0.4 | |
frame = cv2.resize(frame,(-1,-1), fx=scaleFactor, fy=scaleFactor) | |
# tmp = getFaceXYHWAndEyeXYHW(frame) # Haar outputs [x, y, w, h] format | |
face = getFace(frame) | |
if face is not None: | |
# if tmp != None: | |
# face, eye1, eye2 = tmp | |
# eyeleft, headTop, eyeright, eyeTop\ | |
# tmpHeadbox = getHeadbox(face, eye1, eye2) | |
tmpHeadbox = getHeadboxFromHead(face) | |
a = .4 | |
eyeleft = int(tmpHeadbox[0]*a + (1-a)*eyeleft) | |
headTop = int(tmpHeadbox[1]*a + (1-a)*headTop) | |
eyeright = int(tmpHeadbox[2]*a + (1-a)*eyeright) | |
eyeTop = int(tmpHeadbox[3]*a + (1-a)*eyeTop) | |
ROI = frame[headTop:eyeTop, eyeleft:eyeright, 1] | |
intensity = ROI.mean() | |
# intensity = np.median(ROI) # works, but quite chunky. | |
intensities.append(intensity) | |
# Draw the forehead box: | |
curFrame[0] = cv2.rectangle(frame, (eyeleft, headTop), | |
(eyeright, eyeTop), (0, 255, 0), 1) | |
cropBoxBounds[0] = [headTop + 2, eyeTop - 2, eyeleft + 2, eyeright - 2] | |
if (len(intensities) > dataLen): | |
intensities.pop(0) | |
camTimes.append(time.time() - now) | |
now = time.time() | |
camTimes.pop(0) | |
# Start thread of webcam reading | |
cropBoxBounds = [0] | |
curFrame = [0] | |
t1 = threading.Thread(target = readIntensity, daemon=True, args=(intensities, curFrame, cropBoxBounds)) | |
t1.start() | |
time.sleep(.5) | |
while True: | |
frame = curFrame[0] | |
bb = cropBoxBounds[0] | |
ROI = frame[bb[0]:bb[1], bb[2]:bb[3], 1] | |
adjusted_image = ROI.astype('float') - intensities[-1] | |
adjusted_image[adjusted_image < -10] = -10.0 | |
adjusted_image[adjusted_image > 10] = 10.0 | |
# adjusted_image = adjusted_image + adjusted_image | |
frame[bb[0]:bb[1], bb[2]:bb[3], 2] = frame[bb[0]:bb[1], bb[2]:bb[3], 2] + (adjusted_image - adjusted_image.min())*5 | |
frame[bb[0]:bb[1], bb[2]:bb[3], 0] = 20 | |
frame[bb[0]:bb[1], bb[2]:bb[3], 1] = 20 | |
cv2.imshow("yea", frame) | |
if cv2.waitKey(1) & 0xFF == ord('q'): | |
break | |
On Windows I had to add cap = cv2.VideoCapture(0,cv2.CAP_DSHOW)
option and additional delay (2-3s) after starting t1
thread
I modified my code to remove the matplotlib, which I know throws out a lot of errors. Now this will only show the face and the bounding box, and print out the BPM.
The program will spit out junk data for a few seconds, depending on your webcam capture rate and processor speed.
# Must be run in console to work properly
import numpy as np
import cv2
import time
from scipy import signal
import threading
import scipy.signal as sig
face_cascade = cv2.CascadeClassifier('haarcascade_frontalface_default.xml')
eye_cascade = cv2.CascadeClassifier('haarcascade_eye.xml')
def applyFF(data, sampleFreq):
if sampleFreq > 3:
sos = sig.iirdesign([.66, 3.0], [.5, 4.0], 1.0, 40.0, fs=sampleFreq, output='sos')
return sig.sosfiltfilt(sos, data)
else:
return data
def getFaceXYHWAndEyeXYHW(im):
gray = cv2.cvtColor(im, cv2.COLOR_BGR2GRAY)
# Only take one face, the first
faces = face_cascade.detectMultiScale(gray, 1.3, 5)
if len(faces) < 1:
return None
(x, y, w, h) = faces[0]
# Crop out faces and detect eyes in that window.
roi_gray = gray[y:y + h, x:x + w]
eyes, numDetects = eye_cascade.detectMultiScale2(roi_gray, minNeighbors=10)
if len(numDetects) < 2:
return None
# Change eye coords to be in image coordinates instead of face coordinates
eyes[0][0] += x
eyes[1][0] += x
eyes[0][1] += y
eyes[1][1] += y
return [faces[0], eyes[0], eyes[1]]
def getFace(im):
gray = cv2.cvtColor(im, cv2.COLOR_BGR2GRAY)
# Only take one face, the first
faces = face_cascade.detectMultiScale(gray, 1.3, 5)
if len(faces) < 1:
return None
return faces[0]
dataLen = 120
camTimes = [0]*dataLen
intensities = []
x = list(range(len(intensities)))
def getHR():
fs = 1 / (sum(camTimes) / dataLen)
tmpIntens = sig.detrend(applyFF(intensities, fs))
freqs, pows = signal.welch(tmpIntens, fs=fs, nperseg=256)
bpm = round(freqs[np.argmax(pows)] * 60, 2)
print("output BPM: ", bpm, fs)
return bpm
def getHeadboxFromHead(face):
return face[0] + face[2] // 4, face[1] + face[3] // 2, face[0] + 3 * face[
2] // 4, face[1] + face[3] // 2 + 50
cap = cv2.VideoCapture(0)
def readIntensity(intensities, curFrame, cropBoxBounds):
now = 0
eyeleft = 0
headTop = 0
eyeright = 0
eyeTop = 0
while True:
ret, frame = cap.read()
scaleFactor = 0.4
frame = cv2.resize(frame,(-1,-1), fx=scaleFactor, fy=scaleFactor)
# tmp = getFaceXYHWAndEyeXYHW(frame) # Haar outputs [x, y, w, h] format
face = getFace(frame)
if face is not None:
# if tmp != None:
# face, eye1, eye2 = tmp
# eyeleft, headTop, eyeright, eyeTop\
# tmpHeadbox = getHeadbox(face, eye1, eye2)
tmpHeadbox = getHeadboxFromHead(face)
a = .4
eyeleft = int(tmpHeadbox[0]*a + (1-a)*eyeleft)
headTop = int(tmpHeadbox[1]*a + (1-a)*headTop)
eyeright = int(tmpHeadbox[2]*a + (1-a)*eyeright)
eyeTop = int(tmpHeadbox[3]*a + (1-a)*eyeTop)
ROI = frame[headTop:eyeTop, eyeleft:eyeright, 1]
intensity = ROI.mean()
# intensity = np.median(ROI) # works, but quite chunky.
intensities.append(intensity)
# Draw the forehead box:
curFrame[0] = cv2.rectangle(frame, (eyeleft, headTop),
(eyeright, eyeTop), (0, 255, 0), 1)
cropBoxBounds[0] = [headTop + 2, eyeTop - 2, eyeleft + 2, eyeright - 2]
if (len(intensities) > dataLen):
intensities.pop(0)
camTimes.append(time.time() - now)
now = time.time()
camTimes.pop(0)
cropBoxBounds = [0]
curFrame = [0]
t1 = threading.Thread(target = readIntensity, daemon=True, args=(intensities, curFrame, cropBoxBounds))
t1.start()
time.sleep(1)
while True:
frame = curFrame[0]
bb = cropBoxBounds[0]
ROI = frame[bb[0]:bb[1], bb[2]:bb[3], 1]
getHR()
cv2.imshow("yea", frame)
if cv2.waitKey(1) & 0xFF == ord('q'):
break
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
From the blog post here
If you want to use this gist, you'll need to download the two Haar cascades as .xmls from the Opencv Github, "haarcascade_frontalface_default.xml" and "haarcascade_eye.xml"