计算机视觉

一、图像处理

main

denoise

二、Harris角点检测

三、Hough变换直线检测

四、直方图显著性检测

五、人脸识别

六、kmeans

import

函数

kmeanstext

七、神经网络

常用函数：

imread----------读取图像
imshow---------显示图像
rgb2hsv---------RGB转HSV
hsv2rgb---------HSV转RGB
imhist-----------显示图像直方图
rgb2gray-------彩色图转灰度图
imnoise--------对图像添加噪声
fspecial--------构建均值滤波模板
imfilter--------均值滤波函数
imgaussfilt-----高斯滤波函数

一、图像处理

main

clc
clear all
img = double(imread('Fig1.bmp'));

R = img(:,:,1);
G = img(:,:,2);
B = img(:,:,3);

figure(1);
subplot(2,2,1), imshow(uint8(img)),title('原图');
subplot(2,2,2), imshow(uint8(R));title('R分量图');
subplot(2,2,3), imshow(uint8(G));title('B分量图');
subplot(2,2,4), imshow(uint8(B));title('B分量图');
%rgb2hsv
img_hsv = rgb2hsv(img);
img_hsv(:,:,1) = img_hsv(:,:,1)*3;
img_hsv(:,:,2) = img_hsv(:,:,2);
img_hsv(:,:,3) = img_hsv(:,:,3);
img_rgbh = hsv2rgb(img_hsv);

img_hsv = rgb2hsv(img);
img_hsv(:,:,1) = img_hsv(:,:,1);
img_hsv(:,:,2) = img_hsv(:,:,2)*3;
img_hsv(:,:,3) = img_hsv(:,:,3);
img_rgbs = hsv2rgb(img_hsv);

img_hsv = rgb2hsv(img);
img_hsv(:,:,1) = img_hsv(:,:,1);
img_hsv(:,:,2) = img_hsv(:,:,2);
img_hsv(:,:,3) = img_hsv(:,:,3)*3;
img_rgbv = hsv2rgb(img_hsv);

figure(2);
subplot(2,2,1), imshow(uint8(img)),title('原图');
subplot(2,2,2), imshow(uint8(img_rgbh));title('H增强分量图');
subplot(2,2,3), imshow(uint8(img_rgbs));title('S增强分量图');
subplot(2,2,4), imshow(uint8(img_rgbv));title('V增强分量图');
%直方图

figure(3);
subplot(2,2,1), imshow(uint8(img)),title('原图');
subplot(2,2,2), imhist(uint8(img(:,:,1)));title('R分量直方图');
subplot(2,2,3), imhist(uint8(img(:,:,2)));title('B分量直方图');
subplot(2,2,4), imhist(uint8(img(:,:,3)));title('B分量直方图');

denoise

clc
clear all

img = imread('Fig6.bmp');

figure(1)
subplot(1,1,1),imshow(uint8(img)),title('原图');
img = rgb2gray(img);
%添加噪声
im_noise_g = imnoise(img, 'gaussian',0,0.01);
im_noise_s = imnoise(img, 'salt & pepper',0.03);

figure(2)
subplot(1,2,1),imshow(uint8(im_noise_g)),title('高斯噪声图');
subplot(1,2,2),imshow(uint8(im_noise_s)),title('椒盐噪声图');

%均值滤波
A_3 = fspecial('average',[3,3]);
A_9 = fspecial('average',[9,9]);

im_denoise_g_3 = imfilter(im_noise_g,A_3);
im_denoise_g_9 = imfilter(im_noise_g,A_9);

figure(3)
subplot(1,2,1),imshow(uint8(im_denoise_g_3)),title('高斯去噪图-均值—3*3');
subplot(1,2,2),imshow(uint8(im_denoise_g_9)),title('高斯去噪图-均值—9*9');

%高斯滤波

im_denoise_g_1 = imgaussfilt(im_noise_g,1);
im_denoise_g_5 = imgaussfilt(im_noise_g,5);

figure(4)
subplot(1,2,1),imshow(uint8(im_denoise_g_1)),title('高斯去噪图-高斯-1');
subplot(1,2,2),imshow(uint8(im_denoise_g_5)),title('高斯去噪图-高斯-5');

