求图割算法 graph cut 的matlab代码

function [Ncut] = graphcuts(I,pad,MAXVAL)

% function [Ncut] = graphcuts(I)

% Input: I image

% pad: spatial connectivityeg. 3

% MAXVAL: maximum image value

% Output: Ncut: Binary map 0 or 1 corresponding to image segmentation

I = double(I)[H,W] = size(I)

% Find weights between nodes I1 and I2, w = exp(a*abs(I1-I2))

% Set a to have a weight of 0.01 for diff = MAXVAL

a = log(0.01)/MAXVALx = [0:MAXVAL/100:MAXVAL]'y = exp(a*x)

figureplot(x,y)xlabel('intensity diff')ylabel('weights')title('weights')

ws = 2*pad + 1

if(ws <= 3)

ws = 3

end

%Build the weight matrix

disp('Building Weight Matrix')close alltic

WM = zeros(H*W,H*W)countWM = 0

for kk = 1:W

for jj = 1:H

mask = logical(zeros(H,W))

cs = kk-padce = kk+padrs = jj-padre = jj+pad

if(cs<1)

cs = 1

end

if(ce>W)

ce = W

end

if(rs<1)

rs = 1

end

if(re>H)

re = H

end

mask(rs:re,cs:ce) = 1

idx = find(mask==1)

p = abs(I(idx) - I(jj,kk))p = exp(a*p)

countWM = countWM + 1WM(countWM,idx) = p(:)'

end

end

ttime = tocdisp(sprintf('Time for generating weight matrix = %f',ttime))clear countWM

% Weight between a node and iteself is 0

for jj = 1:H*W

WM(jj,jj) = 0

end

WM = sparse(WM)

% Shi and Malik Algorithm: second smallest eigen vector

disp('Finding Eigen Vector')

d = sum(WM,2)D = diag(d)tic

B = (D-WM)B = (B+B')/2OPTS.disp = 0

[v,d,flag] = eigs(B,D,2,'SA',OPTS)ttime = toc

disp(sprintf('Time for finding eigen vector = %f',ttime))clear OPTS

y = v(:,2)

Ncut = reshape(y,H,W)

Ncut = Ncut >0

image_1=imread('E:\ebook\lena.bmp')%读入图片

image_1=rgb2gray(image_1)%灰度化

[m,n]=size(image_1)%计算图片的像素点个数,行列,n是列数，Gray

num=zeros(1,256)%存放各灰度级出现的次数

p=zeros(1,256)%存放慧铅各灰度级的比率

image_1=double(image_1)%双精度化

for i=1:m

for j=1:n

num(image_1(i,j)+1)=num(image_1(i,j)+1)+1%统计各灰度级的像素点个数

end

end

for i=1:256

p(i)=num(i)/(m*n)%计中碧销算各灰度级出现的比率

end

for i=2:256

if p(i)~=0

st=i+1%实现寻找出现比卖游率不为0的最小灰度值

break

end

end

for i=256:-1:1

if p(i)~=0

nd=i-1%实现找出出现比率不为0的最大灰度值

break

end

end

%以下程序实现利用最小方差和法找出门阈值

w=infth=0

for t=st:nd%最小非零比率灰度值到最大非零比率灰度值

qt1=0qt2=0%前景后景像素点比率

u1=0 u2=0%前景后景均值

v1=0 v2=0%

for i=1:t

qt1=qt1+p(i)

end

for i=1:t

u1=u1+i*p(i)/qt1

end

for i=1:t

v1=v1+((i-u1)^2)*p(i)/qt1

end

for i=t+1:256

qt2=qt2+p(i)

end

for i=t+1:256

u2=u2+i*p(i)/qt2

end

for i=t+1:256

v2=v2+((i-u2)^2)*p(i)/qt2

end

if qt1*v1+qt2*v2<w

th=tw=qt1*v1+qt2*v2

end

end

for i=1:m

for j=1:n

if (image_1(i,j)+1>th)

image_2(i,j)=255

else

image_2(i,j)=0

end

end

end

image_2=uint8(image_2)%读入读出变换

figure,imshow(image_2)%显示二值化后的图片

---