Including Wang's FVP algorithm and input video format

[terbed]

Hi!

First of all thank you for your contribution to the field!

I have concerns about the input video format (.mp4) as these typically apply some compression to the images. It would be more faithful if for example the input video data also would be in .mat (uncompressed raw data). Of course a general large benchmark dataset could be constructed in this way to benchmark the different algorithms that is in line with your codes.

It would be also nice to include Wang's FVP algo.:

function H = FVP(frames, Fs)
    %FVP Full Video Pulse extraction
    %   frames: sequence of RGB images (N, h, w, c)
    %   Fs: Frame rate
    
    N = size(frames, 1);
    K = 6; 
    L1 = Fs; 
    u0 = 1; 
    L2 = 128; 
    B = [6, 24];            % pulse band (calculated for L2!)
    Jt = zeros(4*K, N, 3);
    
    Pt = zeros(4 * K, N);
    Zt = zeros(4*K, N); 
    H = zeros(1, N);
    
    for i=1:N
        % Generate weighting masks
        Id = imresize(squeeze(frames(i,:,:,:)), [20, 20], 'box'); % downsize
        D = reshape(bsxfun(@rdivide, Id, sum(Id, 3)), [size(Id, 1)*size(Id, 2), 3]);
        
        A = pdist2(D, D, 'euclidean');
        
        [u, ~, ~] = svds(A, K);
        
        % correct the eigenvector signs
        u = bsxfun(@times, u, sign(sum(u0.*u, 1))); 
        u0 = u; 
        
        w = [u(:,1:K), -u(:, 1 : K)]'; 
        w = bsxfun(@minus, w, min(w, [], 2));
        w = bsxfun(@rdivide, w, sum(w, 2));
        
        % Weight and combine spatial pixel values
        J = bsxfun(@times, reshape(Id, [1, size(Id, 1) * size(Id, 2), 3]), w);
        Jt(:, i, :) = [mean(J, 2); var(J, [], 2)];
        
        % Extract rPPg-signals (with POS)
        if(i >= L1)
            C = Jt(:, i - L1 + 1 : i, :); 
            Cn = bsxfun(@rdivide, C, mean(C, 2)) - 1;
            X = Cn(:, :, 2) - Cn(:, :, 3); Y = Cn(:, :, 2) + Cn(:, :, 3) - 2 * Cn(:, :, 1);
            P = X + bsxfun(@times, std(X, [], 2)./std(Y, [], 2), Y); 
            Z = Cn(:, :, 1) + Cn(:, :, 2) + Cn(:, :, 3);
            
            Pt(:, i - L1 + 1 : i) = Pt(:, i - L1 + 1 : i) + bsxfun(@rdivide, bsxfun(@minus, P, mean(P, 2)), std(P, [], 2));
            Zt(:, i - L1 + 1 : i) = Zt(:, i - L1 + 1 : i) + bsxfun(@rdivide, bsxfun(@minus, Z, mean(Z, 2)), std(Z, [], 2));
        end
        
        % Combine rPPG signals into pulse
        if(i >= L2)
            P = Pt(:, i - L2 + 1 : i); 
            Z = Zt(:, i - L2 + 1 : i);
            
            Fp = fft(bsxfun(@rdivide, bsxfun(@minus, P, mean(P, 2)), std(P, [], 2)), [], 2)/L2;
            Fz = fft(bsxfun(@rdivide, bsxfun(@minus, Z, mean(Z, 2)), std(Z, [], 2)), [], 2)/L2;
            
            W = abs(Fp.* conj(Fp))./(1 + abs(Fz.* conj(Fz)));
            
            W(:, 1 : B(1) - 1) = 0;
            W(:, B(2) + 1 : end) = 0;
            h = real(ifft(sum(W.* Fp, 1), [], 2));

            H(:, i - L2 + 1 : i) = H(:, i - L2 + 1 : i) + (h - mean(h))/std(h);
        end
    end
end

And an example for the usage of the code:

addpath('~/CODES/MATLAB/algorithms');

N = 6000;                             % Num of frames in the directory
Fs = 20;
L2 = 256;
stride = 128;

result = zeros(1, N);
frames = zeros(1, 500, 500, 3);

