Using pixels to find path of descent for rainwater

So I've been given a variable called "map" which is a matrix of heights in meters. It's a 100*100 matrix. I've used

imagesc(map);axis equal;colormap gray;

to display the map as a shaded image, then I used

[r,c]=size(map) to assign r and c to the number of rows and columns.

My goal is to follow the path of steepest descent for rainwater by finding the direction to the lowest height of the 9 pixels surrounding a given pixel in the matrix of map, including the pixel itself. I need to use the row and column offsets to do this. from my limited knowledge about pixels, I'm gathering that offsets can have values of -1, 0, or 1. Also, The boundary pixels will simply be defined as 0,0 for simplicity's sake.

How do I find these offsets? Also, how can I store row and column offsets for each map point r,c to the smallest neighbor? I need to be able to write a function that returns two matrices of the same size as the original matrix so that pixel map(r,c) has pixel map(r+roffset(r,c),c+coffset(r,c))as the minimum of those 9 heights.

The first of the matrices is supposed to find, for every cell in the original matrix, the row offset of the lowest neighbor to that cell, while the second is supposed to find the column offsets of the lowest neighbor for the same cell.

The example I've been given is that if a given pixel x,y has a row offset of -1 and column offset of 1, then the lowest point neighboring that pixel would be at row x-1 and column y+1.

Answer 1

[Image Analyst]

Well, there's the intuitive, straightforward, brute force approach with 4 nested for loops. Did you try that? It goes just like you would think to do if you explained it to me. Scan every row and every column, and for each pixel at that particular row and column, scan a 3 b y 3 window around it looking for the min value and recording its direction.

for col = 1 : columns
  for row = 1 : rows
    minHere = inf;
    for r = -1 : 1
      for c = -1 : 1
        if grayImage(row+r, col+c) < minHere
          direction(row, col) = .......

and so on

[Johny]

Is the "columns" and "rows" you have at the first two lines supposed to be where I substitute the number of columns and rows in the original matrix respectively? also, r and c are supposed to be the row and column values of our original matrix, they do not represent the offsets themselves, which is what I think you have done by that 4th and 5th line of code.

I don't know how that helps, what is the end result supposed to look like? and how does this loop close?

I'm sorry about all the questions.

[Johny]

I'd like to just illustrate my question with a quick example. let's say my original matrix was

1 3 8 9 5 6

2 4 7 3 1 9

3 8 9 2 7 3

7 9 3 7 1 4

For the numbers along the border, the offset for both rows and columns will be 0. Also if the center pixel is the lowest value its offsets will also be 0.

Considering point (2,2) in this matrix, which is 4. Its lowest neighbor, which is the 1 value in the top left corner, has an offset of 1,1 (let's quickly ignore the fact that it is on the edge of the matrix for convenience's sake. normally it would have an offset of 0,0) because it comes 1 row before as well as 1 column before the central value, which is our 4. I need to repeat this for every point on the map that comes in as input to the function and store all the collected row and column offsets in 2 separate matrices. I hope that helps my explanation.

[Image Analyst]

[rows, columns] = size(grayImage);

I don't know how you want to indicate direction. Maybe you want upper left to be 1 and increase clockwise until the left pixel is 8. So the directions might be like this:

1 2 3
8   4
7 6 5

Whatever, it doesn't really matter - just pick some method. Then convert r and c into one of the 8 directions 1-8. So an r of +1 and a c of +1 would be 5, for example.

[Johny]

how do I make your method part of the loop? I would need to have a way of taking those offset data and placing them into two different matrices, one for row offsets and the other for column offsets.

[Image Analyst]

Again, let's just take the intuitive, straightforward brute force approach. How would you think you might do it? Do you think the following might work?

if r == -1 && c == -1
  theDirection = 1;
elseif r == -1 && c == 0
  theDirection = 2;
elseif r == -1 && c == 1
  theDirection = 3;
elseif r == 0 && c == 1
  theDirection = 4;
elseif r == 1 && c == 1
  theDirection = 5;
elseif r == 1 && c == 0
  theDirection = 6;
elseif r == 1 && c == -1
  theDirection = 7;
elseif r == 0 && c == -1
  theDirection = 8;
end
direction(row, col) = theDirection;

By the way, you might want look at bwdistgeodesic(): http://blogs.mathworks.com/steve/2011/11/01/exploring-shortest-paths-part-1/ or the watershed functions.

Answer 2

[Image Analyst]

You can also solve this with Dynamic Programming http://en.wikipedia.org/wiki/Dynamic_programming#Checkerboard if your starting and ending points are known, or Djikstra's method if you have a starting point and just want to do a greedy depth first approach. Or you could use A* if you want to do an optimized, non-greedy approach. All of these are explained in Wikipedia.