% ss_k_planes( flim_gray, plane_fileprefix, k_max, resample_factor, b_display )
% ----------------------------------------------------------------------------------
%
% Component of Scale Space Edges, version 1.0b.
%
% Preprocessing for an image, including resampling and gaussian blur, resulting
% in a set of scale space sample planes. 
%
% Input:
%
%    flim_gray         - Floating point grayscale image array.	
%
%    plane_fileprefix  - Prefix for output file names, like 'Clown_plane_k'.
%
%    k_max             - Maximum scale index (scale = 2^k).  Gaussian 
%			             smoothed raw images are generated for 
%			             [0, k_max], if the size of the image allows.
%    
%    resample_factor   - Upsampling ratio.
%
%    b_display         - [0 1] Show matlab figures with resampled image and
%                        smoothed images (scale space sample planes).
%    
% Ouput:
%
%     Writes a grayscale image (scale space sample plane), along with processing 
%     parameters, for each scale index k = [ k_min, k_max ] in matlab files named 
%     [plane_fileprefix][k_plane #].mat.
%
%
% Mark Dow, last modified July 2003
% University of Oregon
% Brain Development Lab / Lewis Center for NeuroImaging
% dow@braindev.uoregon.edu
%

 

function ss_k_planes( flim_gray, plane_fileprefix, k_max, resample_factor, b_display );
    
    sz_Image = size( flim_gray );
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Interpolate grayscale image.
    %

    if resample_factor ~= 1
    
        [xi, yi] = meshgrid( 1:1/resample_factor:sz_Image(2), ...
                             1:1/resample_factor:sz_Image(1) );
        flim_gray = interp2( flim_gray, xi, yi, 'cubic'   );
    
        sz_Image = size( flim_gray );
        
        msg = ['    Size of resampled image: ' num2str( sz_Image(2) ) ' w x ' num2str( sz_Image(1) ) ' h\n\n' ];
        fprintf( msg );
    end
    
    %
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Display resampled image.
    
        if b_display == 1
     
            % Draw full resampled image.
            figure;
            imagesc( flim_gray );
            colormap(gray)
            zoom
            end
        %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    % Check k_max compared to image size -- half width of gaussian 
    % filter must be less than the minimum image dimension.

    k_max_limit  = floor( log2( min( sz_Image(1), sz_Image(2) )/4  ) ); 

    if k_max > k_max_limit

        k_max = k_max_limit; 
    
        msg = ['\n    Maximum k reset to ' num2str(k_max_limit) ', based on minimum image dimension. \n\n'];
        fprintf( msg );
    end        


    % For each base scale:
    
    nth = 0;
    for k = 0 : k_max

        nth = nth + 1;
        msg = ['\n    Working on ' num2str(nth) ' of ' num2str( k_max+1 ) ' k planes ....'];
        fprintf( msg );

        % Calculate filter parameters
%        sigma = 2^k;
%        filter_size = 8*(2^k) + 1;    % odd
        if k == 0
            sigma = 1;
        else 
            sigma = sqrt( 2^(2*k) - 2^(2*(k-1)) );
        end 
        filter_size = 12*(ceil(sigma)) + 1;    % odd
%         filter_size = 27;    % odd
        
        e = ( filter_size - 1 )/2; 
        s1 = sz_Image(1); 
        s2 = sz_Image(2);
         
        % Reflection padded copy:
        flim_gray_padded = zeros( [ ( s1 + 2*e ) (s2 + 2*e ) ] );
        % Copy original to center of padded.
        flim_gray_padded(  e+1:e+s1,  e+1:e+s2 ) = flim_gray;
        
        % top
        flim_gray_padded( 1:e, e+1:e+s2 ) ...
                =  flipud( flim_gray( 1:e, : ) );
        % bottom
        flim_gray_padded( e+s1+1:2*e+s1, e+1:e+s2 ) ...
                =  flipud( flim_gray( s1-e+1:s1, : ) );
        % left
        flim_gray_padded( e+1:e+s1, 1:e ) ...
                =  fliplr( flim_gray( : , 1:e ) );
        % right
        flim_gray_padded( e+1:e+s1, e+s2+1:2*e+s2 ) ...
                =  fliplr( flim_gray( : , s2-e+1:s2 ) );
                
        % top/left
        flim_gray_padded( 1:e, 1:e ) = flipud( fliplr( flim_gray( 1:e, 1:e ) ) );
            
        % bottom/left
        flim_gray_padded( e+s1+1:2*e+s1, 1:e ) = ...
                flipud( fliplr( flim_gray( s1-e+1:s1, 1:e ) ) );
            
        % top/right
        flim_gray_padded( 1:e,  e+s2:2*e+s2 ) = ...
                flipud( fliplr( flim_gray( 1:e, s2-e:s2 ) ) );
            
        % bottom/right
        flim_gray_padded( e+s1+1:2*e+s1, e+s2:2*e+s2 ) = ...
                flipud( fliplr( flim_gray( s1-e+1:s1, s2-e:s2 ) ) );
              
%% Draw padded image.
%figure;
%imagesc( flim_gray_padded );
%colormap(gray)
%zoom

       % Apply gaussian smoothing filter.
       h = fspecial( 'gaussian', filter_size, sigma );
       flim_gray = filter2( h, flim_gray_padded, 'valid' );
 
 
%         h1 = ss_DoGauss( 0, sigma, sigma, 0 );       
%         for i = 1:length(h1)
%             for x = e+1 : s1+e
%                 for y = e+1 : s2+e
% 
%                     flim_gray( x-e, y-e) = h1(i,3)*flim_gray_padded( x + h1(i,1), y + h1(i,2) ) + h1(i,3)*flim_gray_padded( x - h1(i,1), y + h1(i,2) );
%                 end
%             end
%         end
      
        clear flim_gray_padded;
        
        % Save.
        fullfilestem = [ plane_fileprefix num2str(k) ];
        save( fullfilestem, 'flim_gray' );

 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Display k plane.
    
        if b_display == 1
     
            % Draw k plane image.
            figure;
%            imagesc( flim_gray );
            image( floor( (63*flim_gray)+1.5) );
            colormap(gray)
            zoom
        end
        %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    end
    
    filestem = [ plane_fileprefix '#' ];    
    msg = ['\n\n    k planes stored in:  ' filestem '.mat\n'];
    fprintf( msg );