How to arrange pixels for appropriate threshold

I have two mammography image .

I want to segment tumor area and for this using threshold method.

But image intensities are different for every image so it affect my result badly.

I tried to calculate mean intensity of image with using "meanIntensity = mean(grayImage(:)) " and if condition statement according this value for threshold, but it doesn't work again in a right way.

I mean mdb005 image I am using threshold > 155 to get good result but for mdb013 I have to choose >210 or other values. (İf I choose >155 it shows huge for mdb013)

But if I choose >210 this time mdb005 tumor dissappear.

What to do for these two image to segment exact tumor automatically.

I used @Image Analystcodes in the past questions and combine with other codes and created this codes.

images attached

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 = 25;

%===============================================================================
% Get the name of the image the user wants to use.
% baseFileName = 'pat00002_1-1-1-1-2-10-1.mp4.cover.png';
baseFileName = 'mdb005.png';
% Get the full filename, with path prepended.
folder = ['C:\Users\Ali\Desktop\Miaspng']; % 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();

% 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', 'NumberTitle', 'Off')
drawnow;


meanIntensity = mean(grayImage(:))

if meanIntensity<35
    mask = grayImage > 190;
elseif meanIntensity<40
    mask = grayImage > 185;
elseif meanIntensity<45
    mask = grayImage > 182;
elseif meanIntensity<50
    mask = grayImage > 175;
elseif meanIntensity<55
    mask = grayImage > 170;
elseif meanIntensity<60
    mask = grayImage > 165;
elseif meanIntensity<65
    mask = grayImage > 155;
elseif meanIntensity<70
    mask = grayImage > 150; 
end


% Display the mask image.
subplot(2, 3, 2);
imshow(mask);
axis on;
axis image; % Make sure image is not artificially stretched because of screen's aspect ratio.
title('Binary Image Mask', 'fontSize', fontSize);
drawnow;

% Find the areas
labeledImage = bwlabel(mask);
props = regionprops(labeledImage, 'Area', 'Centroid');
allAreas = sort([props.Area], 'descend');
allCentroids = [props.Centroid];
centroidsX = allCentroids(1:2:end);
centroidsY = allCentroids(2:2:end);

% Make a margin of 10% of the image size.
marginx = 0.34 * columns;
marginy = 0.15 * rows;
keepers = (centroidsX > marginx & centroidsX < (columns - marginx)) & ...
          (centroidsY > marginy & centroidsY < (rows - marginy));
indexes = find(keepers);
% Get a mask with only the keepers in it
newMask = ismember(labeledImage, indexes);
% Display the mask image.
subplot(2, 3, 3);
imshow(newMask);
axis on;
axis image; % Make sure image is not artificially stretched because of screen's aspect ratio.
title('Tags Remove', 'fontSize', fontSize);
drawnow;

% Mask the gray scale image
maskedGrayImage = grayImage; % Initialize.
maskedGrayImage(~newMask) = 0;
% Display the masked grayscale image.
subplot(2, 3, 4);
imshow(maskedGrayImage);
axis on;
axis image; % Make sure image is not artificially stretched because of screen's aspect ratio.
title('Masked Gray Scale Image', 'fontSize', fontSize);
drawnow;


IM_cb = imclearborder(maskedGrayImage);
 BW2 = bwareaopen(IM_cb, 150);
BW_filled = imfill(BW2, 'holes');
 subplot(2, 3, 5);
%  figure,
 imshow(BW_filled);
 axis on;
axis image; % Make sure image is not artificially stretched because of screen's aspect ratio.
title('Small Lines Eliminated', 'fontSize', fontSize);
drawnow;
 
 labeledImage = bwlabel(BW_filled);
measurements = regionprops(labeledImage, 'Area');
% Get all the areas
allAreas = [measurements.Area]; 
[biggestArea, indexOfBiggest] = sort(allAreas, 'descend');
% Extract biggest
biggestBlob = ismember(labeledImage, indexOfBiggest(1));
% Convert to binary
biggestBlob = biggestBlob > 0;
 subplot(2, 3, 6);
