%
% L_system_tiling( strOutFilestemOpt, nGenerations, nSymbols, nDecimate,  
%                     nOffset, strMotifFilename, arrayInit, varargin )
% -------------------------------------------------------------------------
%
%   Creates a 2-D block replacement (D0L) L-system image, with symbols 
%   represented by grayscale pixels or replaced by a set of motifs.
%
% USAGE:    imOut = L_system_tiling( 'Thue-Morse',  8, 2, 1, 0, 'motif_image.png', 0, [0 1; 1 0], [1 0; 0 1] );
%        OR imOut = L_system_tiling( 'Thue-Morse',  8, 2, 1, 0, '',                0, [0 1; 1 0], [1 0; 0 1] );
%
%   To genetate a random tesselation, with a fixed seed:
%
%        imOut = L_system_tiling( 'output_name', 8, 2, 1, 0, '', 0, 0 ); 
%
%   To perform a simple replacement of a 2-D array of symbols with a corresponding
%   set of motifs:
%
%        imOut = L_system_tiling( 'output_name', 8, 2, 1, 0, 'motif_image.png', 'BW.png', 0, 1 );
%
% To Do: Suppress the diagram generation appropriately with symbol replacement. 
% Does it support more than 2 symbols?
%
% INPUT:
%
%   strOutFilestemOpt - Optional filestem for resulting PNG image. The
%                       default filestem is 'LS'. All stems will have
%                       L-system parameters appended.
%
%   nGenerations      - [positive integer] Number of recursive iterations 
%                       for which the rules are applied to the initial 
%                       symbol(s). The tesselation dimension will be
%                       (dimension of arrayStart)*(dimension of rule )^nGenerations.
%
%   nSymbols          - Number of symbols (alphabet) in system. This must 
%                       match thenumber of unique symbols in the rule 
%                       arrays (see varargin).
%
%   nDecimate         - [1,n] Decimate system by only including every nth 
%                       symbol,where n is nDecimate.
%
%   nOffset           - Offset to the first item included in the decimation.
%
%   strMotifFilename  - Filename of an image containing tiling motifs.
%
%   arrayStart        - The initial state of the system (the axiom or
%                       initiator), an array composed of symbols. 
%
%   varargin          - Block replacement (production) rules. The number 
%                       of rules mustmatch the number of symbols
%                       (nSymbols), and all be of the same dimension. The 
%                       (n-1)th varargin is the predecessor of the symbol 
%                       (n-1) -- the first is the rule for the symbol '0'.
%
%
% HARDCODED (immediately below function declaration):
%
%   bShowResults:      - [0, 1] Show progress, display final image and 
%                              algorithm graphic in figure.
%
%   bWriteResults:    - [0, 1] Write final image and algorithm graphic to  
%                              an image file.
%
%   bCreateHistogram: - [0, 1] Write a histogram of rule values numbered
%                              from 1 to nSymbols.
%                      
%   nByGenOversample  - [positive integer] Oversampling ratio for 
%                              generation frames graphic.
%
%   nByGenWidthMax    - [positive integer] Maximum width for 
%                              a generation graphics frame.
%
%   bCreateGAFrames   - [positive integer]Not currently used.
%
%   strOutFileName    - Convention for default (if strOutFilename == '') 
%                       output file name and format.
%                       Default naming is:
%                       TM[nNbrRule]_[nGenerations]_[strMotifFilename].png
%
% OUTPUT:
%
%   Returns an RGB image array. Depending on hardcoded options,the program 
%   will display and/or write a variety of auxilliary graphics. For the
%   simplest systems (no motif as image,a single element initial initial 
%   condition, 2 x 2 rules) the dimensions of the final image will be 
%   2^nGenerations x 2^nGenerations.
%
%   If hardcoded bWriteResults == 1 creates an 
%   image file called [strOutFilestemOpt]-[size of rule arrays]-[strMotifFilename]-
%                     [rule rotate and mirror codes]_[rule codes]-[nGenerations 'g'].png, 
%                     with a default name of
%                     LS-[size of rule arrays]-[rule rotate and mirror codes]_
%                     [rule codes]-[nGenerations]g.png 
%   and if bShowResults == 1 displays progress and a Matlab figure of the 
%   image.
%
%
%   FUNCTION CALLS and DEPENDENCIES: 
%       
%       -> make_image_pyramid.m
%           -> interp2_fa.m
%               -> interp1_fa.m
%       -> get_number_pairs.m
%       -> add_alternating_pairs.m
%       -> get_rule_symmetries.m
%       -> interp2_n.m
%       -> interp2_fa.m
%           -> interp1_fa.m
%       -> dft_of_image.m
%           -> interp2_n.m
%               -> hl.m
%                   -> hl_color_cycle.m
%
%   Requires these image resources in same directory:
%
%       -> LS_rule_graphic_1.png
%       -> LS_rule_graphic_2.png
%       -> LS_rule_graphic_3.png
%       -> LS_rule_graphic_4.png
%       -> Scale_symmetry_L_8.png
%       -> Scale_symmetry_U_8.png
%       -> Scale_symmetry_C_8.png
%       -> Scale_symmetry_F_8.png
%       -> Scale_symmetry_O_8.png
%       -> Scale_symmetry_D_8.png
%       -> Scale_symmetry_D-O_8.png
%
% Mark Dow,             January 1, 2007
% Mark Dow, modified       June 5, 2007     (derived from Thue_Morse_tiling.m) 
% Mark Dow, modified    August 23, 2007     (motif usage modifications) 
% Mark Dow, modified September 15, 2007     (handle images as rules) 
% Mark Dow, modified September 22, 2007     (output rule set graphic) 
% Mark Dow, modified   October 16, 2007     (by generation frames graphic) 
% Mark Dow, modified   December 9, 2007     (coding of parameters in output filestem) 
% Mark Dow, modified   January 24, 2008     (bug fix for non-square motif files) 
% Mark Dow, modified  February 12, 2008     (output naming convention, fixed random option) 
% Mark Dow, modified      March 6, 2008     (input pattern and motifs displayed in the algorithm graphic) 
% Mark Dow, modified       July 3, 2008     (adding rule rotation decoding and implementation) 
% Mark Dow, modified       July 3, 2008     (rule rotation supported,along with altered naming conventions) 
% Mark Dow, modified      July 25, 2008     (added rule mirroring support) 
% Mark Dow, modified      July 31, 2008     (started adding symmetry graphic construction) 
% Mark Dow, modified    August 20, 2008     (adding symmetry graphic construction) 
% Mark Dow, modified    January 8, 2009     (random option bug fixes) 
% Mark Dow, modified      March 8, 2009     (removed symmetry graphic for > 2x2 rules to prevent errors) 
% Mark Dow, modified     March 10, 2009     (allows non-square initial symbol array [arrayStart]) 
% Mark Dow, modified     March 17, 2009     (symbol,rule reality check and rules from image bug fix)
% Mark Dow, modified     March 18, 2009     (adding average of generations graphic)
% Mark Dow, modified     March 28, 2009     (adding DFT graphics)
% Mark Dow, modified      April 7, 2009     (Compatibility with dft_of_image.m)
% Mark Dow, modified        May 17, 2009    (fixing symmetry diagrams, DFT phase drifts, details)
%
%   TERMS FOR USE: There are no restrictions on the use of this code, 
%                  auxilliary code and other required resources. Claiming 
%                  to be the originator, explicitly or implicitly, is bad  
%                  karma. A link (if appropriate), a note to  
%                  dow[at]uoregon.edu, and credit are appreciated but not 
%                  required.
%

% To Do: If not writing or displaying the results( bShowResults and bWriteResults = 0)
% don't generate algorithm, by generation, and symmetry graphics.



function imOut = L_system_tiling( strOutFilestemOpt, nGenerations, nSymbols, nDecimate, nOffset, strMotifFilename, arrayInit, varargin )

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% hardcoded parameters

bShowResults        = 1;    % [true,false]
bWriteResults       = 0;    % [true,false]

bWriteTilingOnly    = 0;    % [true,false]
bDrawAverage        = 1;    % [true,false]
bDrawRulesAsImage   = 0;    % [true,false]
bDrawSmallAlgorithm = 1;    % [true,false]
bShowDFTtoo         = 1;    % [true,false]
bWriteDFTtoo        = 1;    % [true,false]

nOversampleOut      = 1;    % integer oversample some pattern graphics

    flPhaseDriftX           = 1*pi;
    flPhaseDriftY           = 1*pi;
    flPhaseDriftOffset      = 0*pi;
    flAvgPhaseDriftX        = 1*pi;
    flAvgPhaseDriftY        = 1*pi;
    flAvgPhaseDriftOffset   = 0*pi;

bCreateHistogram    = 0;    % [true,false]
nByGenOversample    = 3;
nByGenWidthMax      = 300;
bCreateGAFrames     = 0;    % [true,false] % not used

% hardcoded parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

bRandom             = 0;

arrayRules = [];
arrayRotes = [];
arrayMirrs = [];
szRules    =  0;
imMotifs   = -1;
szMotifs   = -1;
imRules    = -1;

% Sanity check for number of symbols (nSymbols) matching the number of 
% rules given.
if nargin - 7 ~= nSymbols && ~isa( varargin{1}, 'char' )
        fprintf( [ '\nError: The number of symbols you indicated, ' num2str(nSymbols) ', \n' ] );
        fprintf( [ '       does not match with the number of rules given (the last ' num2str( nargin - 7 ) ' parameters). \n' ] );
        imOut = -1;
        return;
end

% If random (2 symbol):
% To Do: 2 symbol reality check.
if nargin == 8 && ~isa( varargin{1}, 'char' ) && ( varargin{1} == 0 )
    
    bRandom = 1;
    
    % 2x2 hardcoded
    szRules( 1:2 ) = 2;
    arrayRules( 1:szRules(1), 1:szRules(2), 1:2 ) = 0;
    
    for i = 1 : nSymbols 
        SymbolColor( i, 1:3 ) = i;
    end
end

% Get an image if the initial condition is given as an image filename.
strInit = '';
if isa( arrayInit, 'char' )
    
    strInit = arrayInit;
    imInit = imread( arrayInit );
    imInit = floor( ((nSymbols-1)/255)*imInit);
    arrayInit = imInit(:,:,1);

end

% % If image array replacement:
% % To Do: Image array replacement reality check.
% if nargin == 8 & ~isa( varargin{1}, 'char' )
%     
%     % 2x2 hardcoded
%     image = varargin{1};
%     szArray = size( varargin{1} );
%     szRules( 1:2 ) = szArray(1:2);
%     arrayRules( 1:szRules(1), 1:szRules(2), 1 ) = image( 1:szRules(1), 1:szRules(2), 1 );
%     
%     for i = 1 : nSymbols 
%         SymbolColor( i, 1:3 ) = i;
%     end
% end

% Assemble an array of rules from function arguements.
if nargin == 8 && bRandom == 0 && isa( varargin{1}, 'char' )
        
    % Extract rules from the eighth arguement's image file:
    strRulesFilename = varargin{1};
    
    % If the rules image is color, assign a numeric rule to each 
    % unique color:
    %imRules = double( imread( strRulesFilename ) );

    % Load the rules image into a 24 bit color array:
    infRules = imfinfo( strRulesFilename );
    if strcmp( infRules.ColorType, 'truecolor' )
        if infRules.Height == 1
            imRules(1,:,:) = imread( strRulesFilename );
        else
            imRules = imread( strRulesFilename );
        end       
    else if strcmp( infRules.ColorType, 'grayscale' )
        imGrayRules = imread( strRulesFilename );
        szGrayRules  = size( imGrayRules );
        imRules( 1:szGrayRules(1), 1:szGrayRules(2), 1:3 ) = 0;
        imRules(:,:,1) = imGrayRules;
        imRules(:,:,2) = imGrayRules;
        imRules(:,:,3) = imGrayRules;
        
        else if strcmp( infRules.ColorType, 'indexed' )
                fprintf( [ 'Indexed image format for the rules is not currently supported.' ] );
                return;
            end
        end
    end

    szRules = size( imRules );
    
    % Find the unique colors of the rule image and initial conditions.
    SymbolColor( 1:nSymbols, 1:3 ) = -1;
    imRulesCopy = imRules;
    
    % Count unique colors in rules before scaling based on total
    % number of rules:
    nClrsInRules = 0;
    imRulesTemp = imRules;
    for i = 1 : szRules(2) 
    for j = 1 : szRules(1) 

        clr = imRulesTemp( j, i, : );
        if clr(1) ~= 999;
            
            imRulesTemp;
            nClrsInRules = nClrsInRules + 1;
            for ii = 1 : szRules(2) 
            for jj = 1 : szRules(1)
                
                if imRulesTemp( jj, ii, : ) == clr;

                    imRulesTemp( jj, ii, : ) = 999;
                end
            end
            end
        end
    end
    end

    imRulesCopyTemp = imRulesCopy;
    for i = 0 : nClrsInRules - 1
        iLowest = find( imRulesCopy == min(min(min(imRulesCopy))) );
        imRulesCopyTemp( iLowest ) = i;
        imRulesCopy( iLowest ) = 99999;
        clear iLowest
    end
    imRulesCopy = imRulesCopyTemp;
   
    imRulesCopy = reshape( imRulesCopy, szRules(1)*szRules(2), 3 );
    szRulesOnly = size(imRulesCopy);
    imInitCopy = arrayInit;
    szInit = size(arrayInit);
    %reshape( arrayInit, szInit(1)*szInit(2), 1 )
    % Append the initial array to the rules.
    for i = 1 : szInit(1)
    for j = 1 : szInit(2)
        imRulesCopy( szRulesOnly(1) + j + szInit(2)*(i-1), 1 ) = imInitCopy( i, j );
        imRulesCopy( szRulesOnly(1) + j + szInit(2)*(i-1), 2 ) = imInitCopy( i, j );
        imRulesCopy( szRulesOnly(1) + j + szInit(2)*(i-1), 3 ) = imInitCopy( i, j );
    end
    end
    szRulesCopy = size( imRulesCopy );
    imRulesGray = double( 2*imRulesCopy(:,1) ) + 3*double( imRulesCopy(:,2) ) + 5*double( imRulesCopy(:,3) );
    for i = 1 : nSymbols 

        iLeastBright = -1;
        nLeastBright = min(min(imRulesGray));
        iLeastBright = find( abs(imRulesGray - nLeastBright) < 4 );
        if iLeastBright(1) <= szRulesOnly(1)
            
            SymbolColor( i, 1 ) = imRules( iLeastBright(1) +   0              );
            SymbolColor( i, 2 ) = imRules( iLeastBright(1) +   szRulesOnly(1) );
            SymbolColor( i, 3 ) = imRules( iLeastBright(1) + 2*szRulesOnly(1) );
        else
            SymbolColor( i, 1 ) = imRulesCopy( iLeastBright(1) +   0              )*255*nSymbols/(nClrsInRules*nClrsInRules);
            SymbolColor( i, 2 ) = imRulesCopy( iLeastBright(1) +   szRulesCopy(1) )*255*nSymbols/(nClrsInRules*nClrsInRules);
            SymbolColor( i, 3 ) = imRulesCopy( iLeastBright(1) + 2*szRulesCopy(1) )*255*nSymbols/(nClrsInRules*nClrsInRules);
        end
        
        imRulesGray( iLeastBright ) = 999;
    end

    SymbolColor = SymbolColor*nClrsInRules/nSymbols;
    imRules     = imRules*nClrsInRules/nSymbols;
% ToDo: Check for more colors that nSymbols
    
     szRules(2) = szRules(2)/nSymbols;
    % To Do: sanity check on size of input tile grid

    arrayRules = -1;
    % Map colors to integer array rules, in the range [0 nSymbols-1].
    for r = 1 : nSymbols
    for j = 1 : szRules(2)
    for k = 1 : szRules(1)
%         if imRules( 1:szRules(1), (i-1)*szRules(2) + 1 : i*szRules(2), : ) == SymbolColor( i, : )
%         arrayRules( 1:szRules(1), 1:szRules(2), i ) = imRules( 1:szRules(1), (i-1)*szRules(2) + 1 : i*szRules(2) );
        for i = 1 : nSymbols
            if ( imRules( k, j + (r-1)*szRules(2), 1 ) == SymbolColor(i,1) && imRules( k, j + (r-1)*szRules(2), 2 ) == SymbolColor(i,2) & imRules( k, j + (r-1)*szRules(2), 3 ) == SymbolColor(i,3) )
                arrayRules( k, j, r ) = i - 1;
            end
        end
    end
    end 
    end
        
    if isempty( strMotifFilename )
        
        for i = 1:nSymbols
            imMotifs( 1:1, i, 1 ) = SymbolColor( i, 1 );
            imMotifs( 1:1, i, 2 ) = SymbolColor( i, 2 );
            imMotifs( 1:1, i, 3 ) = SymbolColor( i, 3 );
        end
    end
    szMotifs = size( imMotifs );

else
    
    if bRandom == 0 && ~isa( varargin{1}, 'char' ) 
        
        % Load explicit rules from arguements:
        for i = 1 : nargin - 7

            szRules = size( varargin{i} );
            arrayRules( 1:szRules(1), 1:szRules(2), i ) = double( varargin{i} );
        end
        
        for i = 1 : nSymbols 
            SymbolColor( i, 1:3 ) = i;
        end
    end
end



