How to do correctly this for loop for segmentation?

Dear all,

I have code for segmentation specific parts. Output from this code is this image:

I would like to gradually in this manner do segmentation the other parts. I think, that another for loop could help. I tried this one:

*for k = 1:16*
cnt=0;
for i=1:numel(tmp)
    if numel(tmp{i})>4
    cnt=cnt+1;
    tmp2{cnt}=tmp{i};
   end
end
tmp = tmp2;
*end*

But it didn´t solve my problem. I would like to for loop, which works like this: for k = 1 this output:

for k = 2 This output:

and so on. I hope that you understand me. I attach my code and original image. Thank you for your answers.

[Image Analyst]

I don't really understand. What regions do you want? So many times in the past I've done things like finding lungs, skulls, tumors, etc. Search my old posts.

[Veronika]

I need segmentation like this

but bloob after bloob. Because I need to filter found objects by size.

Answer 1

[Image Analyst]

Here's the plan (which I can't do because you didn't post the original, non-annotated image):

1. Threshold to find body.
2. Use bwareafilt() to extract largest blob.
3. Use that to mask out non-body stuff in the gray level image.
4. Threshold again at a higher gray level
5. Use bwareaopen() ror bwareafilt() to get only the blobs of the size you want.
6. Call bwboundaries() to get the perimeters of the blobs.
7. Use a for loop and plot() to plot red outlines around the blobs.

Good luck. If you need more help, insert your original image and post your code.

[Veronika]

Okay, I attach original image and code. I think, that I could use this loop

cnt=0;
for i=1:numel(tmp)
    if numel(tmp{i})>4
    cnt=cnt+1;
    tmp2{cnt}=tmp{i};
   end
end
tmp = tmp2;

for my output. Because this for loop I used for one yellow segment i my output, so I thought, that this could be the right way. Thank you for your help.

Answer 2

[Image Analyst]

I didn't understand your comments, so I just wrote the whole thing myself.

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 15;

%===============================================================================
% Get the name of the image the user wants to use.
baseFileName = 'thorax-mdl.jpg'; % Assign the one on the button that they clicked on.
% Get the full filename, with path prepended.
folder = []; % Determine where demo folder is (works with all versions).
fullFileName = fullfile(folder, baseFileName);

%===============================================================================
% Read in a demo image.
grayImage = imread(fullFileName);
% Get the dimensions of the image.  
% numberOfColorChannels should be = 1 for a gray scale image, and 3 for an RGB color image.
[rows, columns, numberOfColorChannels] = size(grayImage);
if numberOfColorChannels > 1
  % It's not really gray scale like we expected - it's color.
  % Use weighted sum of ALL channels to create a gray scale image.
  grayImage = rgb2gray(grayImage); 
  % ALTERNATE METHOD: Convert it to gray scale by taking only the green channel,
  % which in a typical snapshot will be the least noisy channel.
  % grayImage = grayImage(:, :, 2); % Take green channel.
end
% Display the image.
subplot(2, 3, 1);
imshow(grayImage, []);
axis on;
caption = sprintf('Original Gray Scale Image');
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
drawnow;
hp = impixelinfo();

% Let's compute and display the histogram.
[pixelCount, grayLevels] = imhist(grayImage);
subplot(2, 3, 2); 
bar(grayLevels, pixelCount); % Plot it as a bar chart.
grid on;
title('Histogram of original image', 'FontSize', fontSize, 'Interpreter', 'None');
xlabel('Gray Level', 'FontSize', fontSize);
ylabel('Pixel Count', 'FontSize', fontSize);
xlim([0 grayLevels(end)]); % Scale x axis manually.

% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Get rid of tool bar and pulldown menus that are along top of figure.
set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off') 
drawnow;

% Binarize the image by thresholding.
bodyMask = grayImage > 128;
% Get rid of white surround.
bodyMask = imclearborder(bodyMask);
% Extract the largest blob only.
bodyMask = bwareafilt(bodyMask, 1);
subplot(2, 3, 3);
imshow(bodyMask);
axis on;
title('Binary Image of Whole Body', 'fontSize', fontSize);
drawnow;

