function bw1=pancake(H,I,bw1,K,L1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% detects roundedness of amydales to remove non-amydgale detection
%
%http://www.mathworks.com/products/image/demos.html?file=/products/demos/sh
%ipping/images/ipexroundness.html
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
plotting=0;
% clear all;close all;clc
% bw1=imread('messy.bmp');
if plotting == 1
    figure(100),subplot(131),imshow(bw1);
end

F=bw1; %save original in background


% UPthreshold = 1;
% LOthreshold = 0;

    [B,L] = bwboundaries(bw1,'noholes');

    if plotting == 1
        %Display the label matrix
        testim=label2rgb(L, @jet, [.5 .5 .5]);
        figure(100),subplot(132), imshow(testim);
%          figure, imshow(testim)
        hold on
        %  and draw each boundary
        for k = 1:length(B)
          boundary = B{k};
          plot(boundary(:,2), boundary(:,1), 'w', 'LineWidth', 2)
        end
    end

    stats = regionprops(L,'Area','Centroid');

    % loop over the boundaries
    newloc=[]; %location of rounded amydales
for k = 1:length(B)

%       % obtain (X,Y) boundary coordinates corresponding to label 'k'
%       boundary = B{k};
% 
%       % compute a simple estimate of the object's perimeter
%       delta_sq = diff(boundary).^2;
%       perimeter = sum(sqrt(sum(delta_sq,2)));
% 
%       % obtain the area calculation corresponding to label 'k'
%       area = stats(k).Area;
% 
%       % compute the roundness metric
%       metric = 4*pi*area/perimeter^2;
%       
%       if plotting==1
%           %display the results
%           metric_string = sprintf('%2.2f',metric);
%       end

      %select amygdales
%       if (metric < UPthreshold) & (metric > LOthreshold)
            % mark objects above the threshold with a black circle
            centroid = stats(k).Centroid;
            locations = [int2str(centroid(1)) int2str(centroid(2))];
            %check if centroid is within "white area"
            if bw1(str2num(locations)) == 1
                newloc=[newloc; str2num(locations)];
            end
        
            if plotting ==1
                plot(centroid(1),centroid(2),'ko');
            end
%       end

%       if plotting ==1
%             text(boundary(1,2)-35,boundary(1,1)+13,metric_string,'Color','y',...
%            'FontSize',8,'FontWeight','bold');
%       end
end

                %This is where the magic happens..
		above=bw1(newloc)-I(newloc);
	        below=bw1(newloc)-K(newloc);
		%identify single-layer thickness
		sing=find(((bw1(newloc)-I(newloc))==1) & (bw1(newloc)-K(newloc))==1);
		if length(sing) >= 1
		  bw1=~bw1;
		  bw2=imfill(bw1,newloc(sing));
		  bw1=~bw2;	
		end
		
		%identify 1st dual-layer thickness
		dualloc=find((bw1(newloc)-K(newloc))==0);
	        if length(dualloc) >= 1
		   if L(newloc(dualloc))==0
			bw1=~bw1;
			K=~K;
			bw2=imfill(bw1,newloc(dualloc));
			K2=imfill(K,newloc(dualloc));
			bw1=~bw2;
			K=~K2;
		   end
		end

		dualloc=find((bw1(newloc)-I(newloc))==0);
                if length(dualloc) >= 1
                   if H(newloc(dualloc))==0
                        bw1=~bw1;
                        I=~I;
                        bw2=imfill(bw1,newloc(dualloc));
                        I2=imfill(I,newloc(dualloc));
                        bw1=~bw2;
                        I=~I2;
                   end
                end	
		
%                for k=1:length(newloc)
%                    loc=newloc(k);
%                    if (I(loc)==0) & (K(loc)==0)
%                        bw1=~bw1;
%                        bw2=imfill(bw1,loc);
%                        bw1=~bw2;
%		    elseif (H(loc)==0) & (I(loc)==1) & (K(loc)==0)
%			bw1=~bw1;
%                        bw2=imfill(bw1,loc);
%                        bw1=~bw2;
%			I=~I;
%			I2=imfill(I,loc);
%			I=~I2;
%		    elseif (I(loc)==0) & (K(loc)==1) & (L1(loc)==0)
%			bw1=~bw1;
%                        bw2=imfill(bw1,loc);
%                       bw1=~bw2;
%			K=~K;
%			K2=imfill(K,loc);
%			K=~K2;
%                    end
%                end

if plotting ==1
    figure(100), subplot(133),imshow(bw1);
end
