%
% nimOut = Edge_contrast_circle( flRadius, imDim, flAngle, ndegPhase, flShiftMax )
% -------------------------------------------------------------------------
%
%   Creates circle image with a center-surround aura and edge at a given angle. 
%   This image is used as a frame the "directional edge contrast motion" illusion.
%
% USAGE:    nimOut = Edge_contrast_circle( 25, 76, pi/4, 30, 5 );
%
% INPUT:
%
%   flRadius    - Radius of the base circle, in pixels.
%   imDim       - One dimension of the resulting square image, in
%                 pixels, unless this dimension will clip the circle image.
%   ndegPhase   - Phase angle, in degrees
%   flShiftMax  - maximum shift, in pixels
%
% HARDCODED (immediately below function declaration):
%
%   strFilename:      - Filename and extension of output file.
%   nCenterMode:      - [ 0-constant (gray), 1-linear, 2-sinusoidal ]
%                       Circle luminance cycling.
%   nOverSample:      - [1, 16] Flat aperture oversampling rate.
%   bShowResult:      - [0, 1] Show progress, display final image.
%   bWriteResult:     - [0, 1] Write image to file.
%
% RETURN VALUES:
%
%   nimOut = grayscale image array of 8-bit integers (uint8)
%
% OUTPUT:
%
%	If bWriteResult == 1, writes grayscale image to current directory
%	with the filename and format specified in the hardcoded section
%	(immediately below function declaration).
%
% FUNCTION CALLS and DEPENDENCIES: 
%
%       None.
%
% Mark Dow,             May 5, 2008
% Mark Dow, modified    May 6, 2008 (added phase and maximum shift)
% Mark Dow, modified    May 6, 2008 (modified to linear changes in time, parabolic
%                                    spatial surround)
% Mark Dow, modified    May 7, 2008 (modified for any angle)
%
%
%   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.
%

function nimOut = Edge_contrast_circle( flRadius, imDim, flAngle, ndegPhase, flShiftMax )

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

strFilename = 'test_edge_circle.png';

nCenterMode = 0;    % [ 0-constant (gray), 1-linear, 2-sinusoidal ]
flMaskWidth = 3;

nOverSample = 6;

bShowResult         = 0;    % [0,1]
bWriteResult        = 0;    % [0,1]

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

imDimMin = 2*( 2*flShiftMax + floor(flRadius) );
if imDim < imDimMin
    imDim = imDimMin;
end

% Each quadrant of phase angles is a simple transform of the first quadrant;
ndegPhase = mod( ndegPhase, 360 );
bRotate = false;
bInvert = false;
if ndegPhase <= 90
    flradPhase = ndegPhase*pi/180;
    bRotate = false;
    bInvert = false;
end
if ndegPhase >  90 && ndegPhase <= 180
    flradPhase = (180 - ndegPhase)*pi/180;
    bRotate = true;
    bInvert = false;
end
if ndegPhase > 180 && ndegPhase <= 270
    flradPhase = (ndegPhase - 180)*pi/180;
    bRotate = false;
    bInvert = true;
end
if ndegPhase > 270 && ndegPhase <= 360
    flradPhase = (360 - ndegPhase)*pi/180;
    bRotate = true;
    bInvert = true;
end

% Initialize all gray image.
if ~bInvert
    imBackground( 1:imDim, 1:imDim ) = 127; 
    nimOut( 1:imDim, 1:imDim ) = 127;
else
    imBackground( 1:imDim, 1:imDim ) = 128; 
    nimOut( 1:imDim, 1:imDim ) = 128;
end

flAngle = flAngle + 2*pi;

flSamplingIncrement = 1/nOverSample;

for i = 1 + flSamplingIncrement/2 : 1/nOverSample : imDim + 1.0001 - flSamplingIncrement/2
for j = 1 + flSamplingIncrement/2 : 1/nOverSample : imDim + 1.0001 - flSamplingIncrement/2
    
    rc = sqrt( ( i - (imDim+2)/2 )^2 + ( j - (imDim+2)/2 )^2 );
    
    % Draw center constant shade.
    if nCenterMode > 0 && rc <= flRadius
        
        if nCenterMode == 1 
            % linear time
            imBackground( floor(i), floor(j) ) = imBackground( floor(i), floor(j) ) - flradPhase/(pi/2)*( 128/(nOverSample^2) );
            nimOut( floor(i), floor(j) ) = imBackground( floor(i), floor(j) );
        end
        
        if nCenterMode == 2 
            % sinusoidal time
            imBackground( floor(i), floor(j) ) = imBackground( floor(i), floor(j) ) - sin( flradPhase )*( 128/(nOverSample^2) );
            nimOut( floor(i), floor(j) ) = imBackground( floor(i), floor(j) );
        end
    end
    
    % Draw surround as a gradiant from center to gray.
    if rc > flRadius && rc < flRadius + flMaskWidth*flShiftMax
%        % sinusoidal time, linear space
%        imBackground( floor(i), floor(j) ) = imBackground( floor(i), floor(j) ) + sin( flradPhase )*( ( imDim/2 - rc )/25.5 )*128/(nOverSample^2);
%        % linear time, linear space
%        imBackground( floor(i), floor(j) ) = imBackground( floor(i), floor(j) ) + flradPhase/(pi/2)*( ( imDim/2 - rc )/25.5 )*128/(nOverSample^2);
       % linear time, parabolic space
       imBackground( floor(i), floor(j) ) = imBackground( floor(i), floor(j) ) + flradPhase/(pi/2)*( (( flMaskWidth*flShiftMax + flRadius - rc )*( flMaskWidth*flShiftMax + flRadius - rc ))/(( flMaskWidth*flShiftMax )^2) )*128/(nOverSample^2);
       nimOut( floor(i), floor(j) ) = imBackground( floor(i), floor(j) );
    end
end 
end
    
rOffsetP = flShiftMax;
rOffsetM = flShiftMax* ( pi/2 - flradPhase )/(pi/2);

for i = 1 + flSamplingIncrement/2 : 1/nOverSample : imDim + 1.0001 - flSamplingIncrement/2
for j = 1 + flSamplingIncrement/2 : 1/nOverSample : imDim + 1.0001 - flSamplingIncrement/2
    
    rc = sqrt( ( i - (imDim+2)/2 )^2 + ( j - (imDim+2)/2 )^2 );   
    
% To Do: Sanity check.
%    if flradPhase > pi/2

    rp = sqrt( ( ( i + rOffsetP*sin(flAngle) ) - (imDim+2)/2 )^2 + ( ( j + rOffsetP*cos(flAngle) ) - (imDim+2)/2 )^2 );
    rm = sqrt( ( ( i - rOffsetM*sin(flAngle) ) - (imDim+2)/2 )^2 + ( ( j - rOffsetM*cos(flAngle) ) - (imDim+2)/2 )^2 );
     
    flAnglePt = atan2( (  -j + (imDim+2)/2),  i - (imDim+2)/2 ) + 0*pi/2;
    flAnglePt = mod( flAnglePt + 4*pi, 2*pi );
    flAngle = mod( flAngle + 4*pi, 2*pi );
%     && ( abs( flAnglePt - flAngle ) < pi/2 )
%   &&  ( abs( rc - rm ) < .5 )
    
% if floor(j) == 19 && floor(i) == 23
%   flAnglePt  
%   flAngle
% end
if (flAnglePt - flAngle) > pi
  flAnglePt = flAnglePt - 2*pi; 
end
if (flAnglePt - flAngle) < -pi
  flAnglePt = flAnglePt + 2*pi; 
end
% if floor(j) == 19 && floor(i) == 23
%   flAnglePt  
%   flAngle
% end
    % Draw white (plus) edge
    if rc > flRadius && rp <= flRadius
            nimOut( floor(i), floor(j) ) = nimOut( floor(i), floor(j) ) + ( 1 - flradPhase/(pi/2) )*( 255-imBackground( floor(i), floor(j) ) )/(nOverSample^2);
    end
    
    % Draw minus (black) edge
    % Note: The trailing logic blackens between cresent and circle, except for
    % the last frame (flradPhase == pi/2) where double darkening makes
    % circle assymetric. This is a satisfactory approximation, but not
    % exact for very low flShiftMaxs or very slow shifts (dense frames near
    % pi/2)
%    if ( rm <= flRadius ) &&  ( ( rc > flRadius ) ||  ( ( rc > flRadius - .5  )  && flradPhase < .99*(pi/2) ) )
%    if ( rm <= flRadius ) &&  ( ( rc > flRadius ) ||  ( ( rc > flRadius - .5  ) && ( abs( flAnglePt - flAngle ) < pi/2 ) ) && flradPhase < .99*(pi/2) )
    if ( rm <= flRadius ) &&  ( ( rc > flRadius )  )
        nimOut( floor(i), floor(j) ) = nimOut( floor(i), floor(j) ) - imBackground( floor(i), floor(j) )/(nOverSample^2);
    end
    

end
end


if bRotate == true
    nimOut = rot90( nimOut, 2 );
end

if bInvert == true
    nimOut = 255 - nimOut;
end


nimOut = uint8( nimOut );

if bShowResult

    figure
    imshow( nimOut )
end

if bWriteResult

    imwrite( nimOut, strFilename )
end