% Segregate rule numbers, rotation codes and mirror codes.
arrayRotes = zeros( size(arrayRules) );
arrayMirrs = zeros( size(arrayRules) );

arrayMirrs( arrayRules >= 4*nSymbols ) = 1;
arrayRules( arrayRules >= 4*nSymbols ) = arrayRules( arrayRules >= 4*nSymbols ) - 4*nSymbols;

%        arrayRotes( arrayRules >= nSymbols ) = ceil( ( arrayRules( arrayRules >= nSymbols ) - nSymbols )/nSymbols )
arrayRotes( arrayRules >= nSymbols ) = floor( arrayRules( arrayRules >= nSymbols )/nSymbols );
arrayRules( arrayRules >= nSymbols ) = mod( arrayRules( arrayRotes > 0 ), nSymbols );

% Segregate initial rule set, rotation codes and mirror codes.
arrayInitAll = arrayInit;   % file name encoding uses original composite numbering

arrayInitRotes = zeros( size(arrayInit) );
arrayInitMirrs = zeros( size(arrayInit) );

arrayInitMirrs( arrayInit >= 4*nSymbols ) = 1;
arrayInit( arrayInit >= 4*nSymbols ) = arrayInit( arrayInit >= 4*nSymbols ) - 4*nSymbols;

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

if ~isempty( arrayInit( arrayInit >= nSymbols ) )
    arrayInitRotes( arrayInit >= nSymbols ) = floor( arrayInit( arrayInit >= nSymbols )/nSymbols );
end
if ~isempty( arrayInit( arrayInitRotes > 0 ) )
    arrayInit( arrayInit >= nSymbols ) = mod( arrayInit( arrayInitRotes > 0 ), nSymbols );
end

        
        
% Sanity check that the number of rules matches that of the arguement:
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

% Construct output's filestem:

strOutFilestem = [ num2str(szRules(2)) 'x' num2str(szRules(1)) '-' num2str(nSymbols) 's' ];
if ~isempty( find( arrayRotes , 1) )
  strOutFilestem = [ strOutFilestem 'r' ];
end
if ~isempty( find( arrayMirrs , 1) )
  strOutFilestem = [ strOutFilestem 'm' ];
end
strOutFilestem = [ strOutFilestem '-' ];

% Prepend:    
if isempty( strOutFilestemOpt )
    strOutFilestem = [ 'LS-' strOutFilestem ];
else
    strOutFilestem = [ strOutFilestemOpt '-' strOutFilestem ];
end

% Append initial symbol(s) code.
szInit = size( arrayInit );
if max( szInit ) == 1
        strOutFilestem = [ strOutFilestem 'i' num2str( arrayInit ) ];    
else
    if max( szInit ) <= max( szRules )
        
        strOutFilestem = [ strOutFilestem 'i' num2str( szInit(1) ) 'x' num2str( szInit(2) ) '-' ];
        
        for i = 1 : szInit(1)
        for j = 1 : szInit(2)
            if i == 1 && j == 1
                strOutFilestem = [ strOutFilestem num2str( arrayInitAll(i,j) ) ];
            else
                strOutFilestem = [ strOutFilestem  '_' num2str( arrayInitAll(i,j) ) ];
            end
        end
        end
    end
end

% Prepend rotation and mirror codes.
if bRandom == 0 
    if nSymbols <= 7 && ~isempty( find( nCodes(:,2) , 1) )
        for i = 1 : nSymbols
            if i == 1
                strOutFilestem = [ strOutFilestem '-' num2str( nCodes(i,2) ) ];
            else
                strOutFilestem = [ strOutFilestem '_' num2str( nCodes(i,2) ) ];
            end
        end
    end

    % Prepend rule codes.
    if nSymbols <= 7
        for i = 1 : nSymbols
            if i == 1
                strOutFilestem = [ strOutFilestem '-' num2str( nCodes(i,1) ) ];
            else
                strOutFilestem = [ strOutFilestem '_' num2str( nCodes(i,1) ) ];
            end
        end
    end
end

if nDecimate > 1 || nOffset > 0
    strOutFilestem = [ strOutFilestem '-' num2str( nDecimate ) 'd_' num2str( nOffset ) 'o' ];
end

strOutFilestem = [ strOutFilestem ];

if ~isempty( strMotifFilename )
    strOutFilestem = [ strOutFilestem '-' strMotifFilename(1:length(strMotifFilename)-4) ];
end

if isa( arrayInit, 'char' )
    strOutFilestem = [ strOutFilestem '-' strInit ];
end

%strOutFilestem = [ strOutFilestem '.png' ];

szInit = size( arrayInit );

szTiling = [ floor( 1 + ( szInit(1)*(szRules(1)^nGenerations) - (nOffset + 1) )/nDecimate )  floor( 1 + ( szInit(1)*(szRules(1)^nGenerations) - (nOffset + 1) )/nDecimate ) ];

if bShowResults
    fprintf( ['\n'] );
end

% Read or initialize motifs.
if ~isempty( strMotifFilename )
    
    % Load any motif image into a 24 bit color array:
    infMotif = imfinfo( strMotifFilename );
    if infMotif.ColorType == 'truecolor'
        if infMotif.Height == 1
            imMotifs( 1, 1:infMotif.Width, 1:3 ) = imread( strMotifFilename );
        else
            imMotifs = imread( strMotifFilename );
        end       
    else if infMotif.ColorType == 'grayscale'
            
        imMotifsGray = imread( strMotifFilename );
        szGrayMotifs  = size( imMotifsGray )
        imMotifs( 1:szGrayMotifs(1), 1:szGrayMotifs(2), 1:3 ) = 0;
       
        imMotifs(:,:,1) = imMotifsGray;
        imMotifs(:,:,2) = imMotifsGray;
        imMotifs(:,:,3) = imMotifsGray;
        
        else if infMotif.ColorType == 'indexed'
                fprintf( [ '\n\nIndexed image format for the motifs is not currently supported.\n\n' ] );
                return;
            end
        end
    end
    
    szMotifs = size( imMotifs );
    % To Do: sanity check on size of input tile grid
% szTiling(1)
% szMotifs(1)

    if bShowResults
        fprintf( [ 'Size of image:  ' num2str( nOversampleOut*szTiling(1)*szMotifs(1) ) ' x ' num2str( nOversampleOut*szTiling(2)*szMotifs(1) ) ' px. \n' ] );
        fprintf( [ 'Memory required (final image):  ' num2str( nOversampleOut^2*szTiling(1)*szTiling(2) * szMotifs(1)^2 *3/(2^20) ) ' MB \n' ] );
    end

else
    if nargin ~= 8 || bRandom == 1
% szTiling
        nimPattern = 255*ones( szTiling(1), szTiling(2), 'uint8' );
        szMotifs = [1,1];

        if bShowResults
            fprintf( [ 'Size of image:  ' num2str( nOversampleOut*szTiling(1) ) ' x ' num2str( nOversampleOut*szTiling(2) ) ' px. \n' ] );
            fprintf( [ 'Memory required (final image):  ' num2str( nOversampleOut^2*szTiling(1)*szTiling(2)*3/(2^20) ) ' MB \n' ] );
        end
    end
end

imMotifs = double(imMotifs)/255;

% Initialize output image.
if isempty( strMotifFilename )
    imOut = 255*zeros( szTiling(1)*szMotifs(1), szTiling(2)*szMotifs(2),    'uint8' );
else
    imOut = 255*zeros( szTiling(1)*szMotifs(1), szTiling(2)*szMotifs(2), 3, 'uint8' );
end

% Initialize average image.
if bDrawAverage
%    imAvg = zeros( 2^nGenerations, 2^nGenerations, 'uint16' );
	imAvg = zeros( szTiling(1), szTiling(2), 'uint16' );
end

% Construct an eighth generation sequence, for checking global symmetry:
if szRules(2) > 1  && bRandom == 0 && szRules(1) < 3
    
    nByGenerations = 6;

    vctLS0 = arrayInit;
    szLS = size(vctLS0);
    vctRT0 = arrayInitRotes;
    vctMR0 = arrayInitMirrs;
    vctLS( 1:szLS(1)*szRules(1), 1:szLS(2)*szRules(2) ) = 0;
    vctRT( 1:szLS(1)*szRules(1), 1:szLS(2)*szRules(2) ) = 0;
    vctMR( 1:szLS(1)*szRules(1), 1:szLS(2)*szRules(2) ) = 0;

    for ii = 1 : nByGenerations

        % Construct the next generation of the pattern.
        for j = 1 : szLS(1)
        for k = 1 : szLS(2)

            vctLS( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) = rot90( arrayRules( 1:szRules(1), 1:szRules(2), vctLS0( j, k ) + 1 ), -vctRT0( j, k ) );
            if mod( vctMR0( j, k ), 2 ) == 1
                vctLS( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) = fliplr( vctLS( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) );
            end
            vctRT( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) = vctRT0( j, k );
            vctMR( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) = vctMR0( j, k );
            vctRT( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) = vctRT( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) + rot90( arrayRotes( 1:szRules(1), 1:szRules(2), vctLS0( j, k ) + 1 ), -vctRT0( j, k ) );
            if mod( vctMR0( j, k ), 2 ) == 1
                vctMR( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) = vctMR( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) + fliplr( arrayMirrs( 1:szRules(1), 1:szRules(2), vctLS0( j, k ) + 1 ) );
            else
                vctMR( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) = vctMR( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) + arrayMirrs( 1:szRules(1), 1:szRules(2), vctLS0( j, k ) + 1 );
            end

        end
        end

        vctLS0 = vctLS;
        vctRT0 = vctRT;
        vctMR0 = vctMR;
        szLS = size(vctLS0);          
    end
    
    imEighthGen = vctLS;
end

clear vctLS0;
clear vctRT0;
clear vctMR0;
clear vctLS;
clear vctRT;
clear vctMR;

% Construct the symmetry graphic:
if bRandom == 0  &&  szRules(1) > 1 && szRules(1) < 3 &&  szRules(2) > 1 && szRules(2) < 3 && szInit(1) == 1 && szInit(2) == 1
    
% To Do: Accomodate rules that are larger than 2x2 and initializations
% larger than 1x1.
    nSymGenerations = 8;
    if nGenerations + log2(nOversampleOut) == 7
    	nSymGenerations = 7;
    end
    if nGenerations + log2(nOversampleOut) <= 6
    	nSymGenerations = 6;
    end

    fBlendfrac  = .25;
    fBlendfrac2 = 1 - fBlendfrac;
    imScaleSym  =  0;
    
    imScaleSymGraphic8 = 0;

    % Find symmetry of full pattern.
% To Do: Write this as a separate function.
    
%    imRule1Sym = imRules( :,:, arrayInit + 1 );
%    imRule1Sym( find( imRule1Sym ~= arrayInit ) ) = 0;

    strSymmetry = '';
    
    szEighthGen = size(imEighthGen);
    imEighthGenRot90  = flipdim( permute( imEighthGen, [2 1 3] ), 2 );
    imEighthGenRot180 = flipdim( flipdim( imEighthGen, 1 ), 2 );
    imEighthGenRot270 = flipdim( permute( imEighthGen, [2 1 3] ), 1 );
    bDiff1 = sum(sum(abs( imEighthGen - imEighthGenRot90  )));
    bDiff2 = sum(sum(abs( imEighthGen - imEighthGenRot180 )));
%     figure
%     imshow(abs( imEighthGen - imEighthGenRot180 ))
%     imEighthGen - imEighthGenRot180
    bDiff3 = sum(sum(abs( imEighthGen - imEighthGenRot270 )));

    imEighthGenDiagMirror  = permute( imEighthGen, [2 1 3] );
    imEighthGenDiag2Mirror = flipdim( permute( flipdim( imEighthGen, 2 ), [2 1 3] ), 2 );
    bDiff4 = sum(sum(abs( imEighthGen - imEighthGenDiagMirror  )));
    bDiff5 = sum(sum(abs( imEighthGen - imEighthGenDiag2Mirror )));
    
    imEighthGenMirrorX  = flipdim( imEighthGen, 2 );
    imEighthGenMirrorY  = flipdim( imEighthGen, 1 );
    bDiff6 = sum(sum(abs( imEighthGen - imEighthGenMirrorX  )));
    bDiff7 = sum(sum(abs( imEighthGen - imEighthGenMirrorY )));
   
    bSigma7 = sum(std( imEighthGen + imEighthGenMirrorX ));
    bSigma2 = sum(std( imEighthGen + imEighthGenRot180 ));
    
% figure
% imshow(imEighthGen)
% figure
% imshow(imEighthGenMirrorY)



    bSwap = false;
    strSymmetryGlobal   = [];
    strSymmetryRules    = [];
    
    % Check for identity rule swap:
     listSymbolPairs = [];
     nPairs = [];
     [ nPairs, listSymbolPairs ] = get_number_pairs( arrayRules(:,:,1), arrayRules(:,:,2) );

	if nPairs < nSymbols...
        && any( listSymbolPairs( :, 1 ) ~= listSymbolPairs( :, 2 ) )  % not all equal pairs
   
        % Update a symbol pair list by finding identical pairs and checking if a
        % rule system references them alternately, replacing the identical pair
        % with an alternating pair.
        [ nAlternatingReplacements listSymbolPairsSwapUse ] = add_alternating_pairs( listSymbolPairs, arrayRules );
        nPairsSwap = nPairs;
        strSymmetry = 'n_identity';
        bSwap = true;
    end
    
    % Check for symmetric rule swap across center vertical:
    imEighthGenRight = imEighthGen( :, 1:szEighthGen(2)/2 );
    imEighthGenLeft  = imEighthGen( :, szEighthGen(2)/2 + 1 : szEighthGen(2) );
    
    listSymbolPairs = [];
    nPairs = [];
    [ nPairs, listSymbolPairs ] = get_number_pairs( imEighthGenRight, imEighthGenLeft );
    
    if    nPairs < nSymbols ...
       && any( listSymbolPairs( :, 1 ) ~= listSymbolPairs( :, 2 ) )  % not all equal pairs
   
        % Update a symbol pair list by finding identical pairs and checking if a
        % rule system references them alternately, replacing the identical pair
        % with an alternating pair.
        [ nAlternatingReplacements listSymbolPairsSwapUse ] = add_alternating_pairs( listSymbolPairs, arrayRules );
%        if nAlternatingReplacements
            nPairsSwap = nPairs;
            strSymmetry = 'Swap_n_lr';
            bSwap = true;
%        end
    end

%     % Check for mirror symmetric rule swap across center vertical:
%     listSymbolPairs = [];
%     nPairs = [];
%     [ nPairs, listSymbolPairs ] = get_number_pairs( imEighthGen, imEighthGenMirrorX );
% 
%     if nPairs < nSymbols...
%        && any( listSymbolPairs( :, 1 ) ~= listSymbolPairs( :, 2 ) )  % not all equal pairs
%    
%         % Update a symbol pair list by finding identical pairs and checking if a
%         % rule system references them alternately, replacing the identical pair
%         % with an alternating pair.
%         [ nAlternatingReplacements listSymbolPairsSwapUse ] = add_alternating_pairs( listSymbolPairs, arrayRules );
%         if nAlternatingReplacements
%             nPairsSwap = nPairs;
%             strSymmetry = 'Mirror_n_lr';
%             bSwap = true;
%         end
%     end

    % Check for pi rotation symmetric rule swap:
    listSymbolPairs = [];
    nPairs = [];
    [ nPairs, listSymbolPairs ] = get_number_pairs( imEighthGen, imEighthGenRot180 );
    
	if nPairs < nSymbols...
       && any( listSymbolPairs( :, 1 ) ~= listSymbolPairs( :, 2 ) )  % not all equal pairs
   
        % Update a symbol pair list by finding identical pairs and checking if a
        % rule system references them alternately, replacing the identical pair
        % with an alternating pair.
        [ nAlternatingReplacements listSymbolPairsSwapUse ] = add_alternating_pairs( listSymbolPairs, arrayRules );
%        if nAlternatingReplacements
            nPairsSwap = nPairs;
            strSymmetryGlobal = 'Rot180_n';
            bSwap = true;
%        end
        
        % What is the symmetry of the rule set?
        % What symmetry transform, applied to every rule results in an
        % equivalent set of rules?
        %   swap all pairs (identity after swap)
        %       OR
        %   mirror across one of four axes, and swap all pairs
        %       OR
        %   mirror across a pair of four axes, and swap all pairs 
        %strSymmetryRules = ?
        
        % Test mirror about x-axis:
        arrayRulesTransformed = arrayRules;
        % Make a mirrored copy
        for ir = 1 : nSymbols
            arrayRulesTransformed( :, :, ir ) = flipdim( arrayRules( :, :, ir ), 1 );
        end
        arrayRulesTransformedCopy = arrayRulesTransformed;
        % Swap symbols of mirrored copy
        for ip = 1: nPairsSwap
            arrayRulesTransformed( find( arrayRulesTransformedCopy == listSymbolPairsSwapUse( ip, 1 ) ) ) = listSymbolPairsSwapUse( ip, 2 );
            arrayRulesTransformed( find( arrayRulesTransformedCopy == listSymbolPairsSwapUse( ip, 2 ) ) ) = listSymbolPairsSwapUse( ip, 1 );
        end
        % Swap the positions of the rules.
        arrayRulesTransformedCopy = arrayRulesTransformed;