二、Harris角点检测

clear;
% 读取图像
grayImage = imread('1.bmp');
% 转化为灰度图像
%grayImage=rgb2gray(srcImage);
% 求取图像宽高
[ImageHeight,ImageWidth]=size(grayImage);
% 显示原始灰度图
% imshow(ori_im);
% 方法1：x方向梯度算子(用于Harris角点提取算法)
% fx = [-2 -1 0 1 2];
% 方法2：x方向高斯卷积核(高斯尺度空间的多尺度优化)
%fx = [5 0 -5;8 0 -8;5 0 -5];
fx = [1,2,1;0,0,0;-1,-2,-1];

% x方向滤波微分
Ix = filter2(fx,grayImage);
% 显示x方向滤波图
figure;imshow(Ix);
% 方法1：y方向梯度算子(用于Harris角点提取算法)
% fy = [-2;-1;0;1;2];
% 方法2：y方向高斯卷积核(高斯尺度空间的多尺度优化)
%fy = [5 8 5;0 0 0;-5 -8 -5];
fy = [1,0,-1;2,0,-2;1,0,-1];

% y方向滤波微分
Iy = filter2(fy,grayImage);
% 显示y方向滤波图
figure;imshow(Iy);
%（相关参数说明见harris理论，文中前面有链接）
Ix2 = Ix.^2;
Iy2 = Iy.^2;
Ixy = Ix.*Iy;
% 产生7*7的高斯窗函数，sigma=2,用于窗口的高斯平滑
w= fspecial('gaussian',[5 5],2);
Ix2 = filter2(w,Ix2);
Iy2 = filter2(w,Iy2);
Ixy = filter2(w,Ixy);

% 纪录角点位置，角点处值为1
corner = zeros(ImageHeight,ImageWidth);
% 图像中最大的响应值
Rmax = 0;
% 计算各点的响应值
R = zeros(ImageHeight,ImageWidth);
for i = 1:ImageHeight
for j = 1:ImageWidth
M = [Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)];
R(i,j) = det(M)-0.06*(trace(M))^2;
% if R(i,j) > Rmax
% Rmax = R(i,j);
% end;
end;
end;
% 角点个数
cnt = 0;
% 进行非极大抑制，窗口大小3*3
for i = 2:ImageHeight-1
for j = 2:ImageWidth-1
%if ( R(i,j) > 0.01*Rmax && R(i,j) > R(i-1,j-1) && R(i,j) > R(i-1,j) && R(i,j) > R(i-1,j+1) && R(i,j) > R(i,j-1) && R(i,j) > R(i,j+1) && R(i,j) > R(i+1,j-1) && R(i,j) > R(i+1,j) && R(i,j) > R(i+1,j+1))
if ( R(i,j) > 20 && R(i,j) > R(i-1,j-1) && R(i,j) > R(i-1,j) && R(i,j) > R(i-1,j+1) && R(i,j) > R(i,j-1) && R(i,j) > R(i,j+1) && R(i,j) > R(i+1,j-1) && R(i,j) > R(i+1,j) && R(i,j) > R(i+1,j+1))
corner(i,j) = 1;
cnt = cnt+1;
end;
end;
end;

[upix, vpix] = find(corner == 1);
%角点个数
cnt
%绘制角点
figure;
imshow(grayImage)
hold on;
plot(vpix,upix,'r.');

三、Hough变换直线检测

clear;
% 读取图像
I = imread('1.bmp');

subplot(2,2,1);
imshow(I);
title('原图');

BW = edge(I,'canny');
subplot(2,2,2);
imshow(BW);
title('canny算子边缘检测后图像');

%=====hough函数===============
%BW：m*n维二值图
% RhoResolution：ρ的分辨率，缺省为1
% ThetaResolution：θ的分辨率，θ的范围是[-90,90)
% H：hough变换后的累加器矩阵，行数为m，列数为n
% theta：θ
% rho：ρ theta,rho为计算霍夫变换的角度和半径值

