function clusters = findClusters2D(cfg0, mask)

height = size(mask, 1);
width = size(mask, 2);

N_cluster = 0;
clusters = cell(floor(height*width/2), 1);
clusterid = zeros(height, width);

accounted = ~mask;

while(~all(accounted(:)))
    incluster = zeros(height, width);
    N_cluster = N_cluster + 1;

    [y, x] = find(~accounted, 1);

    while(~isempty(y))
        incluster(y, x) = 1;
        accounted(y, x) = 1;
        
        neighbours_x = (x-1):(x+1);
        neighbours_y = (y-1):(y+1);
        neighbours_x = neighbours_x((neighbours_x > 0) & (neighbours_x <= width));
        neighbours_y = neighbours_y((neighbours_y > 0) & (neighbours_y <= height));

        neighbours = zeros(height, width);
        neighbours(neighbours_y, neighbours_x) = 1;
        incluster = incluster | (neighbours & mask);
                
        [y, x] = find(incluster & ~accounted, 1);
    end
    
    clusters{N_cluster} = incluster;
    
    clusterid(incluster) = N_cluster;
end

clusters = clusters(1:N_cluster);

end