How to find pixel counts from an Histogram for 2D uint16 images

I have double MR image (2D uint16), I have plotted hitogram of my original image and I want to find the Otsu threshold level and display it on the histogram. There is my code : How could I correct it please. Thank you.

subplot(2, 4, 2);
h=histogram(grayImage);
title('Histogram of Original Image');
thresholdLevel = graythresh(grayImage);
y = ylim();
line(thresholdLevel, y,'Color', 'r');
grid on;

Answer 1

[Guillaume]

You made two mistakes:

• because you pass a 1x1 x value and a 1x2 y value to line, it actually plots two points, not a line. One point at (thresholdLevel, y(1)) the other at (thresholdLevel, (y(2)).
• graythresh return a normalised threshold in the range [0 1], so to plot your line on your histogram you need to scale the threshold back to the intensity range of your image

So the fix:

line(repelem(threshold * intmax('uint16'), 2), y, 'Color', 'r')

Answer 2

[Image Analyst]

Try this:

h=histogram(grayImage);
title('Histogram of Original Image');
thresholdLevel = graythresh(grayImage)
% Convert to an actual gray level.
thresholdLevel = thresholdLevel * intmax(class(grayImage))
line([thresholdLevel, thresholdLevel], ylim, 'Color', 'r', 'LineWidth', 2); % Plot vertical line.
grid on;
% Threshold the image.
binaryImage = grayImage > thresholdLevel;

[MMSAAH]

Thank You very much. This was really useful.

I have tried this answer for 2D int16 MR images...It gaves erroneous results. So please how to adjust it to this type of images? Thank You

[Image Analyst]

What do you mean by erroneous? The code works perfectly well for 16 bit images. Here's proof:

% Create 2-D 16 bit (uint16) image.
grayImage = uint16(randi(65535, 480, 640));
h=histogram(grayImage);
title('Histogram of Original Image');
thresholdLevel = graythresh(grayImage)
% Convert to an actual gray level.
thresholdLevel = thresholdLevel * intmax(class(grayImage))
line([thresholdLevel, thresholdLevel], ylim, 'Color', 'r', 'LineWidth', 2); % Plot vertical line.
grid on;
% Threshold the image.
binaryImage = grayImage > thresholdLevel;

[MMSAAH]

For signed 2D image (int16), I got thresholdlevel outside the histogram.

[Image Analyst]

What is the attachment? Inside the zip file is a binary file with no extension so Windows doesn't know what program to open it with.

[MMSAAH]

It is a dicom image. I usually open this type of images by ITKSNAP or Amide.

[Image Analyst]

I don't know why graythresh picks that weird threshold for your image. I'd guess it's because it bins the data into only 256 bins when deciding the threshold and most of your data would only be in the lowest 2 bins. Try using imbinarize() instead:

% Create 2-D 16 bit (uint16) image.
grayImage = dicomread('IM15.dcm');
subplot(2, 2, 1);
imshow(grayImage, []);
impixelinfo();
subplot(2, 2, 2);
h=histogram(grayImage);
grid on;
title('Histogram of Original Image');
% Convert to an actual gray level.
thresholdLevel = thresholdLevel * intmax(class(grayImage))
line([thresholdLevel, thresholdLevel], ylim, 'Color', 'r', 'LineWidth', 2); % Plot vertical line.
grid on;
% Threshold the image.
binaryImage = imbinarize(grayImage);
subplot(2, 2, 3);
imshow(binaryImage, []);

[MMSAAH]

Sorry! but why to use the thresholdLevel while I 'm not going to use it to have the binary image ?

[Guillaume]

Note that there shouldn't be any difference difference between using the threshold returned by graythresh and calling imbinarize with no options, since both end up calling otsuthresh, using 256 bins, after having converted the image to uint8.

The problem with using

thresholdLevel = thresholdLevel * intmax(class(grayImage))

is that it is only valid for unsigned classes (which was the case in the original question, but no longer is for the dicom image given). For signed classes, the formula should be:

classrange = [intmin(class(grayImage)), intmax(class(grayImage))];
thresholdLevel = thresholdLevel*diff(classrange) + classrange(1);

That is, for int16 images the threshold is between -32768 and 32767, not 0 and 32767.

Note 2: you should remove the metadata from the images you post here, particularly patient information, unless it's made up.