Matlab code to find the Wave height

I want to find the water wave height and wave length by analysing the video. I write the code for that :

Is this correct?

videoFile = 'wave_video.mp4'; % Replace with your video file name
video = VideoReader(videoFile);
numFrames = floor(video.Duration * video.FrameRate);
frameRate = video.FrameRate;
displacement = zeros(1, numFrames);
% Define a threshold value for detecting the wave (adjust this as needed)
threshold = 100;
for k = 1:numFrames
    frame = readFrame(video);
    grayFrame = rgb2gray(frame);
    
    % Threshold the grayscale image to create a binary mask
    binaryImage = grayFrame > threshold;
    
    % Sum the binary image to get an estimation of the wave's displacement
    displacement(k) = sum(binaryImage(:));
end
time = (0:numFrames-1) / frameRate;
plot(time, displacement);
xlabel('Time (s)');
ylabel('Displacement');
title('Water Wave Displacement Over Time');

How to find the unit of that height.

[Mathieu NOE]

1/ I am not sure if we deal with vertical or horizontal displacement

2/ nor I understand the relationship between color and displacement

3/ maybe you need a "calibration" file (so with known water displacement measured somehow by another mean you can trust)

[Mathieu NOE]

maybe if you could share some video file (a short one) , someone can work on your problem (unless its obvious there's already an answer somewhere)

[Prabin]

1. I am dealing with verticle displacement ie. dislacement with wall.
2. My requirement is to find water wave properties of the swimming pool. Wavelength, displacement.

Answer 1

[Walter Roberson]

image_height = min(sum(cumprod(~binaryImage, 1),1) + 1);

or possibly

image_height = min(sum(cumprod(~flipud(binaryImage), 1),1) + 1);

depending on whether the image starts from the top or from the bottom of the array.

The logic here assumes that the portions without the wave are binary 0. The cumprod(~binaryImage) fills the leading zero elements of each column of the array, producing a "top hat" sort of effect. sum() of that gives the number of leading zeros in each column. min() of that finds the shortest entry -- corresponding to the place the hat was highest.

Perhaps max() would be needed instead of min(), depending how the data is arranged in the images.

[Prabin]

https://drive.google.com/file/d/1n-xEMy44Pq8Z6S6DutCisZ0k7aAwO86o/view?usp=drive_link

i have this video and i need both minimum and maximum wave height. Still i tried a lot to find the solution . i did not get correct graph.

Answer 2

[Mathieu NOE]

hello again

I am by no means an image processing expert, but sometimes I go outside my comfort zone ... hopefully my suggestion is not completely wrong.

One point, all the results given below are "in pixels" dimnesions, so as I commented in first place , there is still some work to do to convert these results into physical units. A very crude first approach would be to identify a know object in the picture and use its dimensions (known) to apply a conversion rule. a more sophisticated approach would require to estimate the position and orientation of the camera , but Ileave that to real professionnals ....

so what can I offer today ?

1/ I tried first to detect the line between water and border (arrow). So there is a lower and upper threshold used on the RGB triplet to convert the image into a binary matrix (water or not). To ease my work I cropped the image to what is inside the red rectangle

2/ then extract the lower boundary of the "non water" objet (ie the border) , detrend it , remove some outliers , and some smoothing on top , and here we get for every frame a more or less accurate image on the water line (for each frame) - this is plotted in the bottom subplot here - the top suplot is the binarized image with the separation in red

3/ from there you can compute some values like mean frequency (interest ?) and evidetly the max / min (= amplitude) of the wave

    % mean period (in pixels) 
    mean_period_in_pixels(k)= 1/meanfreq(y3ds,1); % Fs spatial = 1 pixel
    
    % max wave amplitude (in pixels) 
    amplitude_in_pixels(k)= max(y3ds) - min(y3ds); 

full code :

videoFile = 'wave_video1.mp4'; % Replace with your video file name
video = VideoReader(videoFile);
numFrames = floor(video.Duration * video.FrameRate);
frameRate = video.FrameRate;
% Define a threshold value for detecting the wave (adjust this as needed)
threshold1 = 160;
threshold2 = 190;

for k = 1:numFrames
    frame = readFrame(video);
    ind1 = (frame(:,:,1)>threshold1) & (frame(:,:,2)>threshold1) & (frame(:,:,3)>threshold1);
    ind2 = (frame(:,:,1)<threshold2) & (frame(:,:,2)<threshold2) & (frame(:,:,3)<threshold2);
    binaryImage = ind1 & ind2; % here we get the white/gray border
    
    % crop the bottom of the pictures (to remove the light reflection spots)
    binaryImage = binaryImage(1:200,150:500);
    
    [y,x] = find(binaryImage); % get x,y coordinates 
    
    [x3,y3,x4,y4] = top_bottom_boundary(x,y); %find the separation line between water and border 
    y3 = filloutliers(y3,'linear','movmedian',50);        % remove some outliers
    y3d = detrend(y3,'linear'); % remove linear trend (due to camera orientation)
    y3ds = smoothdata(y3d,'sgolay',50); % smooth a bit the line

        % some plots for debug and demo 
        if k<3 % plot as many as you want
            figure(k),
            subplot(211),imagesc(binaryImage)
            hold on 
            plot(x3,y3,'r*');
            hold off
            subplot(212),plot(x3,y3d,'r*',x3,y3ds,'b');
        end
    
    % mean period (in pixels) 
    mean_period_in_pixels(k)= 1/meanfreq(y3ds,1); % Fs spatial = 1 pixel
    
    % max wave amplitude (in pixels) 
    amplitude_in_pixels(k)= max(y3ds) - min(y3ds); 
    
end

time = (0:numFrames-1) / frameRate;
figure
plot(time(1:k), amplitude_in_pixels);
xlabel('Time (s)');
ylabel('Displacement (in pixels) ');
title('Water Wave Displacement Over Time');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x3,y3,x4,y4] = top_bottom_boundary(x,y)

