{"id":5781,"date":"2022-10-24T14:57:46","date_gmt":"2022-10-24T18:57:46","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/?p=5781"},"modified":"2022-10-24T14:57:46","modified_gmt":"2022-10-24T18:57:46","slug":"double-thresholding","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2022\/10\/24\/double-thresholding\/","title":{"rendered":"Double Thresholding"},"content":{"rendered":"
Back in June, my answer to a <\/span>question on MATLAB Answers<\/span><\/a> used a technique called <\/span>double thresholding<\/span>, also known as <\/span>hysteresis thresholding<\/span>. Although the question raised was a 1-D problem, I wanted to share with you the general idea in two dimensions, making use of the function <\/span>imreconstruct<\/span>.<\/span><\/div>
Double thresholding can be useful when a single threshold produces unsatisfactory results because it either includes too many pixels in the foreground (with a low threshold), or it includes too few (with a high threshold). The idea is to use two thresholds, one that is more selective and one that is less selective, and then let the more selective result \"grow\" by adding touching pixels from the less selective result.<\/span><\/div>
Here's a simple, contrived example to illustrate what I mean.<\/span><\/div>
A = mat2gray(peaks(200));<\/span><\/span><\/div><\/div>
imshow(A)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
Suppose I want to compute a mask of the peaks that exceed 0.9:<\/span><\/div>
T_high = 0.9;<\/span><\/span><\/div><\/div>
mask_high = A > T_high;<\/span><\/span><\/div><\/div>
imshowpair(A,mask_high)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
That's good, but I'd like for the mask to include more of the peak, so I try a lower threshold.<\/span><\/div>
T_low = 0.65;<\/span><\/span><\/div><\/div>
mask_low = A > T_low;<\/span><\/span><\/div><\/div>
imshowpair(A,mask_low)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
Now the mask includes more of the peak of interest, but it also includes some of the smaller peaks that I don't want.<\/span><\/div>
The double thresholding idea is to use the output from the <\/span>higher<\/span> threshold to select all the touching regions in the output from the <\/span>lower<\/span> threshold. Morphological reconstruction can do exactly that. The Image Processing Toolbox function to use is <\/span>imreconstruct<\/span>.<\/span><\/div>
mask_double = imreconstruct(mask_high,mask_low);<\/span><\/span><\/div><\/div>
imshowpair(A,mask_double)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
Now the mask includes more of the high peak, without also including portions of the lower peaks.<\/span><\/div>
Let's try a slightly less contrived example that uses double thresholding on the gradient magnitude of a grayscale image.<\/span><\/div>
B = imread(<\/span>\"coins.png\"<\/span>);<\/span><\/span><\/div><\/div>
imshow(B)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
Bg = mat2gray(imgradient(B));<\/span><\/span><\/div><\/div>
imshow(Bg)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
Can we just get the enclosing edges for each coin using a simple threshold?<\/span><\/div>
T_high = 0.6;<\/span><\/span><\/div><\/div>
Bg_mask_high = Bg > T_high;<\/span><\/span><\/div><\/div>
imshow(Bg_mask_high)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
Try a lower threshold:<\/span><\/div>
T_low = 0.3;<\/span><\/span><\/div><\/div>
Bg_mask_low = Bg > T_low;<\/span><\/span><\/div><\/div>
imshow(Bg_mask_low)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
Now use a double threshold:<\/span><\/div>
Bg_mask_double = imreconstruct(Bg_mask_high,Bg_mask_low);<\/span><\/span><\/div><\/div>
imshow(Bg_mask_double)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
The <\/span>MATLAB Answers post<\/span><\/a> that I responded to concerned a 1-D example.<\/span><\/div>
y = [1 2 0.6 2 3 4 5 4 3 2 0.4 0.3 0.4 0.4 0.6 0.4 0.5 0.4 0.6 0.5 <\/span>...<\/span><\/span><\/div><\/div>
0.3 0.5 0.6 0.2 1 2 3 0.6 1 2 3 2 3 2 5 4];<\/span><\/span><\/div><\/div>
plot(y)<\/span><\/span><\/div><\/div>
grid <\/span>on<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
The problem was to extract low portions below 0.5, but using some kind of loose threshold that allowed portions of the clump to be a bit higher than 0.6. The double threshold does the trick nicely.<\/span><\/div>
T_low = 0.5;<\/span><\/span><\/div><\/div>
mask_low = y < 0.5;<\/span><\/span><\/div><\/div>
T_high = 0.7;<\/span><\/span><\/div><\/div>
mask_high = y < T_high;<\/span><\/span><\/div><\/div>
plot(y)<\/span><\/span><\/div><\/div>
hold <\/span>on<\/span><\/span><\/div><\/div>
plot(T_low*mask_low, LineWidth = 3, Color = <\/span>\"red\"<\/span>)<\/span><\/span><\/div><\/div>
plot(T_high*mask_high, LineWidth = 3, Color = <\/span>\"green\"<\/span>)<\/span><\/span><\/div><\/div>
hold <\/span>off<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
Compute the double threshold using <\/span>imreconstruct<\/span>:<\/span><\/div>
mask_double = imreconstruct(mask_low,mask_high);<\/span><\/span><\/div><\/div>
plot(y)<\/span><\/span><\/div><\/div>
hold <\/span>on<\/span><\/span><\/div><\/div>
plot(T_low*mask_double, LineWidth = 3, Color = <\/span>\"red\"<\/span>)<\/span><\/span><\/div><\/div>
hold <\/span>off<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>\r\n<\/div>