%
% logarithmic_spiral_fixed_edge( nSquare, dr_dtheta, bShow )
% ----------------------------------------------------------
%
%       Creates a black and white image, in unsigned 8-bit format, of two
%   spirals with alternating color between them. The intersection of the 
%   spiral with the square image edge is fixed at the center of the edge. 
%   Use an even size, with any twist/expansion rate, but the result will 
%   not be anti-aliased. For low expansion rates, the Nyquist limit is 
%   exceeded; post-process with block average undersampling.
%
% USAGE:    nimOut = logarithmic_spiral_fixed_edge( 1024,.5, 1 );
%
% nSquare:      [positive integer, even] Size of image, number of pixels on
%               a side.
% dr_dtheta:    [ ~.001, 10 ] Expansion rate of spiral (inverse of twist 
%               rate).
% bShow:        [0, 1] Show progress, and display final image in figure.
%
% The phase of each of the spiral boundaries is hardcoded (see phase, phase2).
%
% CALLS: none

% Mark Dow, December 2006 (extracted and modified from log_spiral_sample.m)

function nimOut = logarithmic_spiral_fixed_edge( nSquare, dr_dtheta, bShow )

nimOut = -1;
if mod( nSquare, 2 ) == 1
    
    fprintf( 'The function logarithmic_spiral_fixed_edge failed: \n' )
    fprintf( 'The size of image (first parameter, nSquare) must be an even integer. \n' )
    return
end

a       = ( nSquare/2 );
b       = dr_dtheta;
phase   = -pi/2;    % currently MUST be -pi/2
phase2  = pi;       % phase of one spiral boundary relative to the other

nimOut( 1:nSquare, 1:nSquare ) = uint8(0);
length_msg = 0;

for i = 1:nSquare

    if bShow
         
        for k = 1:length_msg
            fprintf( '\b' );
        end
        msg = [ 'column =  ' num2str(i) ' of ' num2str(nSquare) ];
        length_msg = length(msg);
        fprintf( msg )
    end
    
    for j = 1:nSquare
    
        x =  nSquare/2 - i + .5;
        y =  nSquare/2 - j + .5;
        
        % Radius from image center to pixel center.
        rp = sqrt( x^2 + y^2 );
        
        % Angle from image center to pixel center.
        ap = atan(y/x);
        
        if x < 0
        
            ap = ap + pi;
        end
        
        k = floor( log(.5/a)/(b*2*pi) ) - 2;                    

        b_past = 0;
        while b_past == 0
        
                    
            r1 = a*exp( b*( ap + k*2*pi + phase ) );
            
            if rp > nSquare/2
                
                b_past = 1;
                
                if x < 0
                
                    nimOut( i, j )  = 1;
                end
                
                if x > 0
                
                    nimOut( i, j )  = 0;
                end
                                
            end
            
            if r1 > rp
            
                nimOut( i, j )  = 1;
                
                b_past = 1;
                      
            else
            
                r2 = a*exp( b*( ap + k*2*pi + phase + phase2 ) );
            
                if r2 > rp
            
                    nimOut( i, j )  = 0;
                
                    b_past = 1;
                end
            end
            
            
            k = k + 1;
        end
                
            
    end
end


if bShow
    
    fprintf( '\n')
    figure;
    imshow( double(nimOut) );
    colormap(gray);
end