Most efficient way to separate numerically labeled objects that share borders?
Show older comments
Consider the following image mask:
IM = [0 0 0 1 1 1 1
0 0 1 1 1 1 1
0 2 2 2 1 1 1
2 2 2 2 2 1 0
2 2 2 2 2 0 0];
Objects 1 and 2 share a diagonal border. If I converted this image as shown above into a logical and passed it to regionprops, it would detect only 1 (merged) object, not 2.
I'm looking for the most computationally efficient way to separate these objects. The end result would look something like this:
want = [0 0 0 1 1 1 1
0 0 0 0 1 1 1
0 2 2 0 0 1 1
2 2 2 2 0 0 0
2 2 2 2 2 0 0];
...where regionprops(logical(want)) would result in 2 objects being detected.
The reason I'm asking is that my actual masks are 25000 x 25000 in size, with about 100,000 objects in each. Separating them would be more efficient in terms of both storage and downstream processing. Also, I'm aware that separating these objects will erode them a bit, and a secondary goal is to minimize the amount of erosion and preserve their overall structures, as some of the actual objects can be as little as 10 pixels wide (they are narrow cells).
So far I've thought of some ways to do this, including brute-forcing it by systematically checking each pixel and converting it to a 0 if one of its neighbors is anything other than a 0 or itself, e.g:
IM2 = padarray(IM, [1 1], 0, 'both'); % make some space
for i = 2:size(IM2, 1)-1
for j = 2:size(IM2, 2)-1
% get that pixel and its neighbors
w = IM2(i-1:i+1, j-1:j+1);
% if bordering another object, turn it to 0
if any(~ismember(w(:), [0 IM2(i, j)]))
IM2(i, j) = 0;
end
end
end
out = IM2(2:end-1, 2:end-1); % remove padding
While this method technically works, it's time-consuming for a large image (about 6 hours per 25000 x 25000) and I have quite a few of these images...
I've also tried various other methods, including isolating individual objects, doing morphological erosion on each, and then mapping the results back on the original image, but this is turned out to be far worse than the brute-force method in terms of both computation time and structural preservation.
Does anyone here have experience with efficiently separating these kinds of objects, which are numerically labeled and could be adjacent to each other, such that they can be converted and stored as logical images without losing too much of their morphological properties?
2 Comments
Ran Yang
on 10 Apr 2023
Answers (2)
---- Update on 2023-04-10 ----
Mostly figured this out. I just had to isolate each object in its bounding box instead of running imdilate across the entire image N times (code updated below).
load('sample.mat') % IM = sample 1000x1000 image
figure, imshow(label2rgb(IM, 'parula', 'k')); title('Adjacent objects');
IM2 = padarray(IM, [1 1], 0, 'both');
% reorder objects 1,2,3... (can comment this out if they are already in order)
IM2 = uint32(categorical(IM2)); IM2 = IM2 - 1;
b = false(size(IM2)); % allocate master boundary mask
R = regionprops(IM2);
for i = 1:size(R, 1)
% isolate each object in its boundingbox
bb = R(i).BoundingBox;
obj = IM2(floor(bb(2)):ceil(bb(2)+bb(4)), floor(bb(1)):ceil(bb(1)+bb(3)));
% remove other objects in this boundingbox
obj(obj ~= i) = 0;
% determine object perim in its boundingbox
obj = imdilate(obj, strel('disk', 1)) - obj;
% locate corresponding bbox in b and append to it
b(floor(bb(2)):ceil(bb(2)+bb(4)), floor(bb(1)):ceil(bb(1)+bb(3))) = ...
b(floor(bb(2)):ceil(bb(2)+bb(4)), floor(bb(1)):ceil(bb(1)+bb(3))) | obj;
end
out = IM2; out(b) = 0; % remove all boundaries
out = out(2:end-1, 2:end-1); % remove padding
figure, imshow(label2rgb(out, 'parula', 'k')); title('Separated objects');
This method is really fast. On the original 25000 x 25000 image with 120,000 objects, it finishes in ~90 seconds (as opposed to my original method, which took 40+ hours).
The only issue now is what is the best algorithm to preserve as much of each object as possible (or as Image Analyst pointed out in a comment, how to determine boundaries of adjacent objects). The above code uses imdilate to define a border around each object as an efficient approximation. This will erode away some objects if they're too small, though I currently don't see a way around sacrificing some morphology.
---- Original idea ----
After giving this some more thought, I've come up with a faster solution, which is to just calculate and remove the perimeters of each object. Method 1 does this via bwperim, while method 2 uses a dilated object minus the original.
% u = nonzeros(unique(IM(:)));
% b = false(size(IM)); % allocate boundary mask
% for i = 1:size(u, 1)
% % isolate each object
% obj = IM == u(i);
% % method 1: find object inner perimeter and append to b
% b(bwperim(obj)) = true;
% % method 2: find object outer perimeter and append to b
% % b(logical(imdilate(obj, strel('disk', 1)) - obj)) = true;
% end
% out = IM;
% out(b) = 0; % split objects by setting their boundaries to 0
Method 1 takes ~6% of the time of my initial brute-force method, but indescriminately erodes every object by a pixel, which might be an issue for smaller objects.
Method 2 is a little slower (~14%), but is better at preserving small objects.
You don't need bwlabel, like KSSV suggested, because your matrix is already labeled. You can simply use
oneBlob = ismember(IM, numberYouWant);
I don't see why bwperim() would not preserve small objects. I guess I'd need to see an example of that to prove it, like
mask = false(4,4);
mask(2, 2:3) = true
perims = bwperim(mask) % See? All preserved
8 Comments
If two regions are adjacent regionprops will give measurements for each assuming you pass it the labeled image. However when you take the labeled image and cast it to logical, all the labels are lost and you have only 0 and 1 so the distinction between blob #1 and blob #2 is lost. So, DON'T DO THAT! Here is fixed code:
IM = repmat([1 1 1 2 2 2 0 0 3 3 3], 4, 1) % 3 objects in this 'image'
% Wrong way:
props = regionprops(logical(IM)); % DO NOT cast a labeled image to binary
size(props) % detects only 2 objects because 1 and 2 are fused (bad!)
% Correct way:
props = regionprops(IM);
size(props) % detects all 3 objects even though 1 and 2 are fused/
Ran Yang
on 10 Apr 2023
Image Analyst
on 11 Apr 2023
If two regions are adjacent, one with a lower number and one with a higher number, which pixel boundary do you want to set to zero?
Ran Yang
on 11 Apr 2023
Image Analyst
on 11 Apr 2023
Like I said, if you cast your labeled matrix to logical, then all your labeled blobs coalesce into one giant blob and regionprops will only give you measurements on the giant blob, not individual blobs.
I really don't know why you think you need to split apart the blobs to call regionprops rather than just put in IM directly and get the measurements for all the blobs at once. Not really sure what you want to measure but if you want the area and perimeter, for all hundred or so blobs, you can simply do this:
props = regionprops(IM, 'table', 'Area', 'Perimeter')
No need to split them apart and measure each one at a time.
Ran Yang
on 11 Apr 2023
Categories
Find more on Region and Image Properties in Help Center and File Exchange
on 11 Apr 2023
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
