Why this demo isn´t working for my image?

Dear all,

I tried use this demo by Image Analyst, but for my image I wasn´t successful. Instead of cameraman.tif I gave my image thorax-mdl.jpg, which I attached. I have this error message in Command Window:

Error using imageDisplayValidateParams>validateCData (line 121)
If input is logical (binary), it must be two-dimensional.

Error in imageDisplayValidateParams (line 31)
common_args.CData = validateCData(common_args.CData,image_type);

Error in imageDisplayParseInputs (line 79)
common_args = imageDisplayValidateParams(common_args);

Error in imshow (line 198)
  [common_args,specific_args] = ...

Error in DP (line 86)
imshow(binaryImage, []);

Can anyone give me advice? Thank you for your answers.

Answer 1

[Image Analyst]

You forgot to attach your code. You must have modified a demo of mine somehow and made a new script called DP.m. Evidently the code is expecting a binary image - a 2D logical image - and you're not giving it that. Maybe you're giving it a color image or something. But I can't tell because you didn't supply your code modifications. Please attach your image and DP.m and anything else we need to run your code.

[Veronika]

Ok. I attach my code and image. Thank you for your help.

[Image Analyst]

As I suspected, your image is not really gray scale - it's color. When you read it in, convert it to gray scale with this code:

grayImage = imread(fullFileName);
% Dostaneme velikost obrazu. numberOfColorBands by mìlo být = 1.
% 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.
  % 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

[Veronika]

Still I have the same error message...

[Image Analyst]

Veronika, you didn't make the changes I said. You modified them again. Look, nowhere in my code , that I suggested, do I have these lines:

numberOfColorChannels = 1;
[rows, columns, numberOfColorBands] = size(grayImage);

The number of color channels/bands is not even named consistently.

This time, I'm attaching the entire file that definitely runs to completion.

[Veronika]

I don´t know, what I´m doing wrong...But I copied your test3.m and again the same error.

[Image Analyst]

Before the line where it throws the error, put these lines

whos binaryImage
[rows, columns, numberOfColorChannel] = size(binaryImage)

It will report those values to the command window. Tell me what you see in the command window.

[Veronika]

Command Window:

Name               Size                Bytes  Class      Attributes

    binaryImage      512x512x3            786432  logical

[Image Analyst]

Sorry, somehow I think I had attached your original code. This code should run.

[Veronika]

Thank you.

[Veronika]

I wanted to even ask, if I choose only one part from my image (for example one lung)and only this part transfer to binarization image.Thank you for your answer.

[Image Analyst]

I don't understand the question.

[Veronika]

Now, I have for example this segmentation (I attached as plíce cévy okolí.fig) from your test3.m. And I would like to segment for example only one lung.

[Image Analyst]

Several ways. I'd probably threshold for the bright body. Then call imclearborder() to get rid of the air and table surround. Then you'll have two blobs. I'd call imfill() to fill holes. Then I'd use ismember() to extract blob #1 (left) or blob #2 (for the right side). Untested code:

binaryImage = grayImage < someValue;
binaryImage = imclearborder(binaryImage);
binaryImage = imfill(binaryImage, 'holes');
labeledImage = bwlabel(binaryImage);
leftLung = ismember(labeledImage, 1);
rightLung = ismember(labeledImage, 2);

[Veronika]

I tried to incorporate into test3.m, it is better, but still there are both lungs.

I attach code.

[Image Analyst]

Here is code to extract each lung:

% Code to segment out each lung from a CT cross sectional image.
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 = pwd
fullFileName = fullfile(folder, baseFileName);

%===============================================================================
% Read in a demo image.
grayImage = imread(fullFileName);
% Get the dimensions of the image.  
% numberOfColorBands should be = 1.
[rows, columns, numberOfColorChannels] = size(grayImage);
if numberOfColorChannels > 1
  % It's not really gray scale like we expected - it's color.
  % Convert it to gray scale by taking only the green channel.
  grayImage = grayImage(:, :, 2); % Take green channel.
end
% Display the image.
subplot(2, 3, 1);
imshow(grayImage, []);
axis on;
caption = sprintf('Original Grayscale Image, %s', baseFileName);
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
drawnow;
hp = impixelinfo();

