-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathhysthresh.m
100 lines (83 loc) · 3.36 KB
/
hysthresh.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
% HYSTHRESH - Hysteresis thresholding
%
% Usage: bw = hysthresh(im, T1, T2)
%
% Arguments:
% im - image to be thresholded (assumed to be non-negative)
% T1 - upper threshold value
% T2 - lower threshold value
%
% Returns:
% bw - the thresholded image (containing values 0 or 1)
%
% Function performs hysteresis thresholding of an image.
% All pixels with values above threshold T1 are marked as edges
% All pixels that are adjacent to points that have been marked as edges
% and with values above threshold T2 are also marked as edges. Eight
% connectivity is used.
%
% It is assumed that the input image is non-negative
%
% Author: Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% pk @ csse uwa edu au http://www.csse.uwa.edu.au/~pk
%
% December 1996 - Original version
% March 2001 - Speed improvements made (~4x)
%
% A stack (implemented as an array) is used to keep track of all the
% indices of pixels that need to be checked.
% Note: For speed the number of conditional tests have been minimised
% This results in the top and bottom edges of the image being considered to
% be connected. This may cause some stray edges to be propagated further than
% they should be from the top or bottom.
%
function bw = hysthresh(im, T1, T2)
if (T2 > T1 | T2 < 0 | T1 < 0) % Check thesholds are sensible
error('T1 must be >= T2 and both must be >= 0 ');
end
[rows, cols] = size(im); % Precompute some values for speed and convenience.
rc = rows*cols;
rcmr = rc - rows;
rp1 = rows+1;
bw = im(:); % Make image into a column vector
pix = find(bw > T1); % Find indices of all pixels with value > T1
npix = size(pix,1); % Find the number of pixels with value > T1
stack = zeros(rows*cols,1); % Create a stack array (that should never
% overflow!)
stack(1:npix) = pix; % Put all the edge points on the stack
stp = npix; % set stack pointer
for k = 1:npix
bw(pix(k)) = -1; % mark points as edges
end
% Precompute an array, O, of index offset values that correspond to the eight
% surrounding pixels of any point. Note that the image was transformed into
% a column vector, so if we reshape the image back to a square the indices
% surrounding a pixel with index, n, will be:
% n-rows-1 n-1 n+rows-1
%
% n-rows n n+rows
%
% n-rows+1 n+1 n+rows+1
O = [-1, 1, -rows-1, -rows, -rows+1, rows-1, rows, rows+1];
while stp ~= 0 % While the stack is not empty
v = stack(stp); % Pop next index off the stack
stp = stp - 1;
if v > rp1 & v < rcmr % Prevent us from generating illegal indices
% Now look at surrounding pixels to see if they
% should be pushed onto the stack to be
% processed as well.
index = O+v; % Calculate indices of points around this pixel.
for l = 1:8
ind = index(l);
if bw(ind) > T2 % if value > T2,
stp = stp+1; % push index onto the stack.
stack(stp) = ind;
bw(ind) = -1; % mark this as an edge point
end
end
end
end
bw = (bw == -1); % Finally zero out anything that was not an edge
bw = reshape(bw,rows,cols); % and reshape the image