matlab中kmeans算法程序如下 我要做图像分类 主程序改怎么写那？知道的写下 谢谢了

function [mu,mask]=kmeans(ima,k)%k为指定类别数

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

% kmeans image segmentation

%

% Input:

% ima: grey color image灰度图像

% k: Number of classes指定的图像中类别数目

% Output:

% mu: vector of class means 每个类的均值

% mask: clasification image mask分类后的图像掩膜（mask)

%

% Author: Jose Vicente Manjon Herrera

%Email: jmanjon@fis.upv.es

% Date: 27-08-2005

%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% check image

ima=double(ima)

copy=ima% make a copy

ima=ima(:) % vectorize ima将图像向量化，即一维化。

mi=min(ima) % deal with negative

ima=ima-mi+1% and zero values

s=length(ima)%获得图像像素个数

% create image histogram%创建图像直方图

m=max(ima)+1%最大像素值加1

h=zeros(1,m)%直方图，有m个bin

hc=zeros(1,m)%标号矩阵，每个像素点的值为该点所隶属的类别号

for i=1:s%s是图像象素个数，即考查每个像素

if(ima(i)>0) h(ima(i))=h(ima(i))+1end%直方图中对应bin加1

end

ind=find(h)%找到直方图中不为零的那些bin的序号。

hl=length(ind)%直方图中非零bin的个数

% initiate centroids

mu=(1:k)*m/(k+1)%k为指定的类别数，mu为不同类的分割点，相当于坐标轴上的整点

% start process

while(true)

oldmu=mu

% current classification

for i=1:hl

c=abs(ind(i)-mu)%就是相当于考察ind(i)在坐标轴上离哪个整点最近！注意mu总共就k个

cc=find(c==min(c))%cc保留距离ind(i)最近整点的序号，序号为1、2、3...k

hc(ind(i))=cc(1)

end

%recalculation of means 下面的程序用于计算每一类的均值位置

for i=1:k,

a=find(hc==i)

mu(i)=sum(a.*h(a))/sum(h(a))%h为直方图

end

if(mu==oldmu) breakend%循环结束条件

end

% calculate mask

s=size(copy)

mask=zeros(s)

mask1=mask%增加一个显示矩阵

size(mask1)

for i=1:s(1),

for j=1:s(2),

c=abs(copy(i,j)-mu)

a=find(c==min(c))

mask(i,j)=a(1)

end

end

mu=mu+mi-1 % recover real range

for i = 1 : k

p=find(mask==i)

mask1(p)=1/k*i

end

figure,imshow(mask1)

[idx,c]=kmeans(X,k)

其中k是聚类中心个数

X是你存储需要处理的坐标的矩阵

c是一个存储了聚类中心点坐标的矩阵

MATLAB 是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。

MATLAB是matrix&laboratory两个词的组合，意为矩阵工厂(矩阵实验室)。是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。

MATLAB和Mathematica、Maple并称为三大数学软件。

MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且MATLAB也吸收了像Maple等软件的优点，使MATLAB成为一个强大的数学软件。

在新的版本中也加入了对C，FORTRAN，C++，JAVA的支持。

程序需要一个数据文件

格式如下：

5 2 3

2 3 4 5 10 12 5 1 12 10

其中，5表示数据的数量，2表示数据的维度，3表示聚类的数量。

后面每两个实数对应一个数据点。

不过这个程序始数据维数固定为2，后来想改成4以下的任意维度，但没有改完，所以其实只能为2维。

我已经对程序做了注释。

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include <conio.h>

#include <math.h>

#define SUCCESS 1

#define FAILURE 0

#define TRUE 1

#define FALSE 0

#define MAXVECTDIM 4 // 数据最大维数 (看来这个程序写了一半，维数实际上只能为2)

#define MAXPATTERN 1588 // 数据数量最大值

#define MAXCLUSTER 10 // 聚类最大值

// 聚类结构

struct aCluster {

double Center[MAXVECTDIM]// 中心/引力数据对象

int Member[MAXPATTERN] // 该聚类中数据在Pattern中的索引

int NumMembers // 数据的数量

}

struct aVector {

double Center[MAXVECTDIM]

int Size

}

