How can I automatically extract banklines from a binary image of a meandering river?

Hi, I am new to image processing and have an interesting problem that I'm hoping to get guidance on before I wander aimlessly through the many possible solutions. I am attempting to automatically extract banklines from a meandering river in the Amazon from Landsat imagery. I already have an algorithm to classify each pixel as either part of river (1) or not part of river (0). Of course, the algorithm returns a noisy image, and there are natural features (tributaries, connecting channels, point bars, etc.) that complicate the process. I have many Landsat scenes to process, and each scene has 30+ years of data with ~1 realization/year, so obviously I'd like to automate the process as much as possible. The desired ultimate output is a set of two polylines for each image: one tracing the left bank, and one the right bank. (I have tools to compute channel widths and migration rates that require vector inputs.)

Here I will show an attempt I made using Matlab's morphological operations and to give you an idea of the kinds of problems I face. Any advice is much appreciated!

Here is the classified image (i.e., what I'm starting with):

I would like to extract the left and right banklines of the main channel. Resolution is 30m, so I'd like to retain as much fidelity to the original image as possible so I can accurately compute migration rates. You can see some of the problems already: small connecting channels between lobes, chunks of 0's at the apex of the bends, some places where the river splits into two roughly-equal channels before rejoining, tributaries, mid-channel islands, etc. I would be happy to just get the inner and outer banks, although it would be nice to (automatically) detect where tributaries enter, or where a mid-channel island was...

So as a first-cut, I ran the following codes:

load BWo
BW1 = bwareaopen(BWo, 20000); % remove small patches of 1s
BW2 = bwmorph(BW1,'erode',2);
BW3 = bwareaopen(BW2, 20000); % remove small patches of 1s
BWclean = bwmorph(BW3, 'dilate',2);

% The rest of this code is just to add back some of the details lost above
% Now subtract cleaned image from original
BWdiff = logical(imsubtract(double(BWo),double(BWclean)));

% Remove large pixel groups (tribs, connecting channels) from BWdiff
CC = bwconncomp(BWdiff);
numPixels = cellfun(@numel,CC.PixelIdxList);
idx = find(numPixels>250);
for g = 1:numel(idx)
    BWdiff(CC.PixelIdxList{idx(g)}) = 0;
end

% Add the small groups back to cleaned image to better resolve banklines
BW4 = BWclean | BWdiff;

% Remove the little stuff not adjacent to river
BWf = bwareaopen(BW4, 20000);

And this is what resulted:

This result isn't perfect, but it's pretty good (and I have been playing with ways to reinsert the lost features adjacent to the channels). Now I am stuck on how to extract the banklines from this image (it would be nice to have a centerline, too, but not necessary). I know how to get the extents using bwperimm but how can I separate into left and right banks?

I also have a few other general questions about my methodology; any advice from more experienced users would be greatly appreciated!

1. The above code works for that particular scene, but I have to adjust the order and number of operations for other scenes (even for the same scene, different year). Is there a more robust way to go about extracting the main channel from scenes that have different features?

2. (Asked above) After de-noising the image, how can I pull out the left and right banks automatically?

3. Is there better methodology that would preserve more of the original image? For instance, I thought of constructing left and right banks by creating some kind of cost function and finding geodesic curves (one left, one right) that minimize that cost function...

4. Would you take a totally different approach?

Any hints, tips, tricks, or advice would be great!

Answer 1

[Image Analyst]

First of all, get rid of islands and noise by calling imfill

binaryImage = imfill(binaryImage, 'holes');

Then call bwboundaries() to get the right and left banks.

boundaries = bwboundaries(binaryImage);

Results are points in the form (row, column), NOT (x,y). In other words the x values of the k'th boundary are

thisBoundary = boundaries{k};
x = thisBoundary(:, 2); % Columns.
y = thisBoundary(:, 1); % Rows.

You're going to have to figure out where to split it apart by looking for the max and min row value.

yTop = ceil(min(y));
yBottom = ceil(max(y));
indexOfTop = y == yTop;
indexOfBottom = y == yBottom;

The thing is we don't know where point #1 is on the boundary, and that makes it complicated. So what can you do? Well, you can use a trick I employed in my maze-solving application. It's a neat clever trick even if I do say so myself. First of all you have to crop your image so that the river cuts the image in two.

croppedImage = binaryImage(yTop:yBottom, :);

Now label the inverted cropped image. You will have two blobs: one on the left side of the image and one on the right side of the image.

labeledImage = bwlabel(~croppedImage);
leftHalf = ismember(labeledImage, 1);
rightHalf = ismember(labeledImage, 2);

Now you call bwboundaries() on each half, but the boundaries will go along the 3 edges of the image as well as down along the river. So you just delete points where y = 1 or the number of rows, and where x = 1 or the number of columns. Can you do that by yourself? Let's see how far you can get on your own. Post your code if you need more help.

[Jon]

Hey IA,

Thanks for the detailed solution! I was able to pull out the banks (in magenta):

They aren't perfect, but I can play with the image morphology processing orders and numbers--I looked at ~10 scenes and all of them were OK .

I am still left with a major-ish problem that I've highlighted with the yellow box: what to do when the river bifurcates and rejoins? If it's small enough, one bifurcation is removed with the erode command, but when the bifurcations are approximately equal width they both remain. And I want to retain both of them. Do you think there's an automatic way to handle this, or must it be manually processed?

Thanks for your help! Your trick is indeed clever; I certainly never would've thought of it (although I still don't really grasp image labeling).

[Image Analyst]

I'm watching the basketball finals tonight. You can do it with clever use of bwareaopen(). See if you can figure it out. Otherwise you'll have to wait until tomorrow for me to do it for you.

[Jon]

Haha, I appreciate your help and don't want you to do my work for me! I need to move on to computing migration rates from the extracted banklines at the moment, but I will return to this in the near future. If you have any more specific advice for how to go about it, I'd be very appreciative. I haven't thought about it much yet, though...