%
% Thue_Morse_rings(  nBandWidth, nOversample, flGenerationOffset, bShow, bWrite  )
%
%   Map sequential generations of Thue-Morse sequence onto concentric rings.
%
% USAGE: imOut = Thue_Morse_rings( 25, 3, 2, true, false );
%
% ARGUMENTS:
%
%   nBandWidth:     Width of a band, in pixels, representing one generation.
%
%   nOversample:	n x n oversampling of each pixel.
%
%   flGenerationOffset: Thue-Morse sequence generation at center of image.
%
% RETURN VALUES:
%
%   imOut:  The grayscale [0,1] image matrix, nSize x nSize.
%
% HARDCODED:  (see below 'function' statement)
%
%   bShow:      Show resulting image in a figure window.
%   bWrite:     Write to .png file with parameters in filename.
%
%   strFilename: If bWrite == true, the full output file name.
%
%   bDiscrete:  [ true, false ] Discrete (integral) or continuous sampling
%               of the Thue-Morse generation.
%
%   nSize:      Size of image.
%
% CALLS:    (none)
%
%
% Mark Dow,             June 22, 2009
%

function imOut = Thue_Morse_rings( nBandWidth, nOversample, flGenerationOffset, bShow, bWrite )

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

%     strFilename = 'Thue-Morse_rings.png';
% bDiscrete = true;
% nSize       = 2048;
%     flC = nSize/2;
% flGenerationOffset = 1.0;
% nBandWidth  = 100;        % Width of a band, in pixels, representing one generation.
% nOversample = 9;

%     strFilename = 'Thue-Morse_airfoil_2.png';
% bDiscrete = false;
% nSize       = 2048;
%     flC = 1300;             % ~10.2*nSize/16 - 1;
% flGenerationOffset = 2.1;
% nBandWidth  = 100;          % Width of a band, in pixels, representing one generation.
% nOversample = 13;

    strFilename = 'Thue-Morse_continuous_-2_25_5_thumbnail.png';
 bDiscrete = false;
 nSize       = 512;
    flC = nSize/2;         

%
%%%%%%%%%%%%%%%%%%%%%%%%

imOut = zeros( nSize );

for j = 1/(2*nOversample) : 1/nOversample : nSize - 1/(2*nOversample)
for i = 1/(2*nOversample) : 1/nOversample : nSize - 1/(2*nOversample)
    
    % Find angle and fractional position on the sequence.
    angle = mod( atan2( ( j - flC ), ( i - flC ) ) + pi + .00000000, 2*pi );
    fPos = angle/(2*pi);
    
    % Find radius and generation.
    r = sqrt( ( j - flC )*( j - flC ) + ( i - flC )*( i - flC ) );
    
    if bDiscrete
%       nG = ceil( .00006*r^2.5/nBandWidth ) + flGenerationOffset;
        nG = floor( r/nBandWidth ) + flGenerationOffset;
    else
        nG =     ( r/nBandWidth ) + flGenerationOffset;
    end

    lengthTM_at_nG = 2^nG;
    nPos = floor( fPos*lengthTM_at_nG);
    
    % Find parity of sum of base 2 bits representation of integer 
    % position on sequence (the Thue-Morse symbol).
    % Iteratively divide by 2, find mod 2, increment an odd counter
    nOnes = 0;
    nPosRemainder = nPos;
    while nPosRemainder > 0
        
        if mod( nPosRemainder, 2 ) == 1
            nOnes = nOnes + 1;
        end
        
        nPosRemainder = floor( nPosRemainder/2 );
    end
    
    % Color white if odd parity.
    if mod( nOnes, 2 ) == 1
        %imOut( round(j+.5), round(i+.5) ) = imOut( round(j+.5), round(i+.5) ) + 1;
        imOut( floor(j+1), floor(i+1) ) = imOut( floor(j+1), floor(i+1) ) + 1;
    end
end
end

imOut = imOut/(nOversample^2);

if bShow
    figure; imshow( imOut ); colormap( gray );
end

if bWrite
    imwrite( imOut, strFilename  )
end