%  Granulometry demo using the Area open close filter
%
%  1999-02-09:DAV:David Cary <d.cary at ieee.or> simplified.
%  1999-01-??:SA:started by Scott Acton
%
% See also AREAOPENCLOSE.

% Optimization possibilities:
% * rather than cycle through all 255 images (inside areaopenclose)
%   for each scale,
%   maybe you could directly calculate the granulometry
%   with only one pass.
%   (i.e., Find smallest area(s),
%   add to granulometry data and delete from temp image;
%   then find next-smallest areas
%   (make sure that previous deletions didn't make them bigger -- how?),
%   add to granulometry data and delete from temp image;
%   ... until you have no more regions left).

clear;

[file, path]=uigetfile('*.RAW', 'RAW File Input');
height=input('Image Height: ');
width=input('Image Width: ');

fid=sprintf('%s%s', path, file);
x=loadraw(fid,width,height);
orig=x;
E=zeros(size(x));
minscale=input('Minimum Scale: ');
maxscale=input('Maximum Scale: ');

outfile=input('Output file name for granulometry: ','s');
trueoutfile=sprintf('%s%s', path, outfile);
g=zeros(1,maxscale-minscale);
normg=g;
normgscale=g;
for s=minscale:maxscale,
    s
    minarea=s;
    
    opened=areaopenclose(x, minarea);% Do it !
    
    figure(1)
    colormap(gray(256))
    
    image(opened);
    title(num2str(s));
    drawnow;
    
    %compute granulometry value
    if s>minscale,
        diff=abs(x-opened);
        totalsum=sum(sum(diff))
        totalsum/s
        totalsum/(height*width)
        normgscale(s-minscale)=totalsum/s;
        normg(s-minscale)=totalsum/(height*width);
        g(s-minscale)=totalsum;
    end
    
    
    x=opened;
    
end

figure(2)
%AXIS([minscale maxscale 0 255])
stem(g)
title('Granulometry: Cumulative Intensity Change');
figure(3)
stem(normgscale)
title('Granulometry: Normalized by Scale');
figure(4)
stem(normg)
title('Granulometry: Normalized by Image Size');
%get(2,'Children');
%for i=1:length(ans),
%   set(ans(i),'XTick',[]);
%   set(ans(i),'YTick',[]);
%end

%figure
%colormap(gray(256))

%image(orig);
%title('Input image');
%X=writeraw(trueoutfile,opened);
save(trueoutfile,'g','normgscale','normg');