% [ flim_regions_labeled, flim_average_of_region ] = ss_region_average( filestem, b_write_image )
% ---------------------------------------------
%
% Component of Scale Space Edges, version 1.2b
%
%
% Mark Dow, last modified April 2003
% University of Oregon
% Brain Development Lab / Lewis Center for NeuroImaging
% dow@braindev.uoregon.edu
%


function [ flim_regions_labeled, flim_average_of_region, n_regions ] = ss_region_average( filestem, b_write_image )


load( [ filestem '_edge_info' ] )

n_erode_max = 0;

sz_Image = size( flim_gray_cropped );

% To Do: Generate the B/W image directly from edge points.
% Generate boundary image ( rom average gradient of segment info).
flim_boundaries = zeros( sz_Image );
for i = 1:seg_n    
    for n_ep = seg(i,1) : seg(i,2)
        
        if seg(i,5) > flim_boundaries( ep(n_ep, 1), ep(n_ep, 2) )
            flim_boundaries( ep(n_ep, 1), ep(n_ep, 2) ) = seg(i,5);
        end
    end
end

flim_average_of_region = flim_gray_cropped;

bwim = im2bw( flim_boundaries, .0000001);
%bwim = imdilate( bwim, strel('diamond',1) ); 
%bwim = bwmorph( bwim, 'skel', 'Inf' ); 
[ flim_regions_labeled, n_regions ] = bwlabel( ~bwim, 4 );

%  figure
%  rgbim = label2rgb(flim_regions_labeled);
%  imagesc( rgbim );

flim_average_of_region( find( bwim > 0 ) ) = 0;

% For each sub-region:
for i = 1:n_regions
%     i
%     num
    posns = find( flim_regions_labeled == i  );
%    if length(posns) > 1
        flim_average_of_region( posns ) = sum( flim_gray_cropped( posns ) ) / length( posns );
%    else
%        flim_average_of_region( posns ) = 0;
%   end
end

bwim = ~bwim;

% Erode and relabel by average several times:
for n_erode  = 1:n_erode_max
 bwim = imerode( bwim, strel('diamond', 1));
 
 [ Lsk, numsk ] = bwlabel( bwim, 4 );
 
 for i = 1:numsk
%      i
%      numsk
     
     posns_sk = find( Lsk == i  );
     
     if ( length(posns_sk) > 5 )
         
         j = flim_regions_labeled( posns_sk(1) ); 
         posns = find( flim_regions_labeled == j  );
         
         flim_average_of_region( posns ) = sum( flim_gray_cropped( posns_sk ) ) / length( posns_sk );
     end
 end
end
 

flim_average_of_region = imdilate( flim_average_of_region, strel( 'diamond', 1 ) );

save( [ filestem '_edge_info' ], 'ep', 'seg', 'node', 'seg_n', 'flim_gray_cropped', 'flim_seg', 'flim_gradient', 'flim_scale', 'flim_regions_labeled', 'flim_average_of_region', 'n_regions' )

if b_write_image == 1
    
    imwrite( flim_average_of_region, [ filestem '_region_ave.tif' ] );
    
end