% figure,
imshow(biggestBlob, []);
axis on;
axis image; % Make sure image is not artificially stretched because of screen's aspect ratio.
title('Final', 'fontSize', fontSize);
drawnow;

Answer 1

[Image Analyst]

Something that simplistic will rarely work. You need to use more sophisticated algorithm like given here:

http://www.visionbib.com/bibliography/contentsmedical.html#Medical%20Applications,%20CAT,%20MRI,%20Ultrasound,%20Heart%20Models,%20Brain%20Models

Look at sections:

20.7 X-Ray Images, Radiography

• 20.7.1 Breast Cancer, Mammograms, Analysis, Mammography
• 20.7.1.1 Breast Cancer Cell Analysis, Pathology, Nuclei Detection
• 20.7.1.2 Mammography, Microcalcifications, Detection, Analysis
• 20.7.1.3 Mammography, Thermal, Infrared Analysis
• 20.7.1.4 Mammography, Texture Based Techniques, Wavelets
• 20.7.1.5 Mammograms, Image Enhancement, Noise Suppression
• 20.7.1.6 Mammograms, Ultrasound
• 20.7.1.7 Mammograms, Density Issues
• 20.7.1.8 Mammograms, MRI, Magnetic Resonance Imaging
• 20.7.1.9 Mammograms, Three Dimensional Analysis, Registration

Find a paper you like and code it up or ask the author for code.

[Ali Zulfikaroglu]

Thank you for sharing but I read many paper . But actually some details are not so clear. So I couldn't apply codes exactly.

Codes are working in a good way with many images in mias databes however some images are not working well like mdb014. Is it possible to add some codes and arrange threshold according to background intensity or pixel values?

This is the idea that comes to my mind for now to make it better.

[Image Analyst]

If you simply want to threshold an image with a global threshold and analyze the bright blobs for area, etc. then see my Image Segmentation tutorial in my File Exchange:

Image Analyst's File Exchange

[Ali Zulfikaroglu]

I found a way but I couldn't apply. Below explanation:

Optimal thresholding can be implemented iteratively by the isodata (iterative selforganizing

data analysis technique algorithm) method. The steps are as follows:

(i) threshold the image using the mean of the two peaks or the mean pixel value, T0;

(ii) calculate the mean value of the pixels below this threshold, μ1, and the mean of the

pixels above this threshold, μ2;

(iii) threshold the image at a new threshold, Ti = (μ1 + μ2)/2;

(iv) repeat steps (ii)–(iii) until Ti – Ti−1 ≤ Δ (where the change, Δ, can be defined in

several different ways, either by measuring the relative change in threshold value or

by the percentage of pixels that change sides (foreground to background or vice

versa) between iterations).

@Image Analyst How can we apply this method ?

[Image Analyst]

@Ali Zulfikaroglu, not sure it's right, but look this over and try it.

% Demo by Image Analyst
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 = 20;

%--------------------------------------------------------------------------------------------------------
%    READ IN IMAGE
folder = pwd;
baseFileName = 'cameraman.tif';
grayImage = imread(baseFileName);

% 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.
	% Extract the red channel (so the magenta lines will be white).
	grayImage = grayImage(:, :, 1);
end

