%%%%%%%%%%%%%%%%%
% 
% 
% Task: Detect heart rate through video
%
%
%
% Authord: Estevan Biau-Loyer & Darren Gallois
% Date: 19/02/2023
%
%%%%%%%%%%%%%%%%%


% Loop through the images selected

% Initialize a matrix to store RGB values
all_rgb = zeros(900, 3, 61, 101);

% Initialize iteration
a=0;
for i = 100:1000
  % Set the image path
  image_path = imread(sprintf('frame2_%d.jpg', i));
  [rows, columns, channels] = size(image_path);
  %Create the ROI
  roi = image_path(90:150, 300:400, :);
  % Save the images adjusted to the ROI
  imwrite(roi, (sprintf('frames/my_imagev2%d.jpg',i)));
end

% Loop through 900 images
for i = 100:1000
  % Read image
  img = imread(sprintf('frames/my_imagev2%d.jpg', i));

  % Convert to double precision
  img = im2double(img);

  % Separate RGB channels
  red = img(:,:,1);
  green = img(:,:,2);
  blue = img(:,:,3);

  % Store RGB channels in the matrix
  all_rgb(i, 1, :, :) = red;
  all_rgb(i, 2, :, :) = green;
  all_rgb(i, 3, :, :) = blue;
  a=a+1;
  
  % Display iteration
  disp(a);
end

% Compute mean RGB values for each image
mean_rgb = squeeze(mean(mean(all_rgb, 4), 3));
save('mean_rgb.m', 'mean_rgb');

% Load the mean RGB values from the .m file
load('mean_rgb.m');

% Get the number of images
num_images = size(mean_rgb, 1);

% Create a vector of indices for each image
indices = 1:num_images;

% Plot the green channel
figure;
plot(indices, mean_rgb(:, 1), 'r');
hold on;

% Plot the red channel
plot(indices, mean_rgb(:, 2), 'g');

% Plot the blue channel
plot(indices, mean_rgb(:, 3), 'b');

% Add legend and labels
legend({'Red', 'Green', 'Blue'});
xlabel('Image index');
ylabel('Mean RGB value');

% Get the standard deviation of each color
red_std = std(mean_rgb(:, 1));
green_std = std(mean_rgb(:, 2));
blue_std = std(mean_rgb(:, 3));

% Get the average value of each color
red_avg = mean(mean_rgb(:, 1));
green_avg = mean(mean_rgb(:, 2));
blue_avg = mean(mean_rgb(:, 3));

% Normalize the values
red_normalized = (mean_rgb(:, 1) - red_avg)/red_std
green_normalized = (mean_rgb(:, 2) - green_avg)/green_std
blue_normalized = (mean_rgb(:, 3) - blue_avg)/blue_std

% Create the 3 channels fft
red_fft = fft(red_normalized)
green_fft = fft(green_normalized)
blue_fft = fft(blue_normalized)

N = length(red_normalized); % length of the signal
freq = fftshift(fftfreq(N)); % frequency vector in Hz

% plot the FFT result for the three channels
figure;
hold on;
plot(freq, abs(red_fft));

plot(freq, abs(green_fft));

plot(freq, abs(blue_fft));
hold off;

xlabel('Frequency (Hz)');
ylabel('Magnitude');
legend({'Red', 'Green', 'Blue'});

% Find indices in f that correspond to the frequency range of interest
idx_range = find(f >= 0.75 & f <= 4);

% Extract the amplitude values in the range of interest for each channel
red_amp_range = abs(red_fft(idx_range));
green_amp_range = abs(green_fft(idx_range));
blue_amp_range = abs(blue_fft(idx_range));

% Find the index of the maximum amplitude value in the range of interest for each channel
[red_max_amp, red_max_idx] = max(red_amp_range);
[green_max_amp, green_max_idx] = max(green_amp_range);
[blue_max_amp, blue_max_idx] = max(blue_amp_range);

% Convert index of maximum amplitude value to frequency, and then to bpm
red_bpm = f(idx_range(red_max_idx)) * 60;
green_bpm = f(idx_range(green_max_idx)) * 60;
blue_bpm = f(idx_range(blue_max_idx)) * 60;

% Determine heart rate by taking the average of the bpm values for the three channels
heart_rate = mean([red_bpm, green_bpm, blue_bpm]);
disp(heart_rate);