%         for ir = 1 : nSymbols
        for ip = 1: nPairsSwap
            arrayRulesTransformed( :, :, listSymbolPairsSwapUse( ip, 1 ) + 1 ) = arrayRulesTransformedCopy( :, :, listSymbolPairsSwapUse( ip, 2 ) + 1 );
            arrayRulesTransformed( :, :, listSymbolPairsSwapUse( ip, 2 ) + 1 ) = arrayRulesTransformedCopy( :, :, listSymbolPairsSwapUse( ip, 1 ) + 1 );
        end
%         end
        
        bIdentical = true;
        for ir = 1 : nSymbols
%             arrayRulesTransformed( :, :, ir ) = flipdim( arrayRules( :, :, ir ), 1 );
            if ~all( arrayRulesTransformed( :, :, ir ) == arrayRules ( :, :, ir ) )
                
                bIdentical = false;
            end
        end
        
        if bIdentical == true
            strSymmetryRules = 'mirror_x'
        end
    end

    % Check for mirror symmetric rule swap across center vertical:
    if    bDiff4 == 0 ...
       && bDiff5 == 0 ...
       
       strSymmetryGlobal = 'Symmetric_dd';
    end
    
    
    
    
     listSymbolPairs = [];
     nPairs = [];
     [ nPairs, listSymbolPairs ] = get_number_pairs( arrayRules(:,:,1), fliplr( arrayRules(:,:,2) ) );

     if nPairs < nSymbols...
        && any( listSymbolPairs( :, 1 ) ~= listSymbolPairs( :, 2 ) )  % not all equal pairs
   
        % Update a symbol pair list by finding identical pairs and checking if a
        % rule system references them alternately, replacing the identical pair
        % with an alternating pair.
        [ nAlternatingReplacements listSymbolPairsSwapUse ] = add_alternating_pairs( listSymbolPairs, arrayRules );
        nPairsSwap = nPairs;
        strSymmetry = 'Mirror_n_lr';
        bSwap = true;
     end
    
     

    listSymbolPairs = [];
    nPairs = [];
    % Construct a symbol pair list of all unique spatially corresponding
    % pairs of symbols. 
    [ nPairs, listSymbolPairs ] = get_number_pairs( imEighthGen, imEighthGenRot180 );
    
    if nPairs < nSymbols...
       && any( listSymbolPairs( :, 1 ) ~= listSymbolPairs( :, 2 ) )  % not all equal pairs 
        
        %   Update a number pair list by finding identical pairs and checking if a
        %   rule system references them alternately, replacing the identical pair
        %   with an alternating pair.
        [ nAlternatingReplacements listSymbolPairsSwapUse ] = add_alternating_pairs( listSymbolPairs, arrayRules );
        nPairsSwap = nPairs;
        strSymmetryGlobal = 'Rot180_n';
        bSwap = true;
        
        % What is the symmetry of the rule set?
        % What symmetry transform, applied to every rule results in an
        % equivalent set of rules?
        %   swap all pairs (identity after swap)
        %       OR
        %   mirror across one of four axes, and swap all pairs
        %       OR
        %   mirror across a pair of four axes, and swap all pairs 
        %strSymmetryRules = ?
        
        % Test mirror about x-axis:
        arrayRulesTransformed = arrayRules;
        % Make a mirrored copy
        for ir = 1 : nSymbols
            arrayRulesTransformed( :, :, ir ) = flipdim( arrayRules( :, :, ir ), 1 );
        end
        
        % Swap symbols of mirrored copy
        arrayRulesTransformedCopy = arrayRulesTransformed;
        for ip = 1: nPairsSwap
            arrayRulesTransformed( find( arrayRulesTransformedCopy == listSymbolPairsSwapUse( ip, 1 ) ) ) = listSymbolPairsSwapUse( ip, 2 );
            arrayRulesTransformed( find( arrayRulesTransformedCopy == listSymbolPairsSwapUse( ip, 2 ) ) ) = listSymbolPairsSwapUse( ip, 1 );
        end
        
        % Swap the positions of the rules.
        arrayRulesTransformedCopy = arrayRulesTransformed;
%         for ir = 1 : nSymbols
        for ip = 1: nPairsSwap
            arrayRulesTransformed( :, :, listSymbolPairsSwapUse( ip, 1 ) + 1 ) = arrayRulesTransformedCopy( :, :, listSymbolPairsSwapUse( ip, 2 ) + 1 );
            arrayRulesTransformed( :, :, listSymbolPairsSwapUse( ip, 2 ) + 1 ) = arrayRulesTransformedCopy( :, :, listSymbolPairsSwapUse( ip, 1 ) + 1 );
        end
%         end
        
        bIdentical = true;
        for ir = 1 : nSymbols
%             arrayRulesTransformed( :, :, ir ) = flipdim( arrayRules( :, :, ir ), 1 );
            if ~all( arrayRulesTransformed( :, :, ir ) == arrayRules ( :, :, ir ) )
                
                bIdentical = false;
            end
        end
        
        if bIdentical == true
            strSymmetryRules = 'mirror_x'
        end
    end

    
    bSigma2;
    if bSigma2 == 0
            strSymmetryGlobal = 'Rot180';
    end
    if bDiff4 == 0
        strSymmetryGlobal = 'Diag_tl_br';
    end
    if bDiff5 == 0
        strSymmetryGlobal = 'Diag_tr_bl';    
    end
    if bDiff6 == 0
        strSymmetryGlobal = 'Mirror_lr';
    end
    if bDiff7 == 0
        strSymmetryGlobal = 'Mirror_tb';
    end
    
    if    bDiff4 == 0 ...
       && bDiff5 == 0
            if    bDiff4 == 0 ...
       && bDiff5 == 0
        strSymmetryGlobal = 'Both diagonals';    
    end
 
    end

    if    bDiff4 == 0 ...
       && bDiff5 == 0 ...    
       && bDiff6 == 0 ...    
       && bDiff7 == 0 ...    
        strSymmetryGlobal = 'Symmetric';
    end
    
    if    bDiff1  ~= 0 ... 
       && bSigma2  ~= 0 ... 
       && bDiff3  ~= 0 ... 
       && bDiff2  ~= 0 ... 
       && bDiff4  ~= 0 ...  
       && bDiff5  ~= 0 ...
       && bDiff6  ~= 0 ...  
       && bDiff7  ~= 0 ...
       && bSigma7 ~= 0 ...
       && ~strcmp( strSymmetryGlobal, 'Symmetric_dd')   ...
       && ~strcmp( strSymmetryGlobal, 'Swap_n_lr')   ...
       && ~strcmp( strSymmetryGlobal, 'Mirror_n_lr') ...
       && ~strcmp( strSymmetryGlobal, 'Rot180_n')
        strSymmetryGlobal = 'Asymmetric';
    end
    
    if isempty(strSymmetryGlobal)
        strSymmetryGlobal = 'undetermined';
    end
    
    if bShowResults
        fprintf( [ 'The symmetry of the pattern is ' strSymmetryGlobal ] );
    end
    
    if strcmp( strSymmetryGlobal, 'Symmetric_d_n180a' )
    	imScaleSymGraphic8 = imread( 'Scale_symmetry_Od_n180a_8.png' );
    end
    if strcmp( strSymmetryGlobal, 'Rot180' ) || strcmp( strSymmetryGlobal, 'Rot180_n' )
    	imScaleSymGraphic8 = imread( 'Scale_symmetry_FrF_8.png' );
    end
    if strcmp( strSymmetryGlobal, 'Mirror_n_lr' )
    	imScaleSymGraphic8 = imread( 'Scale_symmetry_Uv_8.png' );
    end
    if strcmp( strSymmetryGlobal, 'Mirror_lr' )
    	imScaleSymGraphic8 = imread( 'Scale_symmetry_U_8.png' );
    end
    if strcmp( strSymmetryGlobal, 'Mirror_tb' )
        imScaleSymGraphic8 = imread( 'Scale_symmetry_C_8.png' );
    end

    if strcmp( strSymmetryGlobal, 'Diag_tl_br' ) || strcmp( strSymmetryGlobal, 'Diag_tr_bl' )
        imScaleSymGraphic8 = flipdim( permute( imread( 'Scale_symmetry_L_8.png' ), [ 2 1 3 ] ), 2 );
         
    if strcmp( strSymmetryGlobal, 'Diag_tr_bl' )
        imScaleSymGraphic8 = flipdim( permute( imScaleSymGraphic8, [2 1 3] ), 2 );
    end
    end
%    if strcmp( strSymmetryGlobal, 'Symmetric_dd' )
    if strcmp( strSymmetryGlobal, 'Both diagonals' )
        imScaleSymGraphic8 = imread( 'Scale_symmetry_D_8.png' );
    end
    if strcmp( strSymmetryGlobal, 'Symmetric' )
        imScaleSymGraphic8 = imread( 'Scale_symmetry_O_8.png' );
    end
    if strcmp( strSymmetryGlobal, 'Swap_n_lr' )
        imScaleSymGraphic8 = imread( 'Scale_symmetry_FF_8.png' );
    end
    if    strcmp( strSymmetryGlobal, 'Asymmetric' ) ...
       || strcmp( strSymmetryGlobal, 'undetermined' )
%        imScaleSymGraphic8 = flipdim( permute( imread( 'Scale_symmetry_F_8.png' ), [ 2 1 3 ] ), 2 );
        imScaleSymGraphic8 = imread( 'Scale_symmetry_F_8.png' );
    end
    
    imScaleSymGraphic8 = double(imScaleSymGraphic8)/255.0;
    
    % Construct a scale pyramid of the symmetry graphic.
    impyScaleSymGraphic     = make_image_pyramid( imScaleSymGraphic8, 'linear' );
    impyScaleSymGraphicSwap = impyScaleSymGraphic;
    
    % If there are swapped pairs of symbols(?? define this: 
    % "Systems that are identical after a swapping of at least
    % one pair of symbols, and the rules are not identical?),
    % create a symmetry graphic with black and white swapped:
    if bSwap

        imScaleSymGraphic8Swap = imScaleSymGraphic8;

        % Swap black and white.
        szScaleSymGraphic8 = size( imScaleSymGraphic8 );
        for iSG = 1 : szScaleSymGraphic8(1)
        for jSG = 1 : szScaleSymGraphic8(2)

            % Swap black and white on the scale graphic to
            % indicate there is a symbol swap.
            if    imScaleSymGraphic8( iSG, jSG, 1 )  < .0001 ...
               && imScaleSymGraphic8( iSG, jSG, 2 )  < .0001 ...
               && imScaleSymGraphic8( iSG, jSG, 3 )  < .0001 
                imScaleSymGraphic8Swap( iSG, jSG, : ) = 1;
            end
            if    imScaleSymGraphic8( iSG, jSG, 1 )  > .9999 ...
               && imScaleSymGraphic8( iSG, jSG, 2 )  > .9999 ...
               && imScaleSymGraphic8( iSG, jSG, 3 )  > .9999 
                imScaleSymGraphic8Swap( iSG, jSG, : ) = 0;
            end
        end
        end
        
        impyScaleSymGraphicSwap = make_image_pyramid( imScaleSymGraphic8Swap, 'linear' );
    end

    % Maintain transparencies for nearby interpolated values.
	for issg = 1 : size( impyScaleSymGraphic, 2 )
        
        it = find(   impyScaleSymGraphic{issg}(:,:,1) > 128/255.0 - .1 ...
                   & impyScaleSymGraphic{issg}(:,:,2) > 128/255.0 - .1 ...
                   & impyScaleSymGraphic{issg}(:,:,3) > 128/255.0 - .1 ...
                   & impyScaleSymGraphic{issg}(:,:,1) < 128/255.0 + .1 ...
                   & impyScaleSymGraphic{issg}(:,:,2) < 128/255.0 + .1 ...
                   & impyScaleSymGraphic{issg}(:,:,3) < 128/255.0 + .1    );
        impyScaleSymGraphic{issg}(it) = 128/255.0;
        impyScaleSymGraphic{issg}(it + 0*size(impyScaleSymGraphic{issg},1)*size(impyScaleSymGraphic{issg},2)) = 128/255.0;
        impyScaleSymGraphic{issg}(it + 1*size(impyScaleSymGraphic{issg},1)*size(impyScaleSymGraphic{issg},2)) = 128/255.0;
        impyScaleSymGraphic{issg}(it + 2*size(impyScaleSymGraphic{issg},1)*size(impyScaleSymGraphic{issg},2)) = 128/255.0;
        
        it = find(   impyScaleSymGraphicSwap{issg}(:,:,1) > 128/255.0 - .1 ...
                   & impyScaleSymGraphicSwap{issg}(:,:,2) > 128/255.0 - .1 ...
                   & impyScaleSymGraphicSwap{issg}(:,:,3) > 128/255.0 - .1 ...
                   & impyScaleSymGraphicSwap{issg}(:,:,1) < 128/255.0 + .1 ...
                   & impyScaleSymGraphicSwap{issg}(:,:,2) < 128/255.0 + .1 ...
                   & impyScaleSymGraphicSwap{issg}(:,:,3) < 128/255.0 + .1    );
        impyScaleSymGraphicSwap{issg}(it) = 128/255.0;
        impyScaleSymGraphicSwap{issg}(it + 0*size(impyScaleSymGraphicSwap{issg},1)*size(impyScaleSymGraphicSwap{issg},2)) = 128/255.0;
        impyScaleSymGraphicSwap{issg}(it + 1*size(impyScaleSymGraphicSwap{issg},1)*size(impyScaleSymGraphicSwap{issg},2)) = 128/255.0;
        impyScaleSymGraphicSwap{issg}(it + 2*size(impyScaleSymGraphicSwap{issg},1)*size(impyScaleSymGraphicSwap{issg},2)) = 128/255.0;
        
    end
    
                                


    imScaleSym = impyScaleSymGraphic{nSymGenerations+1};
    szimScaleSym = size(imScaleSym);
    imScaleSymCovered     = zeros( szimScaleSym(1), szimScaleSym(1), nSymGenerations );
    imScaleSymCoveredTemp = zeros( szimScaleSym(1), szimScaleSym(1), nSymGenerations );
    % Decrease (compress) range of full scale symmetry graphic background (decrease contrast).
    imScaleSym = fBlendfrac*impyScaleSymGraphic{nSymGenerations+1} + (1-fBlendfrac)/2;
 
    if bRandom == 0 
         
        vctLS0 = arrayInit;
        szLS = size(vctLS0);
        vctRT0 = arrayInitRotes;
        vctMR0 = arrayInitMirrs;
        vctLS( 1:szLS(1)*szRules(1), 1:szLS(2)*szRules(2) ) = 0;
        vctRT( 1:szLS(1)*szRules(1), 1:szLS(2)*szRules(2) ) = 0;
        vctMR( 1:szLS(1)*szRules(1), 1:szLS(2)*szRules(2) ) = 0;

%        imScaleSym = impyScaleSymGraphic{8+1};
%        size(imScaleSym)

        rotarrayRules       = 0;
        rotarrayRotes       = 0;
        rotarrayMirrs       = 0;
        fliplrrotarrayRotes = 0;
        fliplrrotarrayMirrs = 0;

        for ils = 1 : nSymbols
        for irt = 1 : 4
        	rotarrayRules( 1:szRules(1), 1:szRules(2), ils, irt ) = rot90( arrayRules( 1:szRules(1), 1:szRules(2), ils ), -(irt-1) );
        end
        end
        
        for ils = 1 : nSymbols
        for irt = 1 : 4
            rotarrayRotes( 1:szRules(1), 1:szRules(2), ils, irt ) = rot90( arrayRotes( 1:szRules(1), 1:szRules(2), ils ), -(irt-1) );
            rotarrayMirrs( 1:szRules(1), 1:szRules(2), ils, irt ) = rot90( arrayMirrs( 1:szRules(1), 1:szRules(2), ils ), -(irt-1) );
            fliplrrotarrayRotes( 1:szRules(1), 1:szRules(2), ils, irt ) = fliplr( rotarrayRotes( 1:szRules(1), 1:szRules(2), ils, irt ) );
            fliplrrotarrayMirrs( 1:szRules(1), 1:szRules(2), ils, irt ) = fliplr( rotarrayMirrs( 1:szRules(1), 1:szRules(2), ils, irt ) );
        end
        end

        cellLS = cell( zeros( 1, nGenerations ) );
        cellRT = cell( zeros( 1, nGenerations ) );
        cellLS = cell( zeros( 1, nGenerations ) );

        % Initialize 0th generation (first) cell:
        cellLS{0+1} = arrayInit;
        cellRT{0+1} = arrayInitRotes;
        cellMR{0+1} = arrayInitMirrs;
            
        % For each generation up to the size of the symmetry diagram:
        for ig = 1 : nSymGenerations

            % Construct the next generation of the pattern.
            for j = 1 : szLS(1)
            for k = 1 : szLS(2)
                
                % Replace each symbol in the current symbol array (vctLS0, zero based values) 
                % with the current rule block (arrayRules(:,:,vctLS0 + 1), one based index),
                % and rotate by the current rotation (-vctRT0( j, k )).
                vctLS( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) = rotarrayRules( 1:szRules(1), 1:szRules(2), vctLS0( j, k ) + 1 , vctRT0( j, k ) + 1 );
                % If the current mirror is 1, mirror (vertically) the block.
                if vctMR0( j, k ) == 1
                    vctLS( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) = fliplr( rotarrayRules( 1:szRules(1), 1:szRules(2), vctLS0( j, k ) + 1 , vctRT0( j, k ) + 1 ) );
                end

                % Mirror and/or rotate the mirror and rotation array blocks.
                if vctMR0( j, k ) == 1
                    vctRT( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) = vctRT0( j, k ) - fliplrrotarrayRotes( 1:szRules(1), 1:szRules(2), vctLS0( j, k ) + 1 , vctRT0( j, k ) + 1 );
                    vctMR( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) = vctMR0( j, k ) - fliplrrotarrayMirrs( 1:szRules(1), 1:szRules(2), vctLS0( j, k ) + 1 , vctRT0( j, k ) + 1 );
                else
                    vctRT( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) = vctRT0( j, k ) + rotarrayRotes( 1:szRules(1), 1:szRules(2), vctLS0( j, k ) + 1 , vctRT0( j, k ) + 1 );
                    vctMR( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) = vctMR0( j, k ) + rotarrayMirrs( 1:szRules(1), 1:szRules(2), vctLS0( j, k ) + 1 , vctRT0( j, k ) + 1 );
                end

            end
            end

            vctRT = mod( vctRT, 4 );
            vctMR = mod( vctMR, 2 );
            
            % Save each generation array in a cell:
            cellLS{ig+1} = vctLS;
            cellRT{ig+1} = vctRT;
            cellMR{ig+1} = vctMR;

            vctLS0 = vctLS;
            vctRT0 = vctRT;
            vctMR0 = vctMR;
            szLS = size(vctLS0);
        end    

        % Mark first order uniquely non-overlapping self-symmetries 
        % on symmetry diagram.
        %
        % For each generation up to half the size of the symmetry diagram:
%         for ig = 1 %: nSymGenerations - 1
                    
            % Check for any instance of the first order rule appearing at any location 
            % past the second generation of the pattern (a self-symmetry), 
            % and mark these on the self-similarity diagram.
            
            % To Do: Clean up and document.
            % Currently the algorithm for doing this is very obscure, and
            % some details are incorrect.
            
            % The algorithm operates on the 8th (or 7th, nSymGenerations)
            % generation pattern.
            % A match with the 7th (or 6th) generation pattern is
            % made, but only at integral offsets of the 7th (or 6th)
            % generation's width. 
            % [To Do: Check at integral offsets that are half this width?
            %  a corresponding diagram of these "fortuitous juxtapositions"?] 
            % This is equivalent finding correlations of +- 1 between
            % the second and first generation.
            
%            for nOrder = [ 1 0 ]   % 0 : ig - 1
            % For each generation from the size of the symmetry diagram
            % down to first generation:
            for nGenDelta = 1 : nSymGenerations - 0     % range: [ 1, nSymGenerations - 1 ]
%            if ig >= 1 + nOrder % && ig <= nSymGenerations - 1
            nGenTest = nSymGenerations - nGenDelta;     % range: [ nSymGenerations - 1, 1 ]
% To Do: Deprecate ig.
            ig = nGenTest;    

%                 szLSig1 = size( cellLS{   1 + nOrder} );
%                 szLSig  = size( cellLS{ig+1         } );

%                 hBlock   = szLSig1(1)*2^(nSymGenerations-ig) - 1;
%                 wBlock   = szLSig1(2)*2^(nSymGenerations-ig) - 1;
                hBlock   = szInit(1)*2^nGenTest - 1;
                wBlock   = szInit(1)*2^nGenTest - 1;
%                 nszBlock = nSymGenerations+1-ig + nOrder;
                             
%                 imScaleGraphic = impyScaleSymGraphic{nszBlock};
                imScaleGraphic      = impyScaleSymGraphic{nGenTest+1};
                imScaleGraphicSwap  = impyScaleSymGraphicSwap{nGenTest+1};
                
                % Get pattern block from the full pattern:
                szBlock = [ szInit(1)*szRules(1)^nGenTest szInit(2)*szRules(2)^nGenTest ];
                for iii = 0 : 2^nGenDelta - 1
                for jjj = 0 : 2^nGenDelta - 1
                    
                	yBlock   = iii*szInit(1)*szRules(1)^nGenTest + 1;
                    xBlock   = jjj*szInit(2)*szRules(2)^nGenTest + 1;

                    blockLS = cellLS{ nSymGenerations+1 }( yBlock : yBlock + szBlock(1) - 1, xBlock : xBlock + szBlock(2) - 1 );
                    blockRT = cellRT{ nSymGenerations+1 }( yBlock : yBlock + szBlock(1) - 1, xBlock : xBlock + szBlock(2) - 1 );
                    blockMR = cellMR{ nSymGenerations+1 }( yBlock : yBlock + szBlock(1) - 1, xBlock : xBlock + szBlock(2) - 1 );
                
                    % If there is a symbol swapping symmetry, make a swapped pattern block:
                    if bSwap
                        
                        blockLSSwap   = blockLS;
                        szblockLSSwap = size(blockLSSwap);
                        % Swap block symbols
                        for nb = 1 : szblockLSSwap(1)
                        for mb = 1 : szblockLSSwap(2)
                           
                         	for np = 1 : nPairsSwap
                                
                               	%if strcmp( strSymmetryGlobal, 'Swap_n_lr')
                                    if blockLSSwap( nb, mb ) == listSymbolPairsSwapUse( np, 1 )
                                        blockLSSwap( nb, mb ) = listSymbolPairsSwapUse( np, 2 );
                                    else
                                        if blockLSSwap( nb, mb ) == listSymbolPairsSwapUse( np, 2 )
                                            blockLSSwap( nb, mb ) = listSymbolPairsSwapUse( np, 1 );
                                        end
                                    end
                                %end
                            end
                        end
                        end
                    end  % if symbol swapped symmetry
                        
                    
                    % Translation only self-symmetry marking:
                    bTranslation = false; 
                    % If the test block is exact match on the pattern block:
                    if    ~any(any( blockLS - cellLS{nGenTest+1} )) ...
                       && ~any(any( blockRT - cellRT{nGenTest+1} )) ... 
                       && ~any(any( blockMR - cellMR{nGenTest+1} ))

                        bTranslation = true;
                        
                        % If symmetry diagram is not already fully covered by self-symmetries:
                        if length( find( imScaleSymCovered( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1 ) == 0 ) ) > 0

                            if nGenTest >= nSymGenerations - 8
                                imReplaceSquare = imScaleSym( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 );
                                niNotGray = find( imScaleGraphic ~= 128/255.0 );
                                imReplaceSquare( niNotGray ) = imScaleGraphic( niNotGray );
                                imScaleSym(        yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 ) = imReplaceSquare;
                                imScaleSymCovered( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1 ) = 1;
                            end
                            
                                
%                             % Only mark as covered, don't draw?
%                             if nGenTest <= nSymGenerations - 4
%                                 
%                                 imReplaceSquare = imScaleSym( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 );
%                                 imReplaceSquare = fBlendfrac2*imReplaceSquare + (1-fBlendfrac2)/2;
%                                 niNotGray = find( imScaleGraphic ~= 128/255.0 );
%                                 imReplaceSquare( niNotGray ) = imScaleGraphic( niNotGray );
%                                 %imScaleSym(        yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 ) = imReplaceSquare;
%                                 imScaleSymCovered( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1 ) = 1;
%                             end
%                             
%                             % Mark singular points (red dots) on last generation.
%                             if nGenTest == 1
%                                 imScaleSym( 2*(iii+1), 2*(jjj+1), : ) = imScaleGraphic(1,1,:);
% %                                 imScaleSym( 2*(iii+1), 2*(jjj+1), 2 ) = 0.0;
% %                                 imScaleSym( 2*(iii+1), 2*(jjj+1), 3 ) = 0.0;
%                             end
                            
                            imReplaceSquare = 0;
                            niNotGray       = 0;
                        end
                    end
                                       
                     
                    bSwappedExact = false; 
                    % If block is exact symbol swap of first symbol set:
                    if    bSwap && ~bTranslation ...    
                       && ~any(any( blockLSSwap - cellLS{nGenTest+1} )) ...
                       && ~any(any( blockRT     - cellRT{nGenTest+1} )) ... 
                       && ~any(any( blockMR     - cellMR{nGenTest+1} ))
                   
                        bSwappedExact = true;
                        if length( find( imScaleSymCovered( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1 ) == 0 ) ) > 0

                            if nGenTest >= nSymGenerations - 8
                                
                                imReplaceSquare = imScaleSym( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 );
                                imReplaceSquare = fBlendfrac2*imReplaceSquare + (1-fBlendfrac2)/2;
                                niNotGray = find( imScaleGraphicSwap ~= 128/255.0 );
                                imReplaceSquare( niNotGray ) = imScaleGraphicSwap( niNotGray );
                                imScaleSym(        yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 ) = imReplaceSquare;
                                imScaleSymCovered( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1 ) = 1;
                            end
                                                        
                            imReplaceSquare = 0;
                            niNotGray       = 0;
                        end
                    end
                    
                    bSwappedVertical = false;
                    % If block is swapped vertical mirror of first symbol set:
                    if   bSwap && ~bTranslation && ~bSwappedExact ...    
                       && ~any(any( blockLSSwap -  flipdim( cellLS{nGenTest+1}, 2 )          )) ...
                       && ~any(any( blockRT -      flipdim( cellRT{nGenTest+1}, 2 )          )) ... 
                       && ~any(any( blockMR -      flipdim( cellMR{nGenTest+1}, 2 )          )) 
                   
                        bSwappedVertical = true;
                        
                        if imScaleSymCovered( yBlock, xBlock, 5 ) == 0

                            if nGenTest >= nSymGenerations - 8
                                
                                imReplaceSquare = imScaleSym( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 );
                                imReplaceSquare = fBlendfrac2*imReplaceSquare + (1-fBlendfrac2)/2;
                                imSymOverlay = flipdim( imScaleGraphicSwap, 2 );
                                niNotGray = find( imSymOverlay ~= 128/255.0 );
                                imReplaceSquare( niNotGray ) = imSymOverlay( niNotGray );
                                imScaleSym(        yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 ) = imReplaceSquare;
                                imScaleSymCovered( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 5 ) = 1;
                            end
                                                        
                            imReplaceSquare = 0;
                            niNotGray       = 0;
                        end
                    end


%                     % If block is 90 rotation of first symbol set:
%                     if    ~any(any(      blockLS     -      rot90( cellLS{1+nOrder}, -1 )          )) ...
%                        && ~any(any( mod( blockRT, 4 )- mod( rot90( cellRT{1+nOrder}, -1 ) + 1, 4 ) )) ... 
%                        && ~any(any(      blockMR     -      rot90( cellMR{1+nOrder}, -1 )          ))
%                    
%                         if imScaleSymCovered( yBlock, xBlock, 2 ) == 0
% 
%                             if ig <= nSymGenerations - 4
%                                 
%                                 imReplaceSquare = imScaleSym( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 );
%                                 imSymOverlay = flipdim( permute( imScaleGraphic, [2 1 3] ), 2 );
%                                 niNotGray = find( imSymOverlay ~= 128/255.0 );
%                                 imReplaceSquare( niNotGray ) = imSymOverlay( niNotGray );
%                                 imScaleSym(        yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 ) = imReplaceSquare;
%                                 imScaleSymCovered( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 2 ) = 1;
%                             end
%                             
%                             if ig > nSymGenerations - 4 && ig < nSymGenerations
%                                 imReplaceSquare = imScaleSym( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 );
%                                 imSymOverlay = flipdim( permute( imScaleGraphic, [2 1 3] ), 2 );
%                                 niNotGray = find( imSymOverlay ~= 128/255.0 );
%                                 imReplaceSquare( niNotGray ) = imSymOverlay( niNotGray );
% %                                imScaleSym(        yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 ) = imReplaceSquare;
%                                 imScaleSymCovered( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 2 ) = 1;
%                             end
%                            
%                             % Mark singular points (red dots) on last generation.
%                             if ig == nSymGenerations
%                                 imScaleSym( iii, jjj, 1 ) = 1.0;
%                                 imScaleSym( iii, jjj, 2 ) = 0.0;
%                                 imScaleSym( iii, jjj, 3 ) = 0.0;
%                             end
%                             
%                             imReplaceSquare = 0;
%                             niNotGray       = 0;
%                         end
%                     end
%                     
%                     % If block is 90 rotation and symbol swap of first symbol set:
%                     if    bSwap ...    
%                        && ~any(any(      blockLSSwap     -      rot90( cellLS{1+nOrder}, -1 )          )) ...
%                        && ~any(any( mod( blockRT, 4 )    - mod( rot90( cellRT{1+nOrder}, -1 ) + 1, 4 ) )) ... 
%                        && ~any(any(      blockMR         -      rot90( cellMR{1+nOrder}, -1 )          ))
%                    
%                         if imScaleSymCovered( yBlock, xBlock, 2 ) == 0
% 
%                             if ig <= nSymGenerations - 4
%                                 
%                                 imReplaceSquare = imScaleSym( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 );
%                                 imSymOverlay = flipdim( permute( imScaleGraphicSwap, [2 1 3] ), 2 );
%                                 niNotGray = find( imSymOverlay ~= 128/255.0 );
%                                 imReplaceSquare( niNotGray ) = imSymOverlay( niNotGray );
%                                 imScaleSym(        yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 ) = imReplaceSquare;
%                                 imScaleSymCovered( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 2 ) = 1;
%                             end
%                             
%                             if ig > nSymGenerations - 4 && ig < nSymGenerations
%                                 imReplaceSquare = imScaleSym( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 );
%                                 imSymOverlay = flipdim( permute( imScaleGraphicSwap, [2 1 3] ), 2 );
%                                 niNotGray = find( imSymOverlay ~= 128/255.0 );
%                                 imReplaceSquare( niNotGray ) = imSymOverlay( niNotGray );
% %                                imScaleSym(        yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 ) = imReplaceSquare;
%                                 imScaleSymCovered( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 2 ) = 1;
%                             end
%                            
%                             % Mark singular points (red dots) on last generation.
%                             if ig == nSymGenerations
%                                 imScaleSym( iii, jjj, 1 ) = 1.0;
%                                 imScaleSym( iii, jjj, 2 ) = 0.0;
%                                 imScaleSym( iii, jjj, 3 ) = 0.0;
%                             end
%                             
%                             imReplaceSquare = 0;
%                             niNotGray       = 0;
%                         end
%                     end
%                     
%                     % If block is 180 rotation of first symbol set:
%                     if    ~any(any(      blockLS     -      rot90( cellLS{1+nOrder},  2 )          )) ...
%                        && ~any(any( mod( blockRT, 4 )- mod( rot90( cellRT{1+nOrder},  2 ) + 2, 4 ) )) ... 
%                        && ~any(any(      blockMR     -      rot90( cellMR{1+nOrder},  2 )          ))
%                    
%                         if imScaleSymCovered( yBlock, xBlock, 3 ) == 0
% 
%                             if ig <= nSymGenerations - 4
%                                 
%                                 imReplaceSquare = imScaleSym( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 );
%                                 imSymOverlay = flipdim( flipdim( imScaleGraphic, 1 ), 2 );
%                                 niNotGray = find( imSymOverlay ~= 128/255.0 );
%                                 imReplaceSquare( niNotGray ) = imSymOverlay( niNotGray );
%                                 imScaleSym(        yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 ) = imReplaceSquare;
%                                 imScaleSymCovered( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 3 ) = 1;
%                             end
%                             
%                             if ig > nSymGenerations - 4 && ig < nSymGenerations
%                                 imReplaceSquare = imScaleSym( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 );
%                                 imSymOverlay = flipdim( flipdim( imScaleGraphic, 1 ), 2 );
%                                 niNotGray = find( imSymOverlay ~= 128/255.0 );
%                                 imReplaceSquare( niNotGray ) = imSymOverlay( niNotGray );
% %                                imScaleSym(        yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 ) = imReplaceSquare;
%                                 imScaleSymCovered( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 3 ) = 1;
%                             end
%                            
%                             % Mark singular points (red dots) on last generation.
%                             if ig == nSymGenerations
%                                 imScaleSym( iii, jjj, 1 ) = 1.0;
%                                 imScaleSym( iii, jjj, 2 ) = 0.0;
%                                 imScaleSym( iii, jjj, 3 ) = 0.0;
%                             end
%                             
%                             imReplaceSquare = 0;
%                             niNotGray       = 0;
%                         end
%                     end
%                     
%                     % If block is 180 rotation and symbol swap of first symbol set:
%                     if    bSwap ...
%                        && ~any(any(      blockLSSwap     -      rot90( cellLS{1+nOrder},  2 )          )) ...
%                        && ~any(any( mod( blockRT, 4 )    - mod( rot90( cellRT{1+nOrder},  2 ) + 2, 4 ) )) ... 
%                        && ~any(any(      blockMR         -      rot90( cellMR{1+nOrder},  2 )          ))
%                    
%                         if imScaleSymCovered( yBlock, xBlock, 3 ) == 0
% 
%                             if ig <= nSymGenerations - 4
%                                 
%                                 imReplaceSquare = imScaleSym( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 );
%                                 imSymOverlay = flipdim( flipdim( imScaleGraphicSwap, 1 ), 2 );
%                                 niNotGray = find( imSymOverlay ~= 128/255.0 );
%                                 imReplaceSquare( niNotGray ) = imSymOverlay( niNotGray );
%                                 imScaleSym(        yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 ) = imReplaceSquare;
%                                 imScaleSymCovered( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 3 ) = 1;
%                             end
%                             
%                             if ig > nSymGenerations - 4 && ig < nSymGenerations
%                                 imReplaceSquare = imScaleSym( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 );
%                                 imSymOverlay = flipdim( flipdim( imScaleGraphicSwap, 1 ), 2 );
%                                 niNotGray = find( imSymOverlay ~= 128/255.0 );
%                                 imReplaceSquare( niNotGray ) = imSymOverlay( niNotGray );
% %                                imScaleSym(        yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 ) = imReplaceSquare;
%                                 imScaleSymCovered( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 3 ) = 1;
%                             end
%                            
%                             % Mark singular points (red dots) on last generation.
%                             if ig == nSymGenerations
%                                 imScaleSym( iii, jjj, 1 ) = 1.0;
%                                 imScaleSym( iii, jjj, 2 ) = 0.0;
%                                 imScaleSym( iii, jjj, 3 ) = 0.0;
%                             end
%                             
%                             imReplaceSquare = 0;
%                             niNotGray       = 0;
%                         end
%                     end
%                     
%                     % If block is 270 rotation of first symbol set:
%                     if    ~any(any(      blockLS     -      rot90( cellLS{1+nOrder}, 1 )          )) ...
%                        && ~any(any( mod( blockRT, 4 )- mod( rot90( cellRT{1+nOrder}, 1 ) + 3, 4 ) )) ... 
%                        && ~any(any(      blockMR     -      rot90( cellMR{1+nOrder}, 1 )          ))
%                    
%                          if imScaleSymCovered( yBlock, xBlock, 4 ) == 0
% 
%                             if ig <= nSymGenerations - 4
%                                 
%                                 imReplaceSquare = imScaleSym( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 );
%                                 imSymOverlay = flipdim( permute( imScaleGraphic, [2 1 3] ), 1 );
%                                 niNotGray = find( imSymOverlay ~= 128/255.0 );
%                                 imReplaceSquare( niNotGray ) = imSymOverlay( niNotGray );
%                                 imScaleSym(        yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 ) = imReplaceSquare;
%                                 imScaleSymCovered( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 4 ) = 1;
%                             end
%                             
%                             if ig > nSymGenerations - 4 && ig < nSymGenerations
%                                 imReplaceSquare = imScaleSym( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 );
%                                 imSymOverlay = flipdim( permute( imScaleGraphic, [2 1 3] ), 1 );
%                                 niNotGray = find( imSymOverlay ~= 128/255.0 );
%                                 imReplaceSquare( niNotGray ) = imSymOverlay( niNotGray );
%                                 imScaleSymCovered( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 4 ) = 1;
%                             end
%                           
%                             % Mark singular points (red dots) on last generation.
%                             if ig == nSymGenerations
%                                 imScaleSym( iii, jjj, 1 ) = 1.0;
%                                 imScaleSym( iii, jjj, 2 ) = 0.0;
%                                 imScaleSym( iii, jjj, 3 ) = 0.0;
%                             end
%                             
%                             imReplaceSquare = 0;
%                             niNotGray       = 0;
%                         end
%                     end
%                     
%                     % If block is 270 rotation and symbol swap of first symbol set:
%                     if    bSwap ...
%                        && ~any(any(      blockLSSwap     -      rot90( cellLS{1+nOrder}, 1 )          )) ...
%                        && ~any(any( mod( blockRT, 4 )    - mod( rot90( cellRT{1+nOrder}, 1 ) + 3, 4 ) )) ... 
%                        && ~any(any(      blockMR         -      rot90( cellMR{1+nOrder}, 1 )          ))
%                    
%                          if imScaleSymCovered( yBlock, xBlock, 4 ) == 0
% 
%                             if ig <= nSymGenerations - 4
%                                 
%                                 imReplaceSquare = imScaleSym( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 );
%                                 imSymOverlay = flipdim( permute( imScaleGraphicSwap, [2 1 3] ), 1 );
%                                 niNotGray = find( imSymOverlay ~= 128/255.0 );
%                                 imReplaceSquare( niNotGray ) = imSymOverlay( niNotGray );
%                                 imScaleSym(        yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 ) = imReplaceSquare;
%                                 imScaleSymCovered( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 4 ) = 1;
%                             end
%                             
%                             if ig > nSymGenerations - 4 && ig < nSymGenerations
%                                 imReplaceSquare = imScaleSym( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 );
%                                 imSymOverlay = flipdim( permute( imScaleGraphicSwap, [2 1 3] ), 1 );
%                                 niNotGray = find( imSymOverlay ~= 128/255.0 );
%                                 imReplaceSquare( niNotGray ) = imSymOverlay( niNotGray );
%                                 imScaleSymCovered( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 4 ) = 1;
%                             end
%                           
%                             % Mark singular points (red dots) on last generation.
%                             if ig == nSymGenerations
%                                 imScaleSym( iii, jjj, 1 ) = 1.0;
%                                 imScaleSym( iii, jjj, 2 ) = 0.0;
%                                 imScaleSym( iii, jjj, 3 ) = 0.0;
%                             end
%                             
%                             imReplaceSquare = 0;
%                             niNotGray       = 0;
%                         end
%                     end
                     
                    % If block is vertical mirror of first symbol set:
                    if     ~bTranslation && ~bSwappedExact && ~bSwappedVertical ...
                       &&  ~any(any( blockLS -      flipdim( cellLS{nGenTest+1}, 2 )          )) ...
                       &&  ~any(any( blockRT -      flipdim( cellRT{nGenTest+1}, 2 )          )) ... 
                       &&  ~any(any( blockMR -      flipdim( cellMR{nGenTest+1}, 2 )          )) ... % && ~any(any( blockMR - mod( flipdim( cellMR{1+1}, 2 ) + 1, 2 ) ))
                       &&  ~any(any( blockLS -      flipdim( blockLS, 2 )                     ))     % not mirror of self  
                   
                        if imScaleSymCovered( yBlock, xBlock, 5 ) == 0

                            if ig <= nSymGenerations - 9
                                
                                imReplaceSquare = imScaleSym( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 );
                                imReplaceSquare = fBlendfrac2*imReplaceSquare + (1-fBlendfrac2)/2;
                                imSymOverlay = flipdim( imScaleGraphic, 2 );
                                niNotGray = find( imSymOverlay ~= 128/255.0 );
                                imReplaceSquare( niNotGray ) = imSymOverlay( niNotGray );
                                imScaleSym(        yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 ) = imReplaceSquare;
                                imScaleSymCovered( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 5 ) = 1;
                            end
                            
