matlab 非极大值抑制

function outImg = NonMaxLimit( dx  ,dy )
% dx、dy是图像的X、Y方向的梯度
dxy0 = sqrt(dx.^2 + dy.^2);

dxy = zeros( size(dxy0) + 2 );
dxy(2:end-1,2:end-1) = dxy0;

bw = graythresh(dxy0);

outImg = zeros( size(dx) );
for ii = 1:size(dx,1)
for jj = 1:size(dx,2)
the = atan2( dx(ii,jj) , dy(ii,jj));
absX = abs( dx(ii,jj) );
absY = abs( dy(ii,jj) );
if absX > absY
weight = absY / absX;
else
weight = absX / absY;
end
%             1 8 7
%             2 * 6
%             3 4 5
v1 = dxy(ii,jj);
v2 = dxy(ii+1,jj);
v3 = dxy(ii+2,jj);
v4 = dxy(ii+2,jj+1);
v5 = dxy(ii+2,jj+2);
v6 = dxy(ii+1,jj+2);
v7 = dxy(ii,jj+2);
v8 = dxy(ii,jj+1);

if abs(dxy0(ii,jj) <= bw)
outImg(ii,jj) = 0;
else
if the < -pi*3/4 || the < pi/4
g1 = v7;
g2 = v6;
g3 = v3;
g4 = v2;
elseif (the > -pi*3/4 && pi<=-pi/2) || (the>pi/4 && the <= pi/2)
g1 = v7;
g2 = v8;
g3 = v3;
g4 = v4;
elseif (the > -pi/2 && the<-pi/4) || (the > pi/2 && the <= pi*3/4)
g1 = v1;
g2 = v8;
g3 = v5;
g4 = v4;
elseif (the>-pi/4 && the <= 0) || (the > pi*3/4)
g1 = v1;
g2 = v2;
g3 = v5;
g4 = v6;
end

dTmp1 = weight*g1 + (1-weight)*g2;
dTmp2 = weight*g3 + (1-weight)*g4;
if(  dxy0(ii,jj) >= dTmp1 && dxy0(ii,jj) >= dTmp2 )
outImg(ii,jj) = 1;
end
end
end
end

end