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
Copy link
Author

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
Copy link

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

@kongmunist
Copy link
Author

kongmunist commented Aug 9, 2021

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