%hr = heartRateEstimation(imgDirectory, fps, windowDuration, windowShift)
%%%%%%%%%%%%%%%%%
% hr = heartRateEstimation(imgDirectory, fps, windowDuration, windowShift)
% 
% Task:
%
% Inputs:
%	-imgDirectory
% 	-fps
%	-windowDuration
%	-windowShift
%
% Output:
%	-hr:
%
% Author: Guillaume Gibert
% Date: 06/02/2023
%
% Note: images were extracted with ffmpeg -i "video.mkv" ../img/gg_%04d.png
%%%%%%%%%%%%%%%%%

%temp var
imgDirectory = 'img';
fps = 60;
windowDuration = 30;
windowShift = 1;
roi = [225 750 825 1100];

% global 
redChannel = [];
greenChannel = [];
blueChannel = [];

%list the images
listImg = dir([imgDirectory '/*.png']);

for l_img=1:windowDuration*fps
	% load the current image
	image_original = imread([imgDirectory '/' listImg(l_img).name]);
	%image(image_original);
	
	% crop the current image around the face
	image_face = image_original(roi(1):roi(2), roi(3):roi(4), :);
	%image(image_face);
	
	% spatial average for r, g, b channels
	r_mean = mean(mean(image_face(:,:,1)));
	g_mean = mean(mean(image_face(:,:,2)));
	b_mean = mean(mean(image_face(:,:,3)));
	
	% store the current "average" pixel in global arrays
	redChannel = [redChannel r_mean];
	greenChannel = [greenChannel g_mean];
	blueChannel = [blueChannel b_mean];
end

% estimate temporal average and standard deviation
greenChannel_avg = mean(greenChannel);
greenChannel_std = std(greenChannel);

% normalize your data 
% greenChannel_normalized(i) = (greenChannel(i) - greenChannel_avg)/greenChannel_std

% fft(greenChannel_normalized)

% power spectrum (https://www.mathworks.com/help/matlab/ref/fft.html)

% find the peak in the range ([0.75 4] Hz)
% max -> value, index

% convert the index from Hz to bpm

% determine heart rate