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