[H,T,R] = hough(BW);
subplot(2,2,3);
imshow(H,[],'XData',T,'YData',R,'InitialMagnification','fit');
% imshow(H,[]);
title('Hough变换图');
axis on
axis normal
xlabel('\theta')
ylabel('\rho')
hold on
%找到Hough变换结果，累加容器H中所有大于阈值Threshold的峰值中，最大的2个峰值所在的位置P。
% P为矩阵，记录H中被找到的峰值的横纵坐标。
%再hough矩阵中寻找前40个大于hough矩阵中最大值0.24倍峰值
% P = houghpeaks(H,10,'threshold',ceil(0.3*max(H(:)))); %P是峰值的坐标
P = houghpeaks(H,10); %P是峰值的坐标

x = T(P(:,2));
y = R(P(:,1)); %由行列索引转换成实际坐标

plot(x,y,'s','color','white'); %再hough矩阵图像中标出峰值位置

%二值图BW中对应指定（theta，rho），且容器中数值>峰值P的所有线段lines
% 且线段最小长度大于MinLength
% 线段间距离<FillGap的合并
% lines为结构变量，(lines.point1,lines.point2)记录满足条件的线段端点

lines = houghlines(BW,T,R,P,'FillGap',5,'MinLength',7); %合并距离小于5的线段，丢弃所有长度小于7的线段

subplot(2,2,4);
imshow(I);
title('hough变换检测图');
axis on;
hold on;

% mex_len = 0;
for k=1:length(lines) %以此标注各条直线
xy = [lines(k).point1;lines(k).point2];
plot(xy(:,1),xy(:,2),'LineWidth',2,'Color','red');
hold on
% %确定最长线段的端点
% len = norm(lines(k).point1 - lines(k).point2); %最长线段的长度
% if(len > max_len)
% mex_len = len;
% xy_long = xy;
% end
end

% plot(xy_long(:,1),xy_long(:,2),'LineWidth',2,'Color','cyan');

四、直方图显著性检测

%读取图像
%orpic = imread(img_name);
clear all;
clc;
close all;
orpic = imread('test2.bmp');
figure;imshow(orpic);
%获取图像的长、宽和通道值
[row ,column, nn] = size(orpic);
%生成调色板pallet(color,number),默认调色板大小为图像大小
pallet = zeros(row*column,2);
%设置阶梯值，为了将rgb量化后的结果化成一个12进制的数
w = [12*12,12,1];%RGB通道各自的权重，这里的权重不能设置为一样的，这是为了保证不一样的rgb组合得到的值不同
idx = uint32(1);%调色板像数值个数
SUM = uint32(0);%量化和
RGBSUM = zeros(row,column);%原图像为三维彩图，将RGB通道加权后每个像素的值
for i = 1 : row
for j = 1 : column
%获取第i行第j列的R、G、B通道的值，并量化为12阶后再转化为一个数
a1 = floor((double(orpic(i,j,1))/255)*(12-0.0001))*w(3); %乘以权重不一样就是为了每个通道量化后乘以权重的值不重叠了
a2 = floor((double(orpic(i,j,2))/255)*(12-0.0001))*w(2);
a3 = floor((double(orpic(i,j,3))/255)*(12-0.0001))*w(1);
SUM = int32(a1) + int32(a2) + int32(a3); %三个通道量化后颜色和作为计算相似性的依据
% SUM = 3* int32(a1);
RGBSUM(i,j) = SUM;

%对于每个像数值，将其归类到调色板中，先在已经归类的pallet中查找，如果遍历完毕发现未归类，则将索引+1并归类
%找出处理后的数据有多少种颜色值，以及其频数

