对于C-Means，实质也是一个对样本的分类，最初我以为是对特征的分类，让我吃了很多苦！！

一、C-Means的原理：

二、matlab实现代码

%% C_menans
clear ;clc;
male_data = importdata('male.txt');
male_av = mean(male_data);
female_data = importdata('female.txt');
female_av = mean(female_data);
x = [male_data;female_data];
%plot(x(:,1),x(:,2),'.');
k = 2;
%取初值
m = zeros(k,2);
m = [160,50;176,68];  %c==2
% m = [male_av;female_av];
m_ref = zeros(k,2);
data_ad_la = zeros(length(x),1);
value =zeros(k,1);%记录每个类别样本的值，第一列是值，第二列是类别
w_num = zeros(k,1);%每一类的个数
x_temp =zeros(k,1);
hei = zeros(k,1);
wei = zeros(k,1);

%还要改这里
s1 = [];s2 = [];


while (  norm(m-m_ref)>0.0001  )
    m_ref = m;
    for i=1:length(x)
        for j=1:k
            x_temp(j) = norm(x(i,:) -m(j,:) ); %1*2-1*2
        end
        [A,B] = min(x_temp);
        data_ad_la(i) = B;
        hei(B) = hei(B) +x(i,1);
        wei(B) = wei(B) +x(i,2);
        value(B) = value(B)+A; %每一个类别都有个距离的和
        %类别个数不同，这里要改改
        if(B==1)
            s1 = [s1,i];
        else
            s2 = [s2,i];
        end
    end
    for i = 1:k       
        w_num(i) = sum(data_ad_la==i);
        m(i,:) = [hei(i)/w_num(i),wei(i)/w_num(i)];
    end
    Je = sum( value./w_num );
    value = zeros(k,1);
    hei =zeros(k,1);
    wei = zeros(k,1);
%     s1 = [];s2 = [];
end
figure;
% plot(x(s1,1),x(s1,2),'*',x(s2,1),x(s2,2),'o',x(s3,1),x(s3,2),'+');legend('类1','类2','类3');title('C=3时的分类');
plot(x(s1,1),x(s1,2),'*',x(s2,1),x(s2,2),'o');legend('类1','类2');
% title('C=2时的分类,初值：m = [160,53;178,70];');
fprintf('   Je is %f',Je);

测试结果：

        改变类别数，分别进行两类、三类、四类、五类聚类：

注：这里的这个图像我是通过Je/（聚类个数）得到的，按理说随着聚类个数的增加，Je应该会变小，可能代码有问题，还请大家指正

：

结论：随着聚类个数的增加，其Je会减小，但不会无限制的减小。聚类的结果还收初始聚类中心的影响。

1.绪

        图像复原目的：获得关于原始图像的近似估计f(x,y)。

        由于种种原因(拍摄、传输、保存等)，图像变得模糊，或者被噪声污染，图像退化，退化的模型可以被模型化为一个退化函数和一个加性噪声项，这个时候想得到一个比较清晰完整的图像，就需要图像复原。

        主要处理以噪声引起的退化。

2.噪声

2.1 一些重要噪声的概率密度函数

===================================================================== 

========================================================================= 

 =========================================================================

2.2 添加噪声(实例1) 

                   图1.1

对图1.1分别加入上面各种噪声，并显示出其直方图。

=============================MATLAB代码============================

注意将不是灰度图的转为灰度图，如果已经是灰度图就要将rgb2gray函数去掉，后面的代码同理。

如果想显示出比较高清的图片，可以使用imwrite(fig,'a.jpg');fig为要保存的图片，a.jpg为保存的名字和格式，可以直接保存。

直方图可以点击编辑窗口，进行参数范围，字体大小等的调整，使得呈现的图片更加饱满。

噪声添加，也可以使用imnoise(f,type,parameters);其中，f是输入图像，type选择噪声类型，parameters为可调整的参数(噪声的均值和方差等)。

clc;
clear;
%读入图像，并转换为double型
B=imread('huade.jpg');
I = rgb2gray(B);        %转化为灰度图
I_D=im2double(I);
%获得图像大小
[M,N]=size(I_D);
 
