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.