%记录不同颜色出现的次数和出现的频率 pallet第一列是颜色第二列是该颜色出现的次数
if idx == 1 %第一次循环也就是该颜色第一次出现
pallet(idx,1) = SUM;
pallet(idx,2) = pallet(idx,2) + 1;
idx = idx + 1; %颜色值+1
continue;
end
flag = 0;
for m = 1 : idx-1 %检查有没有重复的颜色如果没有 flag=0 里面已经有重复的颜色了，那flag=1，对应的频率值+1
if pallet(m,1) == SUM
pallet(m,2) = pallet(m,2) + 1;
flag = 1;
break;
end
end
if flag == 0
pallet(idx,1) = SUM;
pallet(idx,2) = pallet(idx,2) + 1;
idx = idx + 1;
end
end
end
%获取图中所有的颜色值
idx = idx -1;
pallet = pallet(1:idx,:);%调色板截断，保留真实调色个数

%生成新的调色板num(number,color)，使得颜色按照频数从大到小排列，且一号位为频数大小，二号位为处理后的颜色值
palletemp = pallet(:,2); %把频数取出来
[b,pos] = sort(palletemp,'descend'); %频数降序排列
num = zeros(idx,2);
for i = 1 : idx
num(i,1) = b(i); %pos是降序的位置 b是对应的频数
num(i,2) = pallet(pos(i),1); %pos位置对应的颜色值
%num第一列频数第二列颜色值
end

%计算如果要保留95%的处理后的颜色值，需要保留调色板的多少位
maxnum = idx; %本来是257个颜色值
maxdropnum = floor(row*column*(1-0.95)); %算出覆盖面积 5% 对应多少个像素maxdropnum
crnt = num(maxnum,1); %颜色频数（从大到小排序）这是最后一个的频数
while(crnt<maxdropnum && maxnum>1) %覆盖面积小并且不是第一个
crnt = crnt + num(maxnum-1,1); %crnt 累加倒数第二个颜色频率（颜色频率，也就是这个颜色有多少个像素）
maxnum = maxnum - 1;
end %直到颜色对应的个数大于5%的像素个数了结束
%保证保留的颜色值个数小于256且尽量不太小
maxnum = min(maxnum,256);
if maxnum < 10
maxnum = min(idx,100);
end

%将处理后的颜色值转化回RGB对应的值，此处会有+-1的误差，影响不大
color3i = zeros(idx,3);
for i = 1:idx %是吧257个颜色返回
color3i(i,1) = round(num(i,2)/w(1));
color3i(i,2) = round(mod(num(i,2),w(1))/w(2));
color3i(i,3) = round(mod(num(i,2),w(2)));
end

%建立色板pallet(number,color,pos)，其中pos为保留颜色的序号
a = zeros(idx,1);
pallet = [num,a]; %num第一列频数第二列颜色值第三列序号
for i = 1:maxnum
pallet(i,3) = i;
end

%将舍弃的颜色归到其最相近的颜色当中
for i = (maxnum+1):idx %idx 257 原来的颜色个数 maxnum去掉95%以后的颜色个数要把后面这么多颜色归到周围的颜色中去
simidx = 99999999;simval = 99999999;
for j = 1:maxnum %看看后面去掉的颜色和前面哪些颜色接近计算距离
d_ij = round((color3i(i,1)-color3i(j,1))^2 + (color3i(i,2)-color3i(j,2))^2 + (color3i(i,3)-color3i(j,3))^2);
if d_ij < simval
simval = d_ij;
simidx = j;
end
end %把最小距离的颜色记录下来，那么就把i这个颜色归类到距离最小颜色
pallet(simidx,1) = pallet(simidx,1) + pallet(i,1); %对应的颜色频数+1
pallet(i,3) = simidx; %保留序号就是标记原来的这个颜色归去第4行了 simidx = 4
end

