目录
练习思路
matlab绘图
噪音
坐标编码
邻接矩阵
最大流最小割算法
对最大流最小割算法求解结果转换为图像
源代码:
运行实例:
TIPS:
最近学习了最大流和最小割算法，可以把图看成是一些点的集合，色彩差值的倒数看做是权重，可以利用最小割实现简单的graph-cut操作。现在做个简单的分割椭圆的割图练习。
先用matlab的矩形色彩绘图工具绘制出椭圆的图形，再添加随机噪音，使椭圆不太明显，接着将色块的坐标编码为一维数字，然后做出邻接矩阵，用matlab的最大流最小割函数求最大流，将分割出的数字解编码为坐标，重新赋值颜色，完成graph-cut操作。
matlab绘制矩阵色块有两种方法：
1.pcolor函数，pcolor(X,Y,C),X为横坐标，Y为纵坐标，C为色值矩阵；
pcolor函数绘图时，色值对应的矩阵块是1-2，2-3之间的，所以C可显示范围要少一行一列。
要对应起来的话，可以换成X和Y可以换成0.5开始。如下：
pcolor有方格，而且横坐标是从左向右排列的，纵坐标从下到上排列的。不过本次没有用pcolor用的是imagesc。
2.imagesc(G)。G为色值矩阵。这个和pcolor不一样，色值和坐标是一一对应的。
G（y , x）y行数,对应图像的纵坐标，x列数,对应图像的横坐标。
imagesc没有方格,纵坐标是从上到下排列的。
对于正常图片可以用rgb2grey做成灰度图，然后用randn生成高斯噪音。
对于本次练习，也可以用randn做正态分布随机数，不过本次练习采用了rand采取随机数。
色块的坐标是二维的（a , b)，但是需要简化成一个数字，因为我设置的坐标最大值为M，所以采用n = M * (x - 1) + x + y来编码。对于一个x来说，当前x下的这一组数上限是M*（x + 1），下一组的数下限是（M + 1）*x + 2,刚好可以分隔开。所以，每一个色块有一个独一无二的n作为它的新坐标，用新坐标来做邻接矩阵。
邻接矩阵可以做成有向图和无向图两种，还有正常矩阵和稀疏矩阵两种。无向图需要是对称矩阵。
本次练习是做的无向图。权重是用坐标上限M的任意次幂除去色值差值的平方，放大权重。
最大流最小割是网络流的重要模型，在一个网络流的图中，设置一个源点和汇点，网络流图中的最大流，即能够运输得到最大流量和最小割，即将图中的点分为两拨，两拨之间连线权重最小的和，问题等价。
matlab计算最大流最小割有两个算法：graphmaxflow和maxflow，两者输出的数值略有差别，输入上，graphmaxflow需要输入无向图的稀疏矩阵，maxflow需要输入有向图。本次练习用的是后者。
将最大流算法求出的最小割所分出来的区域，将有源点的区域的点转换为坐标，作图。
imagesc.m
%用imagesc绘制椭圆
function[map,paint] = imagesc_ellipse(M,a,b)
map = zeros(M);
for x = 1:M
for y = 1:M
if (((x-round(M/2))^2*b^2 + (y-round(M/2))^2*a^2) <= a^2*b^2)
map(y , x) = 100;
end
end
end
paint = imagesc(map);imagesc_encode.m
%图矩阵编码列表
function[point_list,centre_code] = imagesc_encode(M)
point_list = 1:M^2;
t = 0;
for x = 1 : M
for y = 1 : M
n = M * (x - 1) + x + y;
t = t + 1;
point_list(t) = n;
end
end
%中心编码
centre_code = find(point_list == M * round(M/2));imagesc_unencode_matrix.m
%制作邻接矩阵
function[matrix] = imagesc_unencode_matrix(M,map,point_list)
matrix = zeros(M^2);m = 50;
for i = 1:M^2
x = ceil(point_list(i)/(M+1));
y = point_list(i) - M*(x-1)-x;
%与下面的点建立联系
if(y~=M)
yy = y + 1;
n = M * (x - 1) + x + yy;
n = find(point_list == n);
store = (map(y,x)-map(yy,x))^2;
matrix(n,i)=round(M^m/store);
matrix(i,n)=round(M^m/store);
end
%与右边的点建立联系
if(x~=M)
xx = x + 1;
n = M * (xx - 1) + xx + y;
n = find(point_list == n);
store = (map(y,x)-map(y,xx))^2;
matrix(n,i)=round(M^m/store);
matrix(i,n)=round(M^m/store);
end
%给边点建立强联系
if(x == 1||x == M||y == 1||y == M)
matrix(M^2,i) = inf;
matrix(i,M^2) = inf;
end
endimagesc_Result_processing.m
%最后割图结果处理
function[result] = imagesc_Result_processing(map,cs,point_list,M)
MAP =map;
for i = 1:length(cs)
cs(i) = point_list(cs(i));
x = ceil(cs(i)/(M+1));
y = cs(i) - M*(x-1)-x;
MAP(y,x) = 100;
end
figure;
result = imagesc(MAP);
gtsf6.m（主程序）
%最大流最小割图割算法
clear;clc;close all;
%绘制椭圆啦
M = 20;A = 7;B = 4;
[map,] = imagesc_ellipse(M,A,B);
%加点点噪声
map = abs(map - rand(M)*40);
piant = imagesc(map);
%编码变成列表
[point_list,centre_code] = imagesc_encode(M);
%绘制邻接矩阵
matrix = imagesc_unencode_matrix(M,map,point_list);
%绘制有权有向图
paint = digraph(matrix);
%最大流最小割
[~,~,cs,] = maxflow(paint, centre_code, M^2);
%最后结果处理
result = imagesc_Result_processing(map,cs,point_list,M);原图像
graph-cut图像
1.在gtsf6.m的第四行是椭圆的形状大小参数；
2.在gtsf6.m的第七行rand(M)后面的常数是调整噪声的参数，40已经很高了，勉强可以分辨出来椭圆；
3.imagesc_unencode_matrix.m的第三行m是M的阶数，50可以满足7*7-40*40的图像，噪音参数是40时，结果不稳定，小噪音下结果稳定。
4.可以通过调整样本矩和噪音的参数，得到不同图像下稳定的结果。
5.当图像是90*90时，噪音参数为40，M阶数为1，成功割图几率较大。
写的不是很详细啦。不过还是希望各位看官看的满意。瑞斯拜~~~