Detecting water level in an image
Behold the power of shortest path.

The solution roughly follows what I did here: Detecting and isolating part of an image.

In this case we first treat the image by limiting every pixel value above the water-level grey to that same intensity.



Then performing shortest path on the inverse edge image created with a difference of Gaussian filter we get:



Resulting in the following segmentation:



My MATLAB functions that are needed for this solution:

http://imageprocessing.com.au/research/code/polarTransform.m http://imageprocessing.com.au/research/code/circularShortestPath.m http://imageprocessing.com.au/research/code/linearShortestPath.m

Code:

clear all
close all

%% Load image
I = imread('CabQx.jpg');
I = double(rgb2gray(I));
figure
imagesc(I);
colormap gray(256); 
title('Image');

%% Adjust intensity
Ia = I;
Ia(Ia &gt; 49) = 49;
figure
imagesc(Ia);
colormap gray(256); 
title('Adjusted Image');

%% Bandpass Kernel
H1 = fspecial('gaussian',[100,100],8); % *** Changed filter size
H2 = fspecial('gaussian',[100,100],2);  % *** 
H = H2-H1;
H = mat2gray(H);
figure;
imagesc(H);
colorbar;
colormap gray(256); 
title('Bandpass kernel');

%% Edge Kernel
Hedge = imag(hilbert(H));
figure;
imagesc(Hedge);
colormap gray(256); 
title('Edge Kernel');

%% Edge filtered image
G = sqrt(imfilter(Ia,Hedge,'replicate').^2 + imfilter(Ia,Hedge.','replicate').^2);
G = max(G(:)) - G;
imagesc(G);
colormap gray(256); 
title('Inverse edge filtered image');

%% Shortest path
% Mid point
[r,c] = size(G);
r = floor(r/2);
c = floor(c/2);
% Min radius and max radius
min_radius = 100;
max_radius = 1000;
% Shortest path
[ path, energy, lastIndex, pathImage ] = ...
    circularShortestPath(G, 100, [c,r], [min_radius,max_radius], [400,720*2]);
figure;
imagesc(I);
hold on;
plot(path(:,1),path(:,2),'r','LineWidth',2);
hold off;
colormap gray(256); 
title('Shortest path');

%% Segment
figure
imagesc(I);
hold on;
fill(path(:,1),path(:,2),'r');
hold off;
colormap gray(256); 
title('Mask')

%%
J = zeros(size(I));
mask = roipoly(J,path(:,1),path(:,2));
segmentedI = I .* mask;
figure;
imagesc(segmentedI);
colormap gray(256); 
title('Segmented image');