%将图像中属于某一“颜色值”的所有颜色的rgb值累加之后求平均再归一，即调整保留下来颜色值
rgbcolor = zeros(idx,4);
for n =1:idx
for i = 1:row
for j = 1:column
if RGBSUM(i,j) == pallet(n,2) %num第一列频数第二列颜色值第三列序号
rgbcolor(n,1) = int32(rgbcolor(n,1)) + int32(orpic(i,j,1));
rgbcolor(n,2) = int32(rgbcolor(n,2)) + int32(orpic(i,j,2));
rgbcolor(n,3) = int32(rgbcolor(n,3)) + int32(orpic(i,j,3));
rgbcolor(n,4) = double(pallet(n,3));
end
end
end
end
for i = 1:maxnum %maxnum 去掉了5%的颜色个数
for j = (maxnum+1):idx
if rgbcolor(j,4) == i
rgbcolor(i,1) = (rgbcolor(i,1) + rgbcolor(j,1));
rgbcolor(i,2) = (rgbcolor(i,2) + rgbcolor(j,2));
rgbcolor(i,3) = (rgbcolor(i,3) + rgbcolor(j,3));
end
end
rgbcolor(i,1) = (rgbcolor(i,1)/pallet(i,1))/255;
rgbcolor(i,2) = (rgbcolor(i,2)/pallet(i,1))/255;
rgbcolor(i,3) = (rgbcolor(i,3)/pallet(i,1))/255;
end

%将颜色格式由RGB转为Lab
labcolor = zeros(maxnum,3);
for i = 1:maxnum
a = zeros(1,1,3);
a(1,1,1) = rgbcolor(i,1);
a(1,1,2) = rgbcolor(i,2);
a(1,1,3) = rgbcolor(i,3);
b = rgb2lab(a);
labcolor(i,1) = b(1,1,1);
labcolor(i,2) = b(1,1,2);
labcolor(i,3) = b(1,1,3);
end

%存放颜色i与颜色j在Lab空间上的距离
distemp = zeros(maxnum,maxnum,2);
%每行将颜色j与颜色i的在Lab空间上的距离按从小到大排列
dist = zeros(maxnum,maxnum,2);
%存放平滑前的颜色显著值
colorsaltemp = zeros(maxnum,2);
for i = 1:maxnum
for j = 1:maxnum
%计算颜色i与其他所有j个颜色在Lab空间上的距离存下来
distemp(i,j,1) = norm(labcolor(i,:)-labcolor(j,:),2);
distemp(i,j,2) = j;
end
for k = 1:maxnum
if k == i
continue;
end
%颜色i在平滑之前的颜色显著值，为S(c)=ΣP(c)D(c,a)，a为其他的颜色概率相加的显著性值
colorsaltemp(i,1) = colorsaltemp(i,1) + (pallet(k,1)/(row*column))*distemp(i,k,1); %num第一列频数第二列颜色值第三列序号
colorsaltemp(i,2) = i;
end

% [b,pos] = sort(distemp(i,:,1),'ascend');
% for k = 1:maxnum
% dist(i,k,1) = b(k);
% dist(i,k,2) = distemp(i,pos(k),2); % dist里面是按照距离排序存放的 b是距离值 pos是位置
% end
end

% % %对各颜色进行平滑操作
% colorsal = zeros(maxnum,2);
% %计算权值时使用的颜色i与和其最邻近颜色的个数
% n = double(round(maxnum/4));
% %保存与颜色i在Lab空间上最相近的n=maxnum/4个点的距离的和
% totaldist = zeros(maxnum,1);
% for i = 1:maxnum
% %计算与颜色i在Lab空间上最相近的n=maxnum/4个点的距离的和
% for j = 2:n
% totaldist(i,1) = totaldist(i,1) + dist(i,j,1);
% end
% valcrnt = 0;
% for j = 1:n
% valcrnt = valcrnt + colorsaltemp(dist(i,j,2),1)*(totaldist(i,1)-dist(i,j,1)); %距离差乘以距离值作为周围四个颜色值的权重
% end
% colorsal(i,1) = valcrnt/((n-1)*totaldist(i,1));
% colorsal(i,2) = i;
% end

%保存原图像各像素的显著值
sal = zeros(row,column);
for i = 1:row
for j = 1:column
for n = 1:idx
if RGBSUM(i,j) == pallet(n,2)
%sal(i,j) = colorsal(pallet(n,3),1);
sal(i,j) = colorsaltemp(pallet(n,3),1);
end
end
end
end
saliency = (sal - min(min(sal)))/(max(max(sal)) - min(min(sal)));
A=fspecial('gaussian',[3,3],2);
saliency=imfilter(saliency,A,'replicate');
%显示显著性图
subplot(1,2,1),imshow(saliency);
subplot(1,2,2),imshow(sal,[]);

