function C = components_3d ( A )

%*****************************************************************************80
%
%% components_3d() assigns contiguous nonzero pixels to a common component.
%
%  Discussion:
%
%    On input, the A array contains values of 0 or 1.
%
%    The 0 pixels are to be ignored.  The 1 pixels are to be grouped
%    into connected components.
%
%    The pixel A(I,J,K) is "connected" to the pixels:
%
%      A(I-1,J,  K  ),  A(I+1,J,  K  ),
%      A(I,  J-1,K  ),  A(I,  J+1,K  ),
%      A(I,  J,  K-1),  A(I,  J,  K+1),
%
%    so most pixels have 6 neighbors.
%
%    On output, COMPONENT_NUM reports the number of components of nonzero
%    data, and the array C contains the component assignment for
%    each nonzero pixel, and is 0 for zero pixels.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    12 March 2025
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer L, M, N, the order of the array.
%
%    integer A(L,M,N), the pixel array.
%
%  Output:
%
%    integer C(L,M,N), the component array.
%
  [ l, m, n ] = size ( A );
%
%  Initialization.
%
  C = zeros ( l, m, n );
  component_index = 0;

  plist = [];

  for i2 = 1 : l
    for j2 = 1 : m
      for k2 = 1: n

        if ( A(i2,j2,k2) ~= 0 && C(i2,j2,k2) == 0 )

          plist = [ plist; i2, j2, k2 ];
          component_index = component_index + 1;

          while ( ~ isempty ( plist ) )

            i = plist(end,1);
            j = plist(end,2);
            k = plist(end,3);
            plist(end,:) = [];
            C(i,j,k) = component_index;

            if ( 1 <= i-1 && A(i-1,j,k) ~= 0 && C(i-1,j,k) == 0 )
              plist = [ plist; i - 1, j, k ];
            end
            if ( i+1 <= l  && A(i+1,j,k) ~= 0 && C(i+1,j,k) == 0 )
              plist = [ plist; i + 1, j, k ];
            end
            if ( 1 <= j - 1 && A(i,j-1,k) ~= 0 && C(i,j-1,k) == 0 )
              plist = [ plist; i, j - 1, k ];
            end
            if ( j+1 <= m && A(i,j+1,k) ~= 0 && C(i,j+1,k) == 0 )
              plist = [ plist; i, j + 1, k ];
            end
            if ( 1 <= k - 1 && A(i,j,k-1) ~= 0 && C(i,j,k-1) == 0 )
              plist = [ plist; i, j, k - 1 ];
            end
            if ( k+1 <= n && A(i,j,k+1) ~= 0 && C(i,j,k+1) == 0 )
              plist = [ plist; i, j, k + 1 ];
            end
          end
        end
      end
    end
  end

  return
end
