% Hysteresis thresholding: keep all points with m>t connected to 
% points with m>t2 (=2.5t)
function [m1] = Hysteresis(m, t);

[L, numc] = bwlabel(m);  
t2 = 2.5*t;

[I,J] = find((L>0) & (m>t2));  % find labeled points with mag.>t2
comps = zeros(numc,1);  % initialize comps

for i=1:size(I),        % components with at least one point mag.>t2
  comps(L(I(i),J(i))) = 1;
end

[I,J] = find(L>0);         % find non-zero points
for i=1:size(I),       
  L(I(i),J(i)) = comps(L(I(i),J(i)));  %replace connected components labels
end;
   
m1 = zeros(size(m));     % initialize m
I2 = find(L>0);          % L can be treated as a vector
m1(I2) = m(I2);