dur = 0;
fi = 1;
for i=1:N
    fname = sprintf('../PIC_12bit_191111/%06d.png', i-1);
    frames(fi, :,:,:) = imread(fname);
    fi = fi+1;
    %imshow(squeeze(frames(i, :,:,:))/max(frames(i, :,:,:), [], 'all'));
    
    if size(frames, 1) == L2
       disp('Calculating POS')
       result(i-L2+1:i) = result(i-L2+1:i) + FVP(frames, Fs);

       
       frames(1:stride,:,:,:) = [];
       fi = fi - stride; % change frame indexing according deletion
    end
    
    disp(i);
end

This is the original code he submitted in his PHD thesis, but I would apply some modifications on it. To be more specific, it is not advisable to zero out frequencies in the Fourier domain and then transform back, that part could be changed to zero phase filtering (the way you do generally on the implemented algos):

% W(:, 1 : B(1) - 1) = 0; 
%  W(:, B(2) + 1 : end) = 0; 
h = real(ifft(sum(W.* Fp, 1), [], 2));

%% Parameters
LPF = 0.7; %low cutoff frequency (Hz) - 0.8 Hz in reference
HPF = 2.5; %high cutoff frequency (Hz) - both 6.0 Hz and 2.0 Hz used in reference
%% Filter, Normalize
NyquistF = 1/2*FS;
[B,A] = butter(3,[LPF/NyquistF HPF/NyquistF]);%Butterworth 3rd order filter - originally specified in reference with a 4th order butterworth using filtfilt function
BVP_F = filtfilt(B,A,(double(h)-mean(h)));

[danmcduff]

Hi,

Thanks so much for this message. I am glad you find the toolbox useful. I'd be happy to add this in. I am at NeurIPS this week but once I am back I'll do it. Is it OK if I credit you in the file header for the contribution?

Regarding mp4 input, I agree that compression can have a big impact. I guess our goal here was to just provide an example for now of how to use the code. But uncompressed input will give the best results.

Daniel

On Sat, Dec 7, 2019, 2:12 AM Terbe Dániel [email protected] wrote:

Hi!

First of all thank you for your contribution to the field!

I have concerns about the input video format (.mp4) as these typically apply some compression to the images. It would be more faithful if for example the input video data also would be in .mat (uncompressed raw data). Of course a general large benchmark dataset could be constructed in this way to benchmark the different algorithms that is in line with your codes.

It would be also nice to include Wang's FVP algo.:

function H = FVP(frames, Fs) %FVP Full Video Pulse extraction % frames: sequence of RGB images (N, h, w, c) % Fs: Frame rate

N = size(frames, 1);
K = 6;
L1 = Fs;
u0 = 1;
L2 = 128;
B = [6, 24];            % pulse band (calculated for L2!)
Jt = zeros(4*K, N, 3);

Pt = zeros(4 * K, N);
Zt = zeros(4*K, N);
H = zeros(1, N);

