我的（但是，还需要做细化操作）
W_H = [ 1   2  1;                  % the model in the horizon direction
        0   0   0;
       -1  -2  -1];
W_V = [1 0 -1;                     % the model in the vetical direction
       2 0 -2;
       1 0 -1];
T = 0.165;                           % the threshold in the 2-value
I = imread('f:\chinalake.bmp','bmp');  % read the image
%I = imread('f:\rice.tif','tif');
[height,width] = size(I);
I1 = double(I);
I2 = zeros(height+2,width+2);      % put the image's data into a bigger array to void the edge
I2(2:height+1,2:width+1) = I1;
for i=2:height+1                   % move the window and calculate the grads
    for j=2:width+1
        sum1 = 0;                  % the cumulus 
        sum2 = 0;
        for m=-1:1
            for n=-1:1
                sum1 = sum1 + W_H(m + 2,n + 2) * I2(i + m,j + n);
            end
        end
        for m=-1:1
            for n=-1:1
                sum2 = sum2 + W_V(m + 2,n + 2) * I2(i + m,j + n);
            end
        end
        grey = abs(sum1) + abs(sum2);        
        I1(i-1,j-1) = grey;
    end
end
big = max(max(I1));                          % 归一化
small = min(min(I1));
for i=1:height
    for j=1:width
        I1(i,j) = (I1(i,j)-small)/(big - small);  % 归一化
        if(I1(i,j) > T)
            I1(i,j) = 1;                          % 二值化
        else
            I1(i,j) = 0;
        end
    end
end

imshow(I1);

关于细化可以参看形态学

非常感谢，多谢你的代码，很受启发！