Steve on Image Processing

October 9th, 2007

At the beginning of the upslope area series, I posted code showing how to compute the direction of maximum slope for a triangular facet, and how to compute the flow direction for a pixel. This purely scalar code ran fine, but more slowly than I wanted, so I decided to vectorize it before submitting the upslope area code onto MATLAB Central.

Here's the facet code:

s1 = (e0 - e1)/d1;             % eqn (1)
s2 = (e1 - e2)/d2;             % eqn (2)

r = atan2(s2, s1);             % eqn (3)
s = hypot(s1, s2);             % eqn (3)

diagonal_angle    = atan2(d2, d1);
diagonal_distance = hypot(d1, d2);

if r < 0
    % eqn (4)
    r = 0;
    s = s1;

elseif r > diagonal_angle
    % eqn (5)
    r = diagonal_angle;
    s = (e0 - e2)/diagonal_distance;
end

The code that needs to change is the part that deals with special cases (r < 0 and r > diagonal_angle). Instead of using logical tests on the scalar value of r, use logical indexing. Here's the new code fragment, from facet_flow.m in the upslope toolbox:

s1 = (e0 - e1) / d1;           % eqn (1)
s2 = (e1 - e2) / d2;           % eqn (2)

r = atan2(s2, s1);             % eqn (3)
s = hypot(s1, s2);             % eqn (3)

% Constrain flow direction to lie within or along the edges of the facet.
too_far_south = r < 0;
r(too_far_south) = 0;
s(too_far_south) = s1(too_far_south);

diagonal_angle    = atan2(d2, d1);
diagonal_distance = hypot(d1, d2);
too_far_north = r > diagonal_angle;
r(too_far_north) = diagonal_angle;
s(too_far_north) = (e0(too_far_north) - e2(too_far_north)) / ...
    diagonal_distance;

Now the code works when e0, e1, and e2 are matrices.

Here's the original pixel flow code:

% Table 1, page 311
e0 = E(i, j);
e1 = [E(i,j+1) E(i-1,j) E(i-1,j) E(i,j-1) ...
      E(i,j-1) E(i+1,j) E(i+1,j) E(i,j+1)];

e2 = [E(i-1,j+1) E(i-1,j+1) E(i-1,j-1) E(i-1,j-1) ...
      E(i+1,j-1) E(i+1,j-1) E(i+1,j+1) E(i+1,j+1)];

ac = [0  1  1  2  2  3  3  4];
af = [1 -1  1 -1  1 -1  1 -1];

rp = zeros(1, 8);
sp = zeros(1, 8);

for k = 1:8
    [rp(k), sp(k)] = facetFlow(e0, e1(k), e2(k), d1, d2);
end

% Find the facet with the maximum down-slope.
[s, k_max] = max(sp);

if s > 0
    % Maximum down-slope is positive.
    r = (af(k_max) * rp(k_max)) + (ac(k_max) * pi / 2);
else
    % The pixel is in a flat area or is a pit.  Flow direction angle is
    % unresolved.
    r = NaN;
end

This code is a little harder to vectorize. For each pixel, we have to look at each of eight triangular faces. I imagine one could write a version with no loops, but it would require lots of memory and be impossible to understand. So I decided to write a version that looped over the eight face directions. With only eight such directions, loop overhead is insignificant.

The way to compute a particular triangular face for every pixel simultaneously is to use linear indexing to access pixel neighbors. I described this technique, which I called neighbor indexing, in a post last March. The new function pixel_flow uses neighbor indexing to compute triangular facet vertices, E0, E1, and E2, for all pixels at the same time.

Start by computing a table of linear offsets for E1 and E2, for each of the facets:

% Compute linear indices at desired locations.
e0_idx = (j - 1)*M + i;

% Table 1, page 311
% Row and column offsets corresponding to e1 and e2 for each
% table entry:
e1_row_offsets = [0 -1 -1  0  0  1  1  0];
e1_col_offsets = [1  0  0 -1 -1  0  0  1];

e2_row_offsets = [-1 -1 -1 -1  1  1  1  1];
e2_col_offsets = [ 1  1 -1 -1 -1 -1  1  1];


% Linear e1 and e2 offsets.
e1_linear_offsets = e1_col_offsets*M + e1_row_offsets;
e2_linear_offsets = e2_col_offsets*M + e2_row_offsets;

Next, initialize E0, E1, and E2 for the first facet and call the vectorized facet_flow function for that facet:

E0 = E(e0_idx);
E1 = E(e0_idx + e1_linear_offsets(1));
E2 = E(e0_idx + e2_linear_offsets(1));

[R, S] = facet_flow(E0, E1, E2, d1, d2);

Finally, loop over the remaining seven facets. For each new facet, check to see if the downward slope is greater than the one previously recorded. If so, update the facet flow direction.

for k = 2:8
    % Compute Rk and Sk corresponding to the k-th facet. Where Sk is positive
    % and greater than S, replace S and recompute R based on Rk.
    E1 = E(e0_idx + e1_linear_offsets(k));
    E2 = E(e0_idx + e2_linear_offsets(k));
    [Rk, Sk] = facet_flow(E0, E1, E2, d1, d2);
    
    new_R = (Sk > S) & (Sk > 0);
    S = max(S, Sk);
    
    R(new_R) = (af(k) * Rk(new_R)) + (ac(k) * pi / 2);
end

I've skipped over some of the code details. See pixel_flow.m for the full story.

I have one more algorithm post planned in this series, about handling NaNs outside catchments, and then I'll post a final recap of the techniques used and lessons learned.

07:00 UTC | Posted in Upslope area | Permalink |

You can follow any responses to this entry through the RSS 2.0 feed. You can skip to the end and leave a response. Pinging is currently not allowed.

2 Responses to “Upslope area - vectorizing pixel flow calculations”

1. Ionut replied on October 12th, 2007 at 06:37 UTC :
   

   hello,

   I have a question which is not connected with your post, but is about image processing and Matlab :). How is possible to do radial/circular average of 2D FFT. I have an image, I did 2D FFT on it, and now I would like to have amplitude over r (r = radius), I could not find any function in Matlab which can do radial average, to obtain this kind of plot. Thank you

2. Steve replied on October 15th, 2007 at 10:40 UTC :
   

   Ionut—A Google search of “MATLAB radial average” pulls up this web page: http://www.yale.edu/eeb/prum/fourier.htm. You might take a look at it to see how they implemented the power spectrum radial average.

Leave a Reply

<pre class="code">
a = magic(3);
sum(a)
</pre>

These postings are the author's and don't necessarily represent the opinions of The MathWorks.