Image Processing of Water Run-up

I need help in conducting image processing on the attached images, where I am interested in automating the process of measuring the water run-up on the front of the cylinder.

Also I am trying to use the first image (0540) as a the calibration.

Any help on this would be much appreciated.

I have put up only 3 images, but I have many, many more (5000+).

[KSSV]

Where are the images?

[Image Analyst]

I had trouble getting the images from your googledrive. Using the zip extractor app, it got stuck on the "Checking authorization..." step. Please attach your images here, not on some third party web site. Make it easy for us to help you, not hard.

[Doug]

Sorry, the file size due to the number of frames was clearly too much.

I have now uploaded 3 frames, which should be enough.

[Doug]

Hey, just wondering if you had a chance to look at my images?

Answer 1

[Image Analyst]

Doug - sorry, had forgotten about it. Have you tried anything in the last 9 days? Well I threw together this in about half an hour and it looks like it does a decent job.

% Demo to find the waterline of a striped pole.
clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 16;

%===============================================================================
% Get the name of the image the user wants to use.
baseFileName = 'Frame 0540.PNG';
% baseFileName = 'Frame 0600.PNG';
% baseFileName = 'Frame 0650.PNG';
folder = pwd;
fullFileName = fullfile(folder, baseFileName);

% Check if file exists.
if ~exist(fullFileName, 'file')
	% The file doesn't exist -- didn't find it there in that folder.
	% Check the entire search path (other folders) for the file by stripping off the folder.
	fullFileNameOnSearchPath = baseFileName; % No path this time.
	if ~exist(fullFileNameOnSearchPath, 'file')
		% Still didn't find it.  Alert user.
		errorMessage = sprintf('Error: %s does not exist in the search path folders.', fullFileName);
		uiwait(warndlg(errorMessage));
		return;
	end
end

%=======================================================================================
% Read in demo image.
rgbImage = imread(fullFileName);
% Get the dimensions of the image.
[rows, columns, numberOfColorChannels] = size(rgbImage)

% Display image.
subplot(2, 3, 1);
imshow(rgbImage, []);
impixelinfo;
axis on;
caption = sprintf('Original Color Image\n%s', baseFileName);
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
hp = impixelinfo(); % Set up status line to see values when you mouse over the image.

% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0.05 1 0.95]);
% Get rid of tool bar and pulldown menus that are along top of figure.
% set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')
drawnow;

% Get the reference mask.
% Read in demo image.
baseFileName = 'Frame reference.PNG';
folder = pwd;
fullFileName = fullfile(folder, baseFileName);
refImage = imread(fullFileName);

% Display the image.
subplot(2, 3, 2);
imshow(refImage);
caption = sprintf('Reference Mask Image');
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;
axis('on', 'image');
drawnow;