%                             if ig > nSymGenerations - 4 && ig < nSymGenerations
%                                 imReplaceSquare = imScaleSym( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 );
%                                 imSymOverlay = flipdim( imScaleGraphic, 2 );
%                                 niNotGray = find( imSymOverlay ~= 128/255.0 );
%                                 imReplaceSquare( niNotGray ) = imSymOverlay( niNotGray );
% %                                imScaleSym(        yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 1:3 ) = imReplaceSquare;
%                                 imScaleSymCovered( yBlock : yBlock + hBlock, xBlock : xBlock + wBlock, 5 ) = 1;
%                             end
%                            
%                             % Mark singular points (red dots) on last generation.
%                             if ig == nSymGenerations
%                                 imScaleSym( iii, jjj, 1 ) = 1.0;
%                                 imScaleSym( iii, jjj, 2 ) = 0.0;
%                                 imScaleSym( iii, jjj, 3 ) = 0.0;
%                             end
                            
                            imReplaceSquare = 0;
                            niNotGray       = 0;
                        end
                    end
                    
                    % If block is 90 rotation then vertical mirror of first symbol set:
                    if 0%        blockLS      ==      flipdim( rot90( arrayInit,      1 ),     2 )          ...
%                        && mod( blockRT, 4 ) == mod( flipdim( rot90( arrayInitRotes, 1 ) + 1, 2 ),     4 ) ... 
%                        &&      blockMR      == mod( flipdim( rot90( arrayInitMirrs, 1 ),     2 ) + 1, 2 )

                        if imScaleSymCovered( 2^(nSymGenerations-igr)*(iii-1)+1, 2^(nSymGenerations-igr)*(jjj-1)+1, 1:3 ) == 0

                            if igr <= nSymGenerations - 4
                                imReplaceSquare = imScaleSym( 2^(nSymGenerations-igr)*(iii-1)+1 : 2^(nSymGenerations-igr)*(iii-1)+1 + 2^(nSymGenerations-igr) - 1, 2^(nSymGenerations-igr)*(jjj-1)+1 : 2^(nSymGenerations-igr)*(jjj-1)+1 + 2^(nSymGenerations-igr) - 1, 1:3 );
                                imSymOverlay = permute( impyScaleSymGraphic{nSymGenerations+1-igr}, [2 1 3] );
                                niNotGray = find( imSymOverlay ~= 128/255.0 );
                                imReplaceSquare( niNotGray ) = imSymOverlay( niNotGray );
                                imScaleSym( 2^(nSymGenerations-igr)*(iii-1)+1 : 2^(nSymGenerations-igr)*(iii-1)+1 + 2^(nSymGenerations-igr) - 1, 2^(nSymGenerations-igr)*(jjj-1)+1 : 2^(nSymGenerations-igr)*(jjj-1)+1 + 2^(nSymGenerations-igr) - 1, 1:3 ) = imReplaceSquare;
                                imScaleSymCovered( 2^(nSymGenerations-igr)*(iii-1)+1 : 2^(nSymGenerations-igr)*(iii-1)+1 + 2^(nSymGenerations-igr) - 1, 2^(nSymGenerations-igr)*(jjj-1)+1 : 2^(nSymGenerations-igr)*(jjj-1)+1 + 2^(nSymGenerations-igr) - 1, 1:3 ) = 1;
                            end
                            
                            if igr > nSymGenerations - 4 && igr < nSymGenerations
                                imReplaceSquare = imScaleSym( 2^(nSymGenerations-igr)*(iii-1)+1 : 2^(nSymGenerations-igr)*(iii-1)+1 + 2^(nSymGenerations-igr) - 1, 2^(nSymGenerations-igr)*(jjj-1)+1 : 2^(nSymGenerations-igr)*(jjj-1)+1 + 2^(nSymGenerations-igr) - 1, 1:3 );
                                imSymOverlay = permute( impyScaleSymGraphic{nSymGenerations+1-igr}, [2 1 3] );
                                niNotGray = find( imSymOverlay ~= 128/255.0 );
                                imReplaceSquare( niNotGray ) = imSymOverlay( niNotGray );
%                                imScaleSym( 2^(nSymGenerations-igr)*(iii-1)+1 : 2^(nSymGenerations-igr)*(iii-1)+1 + 2^(nSymGenerations-igr) - 1, 2^(nSymGenerations-igr)*(jjj-1)+1 : 2^(nSymGenerations-igr)*(jjj-1)+1 + 2^(nSymGenerations-igr) - 1, 1:3 ) = imReplaceSquare;
                                imScaleSymCovered( 2^(nSymGenerations-igr)*(iii-1)+1 : 2^(nSymGenerations-igr)*(iii-1)+1 + 2^(nSymGenerations-igr) - 1, 2^(nSymGenerations-igr)*(jjj-1)+1 : 2^(nSymGenerations-igr)*(jjj-1)+1 + 2^(nSymGenerations-igr) - 1, 1:3 ) = 1;
                            end
                           
                            % Mark singular points (red dots) on last generation.
                            if igr == nSymGenerations
                                imScaleSym( iii, jjj, 1 ) = 1.0;
                                imScaleSym( iii, jjj, 2 ) = 0.0;
                                imScaleSym( iii, jjj, 3 ) = 0.0;
                            end
                            
                            imReplaceSquare = 0;
                            niNotGray       = 0;
                        end
                    end
                    
                    % If block is 180 rotation then vertical mirror of first symbol set:
                    if 0%        blockLS      ==      flipdim( rot90( arrayInit,      2 ),     2 )          ...
%                        && mod( blockRT, 4 ) == mod( flipdim( rot90( arrayInitRotes, 2 ) + 2, 2 ),     4 ) ... 
%                        &&      blockMR      == mod( flipdim( rot90( arrayInitMirrs, 2 ),     2 ) + 1, 2 )

                        if imScaleSymCovered( 2^(nSymGenerations-igr)*(iii-1)+1, 2^(nSymGenerations-igr)*(jjj-1)+1, 1:3 ) == 0

                            if igr <= nSymGenerations - 4
                                
                                imReplaceSquare = imScaleSym( 2^(nSymGenerations-igr)*(iii-1)+1 : 2^(nSymGenerations-igr)*(iii-1)+1 + 2^(nSymGenerations-igr) - 1, 2^(nSymGenerations-igr)*(jjj-1)+1 : 2^(nSymGenerations-igr)*(jjj-1)+1 + 2^(nSymGenerations-igr) - 1, 1:3 );
                                imSymOverlay = flipdim( impyScaleSymGraphic{nSymGenerations+1-igr}, 1 );
                                niNotGray = find( imSymOverlay ~= 128/255.0 );
                                imReplaceSquare( niNotGray ) = imSymOverlay( niNotGray );
                                imScaleSym( 2^(nSymGenerations-igr)*(iii-1)+1 : 2^(nSymGenerations-igr)*(iii-1)+1 + 2^(nSymGenerations-igr) - 1, 2^(nSymGenerations-igr)*(jjj-1)+1 : 2^(nSymGenerations-igr)*(jjj-1)+1 + 2^(nSymGenerations-igr) - 1, 1:3 ) = imReplaceSquare;
                                imScaleSymCovered( 2^(nSymGenerations-igr)*(iii-1)+1 : 2^(nSymGenerations-igr)*(iii-1)+1 + 2^(nSymGenerations-igr) - 1, 2^(nSymGenerations-igr)*(jjj-1)+1 : 2^(nSymGenerations-igr)*(jjj-1)+1 + 2^(nSymGenerations-igr) - 1, 1:3 ) = 1;
                            end
                            
                            if igr > nSymGenerations - 4 && igr < nSymGenerations
                                imReplaceSquare = imScaleSym( 2^(nSymGenerations-igr)*(iii-1)+1 : 2^(nSymGenerations-igr)*(iii-1)+1 + 2^(nSymGenerations-igr) - 1, 2^(nSymGenerations-igr)*(jjj-1)+1 : 2^(nSymGenerations-igr)*(jjj-1)+1 + 2^(nSymGenerations-igr) - 1, 1:3 );
                                imSymOverlay = flipdim( impyScaleSymGraphic{nSymGenerations+1-igr}, 1 );
                                niNotGray = find( imSymOverlay ~= 128/255.0 );
                                imReplaceSquare( niNotGray ) = imSymOverlay( niNotGray );
%                                imScaleSym( 2^(nSymGenerations-igr)*(iii-1)+1 : 2^(nSymGenerations-igr)*(iii-1)+1 + 2^(nSymGenerations-igr) - 1, 2^(nSymGenerations-igr)*(jjj-1)+1 : 2^(nSymGenerations-igr)*(jjj-1)+1 + 2^(nSymGenerations-igr) - 1, 1:3 ) = imReplaceSquare;
                                imScaleSymCovered( 2^(nSymGenerations-igr)*(iii-1)+1 : 2^(nSymGenerations-igr)*(iii-1)+1 + 2^(nSymGenerations-igr) - 1, 2^(nSymGenerations-igr)*(jjj-1)+1 : 2^(nSymGenerations-igr)*(jjj-1)+1 + 2^(nSymGenerations-igr) - 1, 1:3 ) = 1;
                            end
                           
                            % Mark singular points (red dots) on last generation.
                            if igr == nSymGenerations
                                imScaleSym( iii, jjj, 1 ) = 1.0;
                                imScaleSym( iii, jjj, 2 ) = 0.0;
                                imScaleSym( iii, jjj, 3 ) = 0.0;
                            end
                            
                            imReplaceSquare = 0;
                            niNotGray       = 0;
                        end
                    end

                    % If block is 270 rotation then vertical mirror of first symbol set:
                    if 0%        blockLS      ==      flipdim( rot90( arrayInit,      3 ),     2 )          ...
