每个通道都是独立的,最后的可视化图表对比应该是红蓝绿三个通道前后对比的图片,色差计算也是三通道分别计算色差,输出各自的偏差和总体分析, 第一个表也就是RGB目标值体现的信息量其实就是 target_red = (220, 0, 0) # 红色目标值target_green = (0, 220, 0) # 绿色目标值target_blue = (0, 0, 220) # 蓝色目标值并且是一个64*64的led显示屏,请根据我补充的信息修改原先的错误import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.ndimage import median_filter, percentile_filter
from scipy.cluster.vq import kmeans2
from sklearn.metrics import pairwise_distances
from tqdm import tqdm
import warnings
warnings.filterwarnings("ignore", category=RuntimeWarning)
# ====================== 基础配置 ======================
TARGET_VALUE = 220
D65_WHITEPOINT = np.array([0.95047, 1.0, 1.08883]) # D65白点XYZ值
CELL_SIZE = 64 # 显示模块尺寸
GAMMA = 2.4
# ====================== 颜色空间转换函数 ======================
def gamma_correction(rgb, gamma=GAMMA, inverse=False):
"""Gamma校正处理"""
if inverse:
return np.where(rgb <= 0.04045, rgb/12.92, ((rgb+0.055)/1.055)**gamma)
else:
return np.where(rgb <= 0.0031308, 12.92*rgb, 1.055*(rgb**(1/gamma)) - 0.055)
def rgb_to_xyz(rgb, rgb_space='bt2020'):
"""自定义RGB到XYZ转换"""
# BT.2020转换矩阵
if rgb_space == 'bt2020':
M = np.array([
[0.4123908, 0.35758434, 0.18048079],
[0.21263901, 0.71516868, 0.07219232],
[0.01933082, 0.11919478, 0.95053215]
])
return np.dot(rgb, M.T)
def xyz_to_lab(xyz):
"""自定义XYZ到Lab转换"""
epsilon = 216/24389 # (6/29)^3
kappa = 24389/27 # (29/3)^3
xyz_n = xyz / D65_WHITEPOINT
f_xyz = np.where(xyz_n > epsilon, xyz_n**(1/3), (kappa*xyz_n + 16)/116)
L = 116 * f_xyz[:,1] - 16
a = 500 * (f_xyz[:,0] - f_xyz[:,1])
b = 200 * (f_xyz[:,1] - f_xyz[:,2])
return np.column_stack((L, a, b))
def lab_to_xyz(lab):
"""自定义Lab到XYZ转换"""
epsilon = 216/24389
kappa = 24389/27
y = (lab[:,0] + 16) / 116
x = lab[:,1]/500 + y
z = y - lab[:,2]/200
xyz = np.column_stack((x, y, z))
xyz = np.where(xyz**3 > epsilon, xyz**3, (116*xyz - 16)/kappa)
return xyz * D65_WHITEPOINT
import numpy as np
def delta_e_ciede2000(lab1, lab2, kL=1.0, kC=1.0, kH=1.0):
"""
CIEDE2000色差公式完整实现
参数:
lab1, lab2: 形状为(N,3)的numpy数组,包含CIE Lab颜色值
kL, kC, kH: 权重参数(默认1.0)
返回:
deltaE: 形状为(N,)的色差数组
"""
# 拆分Lab分量
L1, a1, b1 = lab1[:,0], lab1[:,1], lab1[:,2]
L2, a2, b2 = lab2[:,0], lab2[:,1], lab2[:,2]
# 步骤1:计算Cabar
C1 = np.sqrt(a1**2 + b1**2)
C2 = np.sqrt(a2**2 + b2**2)
C_avg = (C1 + C2) / 2
# 步骤2:计算G因子
G = 0.5 * (1 - np.sqrt(C_avg**7 / (C_avg**7 + 25**7)))
# 步骤3:调整a分量
a1_prime = a1 * (1 + G)
a2_prime = a2 * (1 + G)
# 步骤4:重新计算C'
C1_prime = np.sqrt(a1_prime**2 + b1**2)
C2_prime = np.sqrt(a2_prime**2 + b2**2)
C_prime_avg = (C1_prime + C2_prime) / 2
# 步骤5:计算h'
h1_prime = np.degrees(np.arctan2(b1, a1_prime)) % 360
h2_prime = np.degrees(np.arctan2(b2, a2_prime)) % 360
# 步骤6:计算Δh'
dh_prime = h2_prime - h1_prime
dh_prime = np.where((C1_prime == 0) | (C2_prime == 0), 0, dh_prime)
dh_prime = np.where(np.abs(dh_prime) <= 180, dh_prime,
np.where(dh_prime > 180, dh_prime - 360, dh_prime + 360))
# 步骤7:计算ΔH'
dH_prime = 2 * np.sqrt(C1_prime * C2_prime) * np.sin(np.radians(dh_prime)/2)
# 步骤8:计算h'_avg
h_prime_avg = np.where((C1_prime == 0) & (C2_prime == 0), 0,
(h1_prime + h2_prime + np.where(np.abs(h1_prime - h2_prime) > 180, 360, 0)) / 2)
# 步骤9:计算T补偿项
T = (1
- 0.17 * np.cos(np.radians(h_prime_avg - 30))
+ 0.24 * np.cos(np.radians(2 * h_prime_avg))
+ 0.32 * np.cos(np.radians(3 * h_prime_avg + 6))
- 0.20 * np.cos(np.radians(4 * h_prime_avg - 63)))
# 步骤10:计算Δθ和RC
delta_theta = 30 * np.exp(-((h_prime_avg - 275)/25)**2)
R_C = 2 * np.sqrt(C_prime_avg**7 / (C_prime_avg**7 + 25**7))
R_T = -np.sin(np.radians(2 * delta_theta)) * R_C
# 步骤11:计算补偿因子
S_L = 1 + (0.015 * ( (L1+L2)/2 - 50 )**2 ) / np.sqrt(20 + ( (L1+L2)/2 -50 )**2 )
S_C = 1 + 0.045 * C_prime_avg
S_H = 1 + 0.015 * C_prime_avg * T
# 步骤12:计算各分量差异
delta_L = L2 - L1
delta_C = C2_prime - C1_prime
# 最终色差计算
term1 = (delta_L/(kL*S_L))**2
term2 = (delta_C/(kC*S_C))**2
term3 = (dH_prime/(kH*S_H))**2
term4 = R_T * (delta_C/(kC*S_C)) * (dH_prime/(kH*S_H))
deltaE = np.sqrt(term1 + term2 + term3 + term4)
return deltaE
# ====================== 数据加载 ======================
def load_data(file_path):
"""加载Excel数据"""
# 读取目标值
target_df = pd.read_excel(file_path, sheet_name='RGB目标值')
# 添加列名检查
print("目标值表的列名:", target_df.columns.tolist())
targets = {
'red': target_df[target_df['220.1'] == TARGET_VALUE].values[0],
'green': target_df[target_df['220.2'] == TARGET_VALUE].values[0],
'blue': target_df[target_df['220.3'] == TARGET_VALUE].values[0]
}
# 读取测量数据
sheets = ['R_R', 'R_G', 'R_B', 'G_R', 'G_G', 'G_B', 'B_R', 'B_G', 'B_B']
raw_data = {}
for sheet in sheets:
raw_data[sheet] = pd.read_excel(file_path, sheet_name=sheet).values
# 构建三维数据立方体
data_cube = np.zeros((63, CELL_SIZE, 3, 3))
# 填充红色通道测量数据
data_cube[:,:,:,0] = np.stack([raw_data['R_R'], raw_data['R_G'], raw_data['R_B']], axis=2)
# 填充绿色通道测量数据
data_cube[:,:,:,1] = np.stack([raw_data['G_R'], raw_data['G_G'], raw_data['G_B']], axis=2)
# 填充蓝色通道测量数据
data_cube[:,:,:,2] = np.stack([raw_data['B_R'], raw_data['B_G'], raw_data['B_B']], axis=2)
return targets, data_cube
# ====================== 预处理 ======================
def preprocess_data(data_cube):
"""数据预处理"""
# Gamma校正
linear_data = gamma_correction(data_cube/255.0, inverse=True)
# 中值滤波+百分位滤波
filtered_data = np.zeros_like(linear_data)
for i in range(3): # R,G,B
for j in range(3): # R,G,B分量
ch_data = linear_data[:,:,:,i]
filtered_data[:,:,:,i] = median_filter(ch_data, size=(3,3,1))
filtered_data[:,:,:,i] = percentile_filter(filtered_data[:,:,:,i], 80, size=(5,5,1))
return filtered_data
# ====================== 动态分区校正 ======================
def find_optimal_clusters(data, max_k=20):
"""使用肘部法则确定最佳分区数"""
distortions = []
for k in tqdm(range(1, max_k+1), desc="计算最佳分区数"):
centroids, labels = kmeans2(data, k, minit='++')
distortion = np.mean(pairwise_distances(data, centroids, metric='euclidean'))
distortions.append(distortion)
# 计算二阶导数确定拐点
diff = np.diff(distortions, 2)
optimal_k = np.argmax(diff) + 2 # 补偿二阶差分索引偏移
return min(optimal_k, max_k)
def dynamic_partition_correction(lab_data):
"""自适应动态分区校正"""
points = lab_data.reshape(-1, 3)
# 自动确定最佳分区数
optimal_k = find_optimal_clusters(points)
print(f"自动确定最佳分区数: {optimal_k}")
# K-means分区
centroids, labels = kmeans2(points, optimal_k, minit='++')
partition_mask = labels.reshape(CELL_SIZE, CELL_SIZE)
# 计算每个分区的校正矩阵
correction_mats = {}
for i in range(optimal_k):
cluster_points = points[labels == i]
if len(cluster_points) == 0:
continue
target = np.mean(cluster_points, axis=0)
correction_mats[i] = np.linalg.lstsq(cluster_points, target, rcond=None)[0]
return partition_mask, correction_mats
# ====================== 异常像素补偿 ======================
def compensate_outliers(corrected_lab, deltaE_map, threshold=3.0):
"""基于色差的异常像素补偿"""
mask = deltaE_map > threshold
compensated = corrected_lab.copy()
# 使用自适应窗口尺寸
for i in range(CELL_SIZE):
for j in range(CELL_SIZE):
if mask[i,j]:
# 动态确定窗口大小
window_size = 3 if (i%8 < 6 and j%8 <6) else 5
i_min = max(0, i-window_size//2)
i_max = min(CELL_SIZE, i+window_size//2+1)
j_min = max(0, j-window_size//2)
j_max = min(CELL_SIZE, j+window_size//2+1)
# 加权中值滤波
neighborhood = corrected_lab[i_min:i_max, j_min:j_max]
weights = 1.0 / (1e-6 + deltaE_map[i_min:i_max, j_min:j_max])
compensated[i,j] = np.average(neighborhood.reshape(-1,3),
weights=weights.flatten(), axis=0)
return compensated
# ====================== 可视化模块 ======================
def plot_3d_gamut(lab_before, lab_after):
"""三维色域可视化"""
fig = plt.figure(figsize=(12, 6))
# 校正前
ax1 = fig.add_subplot(121, projection='3d')
ax1.scatter(lab_before[...,1], lab_before[...,2], lab_before[...,0],
c=lab_before[...,0], cmap='viridis', alpha=0.5)
ax1.set_xlabel('a*')
ax1.set_ylabel('b*')
ax1.set_zlabel('L*')
ax1.set_title('校正前三维色域')
# 校正后
ax2 = fig.add_subplot(122, projection='3d')
sc = ax2.scatter(lab_after[...,1], lab_after[...,2], lab_after[...,0],
c=lab_after[...,0], cmap='viridis', alpha=0.5)
ax2.set_xlabel('a*')
ax2.set_ylabel('b*')
ax2.set_zlabel('L*')
ax2.set_title('校正后三维色域')
fig.colorbar(sc, ax=[ax1, ax2], label='亮度 L*')
plt.tight_layout()
return fig
# ====================== 主流程 ======================
def main():
# 数据加载与预处理
targets, raw_data = load_data("B题附件:RGB数值.xlsx")
processed_data = preprocess_data(raw_data)
# 颜色空间转换
xyz_data = rgb_to_xyz(processed_data.reshape(-1, 3))
lab_data = xyz_to_lab(xyz_data).reshape(CELL_SIZE, CELL_SIZE, 3)
# 动态分区校正
partition_mask, correction_mats = dynamic_partition_correction(lab_data)
# 应用校正
corrected_lab = np.zeros_like(lab_data)
for i in range(CELL_SIZE):
for j in range(CELL_SIZE):
cluster_id = partition_mask[i,j]
if cluster_id in correction_mats:
corrected_lab[i,j] = np.dot(lab_data[i,j], correction_mats[cluster_id])
# 异常像素补偿
deltaE = delta_e_ciede2000(lab_data.reshape(-1,3), corrected_lab.reshape(-1,3))
deltaE_map = deltaE.reshape(CELL_SIZE, CELL_SIZE)
final_lab = compensate_outliers(corrected_lab, deltaE_map)
# 转换回RGB
corrected_xyz = lab_to_xyz(final_lab.reshape(-1, 3))
corrected_rgb = gamma_correction(xyz_to_rgb(corrected_xyz)).reshape(CELL_SIZE, CELL_SIZE, 3)
# 验证指标
final_deltaE = delta_e_ciede2000(lab_data.reshape(-1,3), final_lab.reshape(-1,3))
print(f"\n校正结果:")
print(f"平均色差: {np.mean(final_deltaE):.4f}")
print(f"最大色差: {np.max(final_deltaE):.4f}")
print(f"色度均匀性: {np.std(final_deltaE.reshape(CELL_SIZE, CELL_SIZE)):.4f}")
# 可视化
plt.figure(figsize=(18, 6))
# 色差热力图
plt.subplot(131)
plt.imshow(deltaE_map, cmap='hot', vmin=0, vmax=5)
plt.colorbar(label='ΔE2000')
plt.title('校正色差分布')
# 色域对比
plt.subplot(132)
plt.scatter(lab_data[...,1], lab_data[...,2], c='r', alpha=0.3, label='校正前')
plt.scatter(final_lab[...,1], final_lab[...,2], c='g', alpha=0.3, label='校正后')
plt.xlabel('a*')
plt.ylabel('b*')
plt.title('二维色域对比')
plt.legend()
# 三维色域显示
plot_3d_gamut(lab_data, final_lab)
plt.tight_layout()
plt.savefig('calibration_result.png', dpi=300, bbox_inches='tight')
if __name__ == '__main__':
main()