function test_ms

k=100; %Threshold

im=imread('flower.tif');
figure
imshow(im);
[M,N]=size(im);
edges=[];

for i=1:M-1
    for j=1:N-1
        if im(i,j)<k && im(i+1,j)<k && im(i+1,j+1)<k && im(i,j+1)<k
            edges=edges;
        elseif im(i,j)<k && im(i+1,j)>k && im(i+1,j+1)<k && im(i,j+1)<k
            edges=[edges;i+1-(k-im(i+1,j))/(im(i,j)-im(i+1,j)),j,i+1,j+(k-im(i+1,j))/(im(i+1,j+1)-im(i+1,j))];
        elseif im(i,j)<k && im(i+1,j)<k && im(i+1,j+1)>k && im(i,j+1)<k
            edges=[edges;i+1-(k-im(i+1,j+1))/(im(i,j+1)-im(i+1,j+1)),j,i+1,j+1-(k-im(i+1,j+1))/(im(i+1,j)-im(i+1,j+1))];
        elseif im(i,j)<k && im(i+1,j)>k && im(i+1,j+1)>k && im(i,j+1)<k
            edges=[edges;i+1-(k-im(i+1,j))/(im(i,j)-im(i+1,j)),j,i+1-(k-im(i+1,j+1))/(im(i,j+1)-im(i+1,j+1)),j+1];
        elseif im(i,j)<k && im(i+1,j)<k && im(i+1,j+1)<k && im(i,j+1)>k
            edges=[edges;i+(k-im(i,j+1))/(im(i+1,j+1)-im(i,j+1)),j+1,i,j+1-(k-im(i,j+1))/(im(i,j)-im(i,j+1))];
        elseif im(i,j)<k && im(i+1,j)>k && im(i+1,j+1)<k && im(i,j+1)>k
            edges=[edges;i+1-(k-im(i+1,j))/(im(i,j)-im(i+1,j)),j,i+1,j+(k-im(i+1,j))/(im(i+1,j+1)-im(i+1,j))];
            edges=[edges;i+(k-im(i,j+1))/(im(i+1,j+1)-im(i,j+1)),j+1,i,j+1-(k-im(i,j+1))/(im(i,j)-im(i,j+1))];
        elseif im(i,j)<k && im(i+1,j)<k && im(i+1,j+1)>k && im(i,j+1)>k
            edges=[edges;i,j+1-(k-im(i,j+1))/(im(i,j)-im(i,j+1)),i+1,j+1-(k-im(i+1,j+1))/(im(i+1,j)-im(i+1,j+1))];  
        elseif im(i,j)<k && im(i+1,j)>k && im(i+1,j+1)>k && im(i,j+1)>k
            edges=[edges;i+(k-im(i,j))/(im(i+1,j)-im(i,j)),j,i,j+(k-im(i,j))/(im(i,j+1)-im(i,j))];
      
        elseif im(i,j)>k && im(i+1,j)<k && im(i+1,j+1)<k && im(i,j+1)<k
            edges=[edges;i+(k-im(i,j))/(im(i+1,j)-im(i,j)),j,i,j+(k-im(i,j))/(im(i,j+1)-im(i,j))];  
        elseif im(i,j)>k && im(i+1,j)>k && im(i+1,j+1)<k && im(i,j+1)<k
            edges=[edges;i,j+1-(k-im(i,j+1))/(im(i,j)-im(i,j+1)),i+1,j+1-(k-im(i+1,j+1))/(im(i+1,j)-im(i+1,j+1))];
        elseif im(i,j)>k && im(i+1,j)<k && im(i+1,j+1)>k && im(i,j+1)<k
            edges=[edges;i+1-(k-im(i+1,j))/(im(i,j)-im(i+1,j)),j,i+1,j+(k-im(i+1,j))/(im(i+1,j+1)-im(i+1,j))];
            edges=[edges;i+(k-im(i,j+1))/(im(i+1,j+1)-im(i,j+1)),j+1,i,j+1-(k-im(i,j+1))/(im(i,j)-im(i,j+1))];
        elseif im(i,j)>k && im(i+1,j)>k && im(i+1,j+1)>k && im(i,j+1)<k
            edges=[edges;i+(k-im(i,j+1))/(im(i+1,j+1)-im(i,j+1)),j+1,i,j+1-(k-im(i,j+1))/(im(i,j)-im(i,j+1))];
        elseif im(i,j)>k && im(i+1,j)<k && im(i+1,j+1)<k && im(i,j+1)>k
            edges=[edges;i+1-(k-im(i+1,j))/(im(i,j)-im(i+1,j)),j,i+1-(k-im(i+1,j+1))/(im(i,j+1)-im(i+1,j+1)),j+1];
        elseif im(i,j)>k && im(i+1,j)>k && im(i+1,j+1)<k && im(i,j+1)>k
            edges=[edges;i+1-(k-im(i+1,j+1))/(im(i,j+1)-im(i+1,j+1)),j,i+1,j+1-(k-im(i+1,j+1))/(im(i+1,j)-im(i+1,j+1))];
        elseif im(i,j)>k && im(i+1,j)<k && im(i+1,j+1)>k && im(i,j+1)>k
            edges=[edges;i+1-(k-im(i+1,j))/(im(i,j)-im(i+1,j)),j,i+1,j+(k-im(i+1,j))/(im(i+1,j+1)-im(i+1,j))];
        elseif im(i,j)>k && im(i+1,j)>k && im(i+1,j+1)>k && im(i,j+1)>k
            edges=edges;
        end
    end
end

figure
hold on
for i=1:size(edges,1)
    plot([edges(i,1),edges(i,3)],[edges(i,2),edges(i,4)]);
end