如何科学地利用高光谱图像合成真彩色RGB图像?

1. 前言

参考链接:

色匹配函数是什么？ - 知乎 (zhihu.com)

23. 颜色知识1-人类的视觉系统与颜色 - 知乎 (zhihu.com)

色彩空间基础 - 知乎 (zhihu.com)

色彩空间表示与转换 - 知乎 (zhihu.com)

CIE XYZ - freshair_cn - 博客园 (cnblogs.com)

CIE 1931 color space - Wikipedia

CIE 1931色彩空间 - 维基百科，自由的百科全书 (wikipedia.org)

RGB/XYZ Matrices (brucelindbloom.com)

调色名词浅析——Gamma（伽玛）校正 - 知乎 (zhihu.com)

前置条件：

• 必须有高光谱图像
• 必须有波段信息，(你需要知道高光谱图像每个通道对应的波长是多少)

对于波段信息，我的波段信息是从软件ENVI中获取的txt文件
本文以128波段的高光谱图像为例，txt文件详细内容在文末给出

先上流程图，随后再逐一解释：

2. 计算色匹配函数

色匹配函数其实就是如下三条曲线(使用python手动绘制):

CIE-XYZ色匹配函数曲线如下:

2.1 什么是色匹配函数

追根溯源的话，得从我们人类的视网膜说起。大部分人类的视网膜上有三种感知颜色的感光细胞，叫做视锥细胞，分别对不同波长的光线敏感，称为 L/M/S 型细胞。三种视锥细胞最敏感的波长分别是橙红色（长波，Long），绿色（中波，Medium），蓝色（短波，Short）。这三种视锥细胞的归一化感光曲线如下图所示

不仅如此，人类眼睛对不同颜色光线混合的反应还是 线性 的。根据 格拉斯曼定律（Grassmann’s Law），两束不同颜色的光 C 1 C_1 C1​和 C 2 C_2 C2​，假设某个视锥细胞对他们的反应分别是 r 1 r_1 r1​ 和 r 2 r_2 r2​，现在将他们按照一个比例混合，得到第三种颜色 C 3 = α C 1 + β C 2 C_3=\alpha C_1 + \beta C_2 C3​=αC1​+βC2​，那么视锥细胞对这个混合颜色的反应也将是前两个反应的线性叠加 r 3 = α r 1 + β r 2 r_3=\alpha r_1 + \beta r_2 r3​=αr1​+βr2​。

格拉斯曼定律是一个实验规律，并没有物理或者生物学上的依据。然而这个规律大大简化了我们对人类彩色视觉系统的建模，并且给我们使用线性代数理论分析人类彩色视觉系统提供了一个前提和基础。

前面已经提到，人类视网膜上有三种感知色彩的视锥细胞，所以理论上我们用三种颜色的光就可以混合出自然界中任何一种颜色来。在 20 世纪 20 年代，David Wright 和 John Guild 各自独立地领导了一些实验，通过三种颜色的光源进行匹配，得到了人眼对于不同颜色光的匹配函数。此后，多名科学家多次进行了类似的实验，加深了我们对人类彩色视觉的认识。

实验过程大致是这样的，把一个屏幕用不透光的挡板分割成两个区域，左边照射某个被测试的颜色的光线，这里记为 C C C（以下用大写字母表明颜色，用小写字母表明分量大小），右边同时用三种颜色的光同时照射，这里记为 R R R， G G G， B B B。然后，调节右边三种颜色光源的强度，直到左右两边的颜色看上去一样为止。假设这个时候三种颜色的光源强度分别为 r r r， g g g， b b b，那么根据光色叠加的线性性质，我们可以写出

C = r R + g G + b B C = rR + gG + bB C=rR+gG+bB

也就是说，只要按照 (r,g,b) 的分量来混合 (R,G,B) 三种颜色的光，就可以得到 C 这个颜色的光。于是在这一系列实验里，科学家们把左边的颜色按着光谱顺序，挨个测试了一遍，得到了纯光谱色的混合叠加的数据，这就是 色匹配函数（Color Matching Function） ，并且在这个基准下定义的色彩空间，就是 CIE RGB 色彩空间。下图是 CIE RGB 的色匹配函数曲线，（图片数据来自CVRL main (ucl.ac.uk)）。浅色的细线代表实验中不同参与者个人的色匹配函数曲线，中间深色的粗线代表数据的平均值。