% Find the ribs mask by thresholding and ANDing with the body mask
% to make sure we don't get anything from inside the lungs.
ribsMask = bodyMask & (grayImage > 220);
% Fill holes.
ribsMask = imfill(ribsMask, 'holes');
% Find areas.
labeledImage = bwlabel(ribsMask);
props = regionprops(labeledImage, 'Area');
allAreas = sort([props.Area])
% Get rid of areas below 100 pixels.
ribsMask = bwareaopen(ribsMask, 100);
% Display the image.
subplot(2, 3, 4);
imshow(ribsMask, []);
axis on;
caption = sprintf('Ribs Mask');
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
drawnow;

% bwboundaries() returns a cell array, where each cell contains the row/column coordinates for an object in the image.
% Plot the borders of all the coins on the original grayscale image using the coordinates returned by bwboundaries.
subplot(2, 3, 5);
imshow(grayImage);
title('Outlines, from bwboundaries()', 'FontSize', fontSize); 
axis image; % Make sure image is not artificially stretched because of screen's aspect ratio.
hold on;
boundaries = bwboundaries(ribsMask);
numberOfBoundaries = size(boundaries, 1);
for k = 1 : numberOfBoundaries
  thisBoundary = boundaries{k};
  plot(thisBoundary(:,2), thisBoundary(:,1), 'r', 'LineWidth', 2);
end
hold off;

[Veronika]

This your segmentation I already have, but it dosn´t solve my problem. I would like to create FEM model from this segmentation. But some blobs are too small (1 pixel) for reconstruction, so I think, that I have to choose only bloobs bigger than 1 pixel. This is my code:

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 15;

%===============================================================================
% Get the name of the image the user wants to use.
baseFileName = 'thorax-mdl.jpg'; % Assign the one on the button that they clicked on.
% Get the full filename, with path prepended.
folder = []; % Determine where demo folder is (works with all versions).
fullFileName = fullfile(folder, baseFileName);

%===============================================================================
% Read in a demo image.
grayImage = imread(fullFileName);
% Get the dimensions of the image.  
% numberOfColorChannels  = 1; 
[rows, columns, numberOfColorChannels] = size(grayImage);
if numberOfColorChannels > 1
  % It's not really gray scale like we expected - it's color.
  % Use weighted sum of ALL channels to create a gray scale image.
grayImage = rgb2gray(grayImage); 
  % ALTERNATE METHOD: Convert it to gray scale by taking only the green channel,
  % which in a typical snapshot will be the least noisy channel.
%   grayImage = grayImage(:, :, 2); % Take green channel.
end
% Display the image.
eq_grayImage = histeq(grayImage);
subplot(2, 3, 1);
imshow(grayImage, []);
axis on;
caption = sprintf('Original Gray Scale Image');
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
drawnow;
hp = impixelinfo();

% Let's compute and display the histogram.
[pixelCount, grayLevels] = imhist(grayImage);
subplot(2, 3, 2); 
bar(grayLevels, pixelCount); % Plot it as a bar chart.
grid on;
title('Histogram of original image', 'FontSize', fontSize, 'Interpreter', 'None');
xlabel('Gray Level', 'FontSize', fontSize);
ylabel('Pixel Count', 'FontSize', fontSize);
xlim([0 grayLevels(end)]); % Scale x axis manually.

% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Get rid of tool bar and pulldown menus that are along top of figure.
set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off') 
drawnow;

% Binarize the image by thresholding.
bodyMask = grayImage > 128;
% Get rid of white surround.
bodyMask = imclearborder(bodyMask);
% Extract the largest blob only.
bodyMask = bwareafilt(bodyMask, 1);
subplot(2, 3, 3);
imshow(bodyMask);
axis on;
title('Binary Image of Whole Body', 'fontSize', fontSize);
drawnow;

% Find the ribs mask by thresholding and ANDing with the body mask
% to make sure we don't get anything from inside the lungs.
ribsMask = bodyMask & (grayImage > 220);
% Fill holes.
ribsMask = imfill(ribsMask, 'holes');
% Find areas.
labeledImage = bwlabel(ribsMask);
props = regionprops(labeledImage, 'Area');
allAreas = sort([props.Area])
% Get rid of areas below 100 pixels.
ribsMask = bwareaopen(ribsMask, 100);
% Display the image.
subplot(2, 3, 4);
imshow(ribsMask, []);
axis on;
caption = sprintf('Ribs Mask');
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
drawnow;

