function [opened] = areaopenclose(photo, minarea)
%  Area open close filter
% 1999-02-09:DAV: David Cary <d.cary at ieee.org> converted to a function.
% 1999-01-??:SA:started by Scott Acton
%
% See also AREAOPENCLOSE_DEMO.

% Optimization possibilities:
% * Once you discover that some component is too small and set it to zero,
%   perhaps it would speed things up to go ahead and set it to zero
%   in *all* "higher" level sets (and similar for setting to one).
% * perhaps, rather than labeling each level set individually,
%   just label each connected component region of "1"s with its actual size
%   (rather than a unique number), leave all "0"s as zero,
%   and mash them in a single pass with
%     region = find( area_label < minarea );
%     levelset(region) = 0;
%   or
%     levelset = ( area_label < minarea ) .* levelset;
%   .
%
% more optimization possibilities:
% do a histogram; not only "heights" with zero pixels in them
% but also heights with less pixels than the given area
% can be skipped. [idea from Scott Acton]
% Perhaps sort all pixels ... since the bwlabel
% is what is eating the time,
% perhaps we could have a optimized version
% that did less work, just enough
% for the areaopenclose.
%

% function [opened] = areaopenclose(photo, minarea)
% Pseudocode for a faster version:
%
% Mark all pixels as undone by clearing a "done" flag for each pixel.
% Pick x = any undone pixel.
% 
% do{
%  -- Find and count
%     all pixels that are in the same connected component region as x.
%     They all have intensity photo(x).
%     While you are doing that,
%     remember ``the'' border pixel b that has the maximum intensity
%     of all of the pixels bordering this region.
%     (optional: if you find a border pixel b that
%     has not already been done, and photo(x) < photo(b),
%     then don't bother finding any more pixels
%     in this region).
%
%  -- ASSERT( (b is not done) OR ( photo(x) < photo(b) ).
% 
%  -- Do appropriate task of one of the 4 cases {
% (a) If photo(b) < photo(x),
%     and the area of this region is too small.
%     "delete" this object by reducing all pixels in this region to
%     photo(b).
%     Continue with same x or with x = border pixel b
%     or any other pixel in this region.
%     (optional: remember where you left off finding and counting pixels).
% (b) If photo(b) < photo(x),
%     and the area of this region is big enough.
%     "keep" this object (plateau) by marking all pixels in it "done".
%     Pick x = border pixel b or any other undone pixel.
%     (You could mark all pixels in this region with
%     the next "basin" number for watershedding).
% (c) If photo(x) < photo(b),
%     and b is still undone, forget about this region.
%     Pick x = b, handle that region. Later you will come back to this region.
%     (This is hill-climbing, reminiscent of watershed algorithms).
% (c) If photo(x) < photo(b),
%     and b has already been done,
%     then this region plus that bordering region has plenty of area.
%     Mark all the pixels in this region as "done".
%     (You could mark all pixels in this region
%     with the same "basin" number as b).
%     Pick x = any undone pixel.
%  } end case
%
%  -- ASSERT( (x is undone) OR ( every pixel in photo has been done).
%
% }while pixel x is still undone.
%
% As always, the tricky bits are at the image boundaries.
% Simplest would be just to toroidal wrap the boundaries Asteroids(tm)-style.
%
% -- David Cary 1999-05-06.

[height, width] = size(photo);

L=zeros(height, width, 255);  % stack of level sets
for l=1:255,
   % disp(l);
   levelset=(photo>=l);
   [label, num] = bwlabel(levelset,4) ; %label regions of ones w/ 4 connectivity
   for i=1:num,
      region=find(label==i);
      count=length(region);  %area of an individual level set connected component
      if count < minarea,
         levelset(region)=0;  %zero level set components less than min. area
      end
   end
   complement=~levelset;
   [label, num] = bwlabel(complement,4) ; %label the (originally zeros) w/ 4 connectivity
   for i=1:num,
      region=find(label==i);
      count=length(region);  %area of an individual level set connected component
      if count < minarea,
         levelset(region)=1;  %zero level set components less than min. area
      end
   end

   L(:,:,l)=levelset;
end
opened=sum(L,3);