%========================生成高斯白噪声==================================
a=0;
b=0.001;
N_Gau=a+sqrt(b)*randn(M,N); 
%将噪声叠加到图像上
J_Gaussian=I_D+N_Gau;
%保存噪声数据
 
%======================生成瑞利噪声======================================
a=0;
b=0.08;
B=1;
N_Ray1=a+b*raylrnd(B,M,N);
%将噪声叠加到图像上
J_rayl=I_D+N_Ray1;
 
%=====================叠加伽马噪声=======================================
a=0;
b=0.08;
A=1;
B=2;
N_Gam=a+b*gamrnd(A,B,[M,N]);
%将噪声叠加到图像上
J_Gamma=I_D+N_Gam;
 
%====================叠加指数噪声=======================================
a=0;
b=0.08;
mu=2;
N_exp=a+b*exprnd(mu,[M,N]);
J_exp=I_D+N_exp;
 
%=====================叠加均匀分布噪声=======================================
a=0;
b=0.08;
A=0;
B=2;
N_unif=a+b*unifrnd(A,B,[M,N]);
J_unif=I_D+N_unif;
 
%=====================叠加椒盐分布噪声=======================================
a=0.25;
J_salt = imnoise(I_D,'salt & pepper',a);
 
 
subplot(2,3,1);
imshow(J_Gaussian);
title('高斯噪声');
subplot(2,3,2);
imshow(J_rayl);
title('瑞利噪声');
subplot(2,3,3);
imshow(J_Gamma);
title('伽马噪声');
subplot(2,3,4);
imshow(J_exp);
title('指数噪声');
subplot(2,3,5);
imshow(J_unif);
title('均匀分布噪声');
subplot(2,3,6);
imshow(J_salt);
title('椒盐噪声');
 
figure;
subplot(2,3,1);
imhist(J_Gaussian);
title('高斯噪声');
subplot(2,3,2);
imhist(J_rayl);
title('瑞利噪声');
subplot(2,3,3);
imhist(J_Gamma);
title('伽马噪声');
subplot(2,3,4);
imhist(J_exp);
title('指数噪声');
subplot(2,3,5);
imhist(J_unif);
title('均匀分布噪声');
subplot(2,3,6);
imhist(J_salt);
title('椒盐噪声');

 运行结果：

2.3 周期噪声

        周期噪声是在图像获取时从电力或机电干扰中产生的。这是唯一的一种空间依赖型噪声。通常有正弦噪声和余弦噪声等。周期噪声通过频域滤波可以显著地减少。

2.4 添加正弦噪声(实例2)

                  图2

向图2加入正弦噪声，并画出其频谱图。

=============================MATLAB代码============================

clear;
clc;
%=======加入正弦噪声========
t=imread('einstein.jpg');
[m,n]=size(t);
t_1=t;
for i=1:m
for j=1:n
t_1(i,j)=t(i,j)+40*sin(40*i)+40*sin(40*j);
end
end
imshow(t),title('原图');
imwrite(t_1,'einstein1.jpg');
%===========谱=============
data = imread('einstein1.jpg');
data = rgb2gray(data);
F = fft2(data);
FC = fftshift(F);
S2 = log(1+abs(FC));
imshow(S2,[]);

 运行结果：

3.均值滤波

3.1 算术均值滤波

最简单的均值滤波器，m*n的矩形子图像窗口的各个值相加的和除以m*n。均值简单地平滑了一幅图像的局部变化，在模糊图像的同时减少了噪声。

3.2 几何均值滤波

每个被复原像素由子图像窗口中像素点的乘积并自乘到1/mn次幂。与算术均值相比，在滤波过程中会更少丢失图像细节。

3.3 均值滤波1(实例3)

                   图3

对图3加入均值为0，方差为0.02(归一化后)的高斯噪声，再分别用同样大小的窗口进行算术均值滤波和几何均值滤波。

=============================MATLAB代码============================

%高斯噪声部分
clc;clear;
I=imread('test1.jpg');
I=im2double(I);
 