%                        && mod( blockRT, 4 ) == mod( flipdim( rot90( arrayInitRotes, 3 ) + 3, 2 ),     4 ) ... 
%                        &&      blockMR      == mod( flipdim( rot90( arrayInitMirrs, 3 ),     2 ) + 1, 2 )

                        if imScaleSymCovered( 2^(nSymGenerations-igr)*(iii-1)+1, 2^(nSymGenerations-igr)*(jjj-1)+1, 1:3 ) == 0

                            if igr <= nSymGenerations - 4
                                
                                imReplaceSquare = imScaleSym( 2^(nSymGenerations-igr)*(iii-1)+1 : 2^(nSymGenerations-igr)*(iii-1)+1 + 2^(nSymGenerations-igr) - 1, 2^(nSymGenerations-igr)*(jjj-1)+1 : 2^(nSymGenerations-igr)*(jjj-1)+1 + 2^(nSymGenerations-igr) - 1, 1:3 );
                                imSymOverlay = flipdim( flipdim( permute( impyScaleSymGraphic{nSymGenerations+1-igr}, [2 1 3] ), 1 ), 2 );
                                niNotGray = find( imSymOverlay ~= 128/255.0 );
                                imReplaceSquare( niNotGray ) = imSymOverlay( niNotGray );
                                imScaleSym( 2^(nSymGenerations-igr)*(iii-1)+1 : 2^(nSymGenerations-igr)*(iii-1)+1 + 2^(nSymGenerations-igr) - 1, 2^(nSymGenerations-igr)*(jjj-1)+1 : 2^(nSymGenerations-igr)*(jjj-1)+1 + 2^(nSymGenerations-igr) - 1, 1:3 ) = imReplaceSquare;
                                imScaleSymCovered( 2^(nSymGenerations-igr)*(iii-1)+1 : 2^(nSymGenerations-igr)*(iii-1)+1 + 2^(nSymGenerations-igr) - 1, 2^(nSymGenerations-igr)*(jjj-1)+1 : 2^(nSymGenerations-igr)*(jjj-1)+1 + 2^(nSymGenerations-igr) - 1, 1:3 ) = 1;
                            end
                            
                            if igr > nSymGenerations - 4 && igr < nSymGenerations
                                imReplaceSquare = imScaleSym( 2^(nSymGenerations-igr)*(iii-1)+1 : 2^(nSymGenerations-igr)*(iii-1)+1 + 2^(nSymGenerations-igr) - 1, 2^(nSymGenerations-igr)*(jjj-1)+1 : 2^(nSymGenerations-igr)*(jjj-1)+1 + 2^(nSymGenerations-igr) - 1, 1:3 );
                                imSymOverlay = flipdim( flipdim( permute( impyScaleSymGraphic{nSymGenerations+1-igr}, [2 1 3] ), 1 ), 2 );
                                niNotGray = find( imSymOverlay ~= 128/255.0 );
                                imReplaceSquare( niNotGray ) = imSymOverlay( niNotGray );
%                                imScaleSym( 2^(nSymGenerations-igr)*(iii-1)+1 : 2^(nSymGenerations-igr)*(iii-1)+1 + 2^(nSymGenerations-igr) - 1, 2^(nSymGenerations-igr)*(jjj-1)+1 : 2^(nSymGenerations-igr)*(jjj-1)+1 + 2^(nSymGenerations-igr) - 1, 1:3 ) = imReplaceSquare;
                                imScaleSymCovered( 2^(nSymGenerations-igr)*(iii-1)+1 : 2^(nSymGenerations-igr)*(iii-1)+1 + 2^(nSymGenerations-igr) - 1, 2^(nSymGenerations-igr)*(jjj-1)+1 : 2^(nSymGenerations-igr)*(jjj-1)+1 + 2^(nSymGenerations-igr) - 1, 1:3 ) = 1;
                            end
                           
                            % Mark singular points (red dots) on last generation.
                            if igr == nSymGenerations
                                imScaleSym( iii, jjj, 1 ) = 1.0;
                                imScaleSym( iii, jjj, 2 ) = 0.0;
                                imScaleSym( iii, jjj, 3 ) = 0.0;
                            end
                            
                            imReplaceSquare = 0;
                            niNotGray       = 0;
                        end
                    end

                end   % iii
                end   % jjj
                
                % Mark as all covered any area covered by any one scale symmetry.
%                 if nOrder == 1
                    
                    for iiii = 1 : szimScaleSym(1)
                    for jjjj = 1 : szimScaleSym(2)

                        if any( imScaleSymCovered( iiii, jjjj, : ) )
                            imScaleSymCovered( iiii, jjjj, : ) = 1;
                        end
                    end
                    end
%                 end
                
            end % Mark self-symmetries.
            %end % each norder
%         end  % For each generation up to half the size of the symmetry diagram.
    end
    
    % Downsample if required to match small out put.
    if nGenerations + log2(nOversampleOut) < 6
        
        dy = size(imScaleSym, 1)/( size(imScaleSym, 1)/2 - 1 );
        dx = size(imScaleSym, 2)/( size(imScaleSym, 2)/2 - 1 );
        yi = 1 -.5 + dy/2 : dy : size(imScaleSym, 1) +.5 - dy/2;
        yi = yi';
        xi = 1 -.5 + dx/2 : dx : size(imScaleSym, 2) +.5 - dx/2;
            
        imScaleSymTemp(:,:, 1) = interp2_fa( imScaleSym(:,:, 1), xi, yi, 'linear' );
        imScaleSymTemp(:,:, 2) = interp2_fa( imScaleSym(:,:, 2), xi, yi, 'linear' );
        imScaleSymTemp(:,:, 3) = interp2_fa( imScaleSym(:,:, 3), xi, yi, 'linear' );
        clear imScaleSym;
        imScaleSym = imScaleSymTemp;
        clear imScaleSymTemp;
    end

end

clear vctLS0;
clear vctRT0;
clear vctMR0;
clear vctLS;
clear vctRT;
clear vctMR;


% Construct scaled sequence of each generation ( for the "by generation" graphic):
if szRules(1) > 1 && szRules(2) > 1
    
%    nByGenerations = floor( log10( nByGenWidthMax/(nByGenOversample*szInit(2)*szRules(2)) )/log10( szInit(2)*szRules(2) ) );
%    nByGenerations = floor( log10( nByGenWidthMax/(nByGenOversample*szInit(2)*szRules(2)) )/log10( nByGenerations ) );
    
    nGen = 2;
    nGenW = 0;
    while nGenW < nByGenWidthMax
        nGen = nGen + 1;
        nGenW = nGen*nByGenOversample*szInit(2)*(szRules(2)^nGen);
    end
    nByGenerations = nGen - 1;
%     while nByGenerations*(nByGenOversample + 1)*szInit(2)*(szRules(2)^nByGenerations) < nByGenWidthMax
%         nByGenOversample = nByGenOversample + 1;
%     end
%     nByGenOversample = nByGenOversample - 1;
    
    if bRandom == 0 %& nByGenerations > 1 & nByGenerations >= 2    %strMotifFilename == 'animate'

        vctLS0 = arrayInit;
        szLS = size(vctLS0);
        vctRT0 = arrayInitRotes;
        vctMR0 = arrayInitMirrs;
        vctLS( 1:szLS(1)*szRules(1), 1:szLS(2)*szRules(2) ) = 0;
        vctRT( 1:szLS(1)*szRules(1), 1:szLS(2)*szRules(2) ) = 0;
        vctMR( 1:szLS(1)*szRules(1), 1:szLS(2)*szRules(2) ) = 0;

        nBGSpacing = 6;
%         2*nBGSpacing + 1 + nByGenOversample*szInit(1)*szRules(1)^nByGenerations
%         ( nByGenerations + 1 )*nBGSpacing + 1 + nByGenerations*nByGenOversample*szInit(2)*szRules(2)^nByGenerations
        nimByGenerationGraphic( 1:2*nBGSpacing + 0 + nByGenOversample*szInit(1)*szRules(1)^nByGenerations, 1:( nByGenerations + 1 )*nBGSpacing + nByGenerations*nByGenOversample*szInit(2)*szRules(2)^nByGenerations, 1:3 ) = (10/13)*(nSymbols-1);

        
        rotarrayRules       = 0;
        rotarrayRotes       = 0;
        rotarrayMirrs       = 0;
        fliplrrotarrayRotes = 0;
        fliplrrotarrayMirrs = 0;
        for ils = 1 : nSymbols
        for irt = 1 : 4
           rotarrayRules( 1:szRules(1), 1:szRules(2), ils, irt ) = rot90( arrayRules( 1:szRules(1), 1:szRules(2), ils ), -(irt-1) );
%             rotarrayRules( 1:szRules(1), 1:szRules(2), ils, irt ) = arrayRules( 1:szRules(1), 1:szRules(2) );
        end
        end

        for ils = 1 : nSymbols
        for irt = 1 : 4
           rotarrayRotes( 1:szRules(1), 1:szRules(2), ils, irt ) = rot90( arrayRotes( 1:szRules(1), 1:szRules(2), ils ), -(irt-1) );
           rotarrayMirrs( 1:szRules(1), 1:szRules(2), ils, irt ) = rot90( arrayMirrs( 1:szRules(1), 1:szRules(2), ils ), -(irt-1) );
           fliplrrotarrayRotes( 1:szRules(1), 1:szRules(2), ils, irt ) = fliplr( rotarrayRotes( 1:szRules(1), 1:szRules(2), ils, irt ) );
           fliplrrotarrayMirrs( 1:szRules(1), 1:szRules(2), ils, irt ) = fliplr( rotarrayMirrs( 1:szRules(1), 1:szRules(2), ils, irt ) );
%            rotarrayRotes( 1:szRules(1), 1:szRules(2), ils, irt ) = arrayRotes( 1:szRules(1), 1:szRules(2), ils );
%            rotarrayMirrs( 1:szRules(1), 1:szRules(2), ils, irt ) = arrayMirrs( 1:szRules(1), 1:szRules(2), ils );
%            fliplrrotarrayRotes( 1:szRules(1), 1:szRules(2), ils, irt ) = rotarrayRotes( 1:szRules(1), 1:szRules(2) );
%            fliplrrotarrayMirrs( 1:szRules(1), 1:szRules(2), ils, irt ) = rotarrayMirrs( 1:szRules(1), 1:szRules(2) );
        end
        end

        for ii = 1 : nByGenerations

            % Construct the next generation of the pattern.
            for j = 1 : szLS(1)
            for k = 1 : szLS(2)

                % Replace each symbol in the current symbol array (vctLS0, zero based values) 
                % with the current rule block (arrayRules(:,:,vctLS0 + 1), one based index),
                % and rotate by the current rotation (-vctRT0( j, k )).
                vctLS( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) = rotarrayRules( 1:szRules(1), 1:szRules(2), vctLS0( j, k ) + 1 , vctRT0( j, k ) + 1 );
                % If the current mirror is 1, mirror (vertically) the block.
                if vctMR0( j, k ) == 1
                    vctLS( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) = fliplr( rotarrayRules( 1:szRules(1), 1:szRules(2), vctLS0( j, k ) + 1 , vctRT0( j, k ) + 1 ) );
                end

                % Mirror and/or rotate the mirror and rotation array blocks.
                if vctMR0( j, k ) == 1
                    vctRT( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) = vctRT0( j, k ) - fliplrrotarrayRotes( 1:szRules(1), 1:szRules(2), vctLS0( j, k ) + 1 , vctRT0( j, k ) + 1 );
                    vctMR( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) = vctMR0( j, k ) - fliplrrotarrayMirrs( 1:szRules(1), 1:szRules(2), vctLS0( j, k ) + 1 , vctRT0( j, k ) + 1 );
                else
                    vctRT( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) = vctRT0( j, k ) + rotarrayRotes( 1:szRules(1), 1:szRules(2), vctLS0( j, k ) + 1 , vctRT0( j, k ) + 1 );
                    vctMR( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) = vctMR0( j, k ) + rotarrayMirrs( 1:szRules(1), 1:szRules(2), vctLS0( j, k ) + 1 , vctRT0( j, k ) + 1 );
                end

            end
            end

            vctRT = mod( vctRT, 4 );
            vctMR = mod( vctMR, 2 );
                
            vctLS0 = vctLS;
            vctRT0 = vctRT;
            vctMR0 = vctMR;
            szLS = size(vctLS0);
            
            % Add dark border/background.
            nimByGenerationGraphic( nBGSpacing : nBGSpacing + nByGenOversample*szInit(1)*( szRules(1)^nByGenerations ) + 1, nBGSpacing + (ii-1)*nBGSpacing + nByGenOversample*szInit(2)*(ii - 1)*( szRules(2)^nByGenerations ) : nBGSpacing + nByGenOversample*szInit(2)*( szRules(2)^nByGenerations ) + (ii-1)*nBGSpacing + nByGenOversample*szInit(2)*(ii-1)*( szRules(2)^nByGenerations ) + 1, 1:3 ) = (5/13)*(nSymbols-1);

            % Upsample each generation to match the size of the last generation.
            nUpsampleX = nByGenOversample*szRules(2)^( nByGenerations - ii );
            nUpsampleY = nByGenOversample*szRules(1)^( nByGenerations - ii );
            for i = 1 : szLS(1)
            for j = 1 : szLS(2)
            for k = 1 : nUpsampleY
            for l = 1 : nUpsampleX
                nimByGenerationGraphic( nBGSpacing + 0 + nUpsampleY*(i-1)+ k, nBGSpacing + 0 + nUpsampleX*(j-1) + l + (ii-1)*nBGSpacing + nByGenOversample*szInit(2)*(ii-1)*( szRules(2)^nByGenerations ), 1 ) = vctLS( i, j ); 
                nimByGenerationGraphic( nBGSpacing + 0 + nUpsampleY*(i-1)+ k, nBGSpacing + 0 + nUpsampleX*(j-1) + l + (ii-1)*nBGSpacing + nByGenOversample*szInit(2)*(ii-1)*( szRules(2)^nByGenerations ), 2 ) = vctLS( i, j ); 
                nimByGenerationGraphic( nBGSpacing + 0 + nUpsampleY*(i-1)+ k, nBGSpacing + 0 + nUpsampleX*(j-1) + l + (ii-1)*nBGSpacing + nByGenOversample*szInit(2)*(ii-1)*( szRules(2)^nByGenerations ), 3 ) = vctLS( i, j ); 
            end
            end
            end
            end

        end

        nimByGenerationGraphic = nimByGenerationGraphic/max(max(max(nimByGenerationGraphic)));

        if bShowResults

            if szRules < 3
                figure; imshow( imScaleSym )
            end

            figure; imshow( nimByGenerationGraphic );        
        end

        % Make an alpha (transparency) map, 1 at all non-backgroud pixels,
        % and 0 (transparent) at all background voxels, where background is
        % determined by the color of the upper-left pixe..
        szBG = size( nimByGenerationGraphic );
        imByGenerationAlpha( 1:szBG(1), 1:szBG(2) ) = 1;
        imByGenerationAlpha( find( nimByGenerationGraphic(:,:,1) == nimByGenerationGraphic(1,1) & nimByGenerationGraphic(:,:,2) == nimByGenerationGraphic(1,1) & nimByGenerationGraphic(:,:,3) == nimByGenerationGraphic(1,1) ) ) = 0;
    end
end

clear vctLS0;
clear vctLS;
   
if bRandom == 0
    
    % Construct 2D sequence:
    vctLS0 = arrayInit;
    szLS = size(vctLS0);        
    vctRT0 = arrayInitRotes;
    vctMR0 = arrayInitMirrs;
    vctLS( 1:szLS(1)*szRules(1), 1:szLS(2)*szRules(2) ) = 0;
    vctRT( 1:szLS(1)*szRules(1), 1:szLS(2)*szRules(2) ) = 0;
    vctMR( 1:szLS(1)*szRules(1), 1:szLS(2)*szRules(2) ) = 0;

    % rotarrayRules = 0;
    % for ils = 1 : nSymbols
    % for irt = 1 : 4
    %     rotarrayRules( 1:szRules(1), 1:szRules(2), ils, irt ) = rot90( arrayRules( 1:szRules(1), 1:szRules(2), ils ), -irt );
    % end
    % end

    rotarrayRules       = 0;
    fliplrrotarrayRules = 0;
    rotarrayRotes       = 0;
    rotarrayMirrs       = 0;
    fliplrrotarrayRotes = 0;
    fliplrrotarrayMirrs = 0;

    for ils = 1 : nSymbols
    for irt = 1 : 4
        rotarrayRules( 1:szRules(1), 1:szRules(2), ils, irt ) = rot90( arrayRules( 1:szRules(1), 1:szRules(2), ils ), -(irt-1) );
        fliplrrotarrayRules( 1:szRules(1), 1:szRules(2), ils, irt ) = flipdim( rotarrayRules( 1:szRules(1), 1:szRules(2), ils, irt ), 2 );
%         rotarrayRules( 1:szRules(1), 1:szRules(2), ils, irt ) = arrayRules( 1:szRules(1), 1:szRules(2) );
%         fliplrrotarrayRules( 1:szRules(1), 1:szRules(2), ils, irt ) = rotarrayRules( 1:szRules(1), 1:szRules(2) );
    end
    end

    for ils = 1 : nSymbols
    for irt = 1 : 4
        rotarrayRotes( 1:szRules(1), 1:szRules(2), ils, irt ) = rot90( arrayRotes( 1:szRules(1), 1:szRules(2), ils ), -(irt-1) );
        rotarrayMirrs( 1:szRules(1), 1:szRules(2), ils, irt ) = rot90( arrayMirrs( 1:szRules(1), 1:szRules(2), ils ), -(irt-1) );
        fliplrrotarrayRotes( 1:szRules(1), 1:szRules(2), ils, irt ) = fliplr( rotarrayRotes( 1:szRules(1), 1:szRules(2), ils, irt ) );
        fliplrrotarrayMirrs( 1:szRules(1), 1:szRules(2), ils, irt ) = fliplr( rotarrayMirrs( 1:szRules(1), 1:szRules(2), ils, irt ) );
