function variants = colour_C_variants ( polyomino, C )

%*****************************************************************************80
%
%% colour_C_variants() creates and plots matrices for a coloured region.
%
%  Discussion:
%
%    Given the binary matrix for a polyomino or region, output the
%    matrices PARI1, PARI2, ..., PARIC  corresponding to the C , (2<=C<=6), 
%    colourings of this polyomino or region, i.e., chrominoes.  The elements
%    of these matrices are chosen from the set {0, 1, 2, ..., C} using the modular
%    function, yielding a generalized checkerboard colouring, where 0 elements
%    are used to denote a 'hole' in the region. 
%
%    The C coloured regions are plotted, where each element is associated 
%    with a fixed colour according to the following RGB triplets:
%
%    Element      RGB triplet        colour
%       0           [1 1 1]          white
%       1           [0 0 0]          black
%       2           [0.7 0.7 0.7]    gray
%       3           [1 0 0]          red
%       4           [0 1 0]          green
%       5           [1 1 0]          yellow
%       6           [0 0 1]          blue
%    
%  Example:
%
%    polyomino =   [   0     1     0     0     0
%                      1     1     1     1     0
%                      0     0     1     1     1
%                      0     0     0     1     0  ]
%
%    C = 6
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    13 February 2023
%
%  Author:
%
%    Marcus Garvie
%
%  Input:
%
%    integer polyomino(rows,cols): a matrix of 0's and 1's which defines
%    the shape of a polyomino.
%
%    integer C: the number of colors to be used.
%
%  Output:
%
%    integer variants(rows,cols,C):
%

%
%  Construct matrices PARI1, PARI2, ..., PARIC and pack them into VARIANTS
%
  [ rows, cols ] = size ( polyomino );

  for page = 1 : C
    variants(:,:,page) = zeros ( rows, cols );
  end
  
  for j = 1 : cols
    for i = 1 : rows
      if ( polyomino(i,j) == 1 )
        for page = 1 : C
          variants(i,j,page) = mod ( i+j+page-1, C ) + 1;
        end
      end
    end
  end
%
%  Output some colour counts.
%
  for page = 1 : C
    fprintf ( '\n' );
    fprintf ( '  Variant   %d \n', page );
    for colour = 1 : C
      num = length ( find ( variants(:,:,page) == colour ) );
      fprintf ( '  Number of squares of colour %d = %d \n', colour, num );
    end
  end
%
%  Store all possible RGB triplets
%  
  index(1,3,:) = zeros ( 1, 3 );

  index(1,1:3,1) = [ 1.0, 1.0, 1.0 ];
  index(1,1:3,2) = [ 0.0, 0.0, 0.0 ];
  index(1,1:3,3) = [ 0.7, 0.7, 0.7 ];
  index(1,1:3,4) = [ 1.0, 0.0, 0.0 ];
  index(1,1:3,5) = [ 0.0, 1.0, 0.0 ];
  index(1,1:3,6) = [ 1.0, 1.0, 0.0 ];
  index(1,1:3,7) = [ 0.0, 0.0, 1.0 ];
%
%  Compute the colormap for each of the C variants and plot them.
%
  clf ( );

  for page = 1 : C

    array = variants(:,:,page);
    unq = unique ( array );
    ind = unq + 1;
    len = length ( ind );
    cmap = [];
    for i = 1 : len
      color = ind(i);
      cmap = [ cmap; index(:,:,color) ];
    end
    [ ~, ~, unq ] = unique ( array );
    piece = [ 1 : len ]; 
    arrayscaled = piece(unq);
    arrayscaled = reshape ( arrayscaled, size(array) );  
    subplot ( 1, C, page );
    hAxes = gca;
    imagesc ( hAxes, arrayscaled );
    colormap ( hAxes, cmap )
    set ( gcf, 'color', 'w' )
    title ( sprintf ( 'Variant   %d', page ) );
    axis ( 'equal', 'tight', 'off' )

  end

  return 
end
