利用几种颜色量化方法提取图片颜色色调

这里主要用3种方法：包含两种k-means方法和一种matlab自带的最小方差量化方法。
为了方便对比，文中选用梵高的星空（The Starry Night）作为提取图片。

1 k-means方法（在RGB空间）

颜色量化是指一个减少在图像中的颜色数量的压缩过程。旨在用尽可能少的颜色去还原图片。
所以，也可以用聚类方法去进行处理，即将所有颜色在色彩空间上分布开，然后为每个集群选择所代表的点来创建调色板。

比如梵高的星空：

把所有像素点的颜色在RGB空间显示出来，可以得到：

这时我们再利用k-means聚类算法把调色板设为4种颜色，即N=4，
得到像素的局聚类表示为

质心颜色的连线图如下：

得到的N=4下的图片为

可以看到有较大的偏差。
当取N=20时，得到的图为：

颜色质量要比只用4个颜色的要高出很多。

而在颜色板的分布这张图上，也可以比较明显的看到梵高在做这幅画时使用的主要3种颜色，和wiki百科中提到的ultramarine 、cobalt blue 和indian yellow比较吻合。

2 k-means方法（Lab颜色空间）

事实上，在相同的方法，采用不同的颜色空间进行聚类，可以发现影响最主要的还是数量的选取。可能有些情况下Lab或者HSV方法比RGB看上去更好，但是影响也不是很大。

下图是利用LAB颜色空间进行聚类分析的结果，N=4，并且为了方便对比，将结果依然显示在RGB的坐标系内。

可以看到结果和第一节的相似。

当N=20的时候，观察和第一节的区别，发现只是在一些细节上划分区域有所差别，基本效果很接近。

3 最小方差量化方法rgb2ind()

matlab中有专门的函数用来进行色彩量化，为rgb2ind()，采用的方法为最小方差量化方法。具体可以参见：
https://ww2.mathworks.cn/help/matlab/ref/rgb2ind.html
这里依然取N=4和N=20进行对比。

可以看到用matlab自带的方法要比kmeans方法更好。

4 其余名画

5 代码

K-Means算法用的是matlab里自带的kmeans函数，可以参见：
http://www.mathworks.com/help/stats/kmeans.html
颜色板分布画的连接图用的算法来源于 Python 平面上的点的最短连线网这篇文章：
https://blog.csdn.net/lnotime/article/details/82313355
LAB和RGB之间的转换用的是色彩空间Lab和RGB互转的Matlab代码这篇文章：
https://blog.csdn.net/wendy260310/article/details/17100497

这里由于篇幅，我就只粘贴Kmeans-LAB算法的程序和rgb2ind程序。

