function [ i2, j2, i3, j3, s ] = polyomino_boundary_step ( i1, j1, i2, j2, s, p )

%*****************************************************************************80
%
%% polyomino_boundary_step() takes one more step along the boundary of a polyomino.
%
%  Discussion:
%
%    This function determines the next letter in the boundary word of a
%    polyomino, by taking one more step along the boundary.  The possible
%    directions are 'r', 'u', 'l' and 'd'.
%
%    P is an M by N array indexed by (I,J), containing 0's and 1's.
%
%    Vertices are M+1 by N+1 array indexed by (i,j).
%
%      (i-1,j)--(i-1,j+1)
%      |                |
%      |     P(I,J)     |
%      |                |
%      (i,  j)--(i,  j+1)
%
%    Search directions are:
%    1: turn 90 degrees clockwise from current direction
%    2: continue in current direction
%    3: turn 90 degrees counterclockwise from current direction.
%
%                  (i-1,j)
%                     ^
%         P(i,  j-1)  |    P(i,j)
%                     V
%     (i,j-1) <--> (i,  j) <--> (i,j+1)
%                     ^
%         P(i+1,j-1)  |    P(i+1,j)
%                     V
%                  (i+1,j)
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    21 June 2021
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer I1, J1: the coordinates of the starting node of the edge.
%
%    integer I2, J2: the coordinates of the end node of the edge.
%
%    character string S: the partial boundary word.
%
%    integer P(M,N):
%
%  Output:
%
%    integer I2, J2: the coordinates of the start node of the next edge.
%
%    integer I3, J3: the coordinates of the end node of the next edge.
%
%    character string S: the partial boundary word, updated with the next edge.
%
  [ m, n ] = size ( p );
%
%  R: We just went (i,j-1) to (i,j)
%
%   1|?
%  ->+-
%   0|?
%
%   1 option (i+1,j  )  D  if inside and P(i,  j  ) = 1
%   2 option (i,  j+1)  R  if inside and P(i-1,j  ) = 1
%   3 option (i-1,j  )  U  if inside, no cost, else error
%
  if ( i1 == i2 & j1 == j2 - 1 )

    if ( i2 + 1 <= m + 1 & j2 <= n & p(i2,j2) == 1 )
      i3 = i2 + 1;
      j3 = j2;
      s = [ s, 'd' ];
    elseif ( j2 <= n & p(i2-1,j2) == 1 )
      i3 = i2;
      j3 = j2 + 1;
      s = [ s, 'r' ];
    else
      i3 = i2 - 1;
      j3 = j2;
      s = [ s, 'u' ];
    end
%
%  D: we just went (i-1,j) to (i,j)
%
%     |
%    0V1
%    -+-
%    ?|?
%
%    1 option (i,  j-1)  L if inside and P(i,j-1) = 1
%    2 option (i+1,j  )  D if inside and P(i,j  ) = 1
%    3 option (i,  j+1)  R if inside no cost, else  error
%
  elseif ( i1 == i2 - 1 & j1 == j2 )

    if ( 1 <= j2 - 1 & p(i2,j2-1) == 1 )
      i3 = i2;
      j3 = j2 - 1;
      s = [ s, 'l' ];
    elseif ( i2 + 1 <= m + 1 & p(i2,j2) == 1 )
      i3 = i2 + 1;
      j3 = j2;
      s = [ s, 'd' ];
    else
      i3 = i2;
      j3 = j2 + 1;
      s = [ s, 'r' ];
    end
%
%  L: we just went (i,j+1) to (i,j)
%
%    ?|0
%    -+<--
%    ?|1
%
%    1 option (i-1,j  )  U if inside and P(i-1,j-1) = 1
%    2 option (i,  j-1)  L if inside and P(i,j-1) = 1
%    3 option (i+1,j  )  D if inside no cost, else  error
%
  elseif ( i1 == i2 & j1 == j2 + 1 )

    if ( 1 <= i2 -1 & 1 <= j2 - 1 & p(i2-1,j2-1) == 1 )
      i3 = i2 - 1;
      j3 = j2;
      s = [ s, 'u' ];
    elseif ( 1 <= j2 - 1 & p(i2,j2-1) == 1 )
      i3 = i2;
      j3 = j2 - 1;
      s = [ s, 'l' ];
    else
      i3 = i2 + 1;
      j3 = j2;
      s = [ s, 'd' ];
    end
%
%  U: we just went (i+1,j) to (i,j)  (U)
%
%    ?|?
%    -+-
%    1^0
%     |
%
%    1 option (i,  j+1)  R if inside and P(i-1,j) = 1
%    2 option (i-1,j  )  U if inside and P(i-1,j-1) = 1
%    3 option (i,  j-1)  L if inside no cost, else  error
%
  elseif ( i1 == i2 + 1 & j1 == j2 )

    if ( 1 <= i2 - 1 & j2 + 1 <= n + 1 & p(i2-1,j2) == 1 )
      i3 = i2;
      j3 = j2 + 1;
      s = [ s, 'r' ];
    elseif ( 1 <= i2 - 1 & 1 <= j2 - 1 & p(i2-1,j2-1) == 1 )
      i3 = i2 - 1;
      j3 = j2;
      s = [ s, 'u' ];
    else
      i3 = i2;
      j3 = j2 - 1;
      s = [ s, 'l' ];
    end

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'polyomino_boundary_step(): Error!\n' );
    fprintf ( 1, '  (i1,j1) and (i2,j2) do not define a step.\n' );
    error ( 'polyomino_boundary_step(): Error!' );

  end

  return
end

