function [Centroids] = blob2(Mask, k)
% Given a mask (1s and 0s) in mask, find 'blobs' in the mask
% by finding connected regions of 1's.  
% return a vector of centers [x1,y1;...;xk,yk] of the largest
% k blobs
% do this iteratively because the simple recursive implementation
% kills matlab's stack on real images!

if nargin == 1
    k = 3;
end

[rows, cols] = size(Mask);
BlobMask = zeros(rows, cols);

curBlob = 1;

for row = 2:rows %ignore edges
    for col = 2:(cols-1) %ignore edges
         valueL = 0; valueTL = 0; valueT = 0; valueTR = 0;
         myValue = 0;
         if Mask(row, col) == 1 %need to be in a blob
            valueL  = BlobMask(row, col-1);
            valueTL = BlobMask(row-1, col-1);
            valueT  = BlobMask(row-1, col);
            valueTR = BlobMask(row-1, col+1);

            if valueL ~= 0
                myValue = valueL;
            end
            
            if valueTL ~= 0
                myValue = valueTL; %guaranteed same as valueL by prev iterations
            end
            
            if valueT ~= 0
                if valueL ~= 0 & valueL ~= valueT %we have a U connection
                    inds = find(BlobMask == valueT);
                    BlobMask(inds) = valueL; %change all valueT's to valueL's
                    
                else
                    myValue = valueT; %first blob found
                end
            end
            
            if valueTR ~= 0
                if valueTL ~= 0 & valueTL ~= valueTR % U connection
                    inds = find(BlobMask == valueTR);
                    BlobMask(inds) = valueTL;
                else
                    myValue = valueTR; %first blob found
                end
            end
 
            if myValue == 0
                myValue = curBlob;
                curBlob = curBlob + 1;
            end   
            BlobMask(row, col) = myValue;
        end
    end
end


BlobSizes = zeros(curBlob-1, 1);
for i = 1:(curBlob-1)
    m = find(BlobMask == i);
    s = size(m);
    
    BlobSizes(i,1) = s(1);
    
end

imagesc(BlobMask);
axis image;

stats = regionprops(BlobMask, 'Centroid', 'BoundingBox');

Centroids = zeros(k, 2);
for i = 1:k
    [num, Ind] = max(BlobSizes);
    BlobSizes(Ind) = -1; %take it out of the running in next iter
    Centroids(i,:) = stats(Ind).Centroid;
    BB = stats(Ind).BoundingBox;
    x = BB(1);
    y = BB(2);
    xWidth = BB(3);
    yWidth = BB(4);
    
    hold on;
    plot([x; x + xWidth; x + xWidth; x; x], [y; y; y+yWidth; y+yWidth; y]);
    hold off;
end