- 实验目的
掌握模糊c-均值(FCM)聚类算法基本原理,并基于C语言编程实现其算法过程,应用于图像的分割实验。
- 实验基本原理及步骤(或方案设计及理论计算)
1)原理
模糊聚类算法FCM目标函数为:
如果 p 表示每一个样本的维数,是一个矩阵;N表示样本数目,通常表示图像像素数;C表示聚类数目;是矢量隶属于第 i 类的隶属度函数 ,满足且;聚类中心是矩阵,和更新等式分别为:
对于每一个模糊隶属度,由控制模糊度的权重系数;为相似性度量。
2)变量说明
P为数据样本维数(灰度图像时为1);
N为像素点数目;
为像素 j 特征(灰度图像时表示灰度值);
C为图像分割类别数;
为像素 j 属于第 i 类的隶属度;
为第 i 类聚类中心;
3)算法步骤
- 设置目标函数精度,模糊指数m(m通常取2),最大迭代次数;
- 初始化聚类中心;
- 由式(1)更新模糊划分矩阵和聚类中心;
- 若或者则聚类结束;否则 并转上一步;
- 由所得得到各像素点分类结果。
- 实验结果分析及回答问题(或测试环境及测试结果)
1)原图
2)聚类结果图
如上图即为最终的聚类分割结果,本次聚类一共设置了三个聚类中心。当迭代结束时,每一个聚类中心的值都会赋给各自的样本,所以结果会有三种不同灰度的颜色。在运行算法前,应注意初始聚类中心的选取,若选取不当,则会导致聚类结果只要两种颜色的结果。
- 实验代码
#include<stdio.h>
#include<cstring>
#include<math.h>
#include<windows.h>
#include <time.h>
#define max_pixel 99999999
#define max_k 10
int bmpHeight, bmpWidth, biBitCount, k_num;
int label[max_pixel],rec[max_k];
int center[max_k][3];
float u[max_pixel][max_k][3];
unsigned char *pBmpBuf;
RGBQUAD *pColorTable;
long long int cacJ(int rgb[][3],int h,int w){
long long int J = 0;
for (int k = 0; k < k_num;k++)
for (int j = 0; j < h*w ;j++)
for (int i = 0; i < 3;i++)
J += u[j][k][i] * pow(center[k][i] - rgb[j][i],2);
return J;
}
void fcm(int rgb[][3],int h,int w,float b){
for (int i = 0; i < h * w;i++)
label[i] = -1;
time_t t;
srand((unsigned) time(&t));
for (int i = 0; i < k_num;i++){
int t = rand() % (h*w);
center[i][0] = rgb[t][0];
center[i][1] = rgb[t][1];
center[i][2] = rgb[t][2];
label[t] = i;
}
long long int J = 99999;
long long int L = J+1000;
unsigned int i, j, k, l;
i = 0; j = 0; k = 0;l = 0;
do{
float tmp = 0;
for(j = 0; j < h*w; j++){
for(i = 0; i < 3; i++){
tmp = 0;
for(k = 0; k < k_num; k++){
for(l = 0; l < k_num; l++){
tmp += pow(rgb[j][i]-center[l][i],(-1)/(b-1));
}
}
u[j][k][i] = pow(rgb[j][i]-center[k][i],(-1)/(b-1)) / tmp;
}
}
L = J;
J = cacJ(rgb, h, w);
float tmp1 = 0, tmp2 = 0;
for(k = 0; k < k_num; k++){
tmp1 = 0; tmp2 = 0;
for(i = 0; i < 3; i++){
for(j = 0; j < h*w; j++){
tmp1 += pow(u[j][k][i],b) * rgb[j][i];
tmp2 += pow(u[j][k][i],b);
}
}
center[k][i] = tmp1 / tmp2;
}
}while (((L-J) > 10));
}
int main(){
const char src[] = "test.bmp";
const char dst[] = "segResultTest.bmp";
k_num = 3;
memset(label, 0, sizeof(label));
memset(center, 0, sizeof(center));
memset(rec, 0, sizeof(rec));
//读取 bmp 文件
FILE *fp = fopen(src, "rb");
if (fp == 0)
return 0;
fseek(fp, sizeof(BITMAPFILEHEADER), 0);
BITMAPINFOHEADER head;
fread(&head, 40, 1, fp);
bmpHeight = head.biHeight;
bmpWidth = head.biWidth;
biBitCount = head.biBitCount;
fseek(fp, sizeof(RGBQUAD), 1);
int LineByte = (bmpWidth*biBitCount / 8 + 3) / 4 * 4;
pBmpBuf = new unsigned char[LineByte*bmpHeight];
fread(pBmpBuf, LineByte*bmpHeight, 1, fp);
fclose(fp);
FILE *fp1 = fopen(dst, "wb");
if (fp1 == 0)
return 0;
int LineByte1 = (bmpWidth * 8 / 8 + 3) / 4 * 4;
BITMAPFILEHEADER bfhead;
bfhead.bfType = 0x4D42;
bfhead.bfSize = 14 + 40 + 256 * sizeof(RGBQUAD)+LineByte1*bmpHeight;
bfhead.bfReserved1 = 0;
bfhead.bfReserved2 = 0;
bfhead.bfOffBits = 14 + 40 + 256 * sizeof(RGBQUAD);
fwrite(&bfhead, 14, 1, fp1);
BITMAPINFOHEADER head1;
head1.biBitCount = 8;
head1.biClrImportant = 0;
head1.biCompression = 0;
head1.biClrUsed = 0;
head1.biHeight = bmpHeight;
head1.biWidth = bmpWidth;
head1.biPlanes = 1;
head1.biSize = 40;
head1.biSizeImage = LineByte1*bmpHeight;//修改位图数据的大小
head1.biXPelsPerMeter = 0;
head1.biYPelsPerMeter = 0;
fwrite(&head1, 40, 1, fp1); //将修改后的信息头存入 fp1;
pColorTable = new RGBQUAD[256];
for (int i = 0; i < 256; i++){
pColorTable[i].rgbRed = i;
pColorTable[i].rgbGreen = i;
pColorTable[i].rgbBlue = i;
}
fwrite(pColorTable, sizeof(RGBQUAD), 256, fp1); //将颜色表写入 fp1
//写位图数据
unsigned char *newBmp;
newBmp = new unsigned char[LineByte1*bmpHeight];
int c = 0;
for (int i = 0; i < bmpHeight * bmpWidth;i++,c++)
label[i] = -1;
int rgbMap[bmpHeight * bmpWidth][3];
for (int i = 0; i < bmpHeight;i++){ //取位图 rgb 三通道数据
for (int j = 0; j < bmpWidth;j++){
unsigned char *tmp;
tmp = pBmpBuf + i * LineByte + j * 3;
rgbMap[i * bmpWidth + j][0] = int(*(tmp));
rgbMap[i * bmpWidth + j][1] = int(*(tmp+1));
rgbMap[i * bmpWidth + j][2] = int(*(tmp+2));
}
}
fcm(rgbMap, bmpHeight, bmpWidth,2);
for (int i = 0; i < bmpHeight * bmpWidth;i++){
rec[label[i]]++;
label[i]=int(center[label[i]][0]*0.299+center[label[i]][0]*0.587+center[label[i]][0]*0.114);
}
for (int i = 0;i<k_num;i++)
printf("%d ", rec[i]);
for (int i = 0; i < bmpHeight; i++){
for (int j = 0; j < bmpWidth; j++){
unsigned char *pb;
pb = newBmp + i * LineByte1 + j;
*pb = label[i * bmpWidth + j];
}
}
fwrite(newBmp, LineByte1*bmpHeight, 1, fp1);
fclose(fp1);
system("pause");
return 0;
}