%
% L_system_block_matrix( nGenerations, mtxLS0, varargin )
%
%   Creates a 2-D block replacement (D0L) L_system matrix, with symbols 
%   represented by integers.
%   
%	Written as a utility for L_system_tiling.m. This function does only
%	basic sanity checks on dimensions, bounds, or number of arguments. 
%   See L_system tiling for other examples of bounds checking.
%
%   USAGE:  % Get 8th generation Thue-Morse pattern:
%           mtxOut = L_system_block_matrix( 8, 0, [0 1; 1 0], [1 0; 0 1] );
%           % A 2-rule system with rotations in the first rule.
%
%           % The 30 codes for a 3*pi/2 rotation of rule 0, and
%           % the 10 codes for a 3*pi/2 rotation of rule 0:
%           mtxOut = L_system_block_matrix( 8, 0, [0 30; 10 1], [1 1; 1 1] );
%
%   ARGUMENTS:
%
%       nGenerations:	Number of generations.
%
%       mtxLS0:         L-system initial 2-D matrix, which follows the same
%                       coding conventions as the replacement rules.
%
%       varargin:       A set of L-system block replacement rules. 
%                       These are coded as a list of nxm integer matrices.
%                       The maximum number of rules supported is 10, [0,9]; 
%                       replace_block matrix.m does not have this restriction. 
%
%   RETURN VALUES:
%
%       mtxOut:         Matrix of integers of the system pattern.
%
%   HARDCODED:  (see below 'function' statement)
%
%       bShow:  Show details for debugging, and show output matrix in a
%               figure window.
%
%
% Mark Dow,     August 23, 2009
%


function mtxOut = L_system_block_matrix( nGenerations, mtxLS0, varargin )

%%%%%%%%%%%%%%%%%%%%%%%%
% Hardcoded information:

bShow = false;
%
%%%%%%%%%%%%%%%%%%%%%%%%


% Sanity check: there must be at least one rule, or one cell array of rules:
if nargin - 2 < 1
	fprintf( [ '\nL-system_block_matrix.m: Must have at least three arguments.\n\n' ] );
    return
end

arrayRules      = [];    % rules
arrayRotes      = [];    % rotations
arrayMirrs      = [];    % mirrorings


if nargin - 2 == 1 && iscell( varargin{1} )
    
    nRules = size( varargin{1},2 );
    
    szRules = size( varargin{1}{1} );
    arrayRules( 1:szRules(1), 1:szRules(2), 1:nRules ) = 0;
    for iR = 1 : nRules
        
        arrayRules( 1:szRules(1), 1:szRules(2), iR ) = varargin{1}{iR};
    end
else
    
    nRules = nargin - 2;
    szRules = size( varargin{1} );
    
    % Load rules into n x m x nRules array from arguments:
    for iR = 1 : nargin - 2
        
        arrayRules( 1:szRules(1), 1:szRules(2), iR ) = varargin{iR};
    end
end

% Parse rule symbols, rotation codes and mirror codes:

arrayMirrs = zeros( size(arrayRules) );
arrayRotes = zeros( size(arrayRules) );

arrayMirrs( arrayRules >= 100 ) = 1;
% Strip hundreds digit.
arrayRules( arrayRules >= 100 ) =  mod( arrayRules( arrayRules >= 100 ), 100 );

arrayRotes( arrayRules >= 10 ) = floor( arrayRules( arrayRules >= 10 )/10 );
% Strip tens digit.
arrayRules( arrayRules >= 10 ) =   mod( arrayRules( arrayRules >= 10 ), 10 );


% Parse initial array symbols, rotation codes and mirror codes.

mtxMR0 = zeros( size(mtxLS0) );
mtxRT0 = zeros( size(mtxLS0) );

mtxMR0( arrayRules >= 100 ) = 1;
% Strip hundreds digit.
mtxLS0(     arrayRules >= 100 ) =   mod( arrayRules( arrayRules >= 100 ), 100 );

mtxRT0( arrayRules >= 10  ) = floor( arrayRules( arrayRules >= 10 )/10 );
% Strip tens digit.
mtxLS0(      arrayRules >= 10  ) =   mod( arrayRules( arrayRules >= 10 ), 10 );


% To Do: If non-square system and rotations (and mirrors?) are non-zero, it
% will fail. Insert reality check and abort with message.


% Find the set of (integer) symbols given in rules:
listSymbolsUsed( 1:nRules ) = 0;
for j = 1 : szRules(1)
for k = 1 : szRules(2)
    
    for iR = 1 : szRules
        
        if arrayRules( j, k, iR ) > nRules
            fprintf( [ '\nL-system_block_matrix.m: Rule contains last digit integer symbol that is larger than the number of rules.\n\n' ] );
            return
        end
        
        listSymbolsUsed( arrayRules( j, k, iR ) + 1 ) = 1;
    end
end
end
% Add to list the set of (integer) symbols used in initial symbols:
szInit = size( mtxLS0 );
for j = 1 : szInit(1)
for k = 1 : szInit(2)
    
    if mtxLS0( j, k ) > nRules
        fprintf( [ 'L-system_block_matrix.m: Initial array contains last digit integer symbol that is larger than the number of rules.' ] );
        return
    end

    listSymbolsUsed( mtxLS0( j, k ) + 1 ) = 1;