%读取图像并显示原图像

将图像的RGB值量化为12阶，生成调色板（pallet）

将调色板中的颜色按照频数从大到小排列，舍弃一部分颜色

将舍弃的颜色归到最相近的颜色中，更新调色板

计算保留下来的每种颜色的平均RGB值，并转换为Lab颜色空间

计算每种颜色之间在Lab空间上的距离，得到颜色的显著性数值

将显著性数值归一化，并进行高斯平滑

显示图像的显著性结果

五、人脸识别

clear all
clc
faceDatabase= imageSet('dataSet','recursive');
%% Face detector
Fdetect= vision.CascadeObjectDetector();

[training, test]=partition(faceDatabase,[0.5,0.5]);

featureCount=1;
for i=1:size(training,2)
for j= 1:training(i).Count
Image= read(training(i),j);
BB= step(Fdetect,Image);
if(size(BB,1)==0)
disp('problem here')
face=Image;
else
figure(1);
imshow(Image);
hold on
for ii=1:size(BB,1)
rectangle('Position',BB(ii,:),'LineWidth',5,'LineStyle','--','EdgeColor','r');
end
[m, p]= max(BB(:,3));
try
face=imcrop(Image,BB(p,:));
catch ME
disp('found one');
end
end
face= imresize(face,[224,224]);
X(featureCount,:)= extractHOGFeatures(face);
y{featureCount}= training(i).Description;
featureCount= featureCount+1;
end
end

disp('training Set image collected')

faceClassifier= fitcecoc(X,y);

pred= predict(faceClassifier, X);
pd=zeros(size(X,1),1);

for i=1:size(X,1)
if(strcmp(pred{i},y{i}))
pd(i)=1;
end
end

fprintf('\nTraining Set Accuracy: %f\n', mean(pd) * 100);

featureCount=1;
for i=1:size(test,2)
for j=1:test(i).Count
Image= read(test(i),j);
BB= step(Fdetect,Image);
if(size(BB,1)==0)
disp('problem here')
face=Image;
else
figure(1);
imshow(Image);
hold on
for ii=1:size(BB,1)
rectangle('Position',BB(ii,:),'LineWidth',5,'LineStyle','--','EdgeColor','r');
end
[m, p]= max(BB(:,3));
try
face=imcrop(Image,BB(p,:));
catch ME
disp('found one');
end
end
face= imresize(face,[224,224]);
Xtest(featureCount,:)= extractHOGFeatures(face);
ytest{featureCount}= test(i).Description;
featureCount= featureCount+1;

end
end

disp('test Set image collected')

testPred= predict(faceClassifier, Xtest);

pdt=zeros(size(Xtest,1),1);
for i=1:size(Xtest,1)
if(strcmp(testPred{i},ytest{i})==1)
pdt(i)=1;
end
end

fprintf('\nTest Set Accuracy: %f\n', mean(pdt) * 100);

六、kmeans

import

close all; %matlab代码
clear all;
clc
k=2;
org = imread('1.jpeg'); %读入图像
figure;
subplot(2,2,1);
imshow(org),title('原始图像'); %显示原图像
% 接下来需要知道图片的尺寸（长和宽），如若直接对RGB图像进行操作，如下
%转化为灰度图
gray=rgb2gray(org);
[m,n]=size(gray); %m,n为所求，p=3为通道数
% 将图像进行RGB——3通道分解 % org(:, :, 1)......分别代表rgb通道
A = reshape(org(:, :, 1), m*n, 1);
B = reshape(org(:, :, 2), m*n, 1);
C = reshape(org(:, :, 3), m*n, 1);
data = [A B C]; % r g b分量组成样本的特征，每个样本有三个属性值，共width*height个样本
% 第二张图
% kmeans第一个参数: N*P的数据矩阵，N为数据个数，P为单个数据维度
res = kmeans(double(data), k);
result = reshape(res, m, n); % 反向转化为图片形式
subplot(2,2,2);
% label2rgb功能是转换标记矩阵到RGB图像
imshow(label2rgb(result)),title(strcat('K=',num2str(k),'时RGB通道分割结果')); %显示分割结果
% 第三张图
res = kmeans(double(data), k+1);
result = reshape(res, m, n);
subplot(2,2,3);
imshow(label2rgb(result)),title(strcat('K=',num2str(k+1),'时RGB通道分割结果'));
% 第四张图
res = kmeans(double(data), k+2);
result = reshape(res, m, n);
subplot(2,2,4);
imshow(label2rgb(result)),title(strcat('K=',num2str(k+2),'时RGB通道分割结果'));

