基于DSP的MRF图像分割

发表于 | 分类于

基于MRF的FCM图像分割

一、目标

利用VisualDSP++ 5.0、仿真器、EBF-561实验平台实现该算法，并通过该聚类算法实现图像的不同区域的聚类结果。

二、实现方案

实验原理

MRF理论提供了建模上下文依赖实体的一种方式，实体包括图像像素和相关特征等。

在MRF中，令F={F1,…..,Fm}是一组定义在集合S上的随机变量，其中Fi代表在标签集L中的一个取值fi，Fi=fi代表取值为fi的事件，（F1=f1,…,Fm=fm）代表联合事件（joint event）。假设f={f1,….fm}，那么刚才叙述的联合事件则可简记为F=f。对于离散的标签集合L，Fi=fi的概率记作P(Fi=fi)，简称P(fi)；联合事件的概率记作P(f)。对于连续的标签集合L，我们可以知道pdf p(Fi=fi)和p(F=f)。

对于一个邻域系统N，如果如下两条性质满足：

  1. P(f) >0 , 对任意的f。（非负性）
  2. P(fi| fs-{i})=P(fi|fNi)。（马尔科夫性）

那么，我们就称F是一个马尔科夫随机场。

上述1.2中，S-{i}是i除i之外的所有点，如果对于图像来讲，i是其中一个像素，那么S-{i}就是除i之外的其他所有像素；fs-{i}代表的是S-{i}中的所有点的标签集合；fNi={fi’|i’∈Ni}表示点i的邻域的点的标签集合。

实验步骤

Step1：给定图像初始分割（通过阈值法或聚类方法）；
Step2：由当前分割更新Pi={Ui,Ei}分别是当前第l类区域的均值和标准方差；
Step3：由当前图像参数和上次迭代的分割结果，并计算每一点最大可能的类别；
Step4：判断是否收敛或达到了最高迭代次数，如果满足则退出；否则返回
Step2，进行下一次迭代。

三、核心代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
unsigned int k,m,temp;
double mean1,mean2,mean3;
double var1,var2,var3;
unsigned int u,d,l,r,lu,ld,ru,rd;
int imageData[84][64] , labelData[84][64];
double likelihood[2][2], likelihood0[2][2], likelihood1[2][2], likelihood2[2][2], prior[2][2], posterior[2][2];

/*************用Kmeans对图像进行初始分类************/

int cen1 = image[1], cen2 = image[100],cen3=image[1000];//初始化聚类中心
int sum01 = cen1, sum02 = cen2, sum03=cen3;
int count1 = 0, count2 = 0, count3=0;
int kk = 10;
int label[118*64];
while (kk--)
{
	for (i = 0; i < 84*64; i++)
	{
	    if(abs(image[i]-cen1) < abs(image[i]-cen2))
		{
	        if(abs(image[i]-cen1) < abs(image[i]-cen3))
		    {
			  sum01 += image[i];
		      count1++;
		    }	
		    else
			{
			  sum03 += image[i];
		      count3++;
			}		
		}
	   else
		{
		    if(abs(image[i]-cen2) < abs(image[i]-cen3))
			{
			    sum02 += image[i];
		    	count2++;
			}		
		    else
			{
			    sum03 += image[i];
		    	count3++;
			}						
		}
	}
	cen1 = sum01/count1;
	cen2 = sum02/count2;
	cen3 = sum03/count3;
	count1 = count2 = count3= 0;
	sum01 = sum02 = sum03= 0;
}
for(i = 0; i < 84*64; i++)
{
	if(abs(image[i]-cen1) < abs(image[i]-cen2))
	{
	    if(abs(image[i]-cen1) < abs(image[i]-cen3))
		{
			image[i]=0;
		    label[i]=1;
		}		
		else  
        {
			image[i]=255;
            label[i]=3;
        }
			
	}
	else
	{
		if(abs(image[i]-cen2) < abs(image[i]-cen3))
		{
			image[i]=128;
			label[i]=2;
		}		
		else   
        {
			image[i]=255;
            label[i]=3;
        }			
	}
}

/************将数据转成2维数组***********/
for(i=0;i<Chang;i++)
{
	for(j=0;j<Kuan;j++)
	{
		imageData[i][j] = image[i*Kuan+j];
		labelData[i][j] = label[i*Kuan+j];
	}
}

/******************求三类图像的均值***********/
for(m=0;m<20;m++)
{
   double sum1=0,sum2=0,sum3=0;
   int a1=0,a2=0,a3=0;
  	for(i=0;i<Chang;i++)
	{
		for(j=0;j<Kuan;j++)
		{
			temp = labelData[i][j];
			if(temp==0)
			{
				sum1+=imageData[i][j];
				a1++;
			}
			else if(temp == 1)
			{
				sum2+=imageData[i][j];
				a2++;
			}
			else
            {
                sum3+=imageData[i][j];
				a3++;
            }
		}
	}
	mean1 = sum1/a1;
	mean2 = sum2/a2;
	mean3 = sum3/a3;
	sum1=0;
	sum2=0;
	sum3=0;
}
	