static double Pattern[MAXPATTERN][MAXVECTDIM + 1]// 所有数据存放在这个数组中

// 所以的东西都搁System类里面了

class System {

private :

aCluster Cluster[MAXCLUSTER]// 聚类数组

int NumPatterns // 输入数据的数量

int SizeVector // 数据的维数

int NumClusters // 数据的聚类数

void DistributeSamples()// 根据中心聚类，重新分配数据到不同的聚类

int CalcNewClustCenters() // 重新计算新的聚类中心

double EucNorm(int, int)// 误差准则

int FindClosestCluster(int) // 查找最接近的聚类

public :

System() {}

int LoadPatterns(char * fname) // 从文件中读取数据

void InitClusters()// 初始化聚类

void RunKMeans() // 执行K-Means算法

void ShowClusters()// 显示聚类

void SaveClusters(char * fname)// 保存聚类到文件中

void ShowCenters() // 显示聚类中心数据

}

void System::ShowCenters() {

int i

printf("Cluster centers:\n")

for (i = 0i <NumClustersi++) {

Cluster[i].Member[0] = i

printf("ClusterCenter[%d]=(%f,%f)\n", i, Cluster[i].Center[0],Cluster[i].Center[1])

int b=0

}

printf("\n")

}

int System::LoadPatterns(char * fname) {

FILE * InFilePtr

int i, j

double x

if ( (InFilePtr = fopen(fname, "r")) == NULL)

return FAILURE

// 读入数据的数量，维度和聚类数的数量

fscanf(InFilePtr, "%d", &NumPatterns)

fscanf(InFilePtr, "%d", &SizeVector)

fscanf(InFilePtr, "%d", &NumClusters)

// 读入数据

for (i = 0i <NumPatternsi++) {

for (j = 0j <SizeVectorj++) {

fscanf(InFilePtr, "%lg", &x)

Pattern[i][j] = x

}

}

// 打印读入的数据

printf("Input patterns:\n")

for (i = 0i <NumPatternsi++) {

printf("Pattern[%d]=(%2.3f,%2.3f,%d,%d)\n", i, Pattern[i][0], Pattern[i][1],Pattern[i][2], Pattern[i][3])

}

printf("\n--------------------\n")

return SUCCESS

}

void System::InitClusters() {

int i, j

SizeVector=2// 看来开始数据维数固定为2，后来想改成4以下的任意维度，但没有改完

printf("Initial cluster centers:\n")

// 把前NumClusters个数据作为NumClusters个聚类的数据中心

for (i = 0i <NumClustersi++) {

Cluster[i].Member[0] = i// 记录这个数据的索引到第i个聚类中

for (j = 0j <SizeVectorj++) {

Cluster[i].Center[j] = Pattern[i][j]// 把这个数据作为数据中心

}

}

// 打印聚类的数据中心

for (i = 0i <NumClustersi++) {

printf("ClusterCenter[%d]=(%f,%f)\n", i, Cluster[i].Center[0], Cluster[i].Center[1])

}

printf("\n")

}

void System::RunKMeans() {

int converged// 是否收敛

int pass // 计算的趟数

pass = 1 // 第一趟

converged = FALSE// 没有收敛

while (converged == FALSE) { // 没有收敛的情况下反复跑

printf("PASS=%d\n", pass++)// 打印趟数

DistributeSamples()// 重新分配数据

converged = CalcNewClustCenters()// 计算新的聚类中心，如果计算结果和上次相同，认为已经收敛

ShowCenters()// 显示聚类中心

}

}

// 第p个数据到第c个聚类中心的误差准则(距离的平方)

double System::EucNorm(int p, int c) {

double dist, x// x可能是调试用，没有实际作用

int i

dist = 0

// 将每个维度的差的平方相加，得到距离的平方

for (i = 0i <SizeVectori++) {

x = (Cluster[c].Center[i] - Pattern[p][i]) *

(Cluster[c].Center[i] - Pattern[p][i])

dist += (Cluster[c].Center[i] - Pattern[p][i]) *

(Cluster[c].Center[i] - Pattern[p][i])

}

return dist

}

// 查找第pat个数据的最近聚类