end
end

        
        
% % Sanity check that all rules can be visited:
% szRules = size(arrayRules);
% if nSymbols ~= szRules(3)
% 
%     fprintf( [ '\n\nThe indicated number of rules and symbols, ' num2str( nSymbols ) ...
%                ',\ndoes not match the number of rules given, '  num2str( szRules(3) ) '.\n\n' ] );
% 
%     return
% end
    

% % Extract rule codes:
% if bRandom == 0
%     for i = 1 : nSymbols
% 
%     % To Do: Document this.
%         nCodes( i, 1:2 ) = 0;
%         iDigit = 0;
%         for j = szRules(1) : -1 : 1
%         for k = szRules(2) : -1 : 1
% 
%             nCodes( i, 1 ) = nCodes( i, 1 ) + arrayRules( j, k, i )*nSymbols^iDigit;
%             nCodes( i, 2 ) = nCodes( i, 2 ) + ( arrayRotes( j, k, i ) + 4*arrayMirrs( j, k, i ) )*8^iDigit;
%             iDigit = iDigit + 1;
%         end
%         end
%     end
% end

% mtxLS0
% mtxRT0
% mtxMR0
% 
% arrayRules
% arrayRotes
% arrayMirrs

% Initialize matrices.
    
szLS = size(mtxLS0);
mtxOut( 1:szInit(1)*(szRules(1)^ nGenerations    ), 1:szInit(2)*(szRules(2)^ nGenerations    ) ) = 0;
mtxRT(  1:szInit(1)*(szRules(1)^(nGenerations-1) ), 1:szInit(2)*(szRules(2)^(nGenerations-1) ) ) = 0;
mtxMR(  1:szInit(1)*(szRules(1)^(nGenerations-1) ), 1:szInit(2)*(szRules(2)^(nGenerations-1) ) ) = 0;

mtxOut( 1:szInit(1), 1:szInit(2) ) = mtxLS0;
mtxRT(  1:szInit(1), 1:szInit(2) ) = mtxRT0;
mtxMR(  1:szInit(1), 1:szInit(2) ) = mtxMR0;

mtxOuti = mtxLS0;
mtxRTi  = mtxRT0;
mtxMRi  = mtxMR0;
    
for iG = 1 : nGenerations
    
    % Get the next generation of the pattern.
    % Note: The replacement in this block is a duplication of
    % L_system_block_matrix.m, avoiding some function overhead when
    % passing large matrices.

    for j = 1 : szLS(1)
    for k = 1 : szLS(2)

        y1 = (j-1)*szRules(1) + 1 : j*szRules(1);
        x1 = (k-1)*szRules(2) + 1 : k*szRules(2);

        % Insert rotated replacement pattern.
        mtxOut( y1, x1 ) = rot90( arrayRules( 1:szRules(1), 1:szRules(2), mtxOuti(j,k) + 1 ), ...
                                 -mtxRTi(j,k) );
        % Mirror replaced pattern.
        if mod( mtxMRi(j,k), 2 ) == 1
            mtxOut( y1, x1 ) = fliplr( mtxOut( y1, x1 ) );
        end

        if iG < nGenerations

            % Copy rotate and mirror rules.
            mtxRT( y1, x1 ) =   mtxRTi(j,k);
            mtxMR( y1, x1 ) =   mtxMRi(j,k);
            
            % Rotate rotate and mirror rules.
            mtxRT( y1, x1 ) =   mtxRT( y1, x1 ) ...
                              + rot90( arrayRotes( 1:szRules(1), 1:szRules(2), mtxOuti(j,k) + 1 ), ...
                                       -mtxRTi(j,k) );
            % Mirror rotate and mirror rules.
            if mod( mtxMRi(j,k), 2 ) == 1
                mtxMR( y1, x1 ) =   mtxMR( y1, x1 ) ...
                                  + fliplr( arrayMirrs( 1:szRules(1), 1:szRules(2), mtxOuti(j,k) + 1 ) );
            else
                mtxMR( y1, x1 ) =   mtxMR( y1, x1 ) ...
                                  + arrayMirrs( 1:szRules(1), 1:szRules(2), mtxOuti(j,k) + 1 );
            end
        end
    end
    end
    
    mtxOuti = mtxOut(  1:szInit(1)*(szRules(1)^iG), 1:szInit(2)*(szRules(2)^iG) );
    if iG < nGenerations
        mtxRTi  = mtxRT(  1:szInit(1)*(szRules(1)^iG), 1:szInit(2)*(szRules(2)^iG) );
        mtxMRi  = mtxMR(  1:szInit(1)*(szRules(1)^iG), 1:szInit(2)*(szRules(2)^iG) );
    end
    
    mtxRT = mod( mtxRT, 4 );
    mtxMR = mod( mtxMR, 2 );

    szLS = [ j*szRules(1) k*szRules(2) ];
end

if bShow
    
    figure; imshow( mtxOut )
end