%         rotarrayRotes( 1:szRules(1), 1:szRules(2), ils, irt ) = arrayRotes( 1:szRules(1), 1:szRules(2) );
%         rotarrayMirrs( 1:szRules(1), 1:szRules(2), ils, irt ) = arrayMirrs( 1:szRules(1), 1:szRules(2) );
%         fliplrrotarrayRotes( 1:szRules(1), 1:szRules(2), ils, irt ) = rotarrayRotes( 1:szRules(1), 1:szRules(2) );
%         fliplrrotarrayMirrs( 1:szRules(1), 1:szRules(2), ils, irt ) = rotarrayMirrs( 1:szRules(1), 1:szRules(2) );
    end
    end

    for i = 1 : nGenerations

        % Construct the next generation of the pattern.
        for j = 1 : szLS(1)
        for k = 1 : szLS(2)

            % Replace each symbol in the current symbol array (vctLS0, zero based values) 
            % with the current rule block (arrayRules(:,:,vctLS0 + 1), one based index),
            % and rotate by the current rotation (-vctRT0( j, k )).
            vctLS( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) = rotarrayRules( 1:szRules(1), 1:szRules(2), vctLS0( j, k ) + 1 , vctRT0( j, k ) + 1 );
            % If the current mirror is 1, mirror (vertically) the block.
            if vctMR0( j, k ) == 1
                vctLS( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) = fliplrrotarrayRules( 1:szRules(1), 1:szRules(2), vctLS0( j, k ) + 1 , vctRT0( j, k ) + 1 );
            end

            % Mirror and/or rotate the mirror and rotation array blocks.
            if vctMR0( j, k ) == 1
                vctRT( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) = vctRT0( j, k ) - fliplrrotarrayRotes( 1:szRules(1), 1:szRules(2), vctLS0( j, k ) + 1 , vctRT0( j, k ) + 1 );
                vctMR( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) = vctMR0( j, k ) - fliplrrotarrayMirrs( 1:szRules(1), 1:szRules(2), vctLS0( j, k ) + 1 , vctRT0( j, k ) + 1 );
            else
                vctRT( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) = vctRT0( j, k ) + rotarrayRotes( 1:szRules(1), 1:szRules(2), vctLS0( j, k ) + 1 , vctRT0( j, k ) + 1 );
                vctMR( (j-1)*szRules(1) + 1 : j*szRules(1), (k-1)*szRules(2) + 1 : k*szRules(2) ) = vctMR0( j, k ) + rotarrayMirrs( 1:szRules(1), 1:szRules(2), vctLS0( j, k ) + 1 , vctRT0( j, k ) + 1 );
            end

        end
        end

        vctRT = mod( vctRT, 4 );
        vctMR = mod( vctMR, 2 );

        vctLS0 = vctLS;
        vctRT0 = vctRT;
        vctMR0 = vctMR;
        szLS = size(vctLS0);
        
        % Add scaled average (sum) image.
        if bDrawAverage
            
%             nimN = double( vctLS0 );
            nimN = uint16( vctLS0 );
    
            % Power of 2 resampling to the maximum generation size.
            for ir = 1 : nGenerations - i

                imNtemp = interp2( nimN, 1, 'nearest' );

                szTemp = size( imNtemp );
                % Duplicate first column and row, because they are
                % clipped by the nearest neighbor interpolation.
                nimN( 2:1+szTemp(1), 2:1+szTemp(2) ) = imNtemp;
                nimN( 1, 2:1+szTemp(2) ) = imNtemp( 1, : );
                nimN( 2:1+szTemp(1), 1 ) = imNtemp( :, 1 );
                nimN( 1, 1 ) = imNtemp( 1, 1 );

                %size( nimN )
            end
            
            if nDecimate > 1

                nimNTemp = downsample( nimN, nDecimate, nOffset );
                nimNTemp = nimNTemp';
                nimNTemp = downsample( nimNTemp, nDecimate, nOffset );
                clear nimN;
                nimN = nimNTemp';
            end
% size(imAvg)
% size(nimN)
            imAvg = imAvg + nimN;
            clear nimN;
        end


    end % for each generation

%    szLS0 = size(vctLS0);
end % not random

if bRandom
%     for i = 1:szLS0(1)
%     for j = 1:szLS0(2)
        rand( 'state', 1 );
%        vctLS0( 1:szLS0(1), 1:szLS0(2) ) = floor(  nSymbols * rand( szLS0(1) ) );
        vctLS0( 1:2^nGenerations, 1:2^nGenerations ) = floor(  nSymbols * rand( 2^nGenerations ) );
%     end
%     end
end


if bCreateHistogram

    vct1D = reshape( vctLS0, 1, [] );
%     length(find(vct1D == 0))
    hFig = figure;
    hist( vct1D + .5, nGenerations + 1 );
    %  N = ( histc( vct1D, -.5 : 1 : nGenerations + .5 ) )
    %  bar( -.5 : 1 : nGenerations + .5, N, 'histc' )
    saveas( hFig,[ strOutFilestem '_hist.png' ] );
end


if max(vctLS0(:)) > 0
    imOut = double(vctLS0)/max(vctLS0(:));
end

clear vctLS;


if ~isempty( strMotifFilename ) || ( nargin == 8 && bRandom == 0 )
    
    % Replace sybolic values with corresponding motifs.
    for i =  1 : szTiling(1)
        for j = 1 : szTiling(2)
            
            imOut( 1 + (i-1)*szMotifs(1):i*szMotifs(1), 1 + (j-1)*szMotifs(2)/nSymbols:j*szMotifs(2)/nSymbols, 1:3 )...
                    = imMotifs(  1 : szMotifs(1), 1 + ( vctLS0( i, j ) - 0)*szMotifs(2)/nSymbols : ( vctLS0( i, j ) + 1 )*szMotifs(2)/nSymbols, 1:3 );
        end
    end
    
    clear vctLS0;  

    imOut = imOut/max(imOut(:));
end

% nimOutTemp = imOut;
% nUpsample = 3;
% for i = 1 : szTiling(1)
% for j = 1 : szTiling(2)
% for k = 1 : nUpsample
% for l = 1 : nUpsample
%     imOut( nUpsample*(i-1)+k+1, nUpsample*(j-1)+l+1 ) = nimOutTemp( i, j ); 
% end
% end
% end
% end

if nDecimate > 1

    nimOutTemp1 = downsample( imOut, nDecimate, nOffset );
    nimOutTemp1 = nimOutTemp1';
    nimOutTemp1 = downsample( nimOutTemp1, nDecimate, nOffset );
    clear imOut;
    imOut = nimOutTemp1';
    
%     nimOutTemp1(:,:) = downsample( imOut(:,:,1), nDecimate, nOffset );
%     nimOutTemp2(:,:) = downsample( imOut(:,:,2), nDecimate, nOffset );
%     nimOutTemp3(:,:) = downsample( imOut(:,:,3), nDecimate, nOffset );
%     nimOutTemp1 = nimOutTemp1';
%     nimOutTemp2 = nimOutTemp2';
%     nimOutTemp3 = nimOutTemp3';
%     nimOutTemp1 = downsample( nimOutTemp1, nDecimate, nOffset );
%     nimOutTemp2 = downsample( nimOutTemp2, nDecimate, nOffset );
%     nimOutTemp3 = downsample( nimOutTemp3, nDecimate, nOffset );
%     clear imOut;
%     imOut(:,:,1) = nimOutTemp1';
%     imOut(:,:,2) = nimOutTemp2';
%     imOut(:,:,3) = nimOutTemp3';
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Draw system graphical representation ("algorithm graphic"):
%

nES = 16;   % rule element size (even)
nSS = 4;    % spacing size (even)
nB  = 2;    % element border size (even)

% Read auxilliary graphic icons (arrows, bars, identity, rotation, etc.).
imRuleGraphic1 = imread( 'LS_rule_graphic_1.png' );
szRulesGraphic1 = size( imRuleGraphic1 );
imRuleGraphic2 = imread( 'LS_rule_graphic_2.png' );
szRulesGraphic2 = size( imRuleGraphic2 );
imRuleGraphic3 = imread( 'LS_rule_graphic_3.png' );
szRulesGraphic3 = size( imRuleGraphic3 );
imRuleGraphic4 = imread( 'LS_rule_graphic_4.png' );
szRulesGraphic4 = size( imRuleGraphic4 );

imRotateGraphic(:,:,:,1) = imread( 'LS_rotate_0_graphic.png' );
imRotateGraphic(:,:,:,2) = imread( 'LS_rotate_1_graphic.png' );
imRotateGraphic(:,:,:,3) = imread( 'LS_rotate_2_graphic.png' );
imRotateGraphic(:,:,:,4) = imread( 'LS_rotate_3_graphic.png' );
imRotateGraphic(:,:,:,5) = imread( 'LS_mirror_rotate_0_graphic.png' );
imRotateGraphic(:,:,:,6) = imread( 'LS_mirror_rotate_1_graphic.png' );
imRotateGraphic(:,:,:,7) = imread( 'LS_mirror_rotate_2_graphic.png' );
imRotateGraphic(:,:,:,8) = imread( 'LS_mirror_rotate_3_graphic.png' );
szRotateGraphic = size( imRotateGraphic );


nMotifTopSpace = 0;
if ~isempty( strMotifFilename )
    nMotifTopSpace = 2*nSS + szRules(2)*nES + szRulesGraphic3(1);
end


szInit = size( arrayInit );
imAlgorithmGraphic = 0;

% Draw background color onto full image size.
nInitWidth = 0;
if szInit(2) <= szRules(2) + 1
    nInitWidth = szInit(2)*( nES + nB ) + nB;
else
    nInitWidth = ( 1 + szRules(2) )*( szInit(2)/szInit(1) )*( nES + nB ) + nB;
end

xRight = nInitWidth + nSS*( 2 + nSymbols ) + nES*( nSymbols*szRules(2) ) + nB*( 2 + (nSymbols-1)*( 1 + szRules(2) ) ) + szRules(2) + szRulesGraphic1(2);
if ~isempty( strMotifFilename )
    imAlgorithmGraphic( 1 : 2*nSS + nES + 2*nB + szRulesGraphic2(1) + szRules(1)*(nES + nB) + nB + nMotifTopSpace, 1 : xRight, 1:3 ) = (10/13)*(nSymbols-1);
else
    imAlgorithmGraphic( 1 : 2*nSS + nES + 2*nB + szRulesGraphic2(1) + szRules(1)*(nES + nB) + nB,                  1 : xRight, 1:3 ) = (10/13)*(nSymbols-1);
end

if ~bRandom
    % Draw initial symbol matrix.

    xo = 1 + nSS;
    yo = 1 + nSS + nMotifTopSpace;
    if szInit(2) <= szRules(2) + 1
        for j = 1:szInit(2)
            for k = 1:szInit(1)

                % Draw initial matrix backgrounds.
                imAlgorithmGraphic( yo + (k-1)*( nES + nB )      : yo + k*( nES + nB ) + nB - 1, xo + (j-1)*( nES + nB )      : xo + j*( nES + nB ) + nB - 1, 1:3 ) = (5/13)*(nSymbols-1);
                % Draw initial matrix foregrounds.
    %            imAlgorithmGraphic( yo + (k-1)*( nES + nB ) + nB : yo + k*( nES + nB )      - 1, xo + (j-1)*( nES + nB ) + nB : xo + j*( nES + nB )      - 1, 1:3 ) = arrayInit(k,j);        
                imAlgorithmGraphic( yo + (k-1)*( nES + nB ) + nB : yo + k*( nES + nB )      - 1, xo + (j-1)*( nES + nB ) + nB : xo + j*( nES + nB )      - 1, 1:3 ) = mod( (1/1)*(nSymbols-1)*double(imRotateGraphic(:,:,:, arrayInitRotes(k,j) + 4*arrayInitMirrs(k,j) + 1 ))/255 + arrayInit(k,j), nSymbols );
                imAlgorithmGraphic( imAlgorithmGraphic > nSymbols ) = imAlgorithmGraphic( imAlgorithmGraphic > nSymbols ) - nSymbols;
            end
        end
    else
    
        imInit = arrayInit;
        ndim = nInitWidth;
        dy = ( szInit(1) - 1 )/ndim;
        dx = ( szInit(2) - 1 )/ndim;
        yi = 1 -.5 + dy/2 : dy : szInit(1) +.5 - dy/2;
        yi = yi';
        xi = 1 -.5 + dx/2 : dx : szInit(2) +.5 - dx/2;
        if szInit(2) >= ndim      

            X =   1:szInit(2);
            Y = ( 1:szInit(1) )';
            