% Get a binary mask.
mask = max(refImage, [], 3) > 0;
mask = imfill(mask, 'holes');
% Display the image.
subplot(2, 3, 3);
imshow(mask);
title('Mask', 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;
axis('on', 'image');

% Subtract the two
diffImage = imabsdiff(rgbImage, refImage);
% Mask it to get rid of the surround.
% Mask the image using bsxfun() function to multiply the mask by each channel individually.
diffImage = bsxfun(@times, diffImage, cast(mask, 'like', diffImage));
drawnow;

% Display the image.
subplot(2, 3, 4);
imshow(diffImage);
caption = sprintf('Difference Image');
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;
axis('on', 'image');

% Now, on to determining the height of the water.
% Get a vertical profile of the difference image.
diffImage = rgb2gray(diffImage); % Convert to gray scale.
% Get vertical profile
verticalProfile = sum(diffImage, 2);
% Display the vertical profile.
subplot(2, 3, 5);
plot(verticalProfile, 'b-', 'LineWidth', 2);
grid on;
xlabel('Row', 'FontSize', fontSize);
ylabel('Sum', 'FontSize', fontSize);
title('Vertical Profile', 'FontSize', fontSize);

% Find regions where the sum exceeds some threshold, like 1000
threshold = 1100;  % May need to adjust this.
highSum = verticalProfile > threshold;
hold on
line(xlim, [threshold, threshold], 'Color', 'r', 'LineWidth', 2);
% This may result in several regions.
% In that case, we want to take the bottom most one.
% Find the last one
props = regionprops(highSum, 'PixelList');
% Optional, put a red cross at every peak, just for fun.
for k = 1 : length(props)
	topRow = props(k).PixelList(1,2);
	bottomRow = props(k).PixelList(end, 2);
	y = verticalProfile(round(topRow));
	% Plot the beginning with a red up arrow
	plot(topRow, y, 'r+', 'MarkerSize', 10, 'LineWidth', 2);
	% Plot the ending with a magenta down arrow
	plot(bottomRow, y, 'mv', 'MarkerSize', 10, 'LineWidth', 2);
end

%=========================================================================
% FIND THE WATER LEVEL.
% Get the top most level of the bottom most spike.
waterLevel = props(end).PixelList(1, 2)

% Put up the original image in slot 6
% Display image.
subplot(2, 3, 6);
imshow(rgbImage, []);
impixelinfo;
axis on;
caption = sprintf('Original Color Image\nWater line found at row %d', waterLevel);
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
hp = impixelinfo(); % Set up status line to see values when you mouse over the image.

% Draw the waterline across it.
hold on;
line([1, columns], [waterLevel, waterLevel], 'Color', 'r', 'LineWidth', 2);

% Tell user the answer.
message = sprintf('The water level is at row %d', waterLevel);
uiwait(helpdlg(message));

Basically I create a mask from your reference image. Then I subtract the test image from teh reference image and erase outside the mask. Then I get a vertical profile to see where they are different, like where the water is turbulent. The top most part of the turbulent area is the water level, and I draw a red line across that level. It seems to work well with the waterlines being at 514, then 505, then 483 for your 3 uploaded images as the waterline rises towards the top of the image.

Now you're going to have to develop some sort of calibration where you know what row number in the image corresponds to what real world distance in meters.

Other things that may complicate things are the lighting and the shadow that seems to come across the top of the striped pole. If you snap a reference image close in time to the other images, then you will probably be okay. However if you have fast moving clouds causing the pole to change appearance a lot, then that could affect where it finds the waterline. A more robust way might be some kind of electronic sensor, like a float, to find waterlevel.

[Doug]

Hey, thanks for the reply. Unfortunately I've been caught up with some other work, so I have not done any further work on this.

What you've done so far is incredible.

I have run it straight away but I dont get the same plots as you.

Yes I agree with the calibration, I know that each stripe is 10 mm thick. Are you able to help with that? Say at the start I choose the zero point on the base image and then it calculates the height from that?

Additonally, how would I get it to loop through 5000 frames while still using the single base file image as the zero?

Would I be able to get your email, so that I may contact you directly?

Doug

[Image Analyst]

You can contact me here. I don't know why you wouldn't get the same plots as me with the same images I used.

Basically you need to make up some look up table that says what mm it is for what line of the image it is. Then when you get the line of the image, you just look up the mm. You can do that manually, or narrow the mask and make some kind of automated way to get the middle of the mask and threshold it. One problem is that the pole is curved and if the height is affected by what's on the edge of the pole, then you have different heights in the middle and the sides of the pole. If you want to concentrate on just the middle and ignore the sides, then narrow the mask a little bit (I used Photoshop to create the mask).

[Doug]

Okay, thanks.

Are you able to help in setting up your latest suggestions?

I am very new to matlab, especially the image processing toolbox.

In terms of the point at which I would like to measure, it is at the middle of the pole (that's where the runup is greatest)

[Image Analyst]

Try my attached spatial calibration demo. Basically draw a line from one stripe to another and tell it how many mm it is. It will give you the calibration factor.

For example if one stripe is at line 222 and 30 stripes down is at line 434, that's 434 -222 = 214 lines per 300 mm. So the mm per pixel is 300/214 = 1.4018691588785 mm per pixel.

Then, knowing that and how far down your "zero point" is, you can find the distance. For example, if you have 1.4 mm per pixel, and the top of your pole ("zero") is at line 63, and the water line is at line 514, then there are 514-63 = 451 lines down from your origin. Well, 451 pixels * 1.4 mm/pixel = 632 mm from the top origin (place that is defined to be zero).

[Doug]

Going back to your original answer.

• You said that you made the mask using photoshop, would this then mean that I require an additional input file being the mask? (are you able to flick through that mask so I can check the code)
• How would I get the code to keep say Frame 0540 as my reference image and then from that run through a further 3000+ frames?

I am quite lost as to how I would apply your calibration demo.

[Image Analyst]

1. Yes, it requires a mask derived from the image where you know the water is at the reference level, like when the ship is not moving.
2. Nothing magic about frame 540. I just used it to create the mask. Once you have the reference mask, you can use it on any frame from the first one to the last one in the video.

To use the calibration demo, replace the cameraman image with a reference from your video, Then follow instructions. Draw the reference line and enter the reference distance. It will tell you how many mm per pixel. Then transfer this number into your measurement code where you take the length in pixels and multiply by the number of mm per pixel that you got from the calibration program. See my last comment. I actually did that, using your comment where one stripe was 10 mm and I got 1.4018691588785 mm per pixel. See if you get the same number from the calibration demo.

[Doug]

Hi again,

Sorry it's been a while since I have looked at this.

I am still unsure about:

• How I would make the reference mask in MATLAB?
• How to keep the mask and then loop through all of the images.
• How do I integrate the calubration function and outputs into the one script?
• How do I change the verticl profile from being for each separate image to a running plot of what all of the images are so far? as in each point represents a different image

Thanks