函数

function [C, label, J] = kmeans(I, k)
%C聚类中心
[m, n, p] = size(I);%图片的大小m*n，p代表RGB三层
X = reshape(double(I), m*n, p);%将I的行列排列成m*n行p列,按照列取数据的
rng('default');%将随机数生成器的设置重置为其默认值，因此会生成相同的随机数
C = X(randperm(m*n, k), :);%随机选三个聚类中心；randperm(m*n, k)-在m*n个数中，随机选择k个;X(randperm(m*n, k), :)-选取X中第k行的数据
J_prev = inf;%inf为无穷大
iter = 0;
J = [];
tol = 1e-11;%容忍度10^(-11)

while true
iter = iter + 1;
dist = sum(X.^2, 2)*ones(1, k) + (sum(C.^2, 2)*ones(1, m*n))' - 2*X*C';%计算图片中各个点到K个聚类中心的距离;sum(X, 2)-以行为单位求和;C'-转置；dist-mn行k列
[~,label] = min(dist, [], 2);%label记录最小值的列数，找到dist中那些最小值的索引位置，并放在向量label中返回。C = min(A,[],dim)-返回A中由dim指定的维数范围中的最小值，dim取2时，该函数返回一个列向量，其第i个元素是A矩阵的第i行上的最小值。label的取值在1-k之间
for i = 1:k
C(i, :) = mean(X(label == i , :)); %取新的k个聚类中心;label=1（2、3）的所有行，求每列的均值；
end
J_cur = sum(sum((X - C(label, :)).^2, 2));%距离之和;X - C(label, :)-对应行相减
J = [J, J_cur];
display(sprintf('#迭代次数: %03d, 目标函数: %f', iter, J_cur));%sprintf-将数据格式化为字符串或字符向量;%03d-左边补0 的等宽格式,比如数字12,%03d出来就是:012;sprintf() = display(sprintf());

if norm(J_cur-J_prev, 'fro') < tol%n=norm(A,p)-A和A'的积的对角线和的平方根，即sqrt(sum(diag(A'*A)))，本次与上次距离之差
break
end
if (iter==10),% A和A'的积的对角线和的平方根，即sqrt(sum(diag(A'*A)))，本次与上次距离之差
break
end
J_prev = J_cur;
end

kmeanstext

I = imread('../1.jpeg');
[m, n, p] = size(I);
figure
subplot(1,2,1);imshow(I);title('原始图像');%显示图片I

k = 2;%聚类数量
[C, label2, J2] = kmeans(I, k);
I_seg2 = reshape(C(label2, :), m, n, p);%转换成图片矩阵的格式
k = 4;
[C, label3, J3] = kmeans(I, k);
I_seg3 = reshape(C(label3, :), m, n, p);
k = 8;
[C, label4, J4] = kmeans(I, k);
I_seg4 = reshape(C(label4, :), m, n, p);
k = 16;
[C, label5, J5] = kmeans(I, k);
I_seg5 = reshape(C(label5, :), m, n, p);
figure
subplot(2, 4, 1); imshow(uint8(I_seg2), []); title('2-聚类');%imshow(I,[min(I(:)) max(I(:))])；uint8的范围是0-255。
subplot(2, 4, 2); imshow(uint8(I_seg3), []); title('4-聚类');
subplot(2, 4, 3); imshow(uint8(I_seg4), []); title('8-聚类');
subplot(2, 4, 4); imshow(uint8(I_seg5), []); title('16-聚类');