function [ v_num, v, v_p ] = polyomino_multihedral_variants ( mp, np, p_num, ...
  p, d )

%*****************************************************************************80
%
%% polyomino_multihedral_variants() finds variants of polyominoes in an array.
%
%  Discussion:
%
%    We are supplied with an array P of polyominoes.  We wish to "tile" a
%    region R, using, for each I, D(I) copies of polyomino P(I).
%    We are asked to return an array containing all the variants of these 
%    polynomials that can be generated by rotation and reflection.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    13 June 2018
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer MP, NP, the rows and columns in the representation of the polyomino P.
%
%    integer P_NUM, the number of polyominoes.
%
%    integer P(MP,NP,P_NUM), for each value K, P(:,:,K) is a matrix of 
%    0's and 1's representing the K-th polyomino.  Each matrix should be 
%    "top-left tight", that is, there should be a 1 in row 1, and a 1 in 
%    column 1.  However, there may be trailing rows and columns of zeros.
%
%    integer D(P_NUM), the number of copies of each polyomino that
%    may be used.  If D(I) <= 0, no variants of polyomino I will be generated.
%
%  Output:
%
%    integer V_NUM, the number of variants found.
%
%    integer V(MP,NP,V_NUM), the variants.
%
%    integer V_P(V_NUM); if V(:,:,I) is a variant of the polyomino
%    P(:,:,L), then V_P(I) is set to L.
%
  v_num = 0;
  v = [];
  v_p = [];

  for l = 1 : p_num
%
%  Only use polyominoes for which the user has set D to at least 1.
%
    if ( 1 <= d(l) )
%
%  Make a condensed version of the polyomino called "Q".
%
      [ mq, nq, q ] = polyomino_condense ( mp, np, p(1:mp,1:np,l) );
%
%  Generate all variants.
%
      for reflect = 0 : 1
        for rotate = 0 : 4

          [ ms, ns, s ] = polyomino_transform ( mq, nq, q, rotate, reflect );
%
%  Only proceed if this variant S will fit in [MP,NP].
%
          if ( ms <= mp & ns <= np )
%
%  Expand S to an [MP,NP] array T.
%
            t = zeros(mp,np);
            t(1:ms,1:ns) = s;
%
%  Compare T to all previous variants.
%
            different = true;
            for k = 1 : v_num
              if ( t(1:mp,1:np) == v(1:mp,1:np,k) )
                different = false;
                break;
              end
            end

            if ( different )
              v_num = v_num + 1;
              v(1:mp,1:np,v_num) = t(1:mp,1:np);
              v_p(v_num) = l;
            end

          end

        end
      end

    end

  end

  return
end
