02_ccm_useDeltaE
原始文件为 .py 代码，本文是转换后的 Markdown 文件。
```.py
'''
这是用最优化 DeltaE 来计算 CCM 的代码
'''
import os
os.chdir(os.path.dirname(os.path.abspath(file)))
import numpy as np
from func.cp_hw2 import read_colorchecker_gm, readHDR, writeHDR
np.seterr(all='raise')
def LAB2RGBlinear(Lab):
L, a, b = Lab[0], Lab[1], Lab[2]
# D50 LAB to D50 XYZ
Xn = 96.422 / 100
Yn = 100 / 100
Zn = 82.521 / 100
finv = lambda x: np.where(x > 6/29, x**3, 3 * (6/29)**2 * ( x - 4/29 ))
X = Xn * finv( 1/116 * ( L + 16 ) + 1/500 * a )
Y = Yn * finv( 1/116 * ( L + 16 ) )
Z = Zn * finv( 1/116 * ( L + 16 ) - 1/200 * b )
# D50 XYZ to D65 RGB Linear
M1 = np.array([
[3.1338561, -1.6168667, -0.4906146],
[-0.9787684, 1.9161415, 0.0334540],
[0.0719453, -0.2289914, 1.4052427]
])
# (3x3) @ (3x1) = (3x1)
rgb = np.dot(M1, np.vstack((X, Y, Z)))
return rgb[:, 0]
def RGBlinear2LAB(rgb):
# D65 RGB Linear to D50 XYZ
M1 = np.array([
[3.1338561, -1.6168667, -0.4906146],
[-0.9787684, 1.9161415, 0.0334540],
[0.0719453, -0.2289914, 1.4052427]
])
M1 = np.linalg.inv(M1)
XYZ = np.dot(M1, rgb)
X, Y, Z = XYZ[0], XYZ[1], XYZ[2]
# D50 XYZ to D50 LAB
Xn = 96.422 / 100
Yn = 100 / 100
Zn = 82.521 / 100
f = lambda x: np.where(x > (6/29)**3, x**(1/3), (1/3) * x * (6/29)**(-2) + 4/29)
L = 116 * f(Y / Yn) - 16
a = 500 * (f(X / Xn) - f(Y / Yn))
b = 200 * (f(Y / Yn) - f(Z / Zn))
return np.array([L, a, b])
该函数用于生成色卡块图，用于检测正确性；同时也顺便生成色卡每个块的平均值
def make_patches(coors, img, onerow=6):
rows, cols = 24//onerow, onerow
# patches: 色卡的图，patches_rgb：色卡的平均图(24*3)
patches = np.zeros((rows*80, cols*80, 3))
patches_rgb = np.zeros((24, 3))
for i, coor in enumerate(coors):
nowpatch = img[coor[1]:coor[3], coor[0]:coor[2]]
x = (i // onerow) * 80
y = (i % onerow) * 80
patches[x:x+80, y:y+80] = nowpatch[0:80, 0:80]
patches_rgb[i] = np.mean(nowpatch[0:80, 0:80].reshape(-1, 3), axis=0)
return patches, patches_rgb
coors = []
with open('../data/24color_coor.txt', 'r') as f:
for i, line in enumerate(f):
coors.append([int(x) for x in line.split()])
保存的色卡坐标，索引有点不对，现在重排一下
coors = np.array(coors).reshape(6, 4, 4).transpose(1, 0, 2).reshape(24, 4)
加载基准 RGB
ok_rgb, ok_Lab = read_colorchecker_gm()
处理图片保存下来参数
os.makedirs('../result/stage-2', exist_ok=True)
for filename in os.listdir('../result/stage-1'):
print(filename)
hdrimg = readHDR(f'../result/stage-1/{filename}')
# 先做 AWB，还是一样，先提取色卡，然后只用白色块，最后一行第一个
now_patches, patches_rgb = make_patches(coors, hdrimg)
wr, wg, wb = patches_rgb[3*6+0]
rgain, bgain = wg / wr, wg / wb
hdrimg[..., 0] *= rgain
hdrimg[..., 2] *= bgain
# 根据coors记录的色卡坐标，生成色卡图，顺便把色卡的各个平均值计算一下
now_patches, patches_rgb = make_patches(coors, hdrimg)
# # 方法一：用最小二乘法计算 CCM
# A = np.hstack((patches_rgb, np.ones((24, 1))))
# u, *_ = np.linalg.lstsq(A, ok_rgb, rcond=None)
# 方法二：用最优化 DeltaE 计算 CCM
# ! np.longdouble 类型，很重要，否则很影响最优化结果！！
patches_rgb = patches_rgb.astype(np.longdouble)
# DeltaE 指标
def DeltaE(x):
# ! 这里也要转成 np.longdouble 类型
nowM = np.array([
[1-x[0]-x[1], x[0], x[1]],
[x[2], 1-x[2]-x[3], x[3]],
[x[4], x[5], 1-x[4]-x[5]]
]).astype(np.longdouble)
sum_deltaE = 0
for i in range(24):
now_rgb = np.dot(nowM, patches_rgb[i])
now_Lab = RGBlinear2LAB(now_rgb)
sum_deltaE += np.square(now_Lab - ok_Lab[i]).sum()
return sum_deltaE
# 初始化参数，开始最优化迭代
u = np.array([0, 0, 0, 0, 0, 0])
import scipy.optimize as optimize
x = optimize.minimize(DeltaE, u, method='SLSQP', options={'maxiter': 2000}).x
u = np.array([
[1-x[0]-x[1], x[0], x[1], 0],
[x[2], 1-x[2]-x[3], x[3], 0],
[x[4], x[5], 1-x[4]-x[5], 0],
]).T
# Make CCM
rows, cols = hdrimg.shape[0:2]
hdrimg = np.hstack((hdrimg.reshape(-1, 3), np.ones((rows*cols, 1))))
hdrimg = np.dot(hdrimg.astype(np.longdouble), u.astype(np.longdouble))
# ! 很重要：否则会有异色！！！
hdrimg[hdrimg < 0] = 0
hdrimg = hdrimg.astype(np.float32).reshape((rows, cols, 3))
# WriteHDR, Good Job.
writeHDR(f'../result/stage-2/{filename}', hdrimg)```