%%
%% Estimate of the pixels of the river
%%
%% Segmentation using the topological relationship on the image structure.

filename='satellite.jpg';

PRINT_FLAG = 0;
FONTSIZE = 16;

% Load image
I = imread(filename);
t = rgb2gray(I); % just to be sure that it is a gray levels image
figure; imshow(t);
if PRINT_FLAG, imwrite(t, 'gray_levels.png'); end;
figure; imhist(t);
axis([0 255 0 1600]);
if PRINT_FLAG, print -dpng 'histogram'; end;

if(0),
    figure; imhist(t);
    axis([0 255 0 1600]);
    hold on;
    plot([100 100], [0 200], 'r-', 'LineWidth', 2);
    if PRINT_FLAG, print -dpng 'histogram1'; end;
    t2 = im2bw(t, 100/255);
    figure; imshow(t2);
    if PRINT_FLAG, imwrite(t2, 'bw1.png'); end;

    figure; imhist(t);
    axis([0 255 0 1600]);
    hold on;
    plot([150 150], [0 200], 'r-', 'LineWidth', 2);
    if PRINT_FLAG, print -dpng 'histogram1'; end;
    t2 = im2bw(t, 150/255);
    figure; imshow(t2);
    if PRINT_FLAG, imwrite(t2, 'bw2.png'); end;
end;

figure; imhist(t);
axis([0 255 0 1600]);
hold on;
plot([120 120], [0 200], 'r-', 'LineWidth', 2);
if PRINT_FLAG, print -dpng 'histogram1'; end;
t2 = im2bw(t, 120/255);
figure; imshow(t2);
if PRINT_FLAG, imwrite(t2, 'bw3.png'); end;

cc4 = bwconncomp(~t2,4);
labeled = labelmatrix(cc4);
RGB_label_4 = label2rgb(labeled, 'jet', 'c', 'shuffle');
figure;
imshow(RGB_label_4,'InitialMagnification','fit');
title(['4-adjacency, ' int2str(cc4.NumObjects) ' elements'],...
      'FontSize', FONTSIZE);
if PRINT_FLAG, imwrite(RGB_label_4, 'comp_conn_4.png'); end;

if(0),
    cc8 = bwconncomp(~t2,8);
    labeled = labelmatrix(cc8);
    RGB_label_8 = label2rgb(labeled, 'jet', 'c', 'shuffle');
    figure;
    imshow(RGB_label_8,'InitialMagnification','fit');
    title(['8-adjacency, ' int2str(cc8.NumObjects) ' elements']);
    if PRINT_FLAG, imwrite(RGB_label_8, 'comp_conn_8.png'); end;
end;

% Compute the distance between the connected components
num_cc = length(cc4.PixelIdxList);
num_r = cc4.ImageSize(1);
num_c = cc4.ImageSize(2);
dist4 = zeros(num_cc);
close all;
max_dist = 10; 
for i = 1:num_cc,
    disp(i);
    m = zeros(cc4.ImageSize);
    dist = zeros(cc4.ImageSize);
    m(cc4.PixelIdxList{i}) = 1;
    % Build the distance matrix
    for d = 1:max_dist,
        n1 = [m(2:end,:); zeros(1,num_c)];
        n2 = [zeros(1,num_c); m(1:end-1,:)];
        n3 = [zeros(num_r,1) m(:,1:end-1)];
        n4 = [m(:,2:end) zeros(num_r,1)];
        neighbors = (n1 | n2 | n3 | n4) & not(m);
        dist = dist + d*neighbors;
        m = m | neighbors;
    end;
    dist = dist + (max_dist+1)*not(m);
    for j = 1:num_cc,
      dist4(i,j) = min(min(dist(cc4.PixelIdxList{j})));  
    end;
end;

% Starting from the larger connected components, add to the list
% of the connected components that compose the river 

dim = zeros(1,num_cc);
for i = 1: num_cc,
    dim(i) = length(cc4.PixelIdxList{i});
end;
% large regions (more than min_dim pixels) should belong to the the river 
min_dim = 20; 
list = find(dim > min_dim);
visited = zeros(1,num_cc); % flag: has the corresponding connected
                           % component already been considered?
river = [];  
min_dist = 4; % it must be strictly smaller than max_dist
while not(isempty(list)),
    i = list(1);
    list(1) = []; % delete the considered component from the list
    if not(visited(i)),
        visited(i) = 1;
        river(end+1) = i;
        ind = find((dist4(i,:)<=min_dist) & (visited==0));
        list = [list ind]; % add the components that are enough
                           % large and are not been considered yet
    end;
end;

% new image for the river
river_map = uint8(zeros(cc4.ImageSize));
for i = 1:length(river),
    river_map(cc4.PixelIdxList{river(i)}) = 255;
end;
figure;
imshow(river_map);
title('Pixels representing the river');