%添加高斯噪声
I_noise=double(imnoise(I,'gaussian',0,0.02));%均值为0
imwrite(I_noise,'a.jpg');
%算术均值滤波
%定义子窗口的尺寸
m=3;
n=3;
%确定要扩展的行列数
len_m=floor(m/2);
len_n=floor(n/2);
%将原始图像进行扩展，这里采用了镜像扩展，以进行图像边缘计算
I_2=im2double(I_noise);
[M2,N2]=size(I_2);
I_2_pad=padarray(I_2,[len_m,len_n],'symmetric');
I_3=fspecial('average',[3,3]);%3*3均值滤波    建立预定义的滤波算子
I_3=imfilter(I_2_pad,I_3);%(待处理矩阵，滤波器)
imwrite(I_3,'b.jpg');

%几何均值滤波 
I_D=im2double(I_noise);
[MM,NN]=size(I_D);
%定义子窗口的尺寸
m=3;
n=3;
%确定要扩展的行列数
len_m=floor(m/2);
len_n=floor(n/2);
%将原始图像进行扩展，这里采用了镜像扩展，以进行图像边缘计算
I_D_pad=padarray(I_D,[len_m,len_n],'symmetric');
%获得扩展后的图像尺寸
[M,N]=size(I_D_pad);
J_Geometric=zeros(MM,NN);
%逐点计算子窗口的几何平均
for i=1+len_m:M-len_m
    for j=1+len_n:N-len_n
        %从扩展图像中取出子图像
        Block=I_D_pad(i-len_m:i+len_m,j-len_n:j+len_n);
        %求子窗口的几何均值
        J_Geometric(i-len_m,j-len_n)=(prod(prod(Block)))^(1/(m*n));
    end
end
imwrite(J_Geometric,'c.jpg');

 运行结果：

 (a)原图        (b)高斯噪声污染    

(c)3*3均值滤波     (d)3*3几何均值滤波

3.4 谐波均值滤波

谐波均值滤波器对“盐”噪声效果更好，但是不适用于“胡椒”噪声，它善于处理像高斯噪声那样的其他噪声。

3.5 逆谐波均值滤波

适合减少或在实际中消除椒盐噪声的影响。

Q称为滤波器的阶数。

当Q值为正数时，滤波器用于消除“胡椒”噪声。

当Q值为负数时，滤波器用于消除“盐粒”噪声。

当Q=0时，为算术均值滤波器。

当Q为-1时，为谐波均值滤波器。

3.6 均值滤波2(实例4)

对lena图分别加入概率为0.2的胡椒噪声和盐粒噪声，再用不同Q值得逆谐波滤波器分别对图像进行滤波。

=============================MATLAB代码============================

%椒盐噪声部分
clc;clear;
I=imread('lena.bmp');
%I=rgb2gray(I);
imshow(I);
I=im2double(I);
%添加胡椒信号
J = imnoise(I,'salt & pepper',0.1);
[M,N]=size(J);
for i=1:M
    for j=1:N
        if(J(i,j)==1)
            J(i,j)=0;
        else
            J(i,j)=J(i,j);
        end
    end
end

figure;imshow(J);
%Q=1.5逆谐波滤波器滤波对胡椒
Q=1.5;
I_mean=imfilter(J.^(Q+1),fspecial('average',3))./imfilter(J.^Q,fspecial('average',3));
figure;imshow(I_mean);
% Salt Noise
H = imnoise(I,'salt & pepper',0.1);
[M,N]=size(H);
for i=1:M
    for j=1:N
        if(H(i,j)==0)
            H(i,j)=1;
        else
            H(i,j)=H(i,j);
        end
    end
end

figure;imshow(H);
%Q=-1.5逆谐波滤波器滤波对盐粒
Q=-1.5;
I_mean=imfilter(H.^(Q+1),fspecial('average',3))./imfilter(H.^Q,fspecial('average',3));
figure;imshow(I_mean);

 运行结果：

如果Q值选取错误，即用正值的Q对“盐粒”进行滤波，用负值的Q对“胡椒”进行滤波，会导致更不好的效果。

=============================MATLAB代码============================

