图像主题色提取算法

备份自：http://blog.rainy.im/2015/11/25/extract-color-themes-from-images/

许多从自然场景中拍摄的图像，其色彩分布上会给人一种和谐、一致的感觉；反过来，在许多界面设计应用中，我们也希望选择的颜色可以达到这样的效果，但对一般人来说却并不那么容易，这属于色彩心理学的范畴（当然不是指某些伪神棍所谓的那种）。从彩色图像中提取其中的主题颜色，不仅可以用于色彩设计（参考网站：Design Seeds），也可用于图像分类、搜索、识别等，本文分别总结并实现图像主题颜色提取的几种算法，包括颜色量化法（Color Quantization）、聚类(Clustering)和颜色建模的方法（颜色建模法仅作总结），源码可见：GitHub: ImageColorTheme。

1. 颜色量化算法

彩色图像一般采用RGB色彩模式，每个像素由RGB三个颜色分量组成。随着硬件的不断升级，彩色图像的存储由最初的8位、16位变成现在的24位、32真彩色。所谓全彩是指每个像素由8位（$2^8$=0~255）表示，红绿蓝三原色组合共有1677万（$256*256*256$）万种颜色，如果将RGB看作是三维空间中的三个坐标，可以得到下面这样一张色彩空间图：

RGB color cube

当然，一张图像不可能包含所有颜色，我们将一张彩色图像所包含的像素投射到色彩空间中，可以更直观地感受图像中颜色的分布：

Image in Color space

因此颜色量化问题可以用所有矢量量化（vector quantization, VQ）算法解决。这里采用开源图像处理库 Leptonica 中用到的两种算法：中位切分法、八叉树算法。

1.1. 中位切分法（Median cut）

GitHub: color-theif 项目采用了 Leptonica 中的用到的（调整）中位切分法，Js 代码比 C 要易读得多。中位切分算法的原理很简单直接，将图像颜色看作是色彩空间中的长方体（VBox），从初始整个图像作为一个长方体开始，将RGB中最长的一边从颜色统计的中位数一切为二，使得到的两个长方体所包含的像素数量相同，重复上述步骤，直到最终切分得到长方体的数量等于主题颜色数量为止。

Leptonica 作者在报告 Median-Cut Color Quantization 中总结了这一算法存在的一些问题，其中主要问题是有可能存在某些条件下 VBox 体积很大但只包含少量像素。解决的方法是，每次进行切分时，并不是对上一次切分得到的所有VBox进行切分，而是通过一个优先级队列进行排序，刚开始时这一队列以VBox仅以VBox所包含的像素数作为优先级考量，当切分次数变多之后，将体积*包含像素数作为优先级。

Python 3 中内置了PriorityQueue：

from queue import PriorityQueue as PQueue

class VBox(object):
  def __init__(self, r1, r2, g1, g2, b1, b2, histo):
    self.vol = calV()
    self.npixs = calN()
    self.priority = self.npixs * -1 # PQueue 是按优先级自小到大排序

boxQueue.put((vbox0.priority, vbox0))

vbox.priority *= vbox.vol
boxQueue.put((vbox0.priority, vbox0))

除此之外，算法中最重要的部分是统计色彩分布直方图。我们需要将三维空间中的任意一点对应到一维坐标中的整数，这样才能以最快地速度定位这一颜色。如果采用全部的24位信息，那么我们用于保存直方图的数组长度至少要是$2^{{24}=16777216$，既然是要提取颜色主题（或是颜色量化），我们可以将颜色由RGB各8位压缩至5位，这样数组长度只有$2}{15}=32768$：

def getColorIndex(self, r, g, b):
    return (r << (2 * self.SIGBITS)) + (g << self.SIGBITS) + b
def getPixHisto(self):
    pixHisto = np.zeros(1 << (3 * self.SIGBITS))
    for y in range(self.h):
        for x in range(self.w):
            r = self.pixData[y, x, 0] >> self.rshift
            g = self.pixData[y, x, 1] >> self.rshift
            b = self.pixData[y, x, 2] >> self.rshift

            pixHisto[self.getColorIndex(r, g, b)] += 1
    return pixHisto

分别对4张图片进行切分、提取：

def testMMCQ(pixDatas, maxColor):
    start  = time.process_time()
    themes = list(map(lambda d: MMCQ(d, maxColor).quantize(), pixDatas))
    print("MMCQ Time cost: {0}".format(time.process_time() - start))
    return themes
imgs = map(lambda i: 'imgs/photo%s.jpg' % i, range(1,5))
pixDatas = list(map(getPixData, imgs))
maxColor = 7

themes = [testMMCQ(pixDatas, maxColor)]
imgPalette(pixDatas, themes, ["MMCQ Palette"])

mmcq

1.2. 八叉树算法（Octree）

八叉树算法的原理可以参考这篇文章：圖片主題色提取算法小結。作者也提供了 Js 实现的代码，虽然与 Leptonica 中 C 实现的方法差别很大，但原理上是一致的。

