# Heart Rate Estimation from Video

Changes to skin color due to blood flow in the skin can be captured on video and used to find an estimate of heart rate. Although such a change is too small to be seen by the human eye, the signal can be nicely extracted and analyzed from the skin pixels in consecutive frames.

Acquire frames of a video and its acquisition times.

 In[2]:= XvideoTs = 0.001*{0, 38, 93, ..., 10203, 10258};

Determine the bounding box of the face in all frames and highlight the face for the first frame.

 In[3]:= XfaceBoxes = Map[First@FindFaces[#] &, video]; HighlightImage[video[[1]], Graphics[{EdgeForm[Orange], FaceForm[], Rectangle @@ faceBoxes[[1]]}]]
 Out[3]=

Determine the regularized motion of the bounding box and trim all frames to the detected faces using a median bounding box size.

 In[4]:= XavgFaceBoxSize = Round[{1/GoldenRatio, 1} Map[Median, Transpose@Map[Differences, faceBoxes][[All, 1]]], 2];
 In[5]:= XavgFaceBoxPositions = Transpose@ Map[GaussianFilter[#, 4] &, Transpose@Map[Mean, faceBoxes]];
 In[6]:= Xfaces = MapThread[ ImageTrim[#1, {#2}, avgFaceBoxSize/2] &, {video, avgFaceBoxPositions}];

Determine the face-shifts with respect to the first frame and stabilize the video by undoing the shifts.

 In[7]:= XfaceShifts = Map[Last@FindGeometricTransform[First@faces, #, Method -> "Linear", TransformationClass -> "Translation"] &, faces];
 In[8]:= XregFaces = MapThread[ ImagePerspectiveTransformation[#1, #2, DataRange -> All, Resampling -> "Cubic"] &, {faces, faceShifts}];

Determine the average skin probability distribution in the face bounding box using a typical skin classifier based on Lab color to detect facial skin.

 In[9]:= XskinClassifier = Compile[{{Lab, _Real, 1}}, Exp[-(Lab[[2]] - 0.25)^2/(2 0.06^2)] Exp[-(Lab[[3]] - 0.16)^2/(2 0.04^2)]];
 In[10]:= XLabFaces = Map[ColorConvert[#, "LAB"] &, regFaces];
 Out[11]=

For each frame, extract from all pixels the average Lab colors weighted by the above skin probability distribution.

 In[12]:= XLabSignal = Map[ImageMeasurements[#, "Mean", Masking -> skinWeight] &, LabFaces];
 In[13]:= XListLinePlot[Rest@Transpose@LabSignal, PlotLegends -> {"a", "b"}]
 Out[13]=

Find the optimal demixing angle α by searching for a signal with the least differential volatility.

 In[14]:= Xopt = Last@ FindMinimum[ Total[Differences[{0, Cos[\[Alpha]], Sin[\[Alpha]]}.Transpose[ LabSignal]]^2], {\[Alpha], 0}]
 Out[14]=

Extract the optimal pulse signal and regularize the pulse signal with a bandpass filter that selects frequencies between 0.5Hz and 3Hz.

 In[15]:= Xsignal = ({0, Cos[\[Alpha]], Sin[\[Alpha]]} /. opt).Transpose[ LabSignal];
 In[16]:= XfilteredSignal = BandpassFilter[signal, 2 \[Pi]/18 {0.5, 3}, 11];
 In[17]:= XListLinePlot[filteredSignal, Axes -> {Automatic, None}]
 Out[17]=

Extract the heart beats using FindPeaks.

 In[18]:= XheartBeats = FindPeaks[filteredSignal, 2, InterpolationOrder -> 3, Padding -> 0.01]
 Out[18]=
 In[19]:= XListLinePlot[filteredSignal, Epilog -> {Gray, Map[Line[{{First[#], 0}, #}] &, heartBeats]}, Axes -> {Automatic, None}]
 Out[19]=

Convert the frame numbers into times and extract heart beat intervals.

 In[20]:= XtimeFunction = ListInterpolation[videoTs, InterpolationOrder -> 1]; Ts = Map[timeFunction, heartBeats[[All, 1]]]; \[CapitalDelta]Ts = Differences[Ts]
 Out[20]=

Median heart rate and its deviation.

 In[21]:= Xm\[CapitalDelta]T = Quantity[1/Median[\[CapitalDelta]Ts], 1/"Seconds"]
 Out[21]=
 In[22]:= XUnitConvert[m\[CapitalDelta]T, 1/"Minutes"]
 Out[22]=
 In[23]:= Xd\[CapitalDelta]T = Quantity[MedianDeviation[\[CapitalDelta]Ts], "Seconds"]
 Out[23]=

## Mathematica + Mathematica Online

Questions? Comments? Contact a Wolfram expert »