Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Webcam Photoplethysmyography
# 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
@kongmunist

This comment has been minimized.

Copy link
Owner Author

@kongmunist kongmunist commented Dec 26, 2020

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"

@ytmytm

This comment has been minimized.

Copy link

@ytmytm ytmytm commented Dec 27, 2020

On Windows I had to add cap = cv2.VideoCapture(0,cv2.CAP_DSHOW) option and additional delay (2-3s) after starting t1 thread

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment