%
% dft_of_image( Image, nOversampleOut, bShow, bWrite, strOutFileStem, 
%               flPhaseDriftX, flPhaseDriftY, flPhaseDriftOffset )
%
%   Generates logarithmically scaled DFT (discrete Fourier transform) 
%   graphics of an image or other matrix, using MATLAB's 2-D fast 
%   Fourier transform(FFT) function, fft2.
%
%   Produces amplitude, phase, and amplitude colored by phase image,
%   optionally centered on highest frequencies or the DC (zero frequency) 
%   component. 
%
%   Currently the image values are rounded to 8-bit precision, 
%   for display purposes only.
%
%   Color image input is reduced to a luminance approximation.
%
%   Currently indexed image input is not supported.
%
% USAGE: imDFT0 = dft_of_image( 'test-2x2-2sm-i0-0_288-6_9-9g.png' );
%                 dft_of_image( 'test-2x2-2sm-i0-0_288-6_9-9g.png', 3 );
%                 dft_of_image( 'test-2x2-2sm-i0-0_288-6_9-9g.png', 1, 1, 0, '', pi, pi );
%
% ARGUMENTS:
%
%       Image:          Either an image matrix or the file name of image, 
%                       with extension.
%
%       nOversampleOut: Integer (nearest neighbor) oversampling of output.
%
%       bShow:          Show all results in figure windows.
%                       Optional; default is to show but not write.
%
%       bWrite:         Write all results to files.
%                       Optional; default is to not write.
%
%       strOutFileStem: strOutFileStem
%                       Optional; default will use the input file stem, 
%                       or ??? if a matrix.
%
%       Specialty parameters for segragating a linear phase drift from 
%       the Fourier phase diagram, generating a separate graphic of this
%       phase drift:
%
%           flPhaseDriftX:      In radians (default = 0).
%           flPhaseDriftY:    	In radians (default = 0).
%
%           flPhaseDriftOffset: In radians (default = 0).
%
%
% RETURN VALUES:
%
%       imDFT0:         The raw (complex) DFT image.
%
% HARDCODED:  (see below 'function' statement)
%
%   bZeroCentered           = true;
%   bShowOriginalToo        = true;
%
%	strOutputType           = '.png';   % Common image file extensions, with the dot prefix:
%                                       % '.png', '.tiff', '.jpg'
%	strOutputPostfix        = 'DFT';    % Prefix for output file names.
% 
%   bVioletGreenLuminanceMap = false;   % Choose one and only one for phase color map. 
%   bRGBMap                 = true;     
%	flRGBRampOverlap = 1.0;             % ( .25, 1 ] *2*pi/3 
%                                       % for .25: R+B+G = 0 at n*pi/3
%                                       % for   1: R+B+G = 255
% 
%                                     
%	flLogScaleMagFT         = 1.0;                                   
%	flLightenPhaseColoredFT = 1.0;
% 
%   bCenterValues           = true;     % Reduce or remove zero frequency 
%                                         byshifting the image such that 
%                                         min = -max value. 
% CALLS:
%       -> interp2_n.m
%       -> phase2hue_cycle
%           -> hl.m
%               -> hl_color_cycle.m
%
%
%
% Mark Dow,             March 19, 2009
% Mark Dow, modified    March 22, 2009  (RGB color wheel phase code)
% Mark Dow, modified    March 26, 2009  (Phase colored magnitude merging)
% Mark Dow, modified    March 29, 2009  (Phase wrapping and phase pattern)
% Mark Dow, modified    April  3, 2009  (Clean up, support integer values)
% Mark Dow, modified    April  7, 2009  (Compatibility with L_system_tiling.m)
% Mark Dow, modified    April  7, 2009  (independent x-y phase drift)
% Mark Dow, modified      May 10, 2009  (cyclic color map with hl)
% Mark Dow, modified      May 14, 2009  (phase drift arguments)
%

function imDFT0 = dft_of_image( Image, nOversampleOut, bShow, bWrite, strOutFileStem, flPhaseDriftX, flPhaseDriftY, flPhaseDriftOffset )

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

bZeroCentered           = false;    % true for low frequencies at center,
                                    % false for high frequencies at center.
bShowOriginalToo        = true;

strOutputType           = '.png';   % Common image file extensions, with the dot prefix:
                                    % '.png', '.tiff', '.jpg'
strOutputPostfix        = 'DFT';    % Postfix for output file names.
     
bRGBMap                 = true;

flLogScaleMagFT         = 3.0;                                   
flLightenPhaseColoredFT = 1.0;

bCenterValues           = true;     % Reduce or remove zero frequency by
                                    % shifting the image such that min =
                                    % -max value.
                                    
% Deprecated, now using hl.m for color map generation.
%flRGBRampOverlap        = 1;      % ( .25, 1 ] *2*pi/3 
%                                    % for .25: R+B+G = 0 at n*pi/3
%                                    % for   1: R+B+G = 255

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


% Parse argument list and set default values:
%if nargin < 8
    
    if nargin < 8
        flPhaseDriftOffset = 0*pi;
    end
    if nargin < 7
        flPhaseDriftY = 0*pi;
    end
       
    if nargin < 6
        flPhaseDriftX = 0*pi; 
    end

    if nargin < 5
        
        strOutFileStem = strFileName;
%         if isa( Image, 'char' )
%             strFileName = Image;
%         else
%             % Add dummy file name.
%             strFileName = 'image.png'
%         end
    end
    
    if nargin < 4
        bWrite = false;
    end
    if nargin < 3
        bShow  = true;
    end
    if nargin < 2
        nOversampleOut  = 1;
    end
%end

bPhaseDrift    = 0;
if flPhaseDriftX ~= 0 || flPhaseDriftY ~= 0

    bPhaseDrift    = 1;
end


% Load the image from file:
if ~isa( Image, 'char' )
    
    imIn = Image;     
else
   
    strFileName = Image;
    infImage = imfinfo( strFileName );

    if nargin == 5
        
    else
        % Strip extension from the file name.
        nLocation_of_dot = find( strFileName == '.' );
        strOutFileStem      = strFileName( 1 : nLocation_of_dot - 1 );
    end

    if strcmp( infImage.ColorType, 'truecolor' )
        if infImage.Height == 1
            imIn24bit( 1, 1:infImage.Width, 1:3 ) = imread( strFileName );
        else
            imIn = imread( strFileName );
        end
    end

    if strcmp( infImage.ColorType, 'grayscale' )

        imInGray = imread( strFileName );
        szGray = size( imInGray );
        imIn( 1:szGray(1), 1:szGray(2) ) = 0;
        imIn = imInGray;

    %     imIn( 1:szGray(1), 1:szGray(2), 1:3 ) = 0;
    %     imIn(:,:,1) = imInGray;
    %     imIn(:,:,2) = imInGray;
    %     imIn(:,:,3) = imInGray;
    end

    if strcmp( infImage.ColorType, 'indexed' )

        fprintf( [ '\n\n    Sorry, the indexed (color palette) image format is not currently supported.\n' ] );
        fprintf( [ '\n    Use an image editor to change the color format to grayscale or 24-bit color.\n\n' ] );
        imDFT = -1;
        return;
    end
end

   
% If input if 3x8-bit color, reduce to approximate luminance.
szIm = size( imIn );        
if length( szIm ) == 3
    
    % Check to see if this is just grayscale in 3-byte format. 
    % If so, retain integer formatting.
    bComparePlanes12 = ( imIn(:,:,1) == imIn(:,:,2) );
    bComparePlanes13 = ( imIn(:,:,1) == imIn(:,:,3) );
    if ~all( bComparePlanes12 ) | ~all( bComparePlanes12 )

        imIn = double( imIn );
        imIn = .30*imIn(:,:,1) + .11*imIn(:,:,2) + .59*imIn(:,:,3);
    else
        imIn = imIn(:,:,1);
    end
end
szIm  = [ szIm(1) szIm(2) ];


% Show (nearest neighbor interpolated) original image.
if bShow && bShowOriginalToo
    
    if nOversampleOut > 1
        imInOversampled = interp2_n( imIn, nOversampleOut );
    else
        imInOversampled = imIn;
    end
    
    figure; imshow( uint8( 255*( ( double(imInOversampled) - double( min(imInOversampled(:)) ) )/double( max(imInOversampled(:)) - min(imInOversampled(:)) ) ) ) )
    
    clear imInOversampled;
end

% Reduce or remove zero frequency by shifting the image 
%such that min = -max value.
if bCenterValues
    
    % imIn = ( ( imIn - min( imIn(:) ) )/( max( imIn(:) ) - min( imIn(:) ) ) ) - .5;
    imIn = double( imIn ) - double( ( max( imIn(:) ) + min( imIn(:) ) ) )/2;
    % max(imIn(:))
    % min(imIn(:))
end

imDFT = fft2( imIn );
clear imIn;


imDFT0 = imDFT;

% %%%%%%%%
% % test inverse DFT
% 
% imDFT = fft2( imIn(:,:,1) );
% 
% % % Middle third only;
% % imDFTCenter = zeros( [ szIm(1)/3 szIm(2)/3 ] );
% % imDFTCenter = imDFT( szIm(1)/3 + 1:2*szIm(1)/3, szIm(2)/3 + 1:2*szIm(2)/3 );
% 
% % % Whole image
% % imDFTCenter = zeros( [ szIm(1) szIm(2) ] );
% % imDFTCenter = imDFT( 1:szIm(1), 1:szIm(2) );
% 
% % Average of non-zero
% imDFTCenter = zeros( [ szIm(1) szIm(2) ] );
% imDFTCenter = imDFT( 1:szIm(1), 1:szIm(2) );
% imDFTCenter( find( imDFT > 0 ) ) = 1/length( find( imDFT > 0 ) );
% 
% % % Zero the middle third only (square high pass filter);
% % imDFTCenter = zeros( [ szIm(1) szIm(2) ] );
% % imDFTCenter = imDFT( 1:szIm(1), 1:szIm(2) );
% % imDFTCenter( szIm(1)/3 + 1:2*szIm(1)/3, szIm(2)/3 + 1:2*szIm(2)/3 ) = 0;
% 
% imCenterReconstructed = ifft2( imDFTCenter, 'symmetric' );
% imCenterReconstructed = uint8(255*imCenterReconstructed/max(imCenterReconstructed(:)));
% imwrite( imCenterReconstructed, 'test.png' );
% figure
% imshow(imCenterReconstructed)
% colormap(gray)
% 
% %%%%%%%%%

% Don't find angles of zero power components (NaN instead of zero).
imDFT( find( imDFT == 0.0 ) ) = NaN;

if bZeroCentered
    imDFT = fftshift( imDFT );
end

imDFTPhase              = angle( imDFT );
%figure; imagesc(imDFTPhase)

% Logritmic transform of DFT values.
% See http://homepages.inf.ed.ac.uk/rbf/HIPR2/pixlog.htm for
% loagaritmic operator form, and
% http://homepages.inf.ed.ac.uk/rbf/HIPR2/pixexp.htm for exponential
% operator form.
imDFT = abs( real( imDFT ) );
imDFT = log( imDFT.*exp( flLogScaleMagFT ) + 1 )/log( 1 + max(imDFT(:) ) );


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Center on DC (zero frequency) component
%
%   Note: This is replaced (above) with the MATLAB function fftshift.
%
% 
%     % Deal with an an odd dimension.
%     szQu1 = [ floor( szIm(1)/2 ) floor( szIm(2)/2 ) ]; % T-L
%     szQu2 = [ floor( szIm(1)/2 ) floor( szIm(2)/2 ) ]; % T-R
%     szQu3 = [ floor( szIm(1)/2 ) floor( szIm(2)/2 ) ]; % B-R
%     szQu4 = [ floor( szIm(1)/2 ) floor( szIm(2)/2 ) ]; % B-L
%     if mod( szIm(1), 2) == 1
%         szQu1(1) = szQu1(1) + 1;
%         szQu2(1) = szQu2(1) + 1;
%     end
%     if mod( szIm(2), 2) == 1
%         szQu1(2) = szQu1(2) + 1;
%         szQu4(2) = szQu4(2) + 1;
%     end
% 
%     if bZeroCentered
% 
%         % Q3 <= Q1
%         imDFT( szQu4(1) + 1 :  szIm(1), szQu2(2) + 1 : szIm(2)  ) = imDFT( 1:szQu1(1), 1:szQu1(2) );
%         % Q2 <= Q4
%         imDFT(            1 : szQu4(1), szQu3(2) + 1 : szIm(2)  ) = imDFT( szQu1(1) + 1 : szIm(1), 1:szQu4(2) );
%         % Q1 <= Q3
%         imDFT(            1 : szQu3(1),            1 : szQu3(2) ) = imDFT( szQu1(1) + 1 : szIm(1), szQu1(2) + 1 : szIm(2) );
%         % Q4 <= Q2
%         imDFT( szQu3(1) + 1 : szIm(1),             1 : szQu3(2) ) = imDFT( 1:szQu2(1), szQu1(2) + 1 : szIm(2) );
% 
%         imDFTPhase( szQu4(1) + 1 :  szIm(1), szQu2(2) + 1 : szIm(2)  ) = imDFTPhase( 1:szQu1(1), 1:szQu1(2) );
%         imDFTPhase(            1 : szQu4(1), szQu3(2) + 1 : szIm(2)  ) = imDFTPhase( szQu1(1) + 1 : szIm(1), 1:szQu4(2) );
%         imDFTPhase(            1 : szQu3(1),            1 : szQu3(2) ) = imDFTPhase( szQu1(1) + 1 : szIm(1), szQu1(2) + 1 : szIm(2) );
%         imDFTPhase( szQu3(1) + 1 : szIm(1),             1 : szQu3(2) ) = imDFTPhase( 1:szQu2(1), szQu1(2) + 1 : szIm(2) );
%         
%         imPhaseWrapped( szQu4(1) + 1 :  szIm(1), szQu2(2) + 1 : szIm(2)  ) = imPhaseWrapped( 1:szQu1(1), 1:szQu1(2) );
%         imPhaseWrapped(            1 : szQu4(1), szQu3(2) + 1 : szIm(2)  ) = imPhaseWrapped( szQu1(1) + 1 : szIm(1), 1:szQu4(2) );
%         imPhaseWrapped(            1 : szQu3(1),            1 : szQu3(2) ) = imPhaseWrapped( szQu1(1) + 1 : szIm(1), szQu1(2) + 1 : szIm(2) );
%         imPhaseWrapped( szQu3(1) + 1 : szIm(1),             1 : szQu3(2) ) = imPhaseWrapped( 1:szQu2(1), szQu1(2) + 1 : szIm(2) );
%     end
%     
%     clear imDFT;
%     clear imDFTPhase;
% 
% %   Center on DC (zero frequency) component
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Supress maximum amplitude to a fixed multiplier of the next largest value.
flMaxPenultRatio = 1.1;
iDFTMax = find( imDFT == max( imDFT(:) ) );
if ~isempty(iDFTMax) && imDFT( iDFTMax(1) ) > 0
    
    vDFTMax = imDFT( iDFTMax(1) );
    imDFT( iDFTMax ) = 0;   % set max location to zero
    iDFTMax2 = find( imDFT == max( imDFT(:) ) );    % get penultimate max
    if ~isempty(iDFTMax2) && imDFT( iDFTMax2(1) ) > 0 && imDFT( iDFTMax2(1) )/imDFT( iDFTMax2(1) ) > flMaxPenultRatio
        imDFT( iDFTMax ) = 1.1*imDFT( iDFTMax2(1) );
        imDFT = imDFT/( flMaxPenultRatio*imDFT( iDFTMax2(1) ) );
    else
        imDFT( iDFTMax ) = vDFTMax; % replace with original value
    end
end

imDFT = imDFT/max(imDFT(:));

% imPhaseWrapped      = imDFTPhase;
imPhaseNoDrift      = imDFTPhase;
imPhaseDrift        = imDFTPhase;
imPhaseSign         = zeros( szIm, 'uint8' );
 
imDFTColored        = zeros(  [ szIm(1), szIm(2), 3 ], 'uint8'  );
imDFTWrappedClr     = zeros(  [ szIm(1), szIm(2), 3 ], 'uint8'  );

imPhaseClr          = zeros( [ szIm(1) szIm(2) 3 ], 'uint8' );
% imPhaseWrappedClr	= zeros( [ szIm(1) szIm(2) 3 ], 'uint8' );
imPhaseDriftClr     = zeros( [ szIm(1) szIm(2) 3 ], 'uint8' );
imPhaseNoDriftClr   = zeros( [ szIm(1) szIm(2) 3 ], 'uint8' );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Extract and remove phase complement pattern.

for ii = 1 : szIm(1)
for jj = 1 : szIm(2)

    if ~isnan( imDFTPhase( ii, jj ) )

%         imPhaseDrift( ii, jj )      = bPhaseDrift*(ii-1+jj-1)/(szIm(1)+szIm(2)) - bPhaseDrift/2;        
%         a      = ( ( flPhaseDriftY*(ii-1)/szIm(1) + flPhaseDriftX*(jj-1)/szIm(2) ) - ( flPhaseDriftY + flPhaseDriftX )/2 );        
        imPhaseDrift( ii, jj )      = mod( ( flPhaseDriftY*(ii-1)/szIm(1) + flPhaseDriftX*(jj-1)/szIm(2) ) - ( flPhaseDriftY + flPhaseDriftX )/2 + flPhaseDriftOffset, 2*pi );        
        while imPhaseDrift( ii, jj ) < -pi
            imPhaseDrift( ii, jj ) = imPhaseDrift( ii, jj ) + 2*pi;
        end
        while imPhaseDrift( ii, jj ) > pi
            imPhaseDrift( ii, jj ) = imPhaseDrift( ii, jj ) - 2*pi;
        end
    end
end
end


imPhaseNoDrift = imDFTPhase - imPhaseDrift;
%figure; hist(imPhaseNoDrift(:),1000)
        
iLow = find( imPhaseNoDrift <= -pi + .0000001 );
imPhaseNoDrift(iLow) = imPhaseNoDrift(iLow) +  2*pi;

iPi = find( imPhaseNoDrift > pi - .0000001 & imPhaseNoDrift < pi + .0000001 );
imPhaseNoDrift(iPi) = pi;
iPz = find( imPhaseNoDrift > -.0000001 & imPhaseNoDrift < .0000001 );
imPhaseNoDrift(iPz) = 0;

iPp2h = find( imPhaseNoDrift >  pi/2 - .0000001 & imPhaseNoDrift <  pi/2 + .0000001 );
imPhaseNoDrift(iPp2h) =  pi/2;
iPp2l = find( imPhaseNoDrift < -pi/2 + .0000001 & imPhaseNoDrift > -pi/2 - .0000001 );
imPhaseNoDrift(iPp2l) = -pi/2;

iHgh = find( imPhaseNoDrift >   pi );
imPhaseNoDrift(iHgh) = imPhaseNoDrift(iHgh) -  2*pi;

% Count the number of unique entries, to decide how to display the phase
% with no drift image.
bPhaseSignOnly = false;
if length( unique( imPhaseNoDrift(:) ) ) <= 5
    bPhaseSignOnly = true; % not currently used
    %flLightenPhaseColoredFT = 3.0*flLightenPhaseColoredFT;
end

% Fill in integer multiple of pi/2 radian map:
% To Do: Replace with find() equivalent.
for ii = 1 : szIm(1)
for jj = 1 : szIm(2)

    if ~isnan( imDFTPhase( ii, jj ) )
        

        if abs( imPhaseNoDrift( ii, jj ) - pi   ) == 0.0 || abs( imPhaseNoDrift( ii, jj ) +  pi   ) == 0.0
            imPhaseSign( ii, jj ) = 255;
        end
        if abs( imPhaseNoDrift( ii, jj ) + pi/2 ) == 0.0
            imPhaseSign( ii, jj ) = 64;
        end
        if abs( imPhaseNoDrift( ii, jj ) -    0 ) == 0.0
            imPhaseSign( ii, jj ) = 128;
        end
        if abs( imPhaseNoDrift( ii, jj ) - pi/2 ) == 0.0 
            imPhaseSign( ii, jj ) = 192;
        end
%         if abs( imPhaseNoDrift( ii, jj ) -   pi ) < .0000001 || abs( imPhaseNoDrift( ii, jj ) +   pi  ) > -.0000001
%             imPhaseSign( ii, jj ) = 255;
%         end

    end
end
end

%imPhaseSign( find( imPhaseSign ==   0 ) ) = 127;
%imPhaseSign( find( imPhaseSign == 128 ) ) = 0;

% imPhaseSignColor = zeros( [ szIm(1), szIm(2), 3 ] );
% 
% imPhaseSignColor( find(imPhaseSign == 254 ) + 0*szIm(1)*szIm(2) ) = 255;
% imPhaseSignColor( find(imPhaseSign == 254 ) + 1*szIm(1)*szIm(2) ) = 128;
% %imPhaseSignColor( find(imPhaseSign ==   0 ) + 1*szIm(1)*sxIm(2) ) = ;
% imPhaseSignColor( find(imPhaseSign == 0 ) + 2*szIm(1)*szIm(2) ) = 255;
% imPhaseSignColor( find(imPhaseSign == 0 ) + 1*szIm(1)*szIm(2) ) = 64;

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

%%%%%%%%%%%%%%% Example code from an earlier abandoned color scheme.
%
% % Recode phase map with red/blue-green for +/- phase and luminance for
% % magnitude.
% if bRedBlueLuminanceMap
%     
%     flPhaseMaxsP = 0;
%     flPhaseP = imDFTPhase( find( imDFTPhase > 0 ) );
%     if ~isempty(flPhaseP)
%         flPhaseMaxsP = max( imDFTPhase( find( imDFTPhase > 0 ) ) );
%     end
% 
%     flPhaseMaxsM = 0;
%     flPhaseM = imDFTPhase( find( imDFTPhase < 0 ) );
%     if ~isempty(flPhaseM)
%         flPhaseMaxsM = max( imDFTPhase( find( imDFTPhase < 0 ) ) );
%     end
% 
%     LogValsPos = log(  imDFTPhase( find( imDFTPhase > 0 ) )/pi + 1 )/log(1+flPhaseMaxsP(1)/pi );
%     LogValsNeg = log( -imDFTPhase( find( imDFTPhase < 0 ) )/pi + 1 )/log(1+flPhaseMaxsM(1)/pi );
%     LogValsPos2 = ( (50).^( imDFTPhase( find( imDFTPhase > 0 ) )/pi - 1 ) );
%     LogValsNeg2 = ( (50).^(-imDFTPhase( find( imDFTPhase < 0 ) )/pi - 1 ) );
% 
%     % size(imDFTPhase( find( imDFTPhase < 0 ) ))
%     % size(imDFTPhase( find( imDFTPhase < 0 ) ))
%     % size(LogValsNeg2)
%     % size(find( imDFTPhase < 0 ))
% 
%     imPhaseClr( find( imDFTPhase > 0 ) + 0*szIm(1)*szIm(2) ) =  255*LogValsPos;
%     imPhaseClr( find( imDFTPhase > 0 ) + 1*szIm(1)*szIm(2) ) =  255*LogValsPos2;
%     imPhaseClr( find( imDFTPhase > 0 ) + 2*szIm(1)*szIm(2) ) =  255*LogValsPos2;
% 
%     imPhaseClr( find( imDFTPhase < 0 ) + 0*szIm(1)*szIm(2) ) = 255*LogValsNeg2;
%     imPhaseClr( find( imDFTPhase < 0 ) + 1*szIm(1)*szIm(2) ) = ((.31 - .1)/.59)*255*LogValsNeg + (1 - ((.31 - .1)/.59))*255*LogValsNeg2;
%     imPhaseClr( find( imDFTPhase < 0 ) + 2*szIm(1)*szIm(2) ) = 255*LogValsNeg;
% end
%%%%%%%%%%%%%%%


if bRGBMap
    
% Deprecated method of color cycle rendering, which effectively duplicated
% MATLAB's "hsv" and "jet" color maps:
%
%     % Color complete phase map:
%     % Red linear ramp from 0 to 2*pi/3 and -2*pi/3 to 2*pi/3
%     i1 = find( imDFTPhase >  2*pi/3 - flRampW & imDFTPhase <= 2*pi/3           );
%     i2 = find( imDFTPhase >  2*pi/3           & imDFTPhase <=   pi             );
%     i3 = find( imDFTPhase >=  -pi             & imDFTPhase < -4*pi/3 + flRampW );
%     
%     imPhaseClr( i1 + 0*szIm(1)*szIm(2) ) = 255*(  imDFTPhase(i1) - ( 2*pi/3 - flRampW )          )/flRampW;
%     imPhaseClr( i2 + 0*szIm(1)*szIm(2) ) = 255*( -imDFTPhase(i2) - ( 2*pi/3 - flRampW ) + 4*pi/3 )/flRampW;
%     imPhaseClr( i3 + 0*szIm(1)*szIm(2) ) = 255*( -imDFTPhase(i3) - ( 4*pi/3 - flRampW )          )/flRampW;
% 
%     % Green linear ramp from 0 to -2*pi/3 and -2*pi/3 to 2*pi/3
%     i1 = find( imDFTPhase < -2*pi/3 + flRampW & imDFTPhase >= -2*pi/3           );
%     i2 = find( imDFTPhase < -2*pi/3           & imDFTPhase >=   -pi             );
%     i3 = find( imDFTPhase <=   pi             & imDFTPhase >   4*pi/3 - flRampW );
%     
%     imPhaseClr( i1 + 1*szIm(1)*szIm(2) ) = 255*( -imDFTPhase(i1) - ( 2*pi/3 - flRampW )          )/flRampW;
%     imPhaseClr( i2 + 1*szIm(1)*szIm(2) ) = 255*(  imDFTPhase(i2) - ( 2*pi/3 - flRampW ) + 4*pi/3 )/flRampW;
%     imPhaseClr( i3 + 1*szIm(1)*szIm(2) ) = 255*(  imDFTPhase(i3) - ( 4*pi/3 - flRampW )          )/flRampW;
%     
%     % Blue linear ramp from 0 to 2*pi/3 and 0 to -2*pi/3
%     i1 = find( imDFTPhase >= 0       & imDFTPhase <=  2*pi/3 );
%     i2 = find( imDFTPhase <  0       & imDFTPhase >= -2*pi/3 );
%     
%     imPhaseClr( i1 + 2*szIm(1)*szIm(2) ) = 255*( -imDFTPhase(i1) + flRampW )/flRampW;
%     imPhaseClr( i2 + 2*szIm(1)*szIm(2) ) = 255*(  imDFTPhase(i2) + flRampW )/flRampW;

    % Map phase values to color cycle rendering:
    imPhaseClr          = phase2hue_cycle( imDFTPhase );
    imPhaseDriftClr     = phase2hue_cycle( imPhaseDrift ); 
    imPhaseNoDriftClr   = phase2hue_cycle( imPhaseNoDrift ); 
     
%     flRampW = flRGBRampOverlap*2*pi/3;  % ( .25, 1 ] *2*pi/3
%  
%     
%     
%     % Color phase drift map:
%     % Red linear ramp from 0 to 2*pi/3 and -2*pi/3 to 2*pi/3
%     i1 = find( imPhaseDrift >  2*pi/3 - flRampW & imPhaseDrift <= 2*pi/3           );
%     i2 = find( imPhaseDrift >  2*pi/3           & imPhaseDrift <=   pi             );
%     i3 = find( imPhaseDrift >=  -pi             & imPhaseDrift < -4*pi/3 + flRampW );
%     
%     imPhaseDriftClr( i1 + 0*szIm(1)*szIm(2) ) = 255*(  imPhaseDrift(i1) - ( 2*pi/3 - flRampW )          )/flRampW;
%     imPhaseDriftClr( i2 + 0*szIm(1)*szIm(2) ) = 255*( -imPhaseDrift(i2) - ( 2*pi/3 - flRampW ) + 4*pi/3 )/flRampW;
%     imPhaseDriftClr( i3 + 0*szIm(1)*szIm(2) ) = 255*( -imPhaseDrift(i3) - ( 4*pi/3 - flRampW )          )/flRampW;
% 
%     % Green linear ramp from 0 to -2*pi/3 and -2*pi/3 to 2*pi/3
%     i1 = find( imPhaseDrift < -2*pi/3 + flRampW & imPhaseDrift >= -2*pi/3           );
%     i2 = find( imPhaseDrift < -2*pi/3           & imPhaseDrift >=   -pi             );
%     i3 = find( imPhaseDrift <=   pi             & imPhaseDrift >   4*pi/3 - flRampW );
%     
%     imPhaseDriftClr( i1 + 1*szIm(1)*szIm(2) ) = 255*( -imPhaseDrift(i1) - ( 2*pi/3 - flRampW )          )/flRampW;
%     imPhaseDriftClr( i2 + 1*szIm(1)*szIm(2) ) = 255*(  imPhaseDrift(i2) - ( 2*pi/3 - flRampW ) + 4*pi/3 )/flRampW;
%     imPhaseDriftClr( i3 + 1*szIm(1)*szIm(2) ) = 255*(  imPhaseDrift(i3) - ( 4*pi/3 - flRampW )          )/flRampW;
%     
%     % Blue linear ramp from 0 to 2*pi/3 and 0 to -2*pi/3
%     i1 = find( imPhaseDrift >= 0       & imPhaseDrift <=  2*pi/3 );
%     i2 = find( imPhaseDrift <  0       & imPhaseDrift >= -2*pi/3 );
%     
%     imPhaseDriftClr( i1 + 2*szIm(1)*szIm(2) ) = 255*( -imPhaseDrift(i1) + flRampW )/flRampW;
%     imPhaseDriftClr( i2 + 2*szIm(1)*szIm(2) ) = 255*(  imPhaseDrift(i2) + flRampW )/flRampW;

    
%     % Color phase minus phase drift map:
%     % Red linear ramp from 0 to 2*pi/3 and -2*pi/3 to 2*pi/3
%     i1 = find( imPhaseNoDrift >  2*pi/3 - flRampW & imPhaseNoDrift <= 2*pi/3           );
%     i2 = find( imPhaseNoDrift >  2*pi/3           & imPhaseNoDrift <=   pi             );
%     i3 = find( imPhaseNoDrift >=  -pi             & imPhaseNoDrift < -4*pi/3 + flRampW );
%     
%     imPhaseNoDriftClr( i1 + 0*szIm(1)*szIm(2) ) = 255*(  imPhaseNoDrift(i1) - ( 2*pi/3 - flRampW )          )/flRampW;
%     imPhaseNoDriftClr( i2 + 0*szIm(1)*szIm(2) ) = 255*( -imPhaseNoDrift(i2) - ( 2*pi/3 - flRampW ) + 4*pi/3 )/flRampW;
%     imPhaseNoDriftClr( i3 + 0*szIm(1)*szIm(2) ) = 255*( -imPhaseNoDrift(i3) - ( 4*pi/3 - flRampW )          )/flRampW;
% 
%     % Green linear ramp from 0 to -2*pi/3 and -2*pi/3 to 2*pi/3
%     i1 = find( imPhaseNoDrift < -2*pi/3 + flRampW & imPhaseNoDrift >= -2*pi/3           );
%     i2 = find( imPhaseNoDrift < -2*pi/3           & imPhaseNoDrift >=   -pi             );
%     i3 = find( imPhaseNoDrift <=   pi             & imPhaseNoDrift >   4*pi/3 - flRampW );
%     
%     imPhaseNoDriftClr( i1 + 1*szIm(1)*szIm(2) ) = 255*( -imPhaseNoDrift(i1) - ( 2*pi/3 - flRampW )          )/flRampW;
%     imPhaseNoDriftClr( i2 + 1*szIm(1)*szIm(2) ) = 255*(  imPhaseNoDrift(i2) - ( 2*pi/3 - flRampW ) + 4*pi/3 )/flRampW;
%     imPhaseNoDriftClr( i3 + 1*szIm(1)*szIm(2) ) = 255*(  imPhaseNoDrift(i3) - ( 4*pi/3 - flRampW )          )/flRampW;
%     
%     % Blue linear ramp from 0 to 2*pi/3 and 0 to -2*pi/3
%     i1 = find( imPhaseNoDrift >= 0       & imPhaseNoDrift <=  2*pi/3 );
%     i2 = find( imPhaseNoDrift <  0       & imPhaseNoDrift >= -2*pi/3 );
%     
%     imPhaseNoDriftClr( i1 + 2*szIm(1)*szIm(2) ) = 255*( -imPhaseNoDrift(i1) + flRampW )/flRampW;
%     imPhaseNoDriftClr( i2 + 2*szIm(1)*szIm(2) ) = 255*(  imPhaseNoDrift(i2) + flRampW )/flRampW;
%     
    
end

% Apply lightening to phase colorings of DFT amplitude.
for i = 1 : szIm(1)
for j = 1 : szIm(2)

%    if bPhaseSignOnly
%         imPhaseClr( i, j, 1 ) = (1+flLightenPhaseColoredFT)*double( imPhaseClr(i,j,1) );
%         imPhaseClr( i, j, 2 ) = (1+flLightenPhaseColoredFT)*double( imPhaseClr(i,j,2) );
%         imPhaseClr( i, j, 3 ) = (1+flLightenPhaseColoredFT)*double( imPhaseClr(i,j,3) );
%    end

    imDFTColored( i, j, 1 ) = (1+flLightenPhaseColoredFT)*double( imDFT(i,j) )*( imPhaseClr(i,j,1) );
    imDFTColored( i, j, 2 ) = (1+flLightenPhaseColoredFT)*double( imDFT(i,j) )*( imPhaseClr(i,j,2) );
    imDFTColored( i, j, 3 ) = (1+flLightenPhaseColoredFT)*double( imDFT(i,j) )*( imPhaseClr(i,j,3) );
%     imPhaseDriftClr( i, j, 1 ) = (1+flLightenPhaseColoredFT)*double( imPhaseDriftClr(i,j,1) );
%     imPhaseDriftClr( i, j, 2 ) = (1+flLightenPhaseColoredFT)*double( imPhaseDriftClr(i,j,2) );
%     imPhaseDriftClr( i, j, 3 ) = (1+flLightenPhaseColoredFT)*double( imPhaseDriftClr(i,j,3) );
    
%     if bPhaseDrift ~= 0
%         imPhaseNoDriftClr( i, j, 1 ) = (1+flLightenPhaseColoredFT)*double( imPhaseNoDrift(i,j,1) );
%         imPhaseNoDriftClr( i, j, 2 ) = (1+flLightenPhaseColoredFT)*double( imPhaseNoDrift(i,j,2) );
%         imPhaseNoDriftClr( i, j, 3 ) = (1+flLightenPhaseColoredFT)*double( imPhaseNoDrift(i,j,3) );
%     end

    if imDFTColored( i, j, 1 ) > 255
        imDFTColored( i, j, 1 ) = 255;
    end
    if imDFTColored( i, j, 1 ) < 0
        imDFTColored( i, j, 1 ) = 0;
    end
    if imDFTColored( i, j, 2 ) > 255
        imDFTColored( i, j, 2 ) = 255;
    end
    if imDFTColored( i, j, 2 ) < 0
        imDFTColored( i, j, 2 ) = 0;
    end
    if imDFTColored( i, j, 3 ) > 255
        imDFTColored( i, j, 2 ) = 255;
    end
    if imDFTColored( i, j, 3 ) < 0
        imDFTColored( i, j, 3 ) = 0;
    end

end
end

imDFT = 255*imDFT;

% Nearest neighbor oversample by an integer.
if nOversampleOut > 1
    
        imDFT               = interp2_n( imDFT,             nOversampleOut );
        imDFTColored        = interp2_n( imDFTColored,      nOversampleOut );
    
    if bPhaseDrift == 0
        imPhaseClr          = interp2_n( imPhaseClr,        nOversampleOut );
    end
    
    if bPhaseDrift ~= 0
        imPhaseDriftClr     = interp2_n( imPhaseDriftClr,	nOversampleOut );
    end
    
    if bPhaseDrift ~= 0
        imPhaseNoDriftClr   = interp2_n( imPhaseNoDriftClr,	nOversampleOut );
    end
    
        imPhaseSign         = interp2_n( imPhaseSign,       nOversampleOut );
end


if bShow
    
        figure; imshow( uint8( imDFTColored ) )
        figure; imshow( uint8( imDFT ) )
    
    if bPhaseDrift == 0
        figure; imshow( uint8( imPhaseClr ) )    
    end

    if bPhaseDrift ~= 0
        figure; imshow( uint8( imPhaseDriftClr ) );
    end
    
    if bPhaseDrift ~= 0

        figure; imshow( uint8( imPhaseNoDriftClr ) );
    end
    
        figure; imshow( uint8( imPhaseSign ) );

end

if bWrite
    
    % Encode some hardcoded parameters into the file stem.
    strStem = [ strOutFileStem '-' strOutputPostfix ];
    strEnd  = [ strOutputType ];
    if nOversampleOut > 1
        strEnd = [ '_x' num2str(nOversampleOut) strEnd ];
    end
    if bZeroCentered
        strEnd = [ '_zc' strEnd ];
    end
    
% Only use some formats with 'BitDepth' specified
    if (   strcmp( strOutputType,'.png'  ) ...
        || strcmp( strOutputType,'.tiff' ) ...
        || strcmp( strOutputType,'.jpg'  ) )
    
            imwrite( uint8(imDFTColored),       [ strStem '_color'          strEnd ], 'BitDepth', 8 );
            imwrite( uint8(imDFT),              [ strStem '_mag'            strEnd ], 'BitDepth', 8 );
    
        if bPhaseDrift == 0
            imwrite( uint8(imPhaseClr),         [ strStem '_phase'          strEnd ], 'BitDepth', 8 );
        end
        
        if bPhaseDrift ~= 0
            imwrite( uint8(imPhaseDriftClr),	[ strStem '_phase_drift'	strEnd ], 'BitDepth', 8 );
        end
        
        if bPhaseDrift ~= 0
            imwrite( uint8(imPhaseNoDriftClr),  [ strStem '_phase_no_drift'	strEnd ], 'BitDepth', 8 );
        end
        
            imwrite( uint8(imPhaseSign),        [ strStem '_phase_sign'     strEnd ], 'BitDepth', 8 );
        
    else
        
            imwrite( uint8(imDFTColored),       [ strStem '_color'          strEnd ] );
            imwrite( uint8(imDFT),              [ strStem '_mag'            strEnd ] );
    
        if bPhaseDrift == 0
            imwrite( uint8(imPhaseClr),         [ strStem '_phase'          strEnd ] );
        end
        
        if bPhaseDrift ~= 0
            imwrite( uint8(imPhaseDriftClr),	[ strStem '_phase_drift'	strEnd ] );
        end
        
        if bPhaseDrift ~= 0
            imwrite( uint8(imPhaseNoDrift),     [ strStem '_phase_no_drift'	strEnd ] );
        end
        
            imwrite( uint8(imPhaseSign),        [ strStem '_phase_sign'     strEnd ] );

    end
end