Asked by gugu
on 1 Jan 2018

I'm using triangle method from https://www.mathworks.com/matlabcentral/fileexchange/28047-gray-image-thresholding-using-the-triangle-method. I want to show the chosen threshold value in the histogram with a red line. The demo from the above link uses the fixed numbers to show the thresholded point. I've attached the histogram image. The threshold value returned from the triangle method is 0.5234375.

Answer by Image Analyst
on 1 Jan 2018

Accepted Answer

You can call line() once you determine the x values.

x = 0.5234375 * 255;

grid on;

hold on;

line([x, x], ylim, 'Color', 'r', 'LineWidth', 2);

hold off;

You might want to look at my code below which basically does the same algorithm with many tweaks, like showing you the plot and lines and how it determined the threshold, and lots of comments and explanations.

% Triangle Thresholding Algorithm

% pixelCounts is the image histogram

% side is either 'R' or 'Right' to find a threshold on the right side of the histogram peak,

% or 'L' or 'Left' to find a threshold on the left side of the histogram peak.

function thresholdLevel = triangle_threshold(pixelCounts, side, showPlot)

% This technique is due to Zack (Zack GW, Rogers WE, Latt SA (1977),

% "Automatic measurement of sister chromatid exchange frequency",

% J. Histochem. Cytochem. 25 (7): 741–53, )

% A line is constructed between the maximum of the histogram at

% (b) and the lowest (or highest depending on context) value (a) in the

% histogram. The distance L normal to the line and between the line and

% the histogram h[b] is computed for all values from a to b. The level

% where the distance between the histogram and the line is maximal is the

% threshold value (level). This technique is particularly effective

% when the object pixels produce skewed unimodal histogram,

% where it will find the "corner" of the histogram.

% Sample call:

% binOfThreshold = triangle_threshold(pixelCount, 'R');

try

side = upper(side(1)); % Get just first letter and capitalize it.

% Find maximum of histogram and its location along the x axis

% Find x and y of the peak of the histogram.

[yPeak, xPeak] = max(pixelCounts);

% Find location of first and last non-zero values.

% Values < (h / 10000) are considered zeros.

nonZeroBins = find(pixelCounts>yPeak / 10000);

firstNonZeroBin = nonZeroBins(1);

lastNonZeroBin = nonZeroBins(end);

% Sometimes, through contrast stretching or whatever, there are gaps in the histogram. Fill those.

gaps = pixelCounts == 0;

% Don't include flat parts on the left or right.

gaps(1:firstNonZeroBin-1) = false;

gaps(lastNonZeroBin+1 : end) = false;

% pixelCounts = regionfill(pixelCounts, gaps); % Doesn't work because regionfill doesn't work with 1-D "images"

xAll = 1 : length(pixelCounts);

xGood = 1 : length(pixelCounts);

xGood(gaps) = []; % Remove gaps

yGood = pixelCounts;

yGood(gaps) = [];

% Fill gaps by interpolating in 1-D across the gaps.

pixelCounts = interp1(xGood, yGood, xAll);

% The max can occur over multiple bins (if histogram hump has a flat top),

% so find the bin that is closest to the furthest non-zero bin on the side the user wants to threshold on.

if side(1) == 'R'

xPeak = find(pixelCounts == yPeak, 1, 'last');

xm = lastNonZeroBin; % Point on far right side, at the end of the tail on the right.

ym = pixelCounts(xm);

else

xPeak = find(pixelCounts == yPeak, 1, 'first');

xm = firstNonZeroBin; % Point on far left side, at the end of the tail on the left.

ym = pixelCounts(xm);

end

% Find the denominator for our point-to-line distance formula.

% Reference: http://mathworld.wolfram.com/Point-LineDistance2-Dimensional.html

denominator = sqrt((xm - xPeak) ^ 2 + (ym - yPeak) ^ 2);

% Get the points along the histogram

xh = 1 : length(pixelCounts);

yh = pixelCounts;

if showPlot

% Plot it

hFig = figure;

