下文摘自本人实习总结,涉及到P3,cieXYZ,sRGB之间的转换矩阵和转换方法,此外还有一些个人遇到的问题和思考,希望给大家带来一点帮助,欢迎各位指正
比较方法
- 将P3图片中的像素RGB值转换到CIEXYZ空间
- XYZ映射到xyY二维色度图上,判断是否超过了sRGB的三角形区域
P3转XYZ
转换矩阵计算
我们可以从apple的P3.icc文件出发来计算P3到XYZ的转换矩阵,具体文件内容在MAC OS下,可以直接从设置->显示器->颜色->色彩同步实用工具中直接查看。具体内容如下:
需要的数据为上图红圈中的四个,他们分别代表:
- wtpt, D65下白点XYZ坐标
- rXYZ, D50下红点XYZ坐标,绿点,蓝点以此类推
- rTRC, 反gamma校正曲线
- chad,颜色适应矩阵,将红点XYZ坐标从D50下映射到D65
颜色适应矩阵的使用如下:
Xs代表D65下的XYZ坐标,Xd代表D50下的XYZ坐标,MA代表上述矩阵
通过上述信息求解转换矩阵M的公式如下:
其中Xw为D65白点XYZ坐标,xr为D65红点XYZ坐标映射到xyY后的x值,公式如下:
更多数学信息可见网页,Welcome to Bruce Lindbloom's Web Site
基于Python的计算脚本
相关环境:
pip install numpy
具体代码:
import numpy as np
#D50下Display P3三原色数据
rXYZd50 = np.array([[0.515], [0.241], [-0.001]])
gXYZd50 = np.array([[0.292], [0.692], [0.042]])
bXYZd50 = np.array([[0.157], [0.067], [0.784]])
#色适应矩阵
BM = np.array([[1.0478112, 0.0228866, -0.0501270],
[0.0295424, 0.9904844, -0.0170491],
[-0.0092345, 0.0150436, 0.7521316]])
BM_inv = np.linalg.inv(BM)
#色适应到D65
rXYZ = np.dot(BM_inv,rXYZd50)
gXYZ = np.dot(BM_inv,gXYZd50)
bXYZ = np.dot(BM_inv,bXYZd50)
#计算色品坐标
rxy = np.array([rXYZ[0,0], rXYZ[1,0]])/rXYZ.sum()
gxy = np.array([gXYZ[0,0], gXYZ[1,0]])/gXYZ.sum()
bxy = np.array([bXYZ[0,0], bXYZ[1,0]])/bXYZ.sum()
#D65 白点坐标
W = np.array([[0.950], [1.000], [1.089]])
wxy = np.array([W[0,0], W[1,0]])/W.sum()
#计算p3 2 xyz转换矩阵M
XYZ = np.array([[ rxy[0]/rxy[1], gxy[0]/gxy[1], bxy[0]/bxy[1]],
[1, 1, 1],
[(1-rxy[0]-rxy[1])/rxy[1], (1-gxy[0]-gxy[1])/gxy[1], (1-bxy[0]-bxy[1])/bxy[1]]])
XYZ_inv = np.linalg.inv(XYZ)
S = np.dot(XYZ_inv, W)
M = S.T*XYZ
print(M)
计算结果:
[[4.86119575e-01 2.65794247e-01 1.98086178e-01]
[2.28623913e-01 6.91647657e-01 7.97284306e-02]
[6.71950077e-05 4.52826221e-02 1.04365018e+00]]
额外补充,sRGB到XYZ转换矩阵计算
与P3的矩阵计算方式一致,相关数据可以用https://www.color.org/profileinspector.xalter中的工具读取Windows下的icc文件,文件路径C:\Windows\System32\spool\drivers\color
具体内容如下:
相关词条定义和P3中完全一致,虽然没有颜色适应矩阵,但该矩阵和P3.icc中的矩阵完全一致,利用上述数据和计算方式,计算sRGB转换矩阵如下:
[[0.41224672 0.35771267 0.18049061]
[0.21255252 0.71525374 0.07219374]
[0.01930836 0.1191795 0.95056214]]
和opencv官方给出的矩阵对比,二者还是相当接近的
转换流程
在获得了转换矩阵后,具体的转换流程如下:
- 读取P3图片像素RGB值(0,255)
- RGB值归一化到(0,1)
- 反gamma变换
- 左乘转换矩阵变为XYZ
- 映射到xyY色度图
归一化就是直接除255,反gamma需要利用到icc文件中读取的函数
x为归一化后的rgb值,f(x)即为反gamma变换后的值
矩阵的使用公式如下,RGB为反gamma后的值,
XYZ映射到xyY的公式和计算矩阵部分提到的一致
这样我们就成功映射像素RGB值到色度图上了
色度图比较
在二维的色度图中,P3和sRGB的范围如下:
每一个P3像素映射到二维点后都可以判断它是否超过了sRGB色域,具体需要利用sRGB三原色转换到色度图的坐标来计算三角形的三条直线方程
问题与思考
性能优化
利用上述判断方法可以计算每个像素是否超过sRGB色域,但对于一个4032*3024的图片,单线程计算时长长达3min,由于针对每个像素的计算任务一致,我们可以采用多进程的方式来加速计算。(注:Python的多线程是伪多线程,使用后不仅没有加速还会因为线程分配开销而增加时长)加速后,一个四核CPU可以提速三倍以上。压缩时间到数十秒。
纯色图
对于如下安卓图:
这个P3图中有两种纯红色,其RGB值分别为(255,0,0)和(237,0,0),一般工程中常用这张图测试显示设备的色彩管理。很多地方的文章会介绍说,(255,0,0)是超过sRGB色域的值,而(237,0,0)是在sRGB色域内的值,如果我们用上面的计算方式测试会发生什么呢?
测试结果是,只要是纯红色,无论是(255.0,0)还是(1,0,0),利用P3的矩阵转换到色度图后都会是P3的红色原点。他们都超出了sRGB色域。因此,“(255,0,0)是超过sRGB色域的值,而(237,0,0)是在sRGB色域内的值”这句话是不正确的。二者其实都超过了sRGB色域。
那么为什么还用这个图来测试色彩管理呢?我们之所以能看到小机器人,是因为P3下255和237的红色存在亮度差异,假设这张P3图片在显示时被转换到sRGB色域,导致两个颜色都变成了(255,0,0),那么我们就看不到小机器人了。
P3到sRGB的色域转换有很多方式,下面给出一种基于截断的最简单的转换方式。
- 读取P3像素RGB值(0,255)
- RGB值归一化(0,1)
- RGB值反gamma变换
- 利用P3toXYZ的矩阵转换到XYZ空间
- 利用sRGBtoXYZ的逆矩阵从XYZ转换到sRGB空间
- RGB值gamma变换
- RGB值反归一化
- 大于255的部分映射到255,小于0的部分映射到0
通过这种方式转换P3图片后。有一个阈值(241,0,0)大于241的纯红色在从P3转换到sRGB时都会映射到(255,0,0)而小于241的纯红色在转换后可以映射到(255,0,0)内的某个值。因此,我们认为P3中大于241的纯红色溢出了,这部分属于超过sRGB且不能正常转换的红色,而小于241的纯红色属于超过sRGB但能正常转换的红色。
在工程上常常想要区分出哪些红色能够正常转换,哪些红色溢出了,如果只在色度图上进行比较,是判断不了的。因此,在计算时需要额外添加了关于纯色情况的判断,这样才能够正常标记出上图中的小机器人部分。
特别注意的是,P3到sRGB的转换不只一种方式,例如使用RGB(sRGB) to P3 Converter - Color Space Converter该网站的自动转换,将得到与上文不一样的阈值。并且,进行了转换并不代表读者当前屏幕能看到小机器人,具体的原因可以详细阅读你真的理解色彩管理吗?安卓苹果色彩管理深度解析 - 少数派中的内容
脚本使用方法
文件名:muti_P3_sRGB.py
相关环境
除了基础的python环境外,还需要opencv和numpy两个包,没有的可以用以下命令安装:
pip install opencv-python pip install numpy
使用方法
必须自行设定的有:
- 图片输入文件夹路径line 144
- 结果输出文件夹路径line 147
可选设定的有:
- 多进程数量line 163
- 图片切片数量line 165(每个切片会输入到一个进程,切片数最好等于进程数)
运行:
- 运行时如果开启了图片显示line 209,当前结果展示后,按任意键处理文件夹中下一张图片
- 运行时没有开启图片显示,会自动处理文件夹内所有图片
结果展示
运行过程中可以同时展示原图和标记超限像素后的图片
标记后的图片会存储到结果输出文件夹中
终端中会输出:
- 图片读取状态
- 多进程处理状态
- 图片所有像素超过sRGB的比例
- 总时间消耗
具体代码
import cv2
import numpy as np
import math
import os
from multiprocessing import Pool
import time
import shutil
# P3 2 sRGB Matrix
M = np.array([[4.86119575e-01, 2.65794247e-01, 1.98086178e-01],
[2.28623913e-01, 6.91647657e-01, 7.97284306e-02],
[6.71950077e-05, 4.52826221e-02, 1.04365018e+00]])
# gamma 校正
def gamma(x):
if(x <= 0.0031308):
return x*12.92
else:
return math.pow((1.055*x), (1/2.4)) - 0.055
# 反gamma 校正
def gamma_inv(x):
if(x <= 0.04045):
return x/12.92
else:
return math.pow(((x+0.055)/1.055), 2.4)
# 针对纯色,判断转换到sRGB时是否会截断
def pure_color(r, g, b):
if(r != 0 and g == 0 and b == 0):
if(r>241):
return True
else:
return False
if(r == 0 and g != 0 and b == 0):
if(g>=250):
return True
else:
return False
if(r == 0 and g == 0 and b != 0):
if(b>=245):
return True
else:
return False
#非纯色情况
return True
# 计算P3像素是否超限sRGB
def pixel_process(indexes):
print('enter process', indexes[0])
start = time.time()
cnt = 0 #当前进程超限像素数量
index = indexes[0]
slice = indexes[1]
height = indexes[2]
width = indexes[3]
img_path = indexes[4]
save_out = indexes[5]
img_name = indexes[6]
img = cv2.imread(img_path)
# 计算当前进程计算范围
h_min = 0
h_max = height
w_min = math.ceil(index*width/slice)
w_max = math.ceil((index+1)*width/slice)
# 裁剪图片
img_crop = img[h_min:h_max, w_min:w_max]
h, w, _ = img_crop.shape
# P3_RGB 2 CIE_XYZ
for row in range(h):
for col in range(w):
#Numpy读取像素 BGR通道
b = img_crop.item(row, col, 0)
g = img_crop.item(row, col, 1)
r = img_crop.item(row, col, 2)
if(not pure_color(r, g, b)):
continue
#归一化
r = r/255
g = g/255
b = b/255
#反gamma变换
r = gamma_inv(r)
g = gamma_inv(g)
b = gamma_inv(b)
#P3 2 XYZ
P3_RGB = np.array([[r], [g], [b]])
P3_XYZ = np.dot(M, P3_RGB)
#跳过黑色
if(P3_XYZ.sum() == 0):
continue
#XYZ 2 xyY
p3_xy = np.array([P3_XYZ[0,0], P3_XYZ[1,0]])/P3_XYZ.sum()
x = p3_xy[0]
y = p3_xy[1]
#判断像素是否处于sRGB范围
flag = True
y1 = 3.6*x-0.48
y2 = (27/49)*x+0.33-(17.28/49)
y3 = (-27/34)*x+0.33+(17.28/34)
if(x>=0.15 and x<=0.3 and y>=y2 and y<=y1):
flag = False
if(x>0.3 and x<=0.64 and y>=y2 and y<=y3):
flag = False
#统计并绘制超过sRGB范围像素
if(flag == True):
cnt += 1
img_crop.itemset((row, col, 0), 255)
img_crop.itemset((row, col, 1), 0)
img_crop.itemset((row, col, 2), 0)
print (f'process {index} cost: {(time.time() - start):.3f} s')
# 保存当前进程图像处理结果
save_out = save_out + 'temp/'
folder = os.path.exists(save_out)
if(not folder):
os.makedirs(save_out)
out_name = img_name.split('.')[0]
out_name_temp = out_name + '_' + str(index)
save_path = save_out + out_name_temp + '.jpg'
cv2.imwrite(save_path, img_crop)
#返回当前进程超限像素数量
return cnt
if __name__ == '__main__':
#待读取P3图片文件夹路径
file_root = "D:/P3sRGB/testinput/"
file_list = os.listdir(file_root)
#结果保存文件夹路径
save_out = "D:/P3sRGB/result/"
for img_name in file_list:
img_path = file_root + img_name
#读取图片
img = cv2.imread(img_path)
print("read image:", img_path)
start = time.time()
height, width, _ = img.shape
#总像素数量和超限像素数量
total = height*width
cnt = 0
#进程池数量
pool = Pool(4)
#图像切片数量
slice = 4
#分配进程
res_list = []
for i in range(slice):
indexes = [i, slice, height, width, img_path, save_out, img_name]
res = pool.apply_async(func=pixel_process, args=(indexes, ))
res_list.append(res)
#关闭进程
pool.close()
pool.join()
#统计每个进程超限像素值
for res in res_list:
cnt += res.get()
print(f'{img_name} 所有像素超过sRGB比例 {cnt/total:.3f}')
print (f'total time consume: {(time.time() - start):.3f} s')
#读取每个进程图像处理结果
file_root_temp = save_out + 'temp/'
file_list_temp = os.listdir(file_root_temp)
#组合每个进程处理结果
count = 0
for img_name_temp in file_list_temp:
img_path_temp = file_root_temp + img_name_temp
count +=1
if(count == 1):
img_out = cv2.imread(img_path_temp)
else:
img_temp =cv2.imread(img_path_temp)
img_out = cv2.hconcat([img_out, img_temp])
#保存最终处理后图像
out_name = img_name.split('.')[0]
save_path = save_out + out_name + '.jpg'
cv2.imwrite(save_path, img_out)
#删除临时文件
shutil.rmtree(file_root_temp)
#结果展示
cv2.imshow("origin", img)
cv2.imshow("outlier", img_out)
cv2.waitKey(0)
cv2.destroyAllWindows()