clear
clcN=20;%设置聚类数目
data = imread('xingkong.JPEG');
dataLAB=single(zeros(size(data)));
[dataLAB(:,:,1),dataLAB(:,:,2),dataLAB(:,:,3)] = RGB2Lab(data(:,:,1),data(:,:,2),data(:,:,3));
hs=size(data);
data2=zeros(hs(1)*hs(2),3);
ii=0;
for j=1:hs(1)for k=1:hs(2)ii=ii+1;data2(ii,:)=dataLAB(j,k,:);end
endXLAB=single(data2);[idx,CLAB] = kmeans(XLAB,N,'Distance','cityblock','Replicates',5);
C=zeros(size(CLAB));
X=zeros(size(XLAB));
[C(:,1),C(:,2),C(:,3)]=Lab2RGB(CLAB(:,1),CLAB(:,2),CLAB(:,3));
[X(:,1),X(:,2),X(:,3)]=Lab2RGB(XLAB(:,1),XLAB(:,2),XLAB(:,3));
X=X*256;
C=C*256;figure%原图
set(gcf,'unit','centimeters','position',[10,10,8,7])
image(data);
%三维绘图
figure
set(gcf,'unit','centimeters','position',[10,10,8,7])
hold on
%'MarkerFaceColor'
for j=1:N%plot(X(idx==j,1),X(idx==j,2),'LineStyle','none','Marker','.','MarkerSize',12)plot3(X(idx==j,1),X(idx==j,2),X(idx==j,3),'LineStyle','none','Marker','.','MarkerSize',12,'MarkerFaceColor',C(j,:)/256,'MarkerEdgeColor',C(j,:)/256);
end
%plot3(C(:,1),C(:,2),C(:,3),'kx','MarkerSize',15,'LineWidth',3)
view(3)
hold off%然后输出这几个基准色对应的图
data2_2=single(data2);j=1:length(idx);
data2_2(j,:)=C(idx(j),:);ii=0;
for j=1:hs(1)for k=1:hs(2)ii=ii+1;data(j,k,:)=data2_2(ii,:);end
end
figure%新图
set(gcf,'unit','centimeters','position',[10,10,8,7])
image(data);figure%点线图
set(gcf,'unit','centimeters','position',[10,10,8,7])
paths=find_min_length_3D(C);%参见：平面上的点的最短连线网
hold on
for j=1:size(paths,1)plot3(paths(j,[1,4]),paths(j,[2,5]),paths(j,[3,6]),'b-')plot3(C(j,1),C(j,2),C(j,3),'LineStyle','none','Marker','.','MarkerSize',40,'MarkerFaceColor',C(j,:)/256,'MarkerEdgeColor',C(j,:)/256)
end
plot3(C(j+1,1),C(j+1,2),C(j+1,3),'LineStyle','none','Marker','.','MarkerSize',40,'MarkerFaceColor',C(j,:)/256,'MarkerEdgeColor',C(j+1,:)/256)
hold off
view(3)

rgb2ind程序

clear
N=20;
%设置图像位数
RGB = imread('Stacks of wheat.jpg');
figure
set(gcf,'unit','centimeters','position',[10,10,8,7])
imagesc(RGB)[IND,map] = rgb2ind(RGB,N);
IND=IND+1;
figure
set(gcf,'unit','centimeters','position',[10,10,8,7])
imagesc(IND)
colormap(map)%三维绘图
hs=size(RGB);
data1=reshape(RGB,hs(1)*hs(2),3);
IND1=reshape(IND,hs(1)*hs(2),1);figure
set(gcf,'unit','centimeters','position',[10,10,8,7])
hold on
for j=1:Nplot3(data1(IND1==j,1),data1(IND1==j,2),data1(IND1==j,3),'LineStyle','none','Marker','.','MarkerSize',12,'MarkerEdgeColor',map(j,:));
end
view(3)
hold offfigure%绘制连线图
set(gcf,'unit','centimeters','position',[10,10,8,7])
paths=find_min_length_3D(map);
hold on
for j=1:size(paths,1)plot3(paths(j,[1,4]),paths(j,[2,5]),paths(j,[3,6]),'b-')plot3(map(j,1),map(j,2),map(j,3),'LineStyle','none','Marker','.','MarkerSize',40,'MarkerEdgeColor',map(j,:))
end
plot3(map(j+1,1),map(j+1,2),map(j+1,3),'LineStyle','none','Marker','.','MarkerSize',40,'MarkerEdgeColor',map(j+1,:))hold off
view(3)

find_min_length_3D函数

function paths=find_min_length_3D(pots)
%输入单，返回最短连接网格
%pots满足N行3列的格式N=size(pots,1);
%初始化
con=pots(1,:);%已经连线的点集，先随便放一个点进去
Ncon=pots(2:end,:);%还没连线的点集
paths=zeros(N-1,6);%所有连线length_total=0;%总长度
%循环
for j=1:N-1potA=con(1,:);potB=Ncon(1,:);lengthAB=sum((potA-potB).^2);for m=1:size(con,1)for n=1:size(Ncon,1)lengthAB_temp=sum((con(m,:)-Ncon(n,:)).^2);if lengthAB_temp<lengthABlengthAB=lengthAB_temp;%如果新线段更短，则采用新的线段potA=con(m,:);potB=Ncon(n,:);endendendpaths(j,:)=[potA,potB];con=[con;potB];Ncon(sum(Ncon==potB,2)==3,:)=[];%这里3D所以是3length_total=length_total+lengthAB;
end
end