利用最大流最小割算法matlab割图

目录

练习思路

matlab绘图

 噪音

坐标编码

邻接矩阵

最大流最小割算法

对最大流最小割算法求解结果转换为图像

源代码:

运行实例:

TIPS:

---

最近学习了最大流和最小割算法，可以把图看成是一些点的集合，色彩差值的倒数看做是权重，可以利用最小割实现简单的graph-cut操作。现在做个简单的分割椭圆的割图练习。

练习思路

先用matlab的矩形色彩绘图工具绘制出椭圆的图形，再添加随机噪音，使椭圆不太明显，接着将色块的坐标编码为一维数字，然后做出邻接矩阵，用matlab的最大流最小割函数求最大流，将分割出的数字解编码为坐标，重新赋值颜色，完成graph-cut操作。

matlab绘图

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
end

imagesc_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图像

TIPS:

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，成功割图几率较大。

写的不是很详细啦。不过还是希望各位看官看的满意。瑞斯拜~~~