int System ::FindClosestCluster(int pat) {

int i, ClustID/*最近聚类的ID*/

double MinDist/*最小误差*/, d

MinDist = 9.9e+99// 初始为一个极大的值

ClustID = -1

// 依次计算3个聚类到第pat个数据的误差，找出最小值

for (i = 0i <NumClustersi++) {

d = EucNorm(pat, i)

printf("Distance from pattern %d to cluster %d is %f\n", pat, i, sqrt(d))

if (d <MinDist) { // 如果小于前最小值，将改值作为当前最小值

MinDist = d

ClustID = i

}

}

if (ClustID <0) { // 没有找到？？！ —— 这个应该不可能发生，如果出现表示出现了不可思议的错误

printf("Aaargh")

exit(0)

}

return ClustID

}

void System ::DistributeSamples() {

int i, pat, Clustid, MemberIndex

// 将上次的记录的该聚类中的数据数量清0，重新开始分配数据

for (i = 0i <NumClustersi++) {

Cluster[i].NumMembers = 0

}

// 找出每个数据的最近聚类数据中心，并将该数据分配到该聚类

for (pat = 0pat <NumPatternspat++) {

Clustid = FindClosestCluster(pat)

printf("patern %d assigned to cluster %d\n\n", pat, Clustid)

MemberIndex = Cluster[Clustid].NumMembers // MemberIndex是当前记录的数据数量，也是新加入数据在数组中的位置

Cluster[Clustid].Member[MemberIndex] = pat// 将数据索引加入Member数组

Cluster[Clustid].NumMembers++ // 聚类中的数据数量加1

}

}

int System::CalcNewClustCenters() {

int ConvFlag, VectID, i, j, k

double tmp[MAXVECTDIM]// 临时记录新的聚类中心，因为需要比较，不能直接记录

char xs[255]

char szBuf[255]

ConvFlag = TRUE

printf("The new cluster centers are now calculated as:\n")

// 逐个计算NumClusters个聚类

for (i = 0i <NumClustersi++) {

strcpy(xs,"")

printf("Cluster Center %d (1/%d):\n",i,Cluster[i].NumMembers)

// tmp所有维度数值清0

for (j = 0j <SizeVectorj++) {

tmp[j] = 0.0

}

// 计算聚类中所有数据的所有维度数值的和，为下一步计算均值做准备

for (j = 0j <Cluster[i].NumMembersj++) {

VectID = Cluster[i].Member[j]

for (k = 0k <SizeVectork++) {

tmp[k] += Pattern[VectID][k]

sprintf(szBuf,"%d ",Pattern[VectID][k])

strcat(xs,szBuf)

}

printf("%s\n",xs)// 输出刚才计算的数据

strcpy(xs,"")

}

// 求出均值，并判断是否和上次相同

for (k = 0k <SizeVectork++) {

if(Cluster[i].NumMembers!=0)

tmp[k] = tmp[k] / Cluster[i].NumMembers

if (tmp[k] != Cluster[i].Center[k]) // 如果不同，则认为没有收敛

ConvFlag = FALSE

Cluster[i].Center[k] = tmp[k]// 用新值代替旧值

}

printf("%s,\n", xs)

}

return ConvFlag// 返回收敛情况

}

// 显示聚类中心数据

void System::ShowClusters() {

int cl

for (cl = 0cl <NumClusterscl++) {

printf("\nCLUSTER %d ==>[%f,%f]\n", cl, Cluster[cl].Center[0],

Cluster[cl].Center[1])

}

}

// 没有实现的函数

void System::SaveClusters(char * fname) {

}

#include <windows.h>

void main(int argc, char * argv[]) {

System kmeans

// 如果没有提供参数，显示用法

if (argc <2) {

printf("USAGE: KMEANS PATTERN_FILE\n")

exit(0)

}

// 参数1 为数据文件名，从中读入数据

if (kmeans.LoadPatterns(argv[1]) == FAILURE) {

printf("UNABLE TO READ PATTERN_FILE:%s\n", argv[1])

exit(0)

}

kmeans.InitClusters()

kmeans.RunKMeans()

kmeans.ShowClusters()

kmeans.ShowCenters()

return

}

---