建立八叉树的原理实际上跟上面提到的统计直方图有些相似，将颜色成分转换成二进制之后，较低位（八叉树中位置较深层）数值将被压缩进较高位（八叉树中较浅层）。八叉树算法应用到主题色提取可能存在的问题是，每次削减掉的叶子数不确定，但是新增加的只有一个，这就导致我们需要的主题色数量并不一定刚好得到满足，例如设定的主题色数量为7，可能上一次叶子时总数还有10个，到了下一次只剩5个了。类似的问题在后面手动实现的KMeans算法中也有出现，为了保证可以得到足够的主题色，不得不强行提高算法中的颜色数量，然后取图像中包含数量较多的作为主题色：

def getColors(self, node):
      if node.isLeaf:
          [r, g, b] = list(map(lambda n: int(n[0] / n[1]), zip([node.r, node.g, node.b], [node.n]*3)))
          self.theme.append([r,g,b, node.n])
      else:
          for i in range(8):
              if node.children[i] is not None:
                  self.getColors(node.children[i])
self.theme = sorted(self.theme, key=lambda c: -1*c[1])
return list(map(lambda l: l[:-1],self.theme[:self.maxColor]))

对比上面两种算法的结果：

def testOQ(pixDatas, maxColor):
    start  = time.process_time()
    themes = list(map(lambda d: OQ(d, maxColor).quantize(), pixDatas))
    print("OQ Time cost: {0}".format(time.process_time() - start))
    return themes
themes = [testMMCQ(pixDatas, maxColor), testOQ(pixDatas, maxColor)]
imgPalette(pixDatas, themes, ["MMCQ Palette", "OQ Palette"])

MMCQ vs OQ

可见八叉树算法可能更适合用于提取调色板，而且两种算法运行时间差异也很明显：

#MMCQ Time cost: 8.238793

#OQ Time cost: 55.173573

除了OQ中采用较多递归以外，未对原图进行抽样处理也是其中原因之一。

2. 聚类

聚类是一种无监督式机器学习算法，我们这里采用K均值算法。虽然说是“机器学习”听起来时髦些，但算法本质上比上面两种更加简单粗暴。

KMeans算法

KMeans算法的原理更加简洁：“物以类聚”。我们目的是将一堆零散的数据（如上面图2）归为k个类别，使得每个类别中的每个数据样本，距离该类别的中心（质心，centroid）距离最小，数学公式为：

$$
\sum_{i=0}^N \min_{ \mu_j \in C} (||x_i - \mu_j||^2)
$$

上文提到八叉树算法可能出现结果与主题色数量不一致的情况，在KMeans算法中，初始的k个类别的质心的选择也可能导致类似的问题。当采用随机选择的方法时，有可能出现在迭代过程中，选择的中心点距离所有其它数据太远而最终导致被孤立。这里分别采用手动实现和scikit-learn的方法实现，根据scikit-learn 提供的API，完成主题色的提取大概只需要几行代码：

from sklearn.cluster import KMeans as KM
import numpy as np

#@pixData      image pixels stored in numpy.ndarray
#@maxColor     theme color number
h, w, d = pixData.shape
data = np.reshape((h*w, d))
km = KM(n_clusters=maxColor)
km.fit(data)
theme = np.array(km.cluster_centers_, dtype=np.uint8)

imgs = map(lambda i: 'imgs/photo%s.jpg' % i, range(1,5))
pixDatas = list(map(getPixData, imgs))
maxColor = 7
themes = [testKmeans(pixDatas, maxColor), testKmeans(pixDatas, maxColor, useSklearn=False)]
imgPalette(pixDatas, themes, ["KMeans Palette", "KMeans DIY"])

测试比较手动实现和scikit-learn的结果如下：

KMeans vs DIY

好吧我承认很惨，耗时方面也是惨不忍睹。

3. 色彩建模

从上面几种算法结果来看，MMCQ和 KMeans在时间和结果上都还算不错，但仍有改进的空间。如果从人类的角度出发，两种算法的策略或者说在解决主题色提取这一问题时采纳的特征（feature）都接近于颜色密度，即相近的颜色凑在一起数量越多，越容易被提取为主题颜色。

最后要提到的算法来自斯坦福可视化组13年的一篇研究：Modeling how people extract color themes from images，实际上比较像一篇心理学研究的套路：建模-找人类被试进行行为实验-调参拟合。文章提取了图像中的79个特征变量并进行多元回归，同时找到普通人类被试和艺术系学生对图像的主题颜色进行选择，结果证明特征+回归能够更好地拟合人类选择的结果。

all color themes

79个特征的多元回归模型，不知道会不会出现过度拟合？另外虽然比前面算法多了很多特征，但仍旧多物理特征。对人类观察者来说，我们看到的并非一堆无意义的色块，虽然有研究表明颜色信息并非场景识别的必要线索，但反过来场景图像中的语义信息却很有可能影响颜色对观察者的意义，这大概就是心理学研究与计算机科学方向上的差异。

总结

以上算法若要应用还需更多优化，例如先抽样再处理，计算密集的地方用C/C++或并行等。另外需要一个对Python每个函数执行时间进行记录的工具，分析运行时间长的部分。

参考

1. Color Quantization
2. Color quantization using modified median cut
3. Median-Cut Color Quantization
4. Wicked Code
5. Clustering - scikit-learn
6. Color Quantization using K-Means
7. Extract Color Themes from Images
8. Lin, S., & Hanrahan, P. (2013). Modeling how people extract color themes from images. Proc of Chi Acm, 3101-3110.