% based on FEX : https://fr.mathworks.com/matlabcentral/answers/299796-tight-boundary-around-a-set-of-points


    %split data into classes and find max/min for each class
    class_label = unique(x);
    upper_boundary = zeros(size(class_label));
    lower_boundary = zeros(size(class_label));
    for idx = 1:numel(class_label)
        class = y(x == class_label(idx));
        upper_boundary(idx) = max(class);
        lower_boundary(idx) = min(class);
    end

%     left_boundary = y(x == class_label(1));
%     right_boundary = y(x == class_label(end));
% 
%     % left_boundary
%     x1 = class_label(1)*ones(size(left_boundary));
%     y1 = left_boundary;
% 
%     % right_boundary
%     x2 = class_label(end)*ones(size(right_boundary));
%     y2 = right_boundary;

    % top boundary
    x3 = class_label;
    y3 = upper_boundary;

    % bottom boundary
    x4 = class_label;
    y4 = lower_boundary;

%     x_out = [x1(:);x2(:);x3(:);x4(:);];
%     y_out = [y1(:);y2(:);y3(:);y4(:);];

end

[Mathieu NOE]

here a slight improvement over my previous code

the line angle correction is done by rotation and not by linear detrending (which is not a geometrical correction)

    % correct tilting angle  (due to camera orientation)
    p = polyfit(x3,y3,1);
    slope = p(1);
    alpha = -atan(slope);
    x3r = x3*cos(alpha) - (y3-y3(1))*sin(alpha);
    y3r = (x3-x3(1))*sin(alpha) + y3*cos(alpha);

hope it helps

videoFile = 'wave_video1.mp4'; % Replace with your video file name
video = VideoReader(videoFile);
numFrames = floor(video.Duration * video.FrameRate);
frameRate = video.FrameRate;
% Define a threshold value for detecting the wave (adjust this as needed)
threshold1 = 160;
threshold2 = 190;

for k = 1:numFrames
    frame = readFrame(video);
    ind1 = (frame(:,:,1)>threshold1) & (frame(:,:,2)>threshold1) & (frame(:,:,3)>threshold1);
    ind2 = (frame(:,:,1)<threshold2) & (frame(:,:,2)<threshold2) & (frame(:,:,3)<threshold2);
    binaryImage = ind1 & ind2; % here we get the white/gray border
    
    % crop the bottom of the pictures (to remove the light reflection spots)
    binaryImage = binaryImage(1:200,150:500);
    
    [y,x] = find(binaryImage); % get x,y coordinates 
    
    [x3,y3,x4,y4] = top_bottom_boundary(x,y); %find the separation line between water and border 
    y3 = filloutliers(y3,'linear','movmedian',50);        % remove some outliers
    
    % correct tilting angle  (due to camera orientation)
    p = polyfit(x3,y3,1);
    slope = p(1);
    alpha = -atan(slope);
    x3r = x3*cos(alpha) - (y3-y3(1))*sin(alpha);
    y3r = (x3-x3(1))*sin(alpha) + y3*cos(alpha);
    
    y3r = y3r - mean(y3r); % remove dc value
    y3rs = smoothdata(y3r,'sgolay',50); % smooth a bit the line

        % some plots for debug and demo 
        if k<3 % plot as many as you want
            figure(k),
            subplot(211),imagesc(binaryImage)
            hold on 
            plot(x3,y3,'r*');
            hold off
            subplot(212),plot(x3r,y3r,'m*',x3,y3rs,'b');
        end
    
    % mean period (in pixels) 
    mean_period_in_pixels(k)= 1/meanfreq(y3rs,1); % Fs spatial = 1 pixel
    
    % max wave amplitude (in pixels) 
    amplitude_in_pixels(k)= max(y3rs) - min(y3rs); 
    
end

time = (0:numFrames-1) / frameRate;
figure
plot(time(1:k), amplitude_in_pixels);
xlabel('Time (s)');
ylabel('Displacement (in pixels) ');
title('Water Wave Displacement Over Time');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x3,y3,x4,y4] = top_bottom_boundary(x,y)

% based on FEX : https://fr.mathworks.com/matlabcentral/answers/299796-tight-boundary-around-a-set-of-points


    %split data into classes and find max/min for each class
    class_label = unique(x);
    upper_boundary = zeros(size(class_label));
    lower_boundary = zeros(size(class_label));
    for idx = 1:numel(class_label)
        class = y(x == class_label(idx));
        upper_boundary(idx) = max(class);
        lower_boundary(idx) = min(class);
    end

%     left_boundary = y(x == class_label(1));
%     right_boundary = y(x == class_label(end));
% 
%     % left_boundary
%     x1 = class_label(1)*ones(size(left_boundary));
%     y1 = left_boundary;
% 
%     % right_boundary
%     x2 = class_label(end)*ones(size(right_boundary));
%     y2 = right_boundary;

    % top boundary
    x3 = class_label;
    y3 = upper_boundary;

    % bottom boundary
    x4 = class_label;
    y4 = lower_boundary;

%     x_out = [x1(:);x2(:);x3(:);x4(:);];
%     y_out = [y1(:);y2(:);y3(:);y4(:);];

end