可以看到，曲线上出现了负数，这是怎么回事？回想一下前面描述的实验过程，左边是被测试的光色，右边是可调节的三色光的混合。如果碰到一种情况，右边三色光无论如何调节比例，都不能混合出左边的颜色，比如某种颜色的光强度已经减小为 0 了，然而看趋势还需要继续减小才能与左边的光色相匹配，怎么办？这时候需要往左边的光色中混入三色光中的一种或者几种，继续调节，直到两边的颜色匹配。在左边（被测试）的色光中添加，那就是相当于在右边的混合光中减去，这就导致了色匹配函数曲线上出现了负数。实际上，这相当于就是光线的「减法」了。

比如，对于 555nm 的黄色光，色匹配函数的值是 (1.30, 0.97, -0.01)，意味着将 1.30 份的红光与 0.97 份的绿光混合放在右边，左边放上 1 份的 555nm 的黄光，以及 0.01 份的蓝光，这样左右两边的光色看上去就一样了。

因为有部分出现了负数，在使用和计算上都有不方便，因此就对这个匹配函数进行了一下线性变换，变换到一个所有分量都是正的空间中。变换后的色彩空间就是 CIE XYZ 色彩空间。 （图片数据来自CVRL main (ucl.ac.uk)）

CIE RGB 色彩空间和 CIE XYZ 色彩空间是完全等价的，两者只是差了一个线性变换。由于允许「减法」的存在，因此 CIE RGB 空间是能够表示所有颜色的；同样的，CIE XYZ 空间也能。

在计算机问世之前，计算带有负值乘法的曲线是很麻烦的；XYZ就是为了把RGB空间转换成另一个没有负值的、方便计算的空间。很明显，XYZ的三色刺激值是想象出来的。

XYZ不用于描述颜色，而用于说明光波如何组合会产生什么样的颜色，因此XYZ是独立于设备的。

对于某些计算任务来说，查表可能变得不切实际。CIE XYZ颜色匹配函数可以通过高斯函数的和来近似。