%椒盐噪声部分
clc;clear;
I=imread('lena.bmp');
%I=rgb2gray(I);
imshow(I);
I=im2double(I);
%添加胡椒信号
J = imnoise(I,'salt & pepper',0.1);
[M,N]=size(J);
for i=1:M
    for j=1:N
        if(J(i,j)==1)
            J(i,j)=0;
        else
            J(i,j)=J(i,j);
        end
    end
end
figure;imshow(J);
%Q=1.5逆谐波滤波器滤波对胡椒
Q=1.5;
I_mean=imfilter(J.^(Q+1),fspecial('average',3))./imfilter(J.^Q,fspecial('average',3));
figure;imshow(I_mean);
% Salt Noise
H = imnoise(I,'salt & pepper',0.1);
[M,N]=size(H);
for i=1:M
    for j=1:N
        if(H(i,j)==0)
            H(i,j)=1;
        else
            H(i,j)=H(i,j);
        end
    end
end
figure;imshow(H);
%Q=-1.5逆谐波滤波器滤波对盐粒
Q=-1.5;
I_mean=imfilter(H.^(Q+1),fspecial('average',3))./imfilter(H.^Q,fspecial('average',3));
figure;imshow(I_mean);

 运行结果：

 (a)Q=-1.5对“胡椒”进行滤波                            (b)Q=1.5对“盐粒”进行滤波

4.统计排序滤波

4.1 中值滤波

最著名的统计排序滤波器就是中值滤波器，就是用该像素的相邻像素的灰度中值来代替该像素的值。相同尺寸下，比线性平滑滤波器引起的模糊更少。对单极或双极脉冲噪声非常有效。

4.2 多次中值滤波(实例5)

对lena图加入0.2的椒盐噪声，再进行中值滤波。

=============================MATLAB代码============================

%椒盐噪声部分
clc;clear;
I=imread('test1.jpg');
I=im2double(I);
I_noise=double(imnoise(I,'salt & pepper',0.2));%salt & pepper注意中间的空格 无空格报错
imshow(I_noise);title('椒盐噪声');
%中值滤波
f1 = medfilt2(I_noise,'symmetric');%symmetric指出图像按照镜像反射方式对称地沿边界扩展
figure;imshow(f1);
%再次滤波
f2 = medfilt2(f1,'symmetric');%symmetric指出图像按照镜像反射方式对称地沿边界扩展
figure;imshow(f2);
%再一次滤波
f3 = medfilt2(f2,'symmetric');%symmetric指出图像按照镜像反射方式对称地沿边界扩展
figure;imshow(f3);

 运行结果：

   (a)被概率为0.2的椒盐噪声污染                     (b)使用3*3的中值滤波对(a)处理

   (c)使用3*3的中值滤波对(b)处理                    (d)使用3*3的中值滤波对(c)处理

4.3 最大值和最小值滤波

 最大值：用该像素的相邻像素的灰度最大值来代替该像素的值，滤除“胡椒”比较好。

 最小值：用该像素的相邻像素的灰度最小值来代替该像素的值，滤除“盐粒”比较好。

4.4 最大最小滤波(实例6)

向lena图分别加入0.5的胡椒噪声和盐粒噪声，然后进行最大值最小值滤波。

=============================MATLAB代码============================

clc;clear;
I=imread('hong.jpg');
I=rgb2gray(I);
I=im2double(I);
%添加胡椒信号
J = imnoise(I,'salt & pepper',0.25);
[M,N]=size(J);
for i=1:M
    for j=1:N
        if(J(i,j)==1)
            J(i,j)=0;
        else
            J(i,j)=J(i,j);
        end
    end
end
% Salt Noise
H = imnoise(I,'salt & pepper',0.25);
[M,N]=size(H);
for i=1:M
    for j=1:N
        if(H(i,j)==0)
            H(i,j)=1;
        else
            H(i,j)=H(i,j);
        end
    end
end
%使用最大值
f3 = ordfilt2(J,9,ones(3,3));
imshow(f3);
%使用最小值
f4 = ordfilt2(H,1,ones(3,3));
figure;imshow(f4);

 运行结果：

4.5 修正后的阿尔法均值滤波器