bar(xh, yh, 'BarWidth', 1);

xlim([1, length(xh)]);

grid on;

% Maximize the window via undocumented Java call.

% Reference: http://undocumentedmatlab.com/blog/minimize-maximize-figure-window

MaximizeFigureWindow;

% Plot a line from peak to end

hold on;

plot([xPeak, xm], [yPeak, ym], 'r-', 'LineWidth', 2);

drawnow;

end

% Find the numerator for our point-to-line distance formula.

numerator = abs((xm-xPeak)*(yPeak-yh) - (xPeak - xh) * (ym - yPeak));

% Compute all the distances

distances = numerator ./ denominator;

% Sometimes because of the shape of the histogram, which bulges up, there are some bins that are above

% the line from point m (the tail) to the peak. This means the distance is being measured from the line

% going upwards, above the hypoteneuse, instead of downwards. This is not what we want.

% We need to find those distances and zero them out (set to -infinity) so they won't be considered.

slope = (yPeak-ym)/(xPeak-xm);

yLine = slope * (xh - xm) + ym;

% plot(xh(1:xPeak), yLine(1:xPeak), 'm-', 'LineWidth', 2)

badDistanceIndexes = yh > yLine;

distances(badDistanceIndexes) = -inf;

% Zero out the distances on the side of the histogram that we are not considering.

if side(1) == 'R'

% Looking on the right (bright) side. Zero out bins to the left of the peak.

startBin = xPeak - 1;

if startBin < 1

startBin = 1;

end

distances(1:startBin) = 0;

% Zero out bins to the right of the right-most point of the hypoteneuse.

startBin = lastNonZeroBin + 1;

if startBin > length(distances)

startBin = length(distances);

end

distances(startBin:end) = 0;

else

% Looking on the left (dark) side. Zero out bins to the right of the peak.

startBin = xPeak + 1;

if startBin > length(distances)

startBin = length(distances);

end

distances(startBin:end) = 0;

% Zero out bins to the left of the left-most point of the hypoteneuse.

startBin = lastNonZeroBin - 1;

if startBin < 1

startBin = 1;

end

distances(startBin:end) = 0;

end

% Find the max distance. Because we set bins above the hypoteneuse to -inf, they won't be found.

[maxDistance, indexOfMaxDistance] = max(distances);

thresholdLevel = indexOfMaxDistance;

if showPlot

% Find the slope and intercept of the peak to end line.

coefficients = polyfit([xPeak, xm], [yPeak, ym], 1);

m = coefficients(1); % Slope.

b = coefficients(2); % Intercept.

x0 = indexOfMaxDistance;

y0 = yh(indexOfMaxDistance);

xOnLine = (-b + m * x0 + y0) / (2*m);

yOnLine = m * xOnLine + b;

% Plot a line from peak to end

hold on;

plot([x0, xOnLine], [y0, yOnLine], 'r-', 'LineWidth', 2);

plot([x0, x0], [0, y0], 'r-', 'LineWidth', 2);

xlabel('Gray Level', 'FontSize', 20)

ylabel('Pixel Count', 'FontSize', 20)

title('Triangle Threshold Determination', 'FontSize', 20)

str = sprintf('Threshold = %.1f Gray Levels', x0);

text(x0 + 9, y0, str, 'FontSize', 20, 'Color', 'r', 'Rotation', 0)

yl = ylim(); % Get range in y direction.

% Put up "Longest Line" text angled over the line, more or less.

text(x0 + 6, y0 + 0.15*yl(2), 'Longest Line', 'FontSize', 20, 'Color', 'r', 'Rotation', 50)

drawnow;

msgboxh('This is the threshold determination');

close(hFig);

end

catch ME

% Some error happened if you get here.

errorMessage = sprintf('Error in program %s, function %s(), at line %d.\n\nError Message:\n%s', ...

mfilename, ME.stack(1).name, ME.stack(1).line, ME.message);

WarnUser(errorMessage);

end

return; % triangle_threshold

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.