设 g ( x ) g(x) g(x)表示分段高斯函数，定义如下:
g ( x ; μ , σ 1 , σ 2 ) = { e x p ( − 1 2 ( x − μ ) 2 / σ 1 2 ) , x < μ e x p ( − 1 2 ( x − μ ) 2 / σ 2 2 ) , x < μ g(x;\mu,\sigma_1,\sigma_2) = \begin{cases} exp(-\frac{1}{2}(x-\mu)^2/\sigma_1^2)&, {x<\mu}\\ exp(-\frac{1}{2}(x-\mu)^2/\sigma_2^2)&, {x<\mu} \end{cases} g(x;μ,σ1​,σ2​)={exp(−21​(x−μ)2/σ12​)exp(−21​(x−μ)2/σ22​)​,x<μ,x<μ​

也就是说， g ( x ) g(x) g(x)类似于钟形曲线，其峰值在 x = μ x=μ x=μ处， σ 1 σ1 σ1的扩散/标准偏差在平均值的左侧， σ 2 σ2 σ2的扩散在平均值右侧。使用以纳米为单位测量的波长 λ λ λ，我们可以近似1931色匹配函数：
x ‾ ( λ ) = 1.056 g ( λ ; 599.8 , 37.9 , 31.0 ) + 0.362 g ( λ ; 442.0 , 16.0 , 26.7 ) − 0.065 g ( λ ; 501.1 , 20.4 , 26.2 ) , y ‾ ( λ ) = 0.821 g ( λ ; 568.8 , 46.9 , 40.5 ) + 0.286 g ( λ ; 530.9 , 16.3 , 31.1 ) , z ‾ ( λ ) = 1.217 g ( λ ; 437.0 , 11.8 , 36.0 ) + 0.681 g ( λ ; 459.0 , 26.0 , 13.8 ) . \begin{aligned} {\overline {x}}(\lambda )&=1.056g(\lambda ;599.8,37.9,31.0)+0.362g(\lambda ;442.0,16.0,26.7)\\[2mu]&\quad -0.065g(\lambda ;501.1,20.4,26.2),\\[5mu] {\overline {y}}(\lambda )&=0.821g(\lambda ;568.8,46.9,40.5)+0.286g(\lambda ;530.9,16.3,31.1),\\[5mu] {\overline {z}}(\lambda )&=1.217g(\lambda ;437.0,11.8,36.0)+0.681g(\lambda ;459.0,26.0,13.8).\end{aligned} x(λ)y​(λ)z(λ)​=1.056g(λ;599.8,37.9,31.0)+0.362g(λ;442.0,16.0,26.7)−0.065g(λ;501.1,20.4,26.2),=0.821g(λ;568.8,46.9,40.5)+0.286g(λ;530.9,16.3,31.1),=1.217g(λ;437.0,11.8,36.0)+0.681g(λ;459.0,26.0,13.8).​

以上数据来自维基百科:CIE 1931 color space - Wikipedia

3. 根据光谱数据计算XYZ空间值

给定这些缩放了颜色匹配函数，带有频谱功率分布 I ( λ ) I(\lambda) I(λ)的一个颜色的RGB 三色刺激值给出为：
R = ∫ 0 ∞ I ( λ )   r ‾ ( λ )   d λ G = ∫ 0 ∞ I ( λ )   g ‾ ( λ )   d λ B = ∫ 0 ∞ I ( λ )   b ‾ ( λ )   d λ R= \int_0^\infty I(\lambda)\,\overline{r}(\lambda)\,d\lambda \\ G= \int_0^\infty I(\lambda)\,\overline{g}(\lambda)\,d\lambda \\ B= \int_0^\infty I(\lambda)\,\overline{b}(\lambda)\,d\lambda R=∫0∞​I(λ)r(λ)dλG=∫0∞​I(λ)g​(λ)dλB=∫0∞​I(λ)b(λ)dλ

实际上计算时多采用XYZ色彩空间，故对频谱分布函数积分时一般使用如下公式:
X = ∫ 0 ∞ I ( λ )   x ‾ ( λ )   d λ Y = ∫ 0 ∞ I ( λ )   y ‾ ( λ )   d λ Z = ∫ 0 ∞ I ( λ )   z ‾ ( λ )   d λ X= \int_0^\infty I(\lambda)\,\overline{x}(\lambda)\,d\lambda \\ Y= \int_0^\infty I(\lambda)\,\overline{y}(\lambda)\,d\lambda \\ Z= \int_0^\infty I(\lambda)\,\overline{z}(\lambda)\,d\lambda X=∫0∞​I(λ)x(λ)dλY=∫0∞​I(λ)y​(λ)dλZ=∫0∞​I(λ)z(λ)dλ
其中 λ \lambda λ代表波长， I ( λ ) I(\lambda) I(λ)对应了该波段下的能量强度值(其实就是灰度值) 。

x ‾ ( λ ) , y ‾ ( λ ) , z ‾ ( λ ) \overline{x}(\lambda),\overline{y}(\lambda),\overline{z}(\lambda) x(λ),y​(λ),z(λ) 分别代表了刚才说的色匹配函数。

4. 将XYZ色彩空间变换到RGB色彩空间

由于 CIE RGB 和 CIE XYZ 两者其实是同一个线性空间的不同表达，因此两者的转换可以通过转换矩阵实现。

由于 CIE XYZ 空间是一个很方便的线性空间，与具体设备无关，因此常用来做各种颜色空间转换的中间媒介。设想某个颜色的光，经过色匹配函数的计算，得到了三个 XYZ 的值，如果直接将这三个值作为 RGB 颜色显示到屏幕上，显然是不对的。我们必须把 XYZ 的值转换到屏幕的 RGB 空间中的值。
[ R l i n G l i n B l i n ] = M [ X Y Z ] \begin{bmatrix}R_{lin}\\G_{lin}\\B_{lin}\end{bmatrix} = M\begin{bmatrix}X\\Y\\Z\end{bmatrix} ​Rlin​Glin​Blin​​ ​=M ​XYZ​ ​
由于不同的显示效果的需求，存在不同RGB色彩空间。不同的RGB色彩空间色域存在差异。

不同RGB色彩空间的转换矩阵在如下链接中可以找到:

RGB/XYZ Matrices (brucelindbloom.com)

因为白色参考点相关的内容确实还没有搞懂，但是对于本文要实现的目标不影响，所以挑选不同的转换矩阵尝试就好了，本文使用D50白色参考点，sRGB色域:

5. 伽马校正

在转换完之后，我们只是得到了一个线性的 RGB 值，要被显示器正确显示的话，还需要进行 gamma 校正。每个 RGB 空间对应的 gamma 校正公式不完全一致，大多数情况下都是 C = C l i n 1 / γ C=C_{lin}^{1/\gamma} C=Clin1/γ​，其中 γ = 2.2 γ=2.2 γ=2.2，少数情况，比如 Apple RGB 采用的是 $γ=1.8 $的情况，再比如上面选择的 sRGB 空间，采用了一种分段的非线性函数进行校正：
C = { 12.92 C l i n , C l i n ≤ 0.0031308 1.055 C l i n 1 / 2.4 − 0.055 , C l i n > 0.0031308 C = \begin{cases} 12.92C_{lin} &,C_{lin} \leq 0.0031308 \\ 1.055C_{lin}^{1/2.4} - 0.055 &,C_{lin} \gt 0.0031308 \end{cases} C={12.92Clin​1.055Clin1/2.4​−0.055​,Clin​≤0.0031308,Clin​>0.0031308​
在粗略的场景下，也可以直接使用 γ=2.2 的公式代替。

6. 代码实现

# coding=utf-8
import scipy.io as sio
import numpy as np
import cv2
import matplotlib.pyplot as plt
import os
import re


class GeneratorRGB:

    def __init__(self, wave_length):
        # 波长
        self.wavelengths = wave_length
        # self.transform_matrix = [[2.36461385, -0.89654057, -0.46807328],
        #                          [-0.51516621, 1.4264081, 0.0887581],
        #                          [0.0052037, -0.01440816, 1.00920446]]
		# 转换矩阵
        self.transform_matrix = [[3.1338561, -1.6168667, -0.4906146],
                                 [-0.9787684, 1.9161415, 0.0334540],
                                 [0.0719453, -0.2289914, 1.4052427]]
	# g(x;u,sigma1,sigma2)
    def piecewise_gaussian(self, x, u, sigma1, sigma2):
        return np.where(x < u, np.exp(-0.5 * (x - u) ** 2 / sigma1 ** 2), np.exp(-0.5 * (x - u) ** 2 / sigma2 ** 2))

    def x_gaussian(self, wavelength, adj=None):
        if not adj:
            adj = [{"u": 0, "sigma1": 0, "sigma2": 0},
                   {"u": 0, "sigma1": 0, "sigma2": 0},
                   {"u": 0, "sigma1": 0, "sigma2": 0}]
        result = 0.362 * self.piecewise_gaussian(wavelength,
                                                 442.0 + adj[0]["u"],
                                                 16.0 + adj[0]["sigma1"],
                                                 26.7 + adj[0]["sigma2"])
        result += -0.065 * self.piecewise_gaussian(wavelength,
                                                   501.1 + adj[1]["u"],
                                                   20.4 + adj[1]["sigma1"],
                                                   26.2 + adj[1]["sigma2"])
        result += 1.056 * self.piecewise_gaussian(wavelength,
                                                  599.8 + adj[2]["u"],
                                                  37.9 + adj[2]["sigma1"],
                                                  31.0 + adj[2]["sigma2"])
        return result

    def y_gaussian(self, wavelength, adj=None):
        if not adj:
            adj = [{"u": 0, "sigma1": 0, "sigma2": 0},
                   {"u": 0, "sigma1": 0, "sigma2": 0}]
        result = 0.286 * self.piecewise_gaussian(wavelength,
                                                 530.9 + adj[0]["u"],
                                                 16.3 + adj[0]["sigma1"],
                                                 31.1 + adj[0]["sigma2"])
        result += 0.821 * self.piecewise_gaussian(wavelength,
                                                  568.8 + adj[1]["u"],
                                                  46.9 + adj[1]["sigma1"],
                                                  40.5 + adj[1]["sigma2"])
        return result

    def z_gaussian(self, wavelength, adj: list = None):
        if not adj:
            adj = [{"u": 0, "sigma1": 0, "sigma2": 0},
                   {"u": 0, "sigma1": 0, "sigma2": 0}]
        result = 0.980 * self.piecewise_gaussian(wavelength,
                                                 437.0 + adj[0]["u"],
                                                 11.8 + adj[0]["sigma1"],
                                                 36.0 + adj[0]["sigma2"])
        result += 0.681 * self.piecewise_gaussian(wavelength,
                                                  459.0 + adj[1]["u"],
                                                  26.0 + adj[1]["sigma1"],
                                                  13.8 + adj[1]["sigma2"])
        return result

        # color matching function

    # 计算色匹配函数值
    def get_CIE_XYZ_weights(self, adj_x: list = None, adj_y: list = None, adj_z: list = None):
        cie_x = self.x_gaussian(self.wavelengths, adj=adj_x)
        cie_y = self.y_gaussian(self.wavelengths, adj=adj_y)
        cie_z = self.z_gaussian(self.wavelengths, adj=adj_z)
        return cie_x, cie_y, cie_z
	
    # 显示CIE-xyz色匹配函数曲线
    def show_weights_curve(self, adj_x: list = None, adj_y: list = None, adj_z: list = None):
        # adj_x = [{"u": 50, "sigma1": 0, "sigma2": 0},
        #          {"u": 50, "sigma1": 0, "sigma2": 0},
        #          {"u": 50, "sigma1": 0, "sigma2": 0}]
        #
        # adj_y = [{"u": -50, "sigma1": 0, "sigma2": 0},
        #          {"u": -50, "sigma1": 0, "sigma2": 0}]
        #
        # adj_z = [{"u": 0, "sigma1": 10, "sigma2": 10},
        #          {"u": 0, "sigma1": 10, "sigma2": 30}]
        # rgb_generator.show_weights_curve(adj_x=adj_x, adj_y=adj_y, adj_z=adj_z)
        cie_x, cie_y, cie_z = self.get_CIE_XYZ_weights(adj_x=adj_x, adj_y=adj_y, adj_z=adj_z)
        plt.figure("xyz_cmf")
        plt.xlabel("wave length/nm")
        plt.ylabel("value")
        plt.plot(self.wavelengths, cie_x, linewidth=1.5, color="r", label="x")
        plt.plot(self.wavelengths, cie_y, linewidth=1.5, color="g", label="y")
        plt.plot(self.wavelengths, cie_z, linewidth=1.5, color="b", label="z")
        plt.legend()
        plt.grid()
        plt.show()
	
    # 伽马校正
    def gamma_correction(self, rgb, gamma=1.5):
        rgb = np.power(rgb, 1 / gamma)
        return rgb
	# 将高光谱图像合成为RGB图像
    # 高光谱图像 宽width,高height,通道数channels
    # hyper_image.shape == (height,width,channels) 
    # adj_x,adj_y,adj_z,用来微调,u和sigma，改变xyz曲线的形状实现特殊效果
    def get_rgb(self, hyper_image, adj_x: list = None, adj_y: list = None, adj_z: list = None, BGR=True,
                gamma_correction=False):

        cie_x, cie_y, cie_z = self.get_CIE_XYZ_weights(adj_x=adj_x, adj_y=adj_y, adj_z=adj_z)
        X = np.sum(hyper_image * cie_x, axis=2)
        Y = np.sum(hyper_image * cie_y, axis=2)
        Z = np.sum(hyper_image * cie_z, axis=2)
        R = X * self.transform_matrix[0][0] + Y * self.transform_matrix[0][1] + Z * self.transform_matrix[0][2]
        G = X * self.transform_matrix[1][0] + Y * self.transform_matrix[1][1] + Z * self.transform_matrix[1][2]
        B = X * self.transform_matrix[2][0] + Y * self.transform_matrix[2][1] + Z * self.transform_matrix[2][2]
        if BGR:
            rgb_img = np.stack((B, G, R), axis=-1)
        else:
            rgb_img = np.stack((R, G, B), axis=-1)
        rgb_img = (rgb_img - np.min(rgb_img)) / (np.max(rgb_img) - np.min(rgb_img))
        if gamma_correction:
            rgb_img = self.gamma_correction(rgb_img)
        return rgb_img

# 根据txt波长信息生成对应的npy波长文件
def envitxt2npy(txt_path, npy_savepath):
    with open(txt_path) as f:
        lines = f.readlines()[3:]
        wavelengths = [float(re.findall(r"\d+\.?\d*", line)[0]) for line in lines]
    wavelengths = np.array(wavelengths)
    np.save(os.path.join(npy_savepath, "wavelengths_" + str(len(wavelengths)) + ".npy"), wavelengths)


if __name__ == "__main__":
    envitxt2npy("envi_plot.txt", ".")
    wavelengths = np.load("wavelengths_128.npy")
    # wavelengths_96 = wavelengths[31:223:2]
    rgb_generator = GeneratorRGB(wavelengths)

    rgb_generator.show_weights_curve()

7. 附录

波段信息txt文件，一共128个波段:

ENVI ASCII Plot File [*** *** 24 11:45:14 1919]
Column 1: Wavelength
Column 2: n_D Class #1~~2##255,0,0
  371.700012   778.943339
  376.100006   770.241712
  380.500000   748.032550
  385.000000   711.031947
  389.399994   706.342978
  393.899994   694.409283
  398.299988   698.135021
  402.799988   710.817360
  407.299988   727.958409
  411.799988   737.475588
  416.299988   739.506932
  420.700012   754.845690
  425.200012   753.913201
  429.799988   766.892706
  434.299988   799.750452
  438.799988   835.210970
  443.299988   893.353225
  447.899994   924.624473
  452.399994   948.472574
  457.000000   921.470766
  461.500000   918.262809
  466.100006   917.420735
  470.700012   922.870404
  475.299988   952.725136
  479.899994   976.078360
  484.399994  1012.474985
  489.100006  1050.302592
  493.700012  1181.345992
  498.299988  1346.992164
  502.899994  1583.062086
  507.500000  1964.079566
  512.200012  2337.871609
  516.799988  2608.419530
  521.500000  2934.751055
  526.200012  3094.039783
  530.799988  3225.940928
  535.500000  3279.060880
  540.200012  3196.213382
  544.900024  3166.150090
  549.599976  3102.237492
  554.299988  3030.236890
  559.000000  2949.843882
  563.700012  2992.290536
  568.500000  2944.970464
  573.200012  2879.842676
  577.900024  2709.165160
  582.700012  2544.176612
  587.400024  2243.330319
  592.200012  2080.393008
  597.000000  1939.819168
  601.799988  1797.521398
  606.500000  1670.897529
  611.299988  1568.100663
  616.099976  1483.772152
  621.000000  1448.660639
  625.799988  1376.030741
  630.599976  1329.355636
  635.400024  1293.602170
  640.299988  1251.915009
  645.099976  1221.805907
  650.000000  1174.126582
  654.799988  1124.404461
  659.700012  1131.710066
  664.599976  1119.727547
  669.400024  1103.228451
  674.299988  1090.617842
  679.200012  1080.449066
  684.099976  1026.330922
  689.099976   993.049427
  694.000000  1033.698614
  698.900024  1022.889090
  703.799988  1008.230862
  708.799988   998.568415
  713.700012   962.066305
  718.700012   888.371308
  723.599976   878.981314
  728.599976   866.410488
  733.599976   888.247137
  738.599976   897.895118
  743.599976   917.882459
  748.599976   911.596745
  753.599976   914.708258
  758.599976   771.996986
  763.599976   730.689572
  768.599976   847.165160
  773.700012   855.122363
  778.700012   836.664858
  783.799988   827.171790
  788.799988   806.916817
  793.900024   790.221820
  798.900024   789.949367
  804.000000   809.836649
  809.099976   807.355033
  814.200012   784.625075
  819.299988   769.038577
  824.400024   786.312839
  829.500000   788.345389
  834.700012   806.656420
  839.799988   826.106088
  844.900024   827.244123
  850.099976   827.270645
  855.200012   816.968053
  860.400024   824.337553
  865.500000   814.603978
  870.700012   798.827004
  875.900024   796.328511
  881.099976   787.040386
  886.299988   772.019289
  891.500000   761.085594
  896.700012   742.569620
  901.900024   697.479204
  907.099976   697.937312
  912.400024   680.795660
  917.599976   674.689572
  922.900024   671.657625
  928.099976   671.411694
  933.400024   649.347197
  938.700012   619.044605
  943.900024   624.757083
  949.200012   619.104882
  954.500000   621.953586
  959.799988   610.510549
  965.099976   612.142857
  970.400024   613.606992
  975.700012   621.637734
  981.099976   613.901145
  986.400024   615.367691
  991.700012   616.359253