计算机视觉

目录

一、图像处理

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-聚类');

七、神经网络

import torch

import numpy as np

from torch import nn, optim

from torch.autograd import Variable

from torch.utils.data import DataLoader

from torchvision import datasets, transforms

import matplotlib.pyplot as plt

class CNN(nn.Module):

   def __init__(self):

       super(CNN, self).__init__()  

       self.layer1 = nn.Sequential(nn.Conv2d(1, 10, kernel_size=5), nn.BatchNorm2d(10), 

           nn.ReLU(), nn.MaxPool2d(kernel_size=2))

       self.layer2 = nn.Sequential(nn.Conv2d(10, 20, kernel_size=5),nn.BatchNorm2d(20),  

           nn.ReLU(),  nn.MaxPool2d(kernel_size=2))

       self.fc = nn.Sequential( nn.Linear(320, 50),  nn.ReLU(),nn.Linear(50, 10))

   def forward(self, x):  

       x = self.layer1(x)  

       x = self.layer2(x) 

       x = x.view(x.size(0), -1)  

       x = self.fc(x)  

       return x

class BP(nn.Module):

   def __init__(self):

       super(BP, self).__init__()

       self.fc = nn.Sequential(nn.Linear(784, 128),nn.ReLU(inplace=True), nn.Linear(128, 64), nn.ReLU(inplace=True),nn.Linear(64, 10))

   def forward(self, x):

       x = x.view(x.size(0), -1)

       x = self.fc(x)

       return x

def train(model, train_loader, optimizer, criterion, device):

   model.train()

   for batch_idx, (data, target) in enumerate(train_loader):

       data, target = data.to(device), target.to(device)

       optimizer.zero_grad()

       output = model(data)

       loss = criterion(output, target)

       loss.backward()

       optimizer.step()

       if batch_idx % 100 == 0:

           print('Train Epoch: {} [{}/{} ({:.0f}%)]\tLoss: {:.6f}'.format(

               epoch, batch_idx * len(data), len(train_loader.dataset),

               100. * batch_idx / len(train_loader), loss.item()))

   return loss.item()

def test(model, test_loader, test_dataset, criterion, device):model.eval()

   test_loss = 0

   correct = 0

   with torch.no_grad():

       for data, target in test_loader:

           data, target = data.to(device), target.to(device)

           output = model(data)

           test_loss += criterion(output, target).item()