for i=1:N
    % Generate weighting masks
    Id = imresize(squeeze(frames(i,:,:,:)), [20, 20], 'box'); % downsize
    D = reshape(bsxfun(@rdivide, Id, sum(Id, 3)), [size(Id, 1)*size(Id, 2), 3]);

    A = pdist2(D, D, 'euclidean');

    [u, ~, ~] = svds(A, K);

    % correct the eigenvector signs
    u = bsxfun(@times, u, sign(sum(u0.*u, 1)));
    u0 = u;

    w = [u(:,1:K), -u(:, 1 : K)]';
    w = bsxfun(@minus, w, min(w, [], 2));
    w = bsxfun(@rdivide, w, sum(w, 2));

    % Weight and combine spatial pixel values
    J = bsxfun(@times, reshape(Id, [1, size(Id, 1) * size(Id, 2), 3]), w);
    Jt(:, i, :) = [mean(J, 2); var(J, [], 2)];

    % Extract rPPg-signals (with POS)
    if(i >= L1)
        C = Jt(:, i - L1 + 1 : i, :);
        Cn = bsxfun(@rdivide, C, mean(C, 2)) - 1;
        X = Cn(:, :, 2) - Cn(:, :, 3); Y = Cn(:, :, 2) + Cn(:, :, 3) - 2 * Cn(:, :, 1);
        P = X + bsxfun(@times, std(X, [], 2)./std(Y, [], 2), Y);
        Z = Cn(:, :, 1) + Cn(:, :, 2) + Cn(:, :, 3);

        Pt(:, i - L1 + 1 : i) = Pt(:, i - L1 + 1 : i) + bsxfun(@rdivide, bsxfun(@minus, P, mean(P, 2)), std(P, [], 2));
        Zt(:, i - L1 + 1 : i) = Zt(:, i - L1 + 1 : i) + bsxfun(@rdivide, bsxfun(@minus, Z, mean(Z, 2)), std(Z, [], 2));
    end

    % Combine rPPG signals into pulse
    if(i >= L2)
        P = Pt(:, i - L2 + 1 : i);
        Z = Zt(:, i - L2 + 1 : i);

        Fp = fft(bsxfun(@rdivide, bsxfun(@minus, P, mean(P, 2)), std(P, [], 2)), [], 2)/L2;
        Fz = fft(bsxfun(@rdivide, bsxfun(@minus, Z, mean(Z, 2)), std(Z, [], 2)), [], 2)/L2;

        W = abs(Fp.* conj(Fp))./(1 + abs(Fz.* conj(Fz)));

        W(:, 1 : B(1) - 1) = 0;
        W(:, B(2) + 1 : end) = 0;
        h = real(ifft(sum(W.* Fp, 1), [], 2));

        H(:, i - L2 + 1 : i) = H(:, i - L2 + 1 : i) + (h - mean(h))/std(h);
    end
end

end

And an example for the usage of the code:

addpath('~/CODES/MATLAB/algorithms');

N = 6000; % Num of frames in the directory Fs = 20; L2 = 256; stride = 128;

result = zeros(1, N); frames = zeros(1, 500, 500, 3);

dur = 0; fi = 1; for i=1:N fname = sprintf('../PIC_12bit_191111/%06d.png', i-1); frames(fi, :,:,:) = imread(fname); fi = fi+1; %imshow(squeeze(frames(i, :,:,:))/max(frames(i, :,:,:), [], 'all'));

if size(frames, 1) == L2
   disp('Calculating POS')
   result(i-L2+1:i) = result(i-L2+1:i) + FVP(frames, Fs);


   frames(1:stride,:,:,:) = [];
   fi = fi - stride; % change frame indexing according deletion
end

disp(i);

end

This is the original code he submitted in his PHD thesis, but I would apply some modifications on it. To be more specific, it is not advisable to zero out frequencies in the Fourier domain and than transform back, that part could be changed to zero phase filtering (the way you do generally on the implemented algos):

% W(:, 1 : B(1) - 1) = 0; % W(:, B(2) + 1 : end) = 0; h = real(ifft(sum(W.* Fp, 1), [], 2));

%% Parameters LPF = 0.7; %low cutoff frequency (Hz) - 0.8 Hz in reference HPF = 2.5; %high cutoff frequency (Hz) - both 6.0 Hz and 2.0 Hz used in reference %% Filter, Normalize NyquistF = 1/2*FS; [B,A] = butter(3,[LPF/NyquistF HPF/NyquistF]);%Butterworth 3rd order filter - originally specified in reference with a 4th order butterworth using filtfilt function BVP_F = filtfilt(B,A,(double(h)-mean(h)));

I'm interested what you think about my suggestions!

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/danmcduff/iphys-toolbox/issues/3?email_source=notifications&email_token=ABLHJUWXHWD57YNUSCQ6OE3QXNZHNA5CNFSM4JXK3GCKYY3PNVWWK3TUL52HS4DFUVEXG43VMWVGG33NNVSW45C7NFSM4H62BCSQ, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABLHJUXJ2CEKUQNMM4AQWWLQXNZHNANCNFSM4JXK3GCA .

[terbed]

Hello,

Have a great conference then! I am grateful if you mention my contribution, but it is not even needed.

By the way, do you plan to continuously maintain this repo? I don't know whether you know about other initiation like this: https://www.idiap.ch/software/bob/docs/bob/docs/stable/bob/bob.rppg.base/doc/index.html#bob-rppg-base

I really support the idea, that the rPPG community should have a common tool to which everyone commits instead of attenuating in distinct directions.

Best wishes, Daniel