% bwboundaries() returns a cell array, where each cell contains the row/column coordinates for an object in the image.
% Plot the borders of all the coins on the original grayscale image using the coordinates returned by bwboundaries.
figure(4)
imshow(grayImage);
title('Outlines, from bwboundaries()', 'FontSize', fontSize); 
axis image; % Make sure image is not artificially stretched because of screen's aspect ratio.
hold on;
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Odstranění panelu nástrojů a rozbalovacího menu.
set(gcf, 'Toolbar', 'none', 'Menu', 'none');
set(gcf, 'Name', 'Segmentace páteře', 'NumberTitle', 'Off') 
boundaries = bwboundaries(ribsMask);
numberOfBoundaries = size(boundaries, 1);
for k = 1 : numberOfBoundaries
  thisBoundary = boundaries{k};
  plot(thisBoundary(:,2), thisBoundary(:,1), 'r', 'LineWidth', 2);
end 
hold off;

%Práh pro vytvoření binárního obrazu okolí
thresholdValue = 180;
binaryImage_okoli = eq_grayImage > thresholdValue;
% Odstranění okolí.
binaryImage_okoli = imclearborder(binaryImage_okoli);
% Vyplnění otvorů.
binaryImage_okoli = imfill(binaryImage_okoli, 'holes');
% Vymazání menších otvorů.
binaryImage_okoli = bwareaopen(binaryImage_okoli, 750);

figure(3)
subplot (2, 3, 1)
imshow(binaryImage_okoli, []);
axis on;
caption = sprintf('Binární obraz okolí');
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');

run C:/Users/ced0005/eidors-v3.8-ng/eidors/startup.m
if 2==exist ('Pater_model.mat')
load ('Pater_model')
disp('Model páteře byl načten.')
else

shape = [];
tmp = bwboundaries(binaryImage_okoli);
shape.thorax = tmp{1};
tmp = bwboundaries(ribsMask);
shape.spine = tmp{1};
size(grayImage)   

% cnt=0;
% for i=1:numel(tmp)
%     if numel(tmp{i})>4
%     cnt=cnt+1;
%     tmp2{cnt}=tmp{i};
%    end
% end
% tmp = tmp2;

% interpolace hraničních bodů, fourier
Npts1 = 35;      % počet bodů interpolovaných hranic
Npts2 = 25;      % počet bodů interpolovaných hranic
Npts3 = 25;       % počet bodů interpolovaných hranic
Ncomps = 20;     % počet fourierových komponentů

shape.thorax = fourier_fit(fourier_fit(shape.thorax, Ncomps), linspace(0, 1, Npts1));
shape.spine = fourier_fit(fourier_fit(shape.spine, Ncomps), linspace(0, 1, Npts3));

figure(7);
imshow(grayImage);
hold on
plot(shape.thorax(:,2), shape.thorax(:,1), 'o-b')
plot(shape.spine(:,2), shape.spine(:,1), 'o-y')
legend('Okolí hrudníku', 'Páteř + žebra');

As you can see, if I put as input ribsMask, the output of function fourier_fit is only one rib, but I need to all ribs and spine. Can you advise me?

[Image Analyst]

This code:

% Get rid of areas below 100 pixels.
ribsMask = bwareaopen(ribsMask, 100);

makes sure you don't have blobs that are 1 pixel. I have no idea how they're surviving that line of code. They should all be removed. You might need to call the Mathworks and ask them why bwareaopen() is not working for you. Alternatively you can use

ribsMask = bwareafilt(ribsMask, [100, inf]);

to extract only those blobs bigger than or equal to 100 pixels.

[Veronika]

Okay, but output is still the same, I think, that this part of code could be helpful for this problem:

% cnt=0;
% for i=1:numel(tmp)
%     if numel(tmp{i})>4
%     cnt=cnt+1;
%     tmp2{cnt}=tmp{i};
%    end
% end
% tmp = tmp2;

But I don´t know, how this code use for all blobs.

[Image Analyst]

I don't know what that code does. What is its purpose? And I don't understand what the code after my code does because I can't read that language. So either convert to English so I can read and understand what's going on, or try yourself by just adjusting the parameters to bwareaopen() or all those other parameters you introduced like Ncomps, etc.

[Veronika]

I tried change the parameters, but without success. I attach english version of my code and original image. I hope, that it will be enough understanding for you.