/******************求方差*****************/
   	for(i=0;i<Chang;i++)
	{
		for(j=0;j<Kuan;j++)
		{
		    temp = labelData[i][j];
				if(temp==0)
				{
					sum1 = sum1+pow(imageData[i][j]-mean1, 2);
				}
				else if(temp ==1)
				{
                    sum2 = sum2+pow(imageData[i][j]-mean2, 2);
				}
				else
                {
                     sum3 = sum3+pow(imageData[i][j]-mean3, 2);
                }
		}
	}
	var1 = sum1/a1;
	var2 = sum2/a2;
	var3 = sum3/a3;
	
/******************求似然值*******************/
for (k = 0;k<N;k++)
{
 	likelihood[0][k] = 1/sqrt(2*PI*var1)*exp(-pow((image[k] - mean1),2)/2/var1);
    likelihood[1][k] = 1/sqrt(2*PI*var2)*exp(-pow((image[k]- mean2),2)/2/var2);
    likelihood[2][k] = 1/sqrt(2*PI*var3)*exp(-pow((image[k] - mean3),2)/2/var3);
}
for(i=0;i<Chang;i++)
{
	for(j=0;j<Kuan;j++)
		{
			likelihood0[i][j] = likelihood[0][i*Kuan+j];
			likelihood1[i][j] = likelihood[1][i*Kuan+j];
			likelihood2[i][j] = likelihood[2][i*Kuan+j];
		}
}

/**************求先验概率****************/
double V1,V2,V3;
for(i=1;i<Chang-1;i++)
{
	for(j=1;j<Kuan-1;j++)
	{
		u = labelData[i][j-1];d = labelData[i][j+1];l = labelData[i-1][j];r = labelData[i+1][j];
        lu = labelData[i-1][j-1];ld = labelData[i-1][j+1];ru = labelData[i+1][j-1];rd = labelData[i+1][j+1];
        temp = 0;
        V1 = (!(temp - u)-1) + (!(temp - d)-1) + (!(temp - l)-1) + (!(temp - r)-1) + (!(temp - lu)-1) + (!(temp - ru)-1) + (!(temp - ld) -1)+(!(temp - rd)-1);
        temp = 1;
        V2 = (!(temp - u)-1) + (!(temp - d)-1) + (!(temp - l)-1) + (!(temp - r)-1) + (!(temp - lu)-1) + (!(temp - ru)-1) + (!(temp - ld) -1)+(!(temp - rd)-1);
        temp = 2;
        V3 = (!(temp - u)-1) + (!(temp - d)-1) + (!(temp - l)-1) + (!(temp - r)-1) + (!(temp - lu)-1) + (!(temp - ru)-1) + (!(temp - ld) -1)+(!(temp - rd)-1);
        k = (i-1)*(Kuan-2)+(j-1);
        prior[0][k] = exp(V1)/(exp(V1)+exp(V2)+exp(V3));
        prior[1][k] = exp(V2)/(exp(V1)+exp(V2)+exp(V3));
        prior[2][k] = exp(V3)/(exp(V1)+exp(V2)+exp(V3));
    }
}

/**************求后验概率并分类***************/
    for (i = 1;i<Chang-1;i++)
    {
        for(j = 1;j<Kuan-1;j++)
        {
            k = (i-1)*(Kuan-2)+(j-1);
            posterior[0][k] = likelihood0[i][j]*prior[0][k];
            posterior[1][k] = likelihood1[i][j]*prior[1][k];
            posterior[2][k] = likelihood2[i][j]*prior[2][k];
        }
    }
    for(i = 1;i<Chang;i++)
    {
        for(j = 1;j <Kuan;j++)
        {
            k = (i-1)*(Kuan-2)+(j-1);
        if (posterior[0][k]<posterior[1][k])
        {
            if  (posterior[1][k]<posterior[2][k]) labelData[i][j] = 2;
            else labelData[i][j] = 1;
        }
        else
        {
            if (posterior[1][k]<posterior[2][k])
            {
                if (posterior[0][k]>posterior[2][k])  labelData[i][j] = 0;
                else labelData[i][j] = 2;
            
            else labelData[i][j] = 0;
        }
    }
     

/******************处理图像***************8***/
    for( i = 0;i<Chang;i++)
        {
            for(j=0;j<Kuan;j++)
            {
                temp = labelData[i][j];
                if (temp == 0)
                    image[i*Kuan+j] = 0;
                else if (temp == 1)
                    image[i*Kuan+j] = 128;
                else
                    image[i*Kuan+j] = 255;
            }
        }