function number = variomino_embed_number ( r, p )

%*****************************************************************************80
%
%% variomino_embed_number() counts variomino embeddings in a region.
%
%  Discusion:
%
%    A polyomino is a shape formed by connecting unit squares edgewise
%
%    A variomino is a polyomino in which each square has been assigned
%    a value between 1 and v.
%
%    A variomino can be represented by an MxN matrix, whose entries are
%    between 1 and v for squares that are part of the variomino, 
%    and 0 otherwise.
%
%    A region R is a variomino represened by an MRxNR matrix.
%
%    The variomino P is represented by an MPxNP matrix.
%
%    An "embedding" of P into R is an offset (MI,NJ) such that,
%    for all nonzero P(I,J), we have:
%      P(I,J) = R(I+MI,J+NJ) 
%      for 1 <= I <= MP, 1 <= J <= NP, and 
%      for 0 <= MI <= MR-MP, 0 <= MJ <= NR-NP.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    16 April 2023
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer r(mr,nr): the matrix representation of the R variomino.
%
%    integer p(mp,np): the matrix representation of the P variomino.
%
%  Output:
%
%    integer number: the number of distinct embeddings of P into R.
%
  number = 0;

  [ mr, nr ] = size ( r );
  [ mp, np ] = size ( p );
%
%  Get the index vector for nonzeros in p.
%
  pvec = reshape ( p, mp * np, 1 );
  p_nonzero_index = find ( pvec ~= 0 );
%
%  For each possible (I,J) coordinate of the upper left corner of a subset of R,
%  see if it matches P.
%
  for mi = 0 : mr - mp
    for nj = 0 : nr - np
      rvec = reshape ( r(1+mi:mp+mi,1+nj:np+nj), mp * np, 1 );
      if ( all ( rvec(p_nonzero_index) == pvec(p_nonzero_index) ) )
        number = number + 1;
      end    
    end
  end

  return
end