%--------------------------------------------------------------------------------------------------------
% Display the image.
subplot(2, 2, 1);
imshow(grayImage, []);
axis('on', 'image');
title('Original Gray Scale Image', 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;
hFig = gcf;
hFig.WindowState = 'maximized'; % May not work in earlier versions of MATLAB.
drawnow;

% Take histogram
subplot(2, 2, 2);
edges = 0 : 256;
histObject = histogram(grayImage, edges)
grid on;
title('Delta From Prior Threshold', 'FontSize', fontSize);
xlabel('Gray Level', 'FontSize', fontSize);
ylabel('Pixel Count', 'FontSize', fontSize);

% Optimal thresholding can be implemented iteratively by the isodata (iterative selforganizing
% data analysis technique algorithm) method. The steps are as follows:
% (i) threshold the image using the mean of the two peaks or the mean pixel value, T0;
% (ii) calculate the mean value of the pixels below this threshold, μ1, and the mean of the
% pixels above this threshold, μ2;
% (iii) threshold the image at a new threshold, Ti = (μ1 + μ2)/2;
% (iv) repeat steps (ii)–(iii) until Ti – Ti−1 ≤ Δ (where the change, Δ, can be defined in
% several different ways, either by measuring the relative change in threshold value or
% by the percentage of pixels that change sides (foreground to background or vice
% versa) between iterations).
counts = histObject.Values;
gls = histObject.BinEdges(1 : end-1);
deltaToQuit = 1;
thisDelta = inf;
loopCounter = 1;
maxIterations = 100;
hLine = 1;
% Set up first.
thisMean(1) = mean2(grayImage)
hLine = xline(thisMean(1), 'Color', 'r', 'LineWidth', 2);
thisDelta(1) = 256;
loopCounter = 2;
% Now iterate
while thisDelta(end) > deltaToQuit && loopCounter < maxIterations
	% Get mean less than threshold.
	threshold = thisMean(loopCounter - 1);
	index = find(gls <= threshold, 1, 'last');
	thisMean(loopCounter) = sum(counts(1:index) .* gls(1:index)) / sum(counts(1:index))
	subplot(2, 2, 2);
	% 	delete(hLine); % get rid of old red line.
	hLine = xline(thisMean(loopCounter), 'Color', 'r', 'LineWidth', 2);
	% Compute delta from prior mean
	if loopCounter >= 2
		thisDelta(loopCounter) = abs(thisMean(loopCounter) - thisMean(loopCounter - 1));
	end
	loopCounter = loopCounter + 1;
	
	% Plot progress.
	subplot(2, 2, 3);
	plot(thisDelta, 'b.-', 'LineWidth', 2, 'MarkerSize', 15);
	grid on;
	title('Delta From Prior Threshold', 'FontSize', fontSize);
	xlabel('Iteration', 'FontSize', fontSize);
	ylabel('Gray Levels', 'FontSize', fontSize);
	
	subplot(2, 2, 4);
	plot(thisMean, 'b.-', 'LineWidth', 2, 'MarkerSize', 15);
	grid on;
	title('Threshold', 'FontSize', fontSize);
	xlabel('Iteration', 'FontSize', fontSize);
	ylabel('Gray Levels', 'FontSize', fontSize);
	drawnow;
end
meanGL = thisMean(end);
message = sprintf('The final threshold is %.2f gray levels.\nIt took %d iterations.', meanGL, loopCounter-1);
uiwait(helpdlg(message));

[Ali Zulfikaroglu]

grayImage = imread('C:\Users\Ali\Desktop\all-mias\mdb005.pgm');

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

if(size(grayImage,3)>1)

grayImage=rgb2gray(grayImage);

end

Bw=grayImage>148;

figure, imshow(Bw);

I segmented tumor with this number for mdb005 image and result is not bad. Threshold value should be something close to this number. But according to your plots what is the best threshold value ? I applied your code for mdb005 image and result of graphs at the attachment . But it doesn't give best threshold value. So I couldn't find which threshold value should I use ?

[Ali Zulfikaroglu]

In this paper, they managed really well.

This is the best tumor segmentation method I have ever seen and result is very successful.

People can read and cite this and can learn .

But I couldn't manage to write codes @Image Analyst

[Image Analyst]

Well I just did the steps you gave me. I didn't think it was very good. Why don't you ask the authors for the code?

[Ali Zulfikaroglu]

I got code from author but it was written by C++ and Actually it is complicated very and I have less knowledge about C++. I will be very glad if you read paper and manage to apply it . Because I searched and read many paper and this is one of the best way . If you manage, believe me you will help many people and even doctors to get good accuracy and give more right decision for patienets.

[Image Analyst]

Sorry I can't put it on my agenda. I simply don't have enough time to code up all the interesting papers I run across. Plus I don't want to rob you of all the glory you said would come your way once you write it.

[Ali Zulfikaroglu]

Thank u @Image Analyst for everything . You helped a lot already . I will figure out some way