%            imInitSmall(:,:) = interp2( X,Y, double( imInit(:,:) ), xi, yi, 'linear' );
            imInitSmall(:,:) = interp2_fa( X,Y, imInit(:,:), xi, yi, 'linear' );
        else
            imInitSmall(:,:) = interp2( 1:szInit(2), 1:szInit(1), imInit(:,:), xi, yi, 'nearest' );
        end

        % Draw motifs.
        imInitSmall = (nSymbols-1)*imInitSmall/max(max(imInitSmall));
        szInitSmall = size( imInitSmall );
        imAlgorithmGraphic( yo : yo + szInitSmall(1) - 1, xo : xo + szInitSmall(2) - 1, 1 ) = imInitSmall;
        imAlgorithmGraphic( yo : yo + szInitSmall(1) - 1, xo : xo + szInitSmall(2) - 1, 2 ) = imInitSmall;
        imAlgorithmGraphic( yo : yo + szInitSmall(1) - 1, xo : xo + szInitSmall(2) - 1, 3 ) = imInitSmall;
    end

    % Draw input symbol and output matrix (each rule).
    for i = 1 : nSymbols

        % Draw input rule background and then foreground.
        xLeft  = 1 + nInitWidth + 2*nSS + (( szRules(2)-1)*(nB + nES))/2 + (i-1)*( nES*szRules(2) + nSS + nB*( szRules(2) + 1 ) ) + szRulesGraphic1(2) - 2 + nB/2;
        xRight = xLeft + 2*nB + nES - 1;
        imAlgorithmGraphic( 1 + nSS + nMotifTopSpace      : nSS + nES + nMotifTopSpace + 2*nB, xLeft      : xRight,      1:3 ) = (5/13)*(nSymbols-1);
        imAlgorithmGraphic( 1 + nSS + nMotifTopSpace + nB : nSS + nES + nMotifTopSpace +   nB, xLeft + nB : xRight - nB, 1:3 ) = i - 1;

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

                % Draw output rule background then foreground.
                yTop = 1 + nSS + nES + 2*nB + szRulesGraphic2(1) + (k-1)*(nES + nB) + nMotifTopSpace;
                yBot =     nSS + nES + 2*nB + szRulesGraphic2(1) +    k *(nES + nB) + nMotifTopSpace + nB;
                xLeft  = 1 + nInitWidth + nSS*( 1 + i ) + nES*( (i-1)*szRules(2) + (j - 1) ) + nB*( (i-1)*( 1 + szRules(2) ) + (j - 1) ) + szRulesGraphic1(2);
                xRight = xLeft + nES + 2*nB - 1;

                % Draw output rule background (border).
                imAlgorithmGraphic( yTop      : yBot,       xLeft      : xRight     , 1:3 ) = (5/13)*(nSymbols-1);
                
                % Add rotation graphic, such that the graphics are luminance
                % opposites of the rule color.
    % Note: This pasting of rotation graphic restricts the size of each symbol
    % graphic to 16 x 16 pxls.
                imAlgorithmGraphic( yTop + nB : yBot - nB,  xLeft + nB : xRight - nB, 1:3 ) = mod( (1/1)*(nSymbols-1)*double(imRotateGraphic(:,:,:, arrayRotes(k,j,i) + 4*arrayMirrs(k,j,i) + 1 ))/255 + arrayRules(k,j,i), nSymbols );
                imAlgorithmGraphic( imAlgorithmGraphic > nSymbols - 1 ) = imAlgorithmGraphic( imAlgorithmGraphic > nSymbols - 1 ) - nSymbols;
                %max(max(imAlgorithmGraphic))
                %imAlgorithmGraphic( yTop + nB : yBot - nB,  xLeft + nB : xRight - nB, 1:3 )               

            end
        end
    end

    if ~isempty( strMotifFilename )

        % Draw motifs graphics at top.
        imRuleGraphic3 = double(imRuleGraphic3)/255;
        imRuleGraphic3( find( imRuleGraphic3 == imRuleGraphic3(1,1) ) ) = imAlgorithmGraphic(1,1,1);
        for i = 1:nSymbols

            ndim = ( szRules(2)*nES + nB ) + nB + 1;
            
            % Motif identity signs:
            xo = 1 + nInitWidth + 2*nSS + (( szRules(2)-1)*(nB + nES))/2 + nB/2 + (i-1)*( nES*szRules(2) + nSS + nB*( szRules(2) + 1 ) ) + szRulesGraphic1(2);
            yo = nSS + ndim + nB;
            imAlgorithmGraphic( yo : yo + szRulesGraphic3(1) - 1, xo : xo + (nB + nES) - 1, 1 ) = imRuleGraphic3;
            imAlgorithmGraphic( yo : yo + szRulesGraphic3(1) - 1, xo : xo + (nB + nES) - 1, 2 ) = imRuleGraphic3;
            imAlgorithmGraphic( yo : yo + szRulesGraphic3(1) - 1, xo : xo + (nB + nES) - 1, 3 ) = imRuleGraphic3;

            ithMotif = imMotifs(  1 : szMotifs(1), 1 + (i-1)*szMotifs(2)/nSymbols : i*szMotifs(2)/nSymbols, : );
            szIthMotifs = size( ithMotif(:,:,1) );

            % Initialize small motif with edge color.
            ithMotifSmall( 1:ndim, 1:ndim, 1:3 ) = (5/13)*(nSymbols-1);
            
            dy = szIthMotifs(1)/(ndim-2);
            dx = szIthMotifs(2)/(ndim-2);

            if szIthMotifs(2) >= ndim
                
                yi = 1 - .5 + dy/2 : dy : szIthMotifs(1) + .5 - dy/2;
                yi = yi';
                xi = 1 - .5 + dx/2 : dy : szIthMotifs(2) + .5 - dx/2;
               	ithMotifSmall( 2:ndim-1, 2:ndim-1, 1 ) = interp2_fa( ithMotif(:,:,1), xi, yi, 'linear' )*(nSymbols-1);
                ithMotifSmall( 2:ndim-1, 2:ndim-1, 2 ) = interp2_fa( ithMotif(:,:,2), xi, yi, 'linear' )*(nSymbols-1);
                ithMotifSmall( 2:ndim-1, 2:ndim-1, 3 ) = interp2_fa( ithMotif(:,:,3), xi, yi, 'linear' )*(nSymbols-1);

                ithMotifSmall( find( ithMotifSmall < 0            ) ) = 0;
                ithMotifSmall( find( ithMotifSmall > (nSymbols-1) ) ) = (nSymbols-1);

            else
                
                yi = .5 : dy : szIthMotifs(1) + .499999;
                yi = yi';
                xi = .5 : dx : szIthMotifs(2) + .499999;
             	ithMotifSmall( 2:ndim-1, 2:ndim-1, 1 ) = interp2( 1:szIthMotifs(1), 1:szIthMotifs(2), ithMotif(:,:,1), xi, yi, 'nearest' )*(nSymbols-1);
                ithMotifSmall( 2:ndim-1, 2:ndim-1, 2 ) = interp2( 1:szIthMotifs(1), 1:szIthMotifs(2), ithMotif(:,:,2), xi, yi, 'nearest' )*(nSymbols-1);
                ithMotifSmall( 2:ndim-1, 2:ndim-1, 3 ) = interp2( 1:szIthMotifs(1), 1:szIthMotifs(2), ithMotif(:,:,3), xi, yi, 'nearest' )*(nSymbols-1);
            end

            % Draw motifs.
            szMotifSmall = size( ithMotifSmall );
            xo = 1 + nInitWidth + 2*nSS + szRulesGraphic1(2) + (i-1)*( nES*szRules(2) + nSS + nB*( szRules(2) + 1 ) );
            yo = 1 + nSS;
            imAlgorithmGraphic( yo : yo + szMotifSmall(1) - 1, xo : xo + szMotifSmall(2) - 1, : ) = ithMotifSmall;

        end
    end

    imAlgorithmGraphic = imAlgorithmGraphic + min( imAlgorithmGraphic(:) );
    imAlgorithmGraphic = imAlgorithmGraphic/max( imAlgorithmGraphic(:) );

    % Insert line graphics:
    % Set feedback graphic background to algorithm background color:
    imRuleGraphic1 = double(imRuleGraphic1)/255;
    imRuleGraphic1( find( imRuleGraphic1 == imRuleGraphic1(1,1) ) ) = imAlgorithmGraphic(1,1);
    % Forward arrow.
    yTop   = 1 + nSS + (nES + 2*nB)/2 + nMotifTopSpace - 9;
    yBot   = yTop + 17;
    xLeft  = 1 + nSS + nInitWidth + nSS/2;
    xRight =  xLeft + szRulesGraphic1(2) - 1;
    imAlgorithmGraphic( yTop : yBot, xLeft : xRight, 1 ) = imRuleGraphic1( 1:18, : );
    imAlgorithmGraphic( yTop : yBot, xLeft : xRight, 2 ) = imRuleGraphic1( 1:18, : );
    imAlgorithmGraphic( yTop : yBot, xLeft : xRight, 3 ) = imRuleGraphic1( 1:18, : );
    % Return arrow.
    yTop   = 1 + nSS + nES + 2*nB + szRulesGraphic2(1) + ( k*(nES + nB) + nB )/2 + nMotifTopSpace - 7;
    yBot   = yTop + 14;
    xLeft  = 1 + nSS + nInitWidth + nSS/2;
    xRight = xLeft + szRulesGraphic1(2) - 1;
    imAlgorithmGraphic( yTop : yBot, xLeft : xRight, 1 ) = imRuleGraphic1( 19:33, : );
    imAlgorithmGraphic( yTop : yBot, xLeft : xRight, 2 ) = imRuleGraphic1( 19:33, : );
    imAlgorithmGraphic( yTop : yBot, xLeft : xRight, 3 ) = imRuleGraphic1( 19:33, : );
    % Vertical bar.
    yTop   = 1 + nSS + (nES + 2*nB)/2 + nMotifTopSpace;
    yBot   = yTop + szRulesGraphic2(1) + ( (k+1)*(nES + nB) + nB )/2;
    xLeft  = 1 + nSS + nInitWidth + 16;
    xRight = xLeft + 1;
    imAlgorithmGraphic( yTop : yBot, xLeft : xRight, 1 ) = 0;
    imAlgorithmGraphic( yTop : yBot, xLeft : xRight, 2 ) = 0;
    imAlgorithmGraphic( yTop : yBot, xLeft : xRight, 3 ) = 0;
    % Rule identity:
    imRuleGraphic2 = double(imRuleGraphic2)/255;
    imRuleGraphic2( find( imRuleGraphic2 == imRuleGraphic2(1,1) ) ) = imAlgorithmGraphic(1,1,1);
    yTop   = 1 + nSS + 2*nB + nES + nMotifTopSpace;
    yBot   = yTop + szRulesGraphic2(1) - 1;
    for i = 1:nSymbols

        xLeft  = 1 + nInitWidth + 2*nSS + (( szRules(2)-1)*(nB + nES))/2 + nB/2 + (i-1)*( nES*szRules(2) + nSS + nB*( szRules(2) + 1 ) ) + szRulesGraphic1(2);
        xRight = xLeft + nB + nES - 1;
        imAlgorithmGraphic( yTop : yBot, xLeft : xRight, 1 ) = imRuleGraphic2;
        imAlgorithmGraphic( yTop : yBot, xLeft : xRight, 2 ) = imRuleGraphic2;
        imAlgorithmGraphic( yTop : yBot, xLeft : xRight, 3 ) = imRuleGraphic2;
    end

    % Set transparent color based on identity with top left corner color?
    szA = size( imAlgorithmGraphic );
    imAlgorithmAlpha( 1:szA(1), 1:szA(2) ) = 1;
    imAlgorithmAlpha( find( imAlgorithmGraphic(:,:,1) == imAlgorithmGraphic(1,1) & imAlgorithmGraphic(:,:,2) == imAlgorithmGraphic(1,1) & imAlgorithmGraphic(:,:,3) == imAlgorithmGraphic(1,1) ) ) = 0;

    if bShowResults
        figure; imshow( imAlgorithmGraphic );
    end

    % end of:
    % Draw system graphical representation ("algorithm graphic"):
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

if     bDrawSmallAlgorithm ... 
   &&  szRules(1) == 2  &&  szRules(2) == 2 ...  
   &&  szInit(1)  == 1  &&  szInit(2) == 1 ...  
   && ~bRandom

    nimSmallAlgorithm(  1:128,  1:128 ) = nimByGenerationGraphic( 1, 1 );
    nimSmallAlgorithm( 74:128, 10:120 ) = nimByGenerationGraphic( 4:58, 57:167 );

    nimSmallAlgorithm(  1:73,  44:128 ) = imAlgorithmGraphic(  2:74, 60:144 );
    nimSmallAlgorithm(  1:26,   1: 26 ) = imAlgorithmGraphic(  2:27,  2: 27 );
    nimSmallAlgorithm(  8:20,  42: 53 ) = imAlgorithmGraphic(  9:21, 50: 61 );
    nimSmallAlgorithm( 10:60,  26: 42 ) = imAlgorithmGraphic( 11:61, 33: 49 );
    szBG = size( nimSmallAlgorithm );
    nimSmallAlgorithmAlpha( 1:szBG(1), 1:szBG(2) ) = 1;
    nimSmallAlgorithmAlpha( find( nimSmallAlgorithm(:,:) == nimSmallAlgorithm(1,1) ) ) = 0;
    
	if bShowResults
        figure; imshow( nimSmallAlgorithm );        
    end
    
    if nGenerations + log2(nOversampleOut) <= 6
        
        nimSmallAlgorithmTemp = nimSmallAlgorithm;
        %nimSmallAlgorithmTemp( 76:125, 11:108 ) = nimSmallAlgorithmTemp( 76:125, 13:110 );
        % Nearest neighbor downsample by x2.
        nimTinyAlgorithm = downsample( rot90( nimSmallAlgorithmTemp, -1 ), 2 );
        nimTinyAlgorithm = downsample( rot90( nimTinyAlgorithm, 1 ), 2 );
        % Erase undersampled arrow regions.
        nimTinyAlgorithm(  1:12, 14:25)	= nimSmallAlgorithm(1,1);
        nimTinyAlgorithm( 13:17, 14:59)	= nimSmallAlgorithm(1,1);
        nimTinyAlgorithm( 17:28, 14:22)	= nimSmallAlgorithm(1,1);
        
        % Fix up details:
        
            % Top row of squares:
            nimTinyAlgorithm(  1:15,  2: 6 ) = nimTinyAlgorithm( 1:15,  1: 5 );
            nimTinyAlgorithm(  7:12,  1:15 ) = nimTinyAlgorithm( 8:13,  1:15 );
            nimTinyAlgorithm(  7:13, 26:59 ) = nimTinyAlgorithm( 8:14, 26:59 );
            nimTinyAlgorithm(  1:13, 27:33 ) = nimTinyAlgorithm( 1:13, 26:32 );
            nimTinyAlgorithm(  1:13, 53:58 ) = nimTinyAlgorithm( 1:13, 54:59 );

            % Rules:
            nimTinyAlgorithm( 16:37, 24:28 ) = nimTinyAlgorithm( 16:37, 23:27 );
            nimTinyAlgorithm( 16:37, 38:42 ) = nimTinyAlgorithm( 16:37, 39:43 );
            nimTinyAlgorithm( 16:37, 49:63 ) = nimTinyAlgorithm( 16:37, 50:64 );
            nimTinyAlgorithm( 16:37, 54:63 ) = nimTinyAlgorithm( 16:37, 55:64 );
            
            nimTinyAlgorithm( 17:36, 25:64 ) = nimTinyAlgorithm( 18:37, 25:64 );
            nimTinyAlgorithm( 21:36, 25:64 ) = nimTinyAlgorithm( 22:37, 25:64 );
            nimTinyAlgorithm( 29:36, 25:64 ) = nimTinyAlgorithm( 30:37, 25:64 );
            
            % By Generations:
            nimTinyAlgorithm( 37:63,  6:59 ) = nimTinyAlgorithm( 38:64,  6:59 );
            nimTinyAlgorithm( 37:63,  6:29 ) = nimTinyAlgorithm( 37:63,  7:30 );

            nimTinyAlgorithm(    37,  6:31 ) = (5/13)*(nSymbols-1);
            nimTinyAlgorithm( 37:62,    31 ) = (5/13)*(nSymbols-1);
            
            nimTinyAlgorithm(    37, 34:59 ) = (5/13)*(nSymbols-1);
            nimTinyAlgorithm( 37:62,    59 ) = (5/13)*(nSymbols-1);
       
        % Set feedback graphic background to algorithm background color:
        imRuleGraphic4 = double(imRuleGraphic4)/255;
        imRuleGraphic4( find( imRuleGraphic4 == imRuleGraphic4(1,1) ) ) = nimSmallAlgorithm(1,1);
        nimTinyAlgorithmTemp( 1:64,  1:64 ) = nimSmallAlgorithm(1,1);
        nimTinyAlgorithmTemp( 1:30, 11:60 ) = imRuleGraphic4;
        iArrow = find( nimTinyAlgorithmTemp ~= nimSmallAlgorithm(1,1) );
        % Insert arrow graphics
        nimTinyAlgorithm( iArrow ) = nimTinyAlgorithmTemp( iArrow );
        
        % Make alpha mask:
        szBG = size( nimTinyAlgorithm );
        nimTinyAlgorithmAlpha( 1:szBG(1), 1:szBG(2) ) = 1;
        nimTinyAlgorithmAlpha( find( nimTinyAlgorithm(:,:) == nimTinyAlgorithm(1,1) ) ) = 0;
        
        if bShowResults
            figure; imshow( nimTinyAlgorithm );        
        end
    end
end

bShowDFT = false;
if bShowResults && bShowDFTtoo
    bShowDFT  = true;
end

bWriteDFT = false;
if bWriteResults && bWriteDFTtoo && ~bWriteTilingOnly
    bWriteDFT = true;
end

bWriteDFTavg = false;
if bWriteResults && bWriteDFTtoo && ~bWriteTilingOnly && bDrawAverage
    bWriteDFTavg = true;
end

if bShowDFTtoo || bWriteDFTtoo
    dft_of_image( imOut, nOversampleOut, bShowDFT, bWriteDFT, [ strOutFilestem '-' num2str(nGenerations) 'g' ], flPhaseDriftX, flPhaseDriftY, flPhaseDriftOffset );
end

imOut = interp2_n( imOut, nOversampleOut );

if bWriteResults
    
    strOx  = '';
    if nOversampleOut > 1
        strOx = [ '-' num2str(nOversampleOut) 'x' ];
    end

    
    if bShowResults
        fprintf( [ '\nWriting image file and graphics to the current directory, \n' ] );
        fprintf( [ 'with the file prefix: ' strOutFilestem '\n' ] );
    end
    
    imwrite( imOut, [ strOutFilestem '-' num2str(nGenerations) 'g' strOx '.png' ] );
end

if bDrawAverage && ( bShowDFTtoo || bWriteDFTtoo )
    dft_of_image( imAvg, nOversampleOut, bShowDFT, bWriteDFTavg, [ strOutFilestem '-' num2str(nGenerations) 'g_avg' ], flAvgPhaseDriftX, flAvgPhaseDriftY, flAvgPhaseDriftOffset );
end

if bDrawAverage
    imAvg = interp2_n( imAvg, nOversampleOut );
end

if bDrawAverage

    imAvg = imAvg/max(imAvg(:));
end

if bShowResults

    % Show tiling and average images in a figure window.
    figure; imshow( imOut );
    
    if bDrawAverage
        
        figure; imshow( imAvg )
    end
end


if bWriteResults
    
    % Write documentation graphics:
    if bWriteTilingOnly == 0

        if ~bRandom
            
            imwrite( uint8(imAlgorithmGraphic*255),     [ strOutFilestem '_algorithm_graphic.png'    ], 'BitDepth', 8, 'Alpha', imAlgorithmAlpha    );
            if szRules(1) > 1 && szRules(2) > 1
                imwrite( uint8(nimByGenerationGraphic*255), [ strOutFilestem '_byGeneration_graphic.png' ], 'BitDepth', 8, 'Alpha', imByGenerationAlpha );
                if     bDrawSmallAlgorithm ... 
                   &&  szInit(1)  == 1  &&  szInit(2) == 1 
                    imwrite( uint8(nimSmallAlgorithm*255), [ strOutFilestem '-7g_algorithm_graphic.png' ], 'BitDepth', 8, 'Alpha', nimSmallAlgorithmAlpha );
                    if nGenerations + log2(nOversampleOut) <= 6
                        imwrite( uint8(nimTinyAlgorithm*255), [ strOutFilestem '-6g_algorithm_graphic.png' ], 'BitDepth', 8, 'Alpha', nimTinyAlgorithmAlpha );
                    end
                end
            end
            if szRules(1) > 1 && szRules(2) > 1 && szRules(1) < 3 && szInit(1) == 1 && szInit(2) == 1
                if nGenerations >= 8
                    imwrite( imScaleSym, [ strOutFilestem '-8g-symmetry_graphic.png' ] );
                else
                    imwrite( imScaleSym, [ strOutFilestem '-' num2str( nGenerations + log2(nOversampleOut) ) 'g-symmetry_graphic.png' ] );
                end
            end
            if bDrawAverage
                imwrite( uint8(imAvg*255),              [ strOutFilestem '-' num2str(nGenerations) 'g_avg' strOx '.png'    ], 'BitDepth', 8 );
            end
        end

        % Construct 1 px/symbol horizontal rule image, and save with a name
        % related to the image and algorithm graphic.
        clear imRules;
        imRules( 1:szRules(1), 1:nSymbols*szRules(2) ) = 0;
        for i = 1:nSymbols
            for j = 1:szRules(2)
                for k = 1:szRules(1)
                    imRules( k, j + (i-1)*szRules(2) ) = arrayRules(k,j,i);
                end
            end
        end

        imRules = imRules/max(max(max(imRules)));
        if bDrawRulesAsImage
            imwrite( imRules, [ strOutFilestem '_rules.png' ] );
        end

        % If a motif (or color mapping of rules) is provided, resave with a 
        % name related to the image and algorithm graphic.
        if ~isempty( strMotifFilename ) || nargin == 8
            imwrite( imMotifs, [ strOutFilestem '_motifs.png' ] );
        end
    end  % End write documentation graphics.
end

if bShowResults
    fprintf( ['\n\n'] );
end