% 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') 

% Display the histogram so we can see what gray level we need to threshold it at.
subplot(2, 3, 2:3);
histogram(grayImage, 256);
grid on;
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');

% Threshold the image to make a binary image.
thresholdValue = 170;
binaryImage = grayImage < thresholdValue;
% Get rid of surround.
binaryImage = imclearborder(binaryImage);
% Fill holes.
binaryImage = imfill(binaryImage, 'holes');
% Display the image.
subplot(2, 3, 4);
imshow(binaryImage, []);
axis on;
caption = sprintf('Binary Image');
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');

labeledImage = bwlabel(binaryImage);
leftLung = ismember(labeledImage, 1);
rightLung = ismember(labeledImage, 2);

% Display the image.
subplot(2, 3, 5);
imshow(leftLung, []);
axis on;
caption = sprintf('Left Lung');
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');

% Display the image.
subplot(2, 3, 6);
imshow(rightLung, []);
axis on;
caption = sprintf('Right Lung');
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');

[Veronika]

Thank you for your previous answer.I would like to remove this part of segmentation.

It´s a pad and it isn´t necessar for my purpose. Thank you so much.

[Image Analyst]

You can call bwareaopen() to remove blobs smaller than a specified size. Or you can call bwareafilt() to extract the largest blob only.

[Veronika]

I tried bwaropen, but nothing changed.

binaryImage_okoli = bwareaopen(binaryImage_okoli, 150);
binaryImage_pater = bwareaopen(binaryImage_pater, 150);

Namely, these binary images: binaryImage_pater and binaryImage_okoli. I attach my code again.

[Veronika]

So it changes (concretely binaryImage_pater), but there are two parts, that I removed too, but I didn´t want to. Is there any chance to segmented this parts without segmentation unwanted part?

Thank you for your answer.

[Image Analyst]

The areas are body=141328, table = 708. So your 150 in bwareaopen() was too small to remove the 708 area blob. Use a bigger number like 2000, or else use bwareafilt() to extract only the largest blob. Find this in your code and replace it with the code I've attached.

%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í otvorù menších jak 50 pixelù.
binaryImage_okoli = bwareaopen(binaryImage_okoli, 2000);

% Find out the areas that we have in these blobs
% Uncomment if you want to see the areas.
% [labeledImage, numBlobs] = bwlabel(binaryImage_okoli);
% props = regionprops(labeledImage, 'Area');
% allAreas = [props.Area]

% Extract the larest blob only.
binaryImage_okoli = bwareafilt(binaryImage_okoli, 1);

[Veronika]

I´m so sorry, I was badly expressed, again...It works for binaryImage_okoli, but not for binaryImage_pater. I attach code only for binaryImage_pater, which I try. You can see, that there are 4 parts (blobs) not segmented in the output image. I can´t use function bwareafilt, because I have Matlab 2011b.

[Image Analyst]

I'm having a little trouble following what you're doing because of the language difference. The code runs and does not throw any errors. And the table at the bottom is not there. I thought removing the table at the bottom was the problem, but it's not there. So what is the problem?

[Veronika]

I have this error:

Undefined function 'bwareafilt' for input arguments of type 'double'.
binaryImage_okoli = bwareafilt(binaryImage_okoli, 1);

Because I have Matlab2011b version. I think I need some other solution.

[per isakson]

bwareafilt was introduced in R2014b

[Image Analyst]

For a replacement to bwareafilt(), see my attached function/demo.

[Veronika]

I try this your demo, but this error appears:

Undefined function 'ExtractNLargestBlobs' for input arguments of type 'double'.

Error in ExtractBiggestBlob (line 94)
biggestBlob = ExtractNLargestBlobs(binaryImage, numberToExtract);

[Image Analyst]

I just ran it to double check. It works fine. How did you modify it? It looks like you just copied the first function in the file and did not copy all the functions. The ExtractNLargestBlobs() function is definitely in there. It's further down the file. I think you didn't save the entire file but just copied a portion of it. Please attach the code you actually rand with the paper clip icon so I can verify that the ExtractNLargestBlobs() function is missing from it. Then copy the whole thing this time and run that.

