How to identify smaller regions sandwiched between larger regions in watershed image?

Hi, first time posting here so if I miss anything please let me know.

I have a binarized waterhsed image to process. Running regionprops, I have found the centroids, pixelList and areas of all regions in the image. I have sorted all 3 properties according to the area of the regions - area < 100 and area >= 100 (thus centroids and pixelLists have also been split into 2 separate lists/cells). In the image at the bottom that you see, the regions with areas < 100 have their centroids marked in red and areas >= 100 are marked in blue. For simplicity, lets call them blue and red regions.

Onto the challenge I am facing - I want to find the red regions that are sharing a boundary with at least 1 blue region. If the red region successfully shares the boundary with the blue region, then we must store the area, pixelList and centroid data for that region in a new array/struct.

I have tried the following (with little success):

1. bwboundaries to find out the pixel indices of each region, use the pixel indices to find actual pixels and then compare the pixels for red regions and blue regions to find shared boundaries and thus, sort the red regions. The problem i faced here is that the bwboundaries command runs on the entire image containing the large and small regions alike, because of which it is impossible to understand which boundaries belong to the red regions and blue regions (probably not impossible but I could not think of any way here).
2. Using knnsearch to identify the nearest 3 neighbors of a red region in the blue regions and then find if these neighbours share a boundary. The problem here, I guess, is that regionprops.PixelList does not record the boundary pixels of regions which renders the approach redundant.

Here is a snippet of the code including the knnsearch attmempt which failed spectacularly (please note that this is all happening inside a large for loop):

B = bwboundaries(Lb);

% Extract properties of the segmented regions
props = regionprops(Lb, 'Centroid', 'PixelList', 'Area');
areas = cat(1, props.Area);
centroids = cat(1, props.Centroid);
num_centroids = size(centroids, 1);

% Indices of large and small fiber regions
idx_b = find(areas >= 100);
idx_f = find(areas < 100);

% Sort the areas according to idx_f and idx_b
area_f = areas(idx_f);
area_b = areas(idx_b);

% Sort the Centroids according to idx_f and idx_b
centroids_f = cell2mat({props(idx_f).Centroid}');
centroids_b = cell2mat({props(idx_b).Centroid}');

% Sort the PixelList according to idx_f and idx_b
pixelList_f = {props(idx_f).PixelList};
pixelList_b = {props(idx_b).PixelList};

%%%%% Code for the failed knnseaerch attempt %%%%%

% Find 3 nearest bundle centroids for the fibril centroids
Idx = knnsearch(centroids_b, centroids_f, "K",3);

% Initialize a cell array to store the areas, centroids, and pixel lists of shared regions
    sharedRegions = {};

    % Loop through each region in PixelList_f
    for k = 1:length(pixelList_f)
        % Get the current region's pixel list
        pixels_f = pixelList_f{k};

        % Get the nearest neighbors from PixelList_b using indices from Idx
        idx_nn1 = Idx(k, 1);
        idx_nn2 = Idx(k, 2);
        idx_nn3 = Idx(k, 3);

        % Retrieve the corresponding pixel lists from PixelList_b
        pixels_nn1 = pixelList_b{idx_nn1};
        pixels_nn2 = pixelList_b{idx_nn2};
        pixels_nn3 = pixelList_b{idx_nn3};

        % Check for shared pixels between PixelList_f and each neighbor in PixelList_b
        sharedPixels1 = intersect(pixels_f, pixels_nn1, 'rows');
        sharedPixels2 = intersect(pixels_f, pixels_nn2, 'rows');
        sharedPixels3 = intersect(pixels_f, pixels_nn3, 'rows');

        % If there are shared pixels, store the details
        if ~isempty(sharedPixels1) || ~isempty(sharedPixels2) || ~isempty(sharedPixels3)

            sharedRegion = struct();
            sharedRegion.Area = area_f(k);
            sharedRegion.Centroid = centroids_f(k, :);
            sharedRegion.PixelList = pixels_f;

            % Add the sharedRegion to the sharedRegions cell array
            sharedRegions{end+1} = sharedRegion;
        end
    end
 

I have also shared Lb figure (as a maltab fig) in case required to run the code.

[Image Analyst]

Can you post the original image here so we can see it? Like the original PNG file? Then I can possibly try something with it. The .fig file doesn't show up here in the web browser.

In the meantime, in general, to find out what blobs are touching other blobs you can use graycomatrix, the gray level cooccurrence matrix, on a labeled image. I've done this before and I think it will work better than your "knnsearch attmempt which failed spectacularly".

[Pranav]

Hi, I have attached the Lb png image below (Lb_t4.png). Also attached the Lb_blue_red image that contains the centroids plotted. I'll give the graycomatrix a shot. Thanks! :D

[Umar]

Hi @Pranav,

I would be honest when it comes to resolving complex issues regarding image analysis, then definitely follow the advice of @Image Analyst, he is a humble professor with years of technical experience with image analysis. You are in good hands.

[Pranav]

Hi Umar and Image Analyst - thanks for your help on this. I am working on the glcm approach but I'm not able to visualize the solution just yet.

@Image Analyst - would it be okay to bother you once again in case I am unable to solve this problem via the glcm approach?

[Image Analyst]

I'm working on a different way but I'm only part way done. Basically these steps.

1. Get a mask for small blobs and large blobs separately.
2. Dilate the large blobs two layers so it will overlap the small blobs.
3. Determine the overlapping pixels.
4. Use the overlapping pixels as markers with imreconstruct to get the small blobs that are next to the big blobs.
5. AND the result with the original small blobs image to get only the small blobs that are next to large blobs (and not the merged blobs).

I'm pretty busy until tomorrow afternoon or evening when I can continue working on it.

Answer 1

[Umar]

Hi @Pranav,

To tackle the problem of identifying red regions that share boundaries with blue regions in a binarized watershed image, we can adopt a more systematic approach. The goal is to determine the shared boundaries and store the relevant properties of the red regions. Below is a detailed explanation of the steps involved, along with an updated code snippet.

First, utilize regionprops to obtain the centroids, pixel lists, and areas of the regions in the image. Then, seperate the regions based on their areas into two categories: red regions (areas < 100) and blue regions (areas >= 100). Use bwboundaries to extract the boundaries of the entire labeled image. For more information on bwboundaries function, please refer to:

https://www.mathworks.com/help/images/ref/bwboundaries.html

This will help in identifying the boundary pixels of both red and blue regions. For each red region, check if any of its boundary pixels coincide with the boundary pixels of any blue region. If a red region shares a boundary with a blue region, store its area, centroid, and pixel list in a structured format. Here is the revised version of your code that implements the above approach:

% Load a binary image
binaryImage = imread('/MATLAB Drive/IMG_7682.PNG'); % Replace with your   
actual image file
% Convert to binary if necessary
binaryImage = imbinarize(rgb2gray(binaryImage)); % Convert to grayscale 
and then binary

% Create a labeled binary image
Lb = bwlabel(binaryImage);

% Extract properties of the segmented regions
props = regionprops(Lb, 'Centroid', 'PixelList', 'Area');
areas = cat(1, props.Area);
centroids = cat(1, props.Centroid);

% Indices of large and small fiber regions
idx_b = find(areas >= 100);
idx_f = find(areas < 100);

% Sort the areas according to idx_f and idx_b
area_f = areas(idx_f);
area_b = areas(idx_b);

% Sort the Centroids according to idx_f and idx_b
centroids_f = cell2mat({props(idx_f).Centroid}');
centroids_b = cell2mat({props(idx_b).Centroid}');

% Sort the PixelList according to idx_f and idx_b
pixelList_f = {props(idx_f).PixelList};
pixelList_b = {props(idx_b).PixelList};

% Initialize a cell array to store the areas, centroids, and pixel lists of   shared. 
 regions
 sharedRegions = {};

% Loop through each red region
for k = 1:length(pixelList_f)
  % Get the current red region's pixel list
  pixels_f = pixelList_f{k};

    % Get the boundary of the current red region
    boundary_f = bwboundaries(reshape(ismember(Lb, idx_f(k)), size(Lb)), 
   'noholes');
    boundary_f = boundary_f{1}; % Extract the first boundary

    % Check against all blue regions
    for j = 1:length(pixelList_b)
        % Get the boundary of the current blue region
        boundary_b = bwboundaries(reshape(ismember(Lb, idx_b(j)), size(Lb)), 
   'noholes');
        boundary_b = boundary_b{1}; % Extract the first boundary

        % Check for shared boundary pixels
        sharedPixels = intersect(boundary_f, boundary_b, 'rows');

        % If there are shared pixels, store the details
        if ~isempty(sharedPixels)
            sharedRegion = struct();
            sharedRegion.Area = area_f(k);
            sharedRegion.Centroid = centroids_f(k, :);
            sharedRegion.PixelList = pixels_f;
            % Add the sharedRegion to the sharedRegions cell array
            sharedRegions{end+1} = sharedRegion;
            break; % Exit the loop once a shared boundary is found
        end
    end
  end

% Final plot to visualize the results
figure;
imshow(Lb, []);
hold on;

% Plot red regions
for k = 1:length(pixelList_f)
  if ismember(k, idx_f)
      plot(centroids_f(k, 1), centroids_f(k, 2), 'ro', 'MarkerSize', 10);
  end
end

% Plot blue regions
for j = 1:length(pixelList_b)
  if ismember(j, idx_b)
      plot(centroids_b(j, 1), centroids_b(j, 2), 'bo', 'MarkerSize', 10);
  end
end

title('Centroids of Red and Blue Regions');
hold off;

Please see attached.

So, in your updated code, the bwboundaries function is used to extract the boundaries of both red and blue regions which allows to focus on the actual boundary pixels rather than the entire pixel list. The intersect function is employed to determine if there are any common pixels between the boundaries of the red and blue regions. If a shared boundary is found, the relevant properties of the red region are stored in a structured format for further analysis.

If you have any further questions or need additional assistance, feel free to ask!

[Pranav]

Hi Umar, thanks for the response. I understood your approach and it does look solid. However, when I ran the script over an isolated image (to test), the script was able to recognize and populate the boundary_f and boundary_b variables, but the shared_regions struct remained empty.

Additionally, running the script over a single 1024x1024 image for the test I ran took nearly 25 minutes. It might make sense when we see that pixelList_f is a 1x1637 cell and pixelList_b is a 1x545 cell - so looping over these variables might take this long - I could be wrong and would defer to someone more experienced in this regard - but if I am correct, then this approach will not be conducive for larger datasets that I have to eventually analyze.

Given that I am pretty new to matlab, I am quite likely not aware of the most efficient way to solve this problem. Would you be aware of any other approaches that might be more efficient than what I proposed for such an analysis?

Thanks!

Edited - grammar and typos.

[Umar]

Hi @Pranav,

Based on your comments, to address the concerns regarding efficiency and the issue of empty shared regions, I did analyze the code and the best approach to optimize it is to minimize the number of loops and utilize vectorized operations wherever possible. Below is a refined version of your code which incorporates these improvements now.

% Load a binary image
binaryImage = imread('/MATLAB Drive/IMG_7682.PNG'); % Replace with your   
actual image file

% Convert to binary if necessary
binaryImage = imbinarize(rgb2gray(binaryImage)); % Convert to grayscale 
and then binary

% Create a labeled binary image
Lb = bwlabel(binaryImage);

   % Extract properties of the segmented regions
   props = regionprops(Lb, 'Centroid', 'PixelList', 'Area');
   areas = cat(1, props.Area);
   centroids = cat(1, props.Centroid);

   % Indices of large and small fiber regions
    idx_b = find(areas >= 100);
   idx_f = find(areas < 100);

    % Sort the areas according to idx_f and idx_b
    area_f = areas(idx_f);
    area_b = areas(idx_b);

    % Sort the Centroids according to idx_f and idx_b
    centroids_f = cell2mat({props(idx_f).Centroid}');
    centroids_b = cell2mat({props(idx_b).Centroid}');

% Initialize a cell array to store the areas, centroids, and pixel lists of shared   
regions
sharedRegions = {};

% Precompute boundaries for all regions
boundaries_f = cellfun(@(x) bwboundaries(reshape(ismember(Lb, x), 
size(Lb)), 'noholes'), num2cell(idx_f), 'UniformOutput', false);
boundaries_b = cellfun(@(x) bwboundaries(reshape(ismember(Lb, x), 
size(Lb)), 'noholes'), num2cell(idx_b), 'UniformOutput', false);

% Loop through each red region
for k = 1:length(boundaries_f)
  boundary_f = boundaries_f{k}{1}; % Extract the first boundary of the red   
region

    % Check against all blue regions
    for j = 1:length(boundaries_b)
        boundary_b = boundaries_b{j}{1}; % Extract the first boundary of the blue   
  region

        % Check for shared boundary pixels
        sharedPixels = intersect(boundary_f, boundary_b, 'rows');

        % If there are shared pixels, store the details
        if ~isempty(sharedPixels)
            sharedRegion = struct();
            sharedRegion.Area = area_f(k);
            sharedRegion.Centroid = centroids_f(k, :);
            sharedRegion.PixelList = props(idx_f(k)).PixelList; % Use the original 
            pixel list
            % Add the sharedRegion to the sharedRegions cell array
            sharedRegions{end+1} = sharedRegion;
            break; % Exit the loop once a shared boundary is found
        end
    end
  end

% Final plot to visualize the results
figure;
imshow(Lb, []);
hold on;

% Plot red regions
plot(centroids_f(:, 1), centroids_f(:, 2), 'ro', 'MarkerSize', 10);

% Plot blue regions
plot(centroids_b(:, 1), centroids_b(:, 2), 'bo', 'MarkerSize', 10);

title('Centroids of Red and Blue Regions');
hold off;

% Display shared regions for verification
disp('Shared Regions:');
disp(sharedRegions);

Summary of changes made to the code

Boundary Precomputation: The boundaries for both red and blue regions are computed once and stored in a cell array. This reduces the need to repeatedly compute boundaries within nested loops.

Vectorization: The use of cellfun allows for more efficient processing of the boundaries, minimizing the overhead of manual loops.

Early Exit: The inner loop breaks as soon as a shared boundary is found, preventing unnecessary checks against remaining blue regions.

Display of Shared Regions: The code now includes a display of shared regions for verification, which can help in debugging and understanding the results.

So, this revised approach should reduce the processing time, especially for larger datasets, while making sure that shared regions are accurately identified.

Hope this should help resolve your problem now.

If you have any further questions or need additional assistance, please feel free to ask!

[Pranav]

Thanks a ton Umar, but even though the boundaries_f and boundaries_p cells are being populated, (subsequent boundary_f and boundary_p variables too get data in them) - the shared_regions variable remains empty at the end of it all. I believe I will have to work with the glcm approach to solve this issue.

[Umar]

Hi @Pranav,

After going through documentation for glcm documentation shared in the link below

https://www.mathworks.com/help/images/ref/graycomatrix.html

It seems that this approach will enrich your analysis by providing texture-based insights, the immediate issue of empty sharedRegions may be better resolved by refining the boundary detection methods or adjusting the segmentation thresholds. By incorporating visual validation and potentially dilating boundaries, you can enhance the accuracy of shared pixel identification. If you decide to proceed with GLCM, it will complement your existing methods by offering a deeper understanding of the textures present in your segmented regions.