[Veronika]

Yes, function ExtractNlargestBlobs() works fine. I forgot copied end of this function. But I have another problem.

Error using iptcheckinput
Function IMHIST expected its first input, I or X, to be two-dimensional.

Error in imhist>parse_inputs (line 277)
iptcheckinput(a, {'double','uint8','int8','logical','uint16','int16','single','uint32', 'int32'}, ...

Error in imhist (line 59)
[a, n, isScaled, top, map] = parse_inputs(varargin{:});

Error in ExtractBiggestBlob (line 41)
[pixelCount, grayLevels] = imhist(grayImage);

I attach my code.

[Image Analyst]

Your thorax_mdl image is color, not gray scale. Put this code in there after imread():

[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

[Veronika]

Ok. But I can´t use bwareafilt function obviously. Because of this error message:

Undefined function 'bwareafilt' for input arguments of type 'double'.

Error in ExtractBiggestBlob (line 101) biggestBlob = bwareafilt(binaryImage, abs(numberToExtract), sizeOption);

[Image Analyst]

Look at the comments and the code:

%---------------------------------------------------------------------------
% Extract the largest area using our custom function ExtractNLargestBlobs().
% This is the meat of the demo!
% biggestBlob = ExtractNLargestBlobs(binaryImage, numberToExtract);
%
% Now, with R2015a release, you can use bwareafilt() instead of my custom ExtractNLargestBlobs() function.
biggestBlob = bwareafilt(binaryImage, abs(numberToExtract), sizeOption);
%---------------------------------------------------------------------------

Since you do not have R2015a or later, you can't use that so you'll have to comment out this line:

biggestBlob = bwareafilt(binaryImage, abs(numberToExtract), sizeOption);

and uncomment this line:

biggestBlob = ExtractNLargestBlobs(binaryImage, numberToExtract);

I'm attaching updated code to make the selection automatic by checking the version number.

[Veronika]

It works. But I would like to combine the smallest and the largest together for the best result. Is it possible?

[Image Analyst]

If they're not connected, how would you like to connect them - the body to the smallest blob (which may be a speck of noise)? Would you like to burn a white line into the image from one centroid to another? And WHY do you want to connect the largest blob to the smallest blob?

[Veronika]

I would like to combine the largest and the smallest blob, because I would like to remove part of segmented image, which I attached previously. And it wasn´t successful, when I used only the largest or only the smallest blob.

[Image Analyst]

I don't know what you mean. Point to the largest blob, and point to the smallest blob, and tell me how you want to connect them, because I'm having trouble visualizing what blobs you're talking about. The largest blob will be the body, and the smallest one will probably be some small speck of noise, like a single pixel, so I'm not sure why you want to keep the smallest one. Or do you want the 2 largest blobs, like the body and the table? But I thought you wanted to remove the table.

So, in addition to the 2 blobs that you want to keep, show me the part that you "would like to remove part of segmented image" - I'm not sure what part that is.

[Veronika]

I would like to send it, but the icone (clip) disappeared and I don´t know why...so nothing I can attach...

[Veronika]

Yes, I want to remove the table, but for different segmentation. This segmentation:

I want to remove parts bordered in blue. I attach also the code and original image. Thank you.

[Image Analyst]

It's always going to be in the same place, so just simply erase it:

grayImage(someRow:end, :) = 0;

You will know what someRow is, and it will stay the same for all images.

[Veronika]

So, for example I choose row number 500, but it doesn´t solve my problem. Code:

      % Kód pro segmentaci jednotlivých orgánů.
      clc;    
      close all;  
      clear;  
      workspace;  
      format long g;
      format compact;
      fontSize = 15;

      % Zvolení souboru, který chceme použít.
      baseFileName = 'thorax-mdl.jpg'; 
      folder = pwd
      fullFileName = fullfile(folder, baseFileName);

      % Načtení obrazu.
      grayImage = imread(fullFileName);
      % Dimenze obrazu. 
      % numberOfColorBands by měl být = 1.
      [rows, columns, numberOfColorChannels] = size(grayImage);
      if numberOfColorChannels > 1
        % Máme barevný obraz, musíme ho převést na černobílý = vybereme zelený kanál
        grayImage = grayImage(:, :, 2); % zelený kanál
          grayImage(500:end, :) = 0;
      end
      eq_grayImage = histeq(grayImage);%ekvalizace pomocí histogramu obrazu
      % Zobrazení obrazu.
      subplot(2, 3, 1);
      imshow(grayImage, []);
      axis on;
      caption = sprintf('Originální černobílý obraz, %s', baseFileName);
      title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
      drawnow;
      hp = impixelinfo();
      % Nastavení obrazu.
      % Zvětšení obrazu na celou obrazovku.
      set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
      % Odstranění panelu nástrojů a rozbalovacího menu.
      set(gcf, 'Toolbar', 'none', 'Menu', 'none');
      % Přidání titulku okna.
      set(gcf, 'Name', 'Segmentace', 'NumberTitle', 'Off') 
      % Zobrazení histogramu pro lepší určení prahu.
      subplot(2, 3, 2)
      imhist(grayImage, 256);
      grid on;
      title(caption, 'FontSize', fontSize, 'Interpreter', 'None');

      % Práh pro vytvoření binárního obrazu páteře.
      thresholdValue = 220;
      binaryImage_pater = grayImage > thresholdValue;
      % Odstranění okolí.
      binaryImage_pater = imclearborder(binaryImage_pater);
      % Vyplnění otvorů
      binaryImage_pater = imfill(binaryImage_pater, 'holes');
      %binaryImage_pater = bwareaopen(binaryImage_pater, 200);
      %Zobrazení segmentace páteře
      imshow(binaryImage_pater, []);
      axis on;
      caption = sprintf('Binární obraz páteře');
      title(caption, 'FontSize', fontSize, 'Interpreter', 'None');

      figure(4)
      imshow(grayImage, []);
      title('Segmentace páteře', 'FontSize', fontSize); 
      axis image; % Ujištění, že se obraz po zvětšení okna nezdeformuje.
      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');
      % Přidání titulku okna.
      set(gcf, 'Name', 'Segmentace páteře', 'NumberTitle', 'Off') 

      %segmentace páteře
      boundaries = bwboundaries(binaryImage_pater);
      numberOfBoundaries = size(boundaries, 1);
      for k = 1 : numberOfBoundaries
        thisBoundary = boundaries{k};
        plot(thisBoundary(:,2), thisBoundary(:,1), 'r', 'LineWidth', 2);
      end

      % Vykreslení značek na hranicích.
      skip = 20;
      for k = 1 : numberOfBoundaries
        thisBoundary = boundaries{k};
        x_pater = thisBoundary(1 : skip : end,2);
        y_pater = thisBoundary(1 : skip : end,1);
        plot(x_pater, y_pater, 'ro', 'LineWidth', 2, 'MarkerSize', 10);
      end

[Image Analyst]

How could it still be finding the table if you erased it???

[Veronika]

I erase the table for segmentation binaryImage_okoli, but not for binaryImage_pater. Here: binaryImage_okoli

That´s right for this segmentation I removed the smallest blob.

binaryImage_pater

But for this segmentation is not so easy. Because, when I removed the table (blob) I also removed some blobs, which I want to. So I used

grayImage(someRow:end, :) = 0;

for example:

grayImage(500:end, :) = 0;

and this is my result:

I think, that it just cropped image. So it doesn´t solve my problem with binaryImage_pater.

[Image Analyst]

Well, erase more of it. Go from 450 or something.

If you still have problems, attach your current code and original image (not the one with red stuff).

[Veronika]

Ok, I attach my code and original image.

[Image Analyst]

I don't see a problem. I don't get the image you got above, with the cluster or red circles along the bottom of the table.

It looks like you could use 425 as the line to erase below. Why do you want those red circles around the boundary? They're distracting.

[Veronika]

It works. You´re right as usual, I don´t need those red circles around the boundary. Thank you so much.

[Veronika]

I still have the very last problem, which I need to resolve. I would like to segment heart from this picture, it means only this part of image:

I tried change threshold value, but I didn´t succes, also I think, that remove the smallest or largest blobs isn´t right way. Because the best segmentation, which I got is this:

I attach code and original image. Any other advice? Thank you for your helpfulness.

[Image Analyst]

It's not clear to me where the heart is (it's gray levels continuously connect with the rest of the body), so you might have to draw it manually with imfreehand. I attach a demo.

[Veronika]

Thank you very much!

[Veronika]

Is there any other option for this segmentation of heart than function imfreehand?

[Image Analyst]

Possibly, though I don't have it. You'd have to come up with some kind of custom algorithm. Check the literature, either Vision Bib or PubMed.

[Veronika]

Dear Image Analyst,

I have another problem...

I have code for segmentation parts of image. And I would like to binarization left lung, but I don´t know, why this code isn´t working. This is my code:

      fontSize = 15;
      baseFileName = '110.jpg'; 
      % baseFileName = 'thorax-mdl.jpg'; 
      folder = pwd
      fullFileName = fullfile(folder, baseFileName);
      % Načtení obrazu.
      grayImage = imread(fullFileName);
      grayImage_pater = imread (fullFileName);
      % Dimenze obrazu. 
      % ***
      [rows, columns, numberOfColorChannels] = size(grayImage);
      if numberOfColorChannels > 1
        % Máme barevný obraz, musíme ho převést na černobílý = vybereme zelený kanál
        grayImage = grayImage(:, :, 2); % zelený kanál
      end
      eq_grayImage = histeq(grayImage);%ekvalizace pomocí histogramu obrazu
      [rows, columns, numberOfColorChannels] = size(grayImage_pater);
      if numberOfColorChannels > 1
        % Máme barevný obraz, musíme ho převést na černobílý = vybereme zelený kanál
        grayImage_pater = grayImage_pater(:, :, 2); % zelený kanál
           %grayImage_pater(425:end, :) = 0;
      end
      eq_grayImage_pater = histeq(grayImage_pater);%ekvalizace pomocí histogramu obrazu
      % Zobrazení obrazu.
      figure(1)
      imshow(grayImage, []);
      axis on;
      caption = sprintf('Originální černobílý obraz, %s', baseFileName);
      title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
      drawnow;
      hp = impixelinfo();
      % Nastavení obrazu.
      % Zvětšení obrazu na celou obrazovku.
      set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
      % Odstranění panelu nástrojů a rozbalovacího menu.
      set(gcf, 'Toolbar', 'none', 'Menu', 'none');
      % Přidání titulku okna.
      set(gcf, 'Name', 'Segmentace', 'NumberTitle', 'Off') 
      thresholdValue = 120;
      binaryImage = grayImage < thresholdValue;
      % Odstranění okolí.
      binaryImage = imclearborder(binaryImage);
      % Vyplnění otvorů.
      binaryImage = imfill(binaryImage, 'holes');
      labeledImage = bwlabel(binaryImage);
      rightLung = ismember(labeledImage, 1);
      leftLung = ismember(labeledImage, 2);
      %Roztažení binárního obrazu pro přesnější segmentaci
      se = strel('line',5,100);
      rightLung= imdilate(rightLung,se);
      %Roztažení binárního obrazu pro přesnější segmentaci
      se = strel('line',5,100);
      leftLung= imdilate(leftLung,se);
      subplot (2, 3, 3)
      imshow(binaryImage, []);
      axis on;
      caption = sprintf('Binární obraz plic');
      title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
      labeledImage = bwlabel(binaryImage);
      rightLung = ismember(labeledImage, 1);
      leftLung = ismember(labeledImage, 2);
      subplot (2, 3, 4)
      imshow(leftLung, []);
      axis on;
      caption = sprintf('Binární obraz levé plíce');
      title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
      subplot (2, 3, 5)
      imshow(rightLung, []);
      axis on;
      caption = sprintf('Binární obraz pravé plíce');
      title(caption, 'FontSize', fontSize, 'Interpreter', 'None');

I attach original image. Can you advise me? Thank you for your answer.