cp_hw2
复制本地路径 | 在线编辑
原始文件为 .py 代码，本文是转换后的 Markdown 文件。
```.py
!/usr/bin/python
""" This is a module for hdr imaging homework (15-463/663/862, Computational Photography, Fall 2020, CMU).
You can import necessary functions into your code as follows:
from cp_hw2 import lRGB2XYZ, XYZ2lRGB, writeEXR, read_colorchecker_gm
Depends on OpenCV to read/write HDR files"""
import numpy as np
import cv2
def read_colorchecker_gm():
"""Returns a 4x6 matrix with sRGB linear values of the Greatg-Macbeth color checker
function uses L*a*b* data under D50 illumination published by Gretag-Macbeth in 2005
(according to http://www.babelcolor.com/main_level/ColorChecker.htm)
data obtained from
Danny Pascale: "RGB coordinates of the Macbeth ColorChecker", page 6
(available from same webpage)
the function performs chromatic adaptation from D50 to D65 (sRGB standard illum.) and a conversion from L*a*b* to linear sRGB values
(c) 200x - 2011; x in {9,10}, Ivo Ihrke
Universitaet des Saarlandes / MPI Informatik
L a* b* data CIE D50 illumination for the Gretag Macbeth color checker
"""
L = [ \
37.986, \
65.711, \
49.927, \
43.139, \
55.112, \
70.719, \
62.661, \
40.020, \
51.124, \
30.325, \
72.532, \
71.941, \
28.778, \
55.261, \
42.101, \
81.733, \
51.935, \
51.038, \
96.539, \
81.257, \
66.766, \
50.867, \
35.656, \
20.461 ]
a = [ \
13.555, \
18.130, \
-4.880, \
-13.095, \
8.844, \
-33.397, \
36.067, \
10.410, \
48.239, \
22.976, \
-23.709, \
19.363, \
14.179, \
-38.342, \
53.378, \
4.039, \
49.986, \
-28.631, \
-0.425, \
-0.638, \
-0.734, \
-0.153, \
-0.421, \
-0.079 ]
b = [ \
14.059, \
17.810, \
-21.925, \
21.905, \
-25.399, \
-0.199, \
57.096, \
-45.964, \
16.248, \
-21.587, \
57.255, \
67.857, \
-50.297, \
31.370, \
28.190, \
79.819, \
-14.574, \
-28.638, \
1.186, \
-0.335, \
-0.504, \
-0.270, \
-1.231, \
-0.973 ]
L = np.reshape(L, (4, 6))
a = np.reshape(a, (4, 6))
b = np.reshape(b, (4, 6))
Lab = np.zeros((4, 6, 3))
Lab[:, :, 0] = L
Lab[:, :, 1] = a
Lab[:, :, 2] = b
# 先转换到 XYZ，其中要注意色卡上的 LAB 值把 D50 作为参考白点
XYZ = Lab_to_XYZ(Lab, 'D50')
# 计算 XYZ 到 RGB 的转换矩阵，不同色彩空间有不同的转换矩阵，以及不同的参考白点！
[ M_XYZ_to_RGB, illuminant ] = XYZ_to_RGB_linear( 'sRGB' )
# 上一步说了不同颜色空间有不同矩阵有不同参考白点，而色卡是 D50 参考白点，所以要把 XYZ 值转换一下
M = chromatic_adaptation_xyz ( 'D50', illuminant, 'Bradford' )
XYZ = apply_color_matrix( XYZ, M )
# 最后一步就是直接用上面算出的矩阵去求 RGB 值
RGB = apply_color_matrix( XYZ, M_XYZ_to_RGB ) # we want linear RGB values for our HDR measurements
r = RGB[:, :, 0]
g = RGB[:, :, 1]
b = RGB[:, :, 2]
rgb = np.vstack((r.flatten(), g.flatten(), b.flatten())).T
return rgb, Lab.reshape(-1, 3)
def Lab_to_XYZ(Lab, illuminant='D65'):
"""© Ivo Ihrke 2011
Universitaet des Saarlandes / MPI Informatik
convert from L*a*b* (CIELAB) to XYZ color space
using one of the CIE standard illuminants
source:
http://en.wikipedia.org/wiki/Lab_color_space#CIE_XYZ_to_CIE_L.2Aa.2Ab.2A_.28CIELAB.29_and_CIELAB_to_CIE_XYZ_conversions
2011-01-06
see also:
illuminant_xyz
"""
def finv(val):
val_out = np.where(val > 6/29, val**3, 3 * (6/29)**2 * ( val - 4/29 ))
return val_out
Xn, Yn, Zn = illuminant_xyz(illuminant)
L = Lab[:, :, 0]
a = Lab[:, :, 1]
b = Lab[:, :, 2]
XYZ = np.zeros_like(Lab)
XYZ[:, :, 0] = Xn * finv( 1/116 * ( L + 16 ) + 1/500 * a )
XYZ[:, :, 1] = Yn * finv( 1/116 * ( L + 16 ) )
XYZ[:, :, 2] = Zn * finv( 1/116 * ( L + 16 ) - 1/200 * b )
return XYZ
def illuminant_xyz(illuminant_in):
"""© Ivo Ihrke 2011
Universitaet des Saarlandes / MPI Informatik
define xyz coordinates of CIE standard illuminants
1931 and 1964
data from
http://en.wikipedia.org/wiki/Standard_illuminant
Name CIE 1931 2° CIE 1964 10° CCT (*) in K Hue Note
x2 y2 x10 y10
A 0.44757, 0.40745, 0.45117, 0.40594, 2856 Incandescent / Tungsten
B 0.34842, 0.35161, 0.34980, 0.35270, 4874 {obsolete} Direct sunlightat noon
C 0.31006, 0.31616, 0.31039, 0.31905, 6774 {obsolete} Average / North sky Daylight
D50 0.34567, 0.35850, 0.34773, 0.35952, 5003 Horizon Light. ICC profile PCS
D55 0.33242, 0.34743, 0.33411, 0.34877, 5503 Mid-morning / Mid-afternoon Daylight
D65 0.31271, 0.32902, 0.31382, 0.33100, 6504 Noon Daylight: Television, sRGB color space
D75 0.29902, 0.31485, 0.29968, 0.31740, 7504 North sky Daylight
E 1/3 , 1/3 , 1/3 , 1/3 , 5454 Equal energy
F1 0.31310, 0.33727, 0.31811, 0.33559, 6430 Daylight Fluorescent
F2 0.37208, 0.37529, 0.37925, 0.36733, 4230 Cool White Fluorescent
F3 0.40910, 0.39430, 0.41761, 0.38324, 3450 White Fluorescent
F4 0.44018, 0.40329, 0.44920, 0.39074, 2940 Warm White Fluorescent
F5 0.31379, 0.34531, 0.31975, 0.34246, 6350 Daylight Fluorescent
F6 0.37790, 0.38835, 0.38660, 0.37847, 4150 Lite White Fluorescent
F7 0.31292, 0.32933, 0.31569, 0.32960, 6500 D65 simulator, Daylight simulator
F8 0.34588, 0.35875, 0.34902, 0.35939, 5000 D50 simulator, Sylvania F40 Design 50
F9 0.37417, 0.37281, 0.37829, 0.37045, 4150 Cool White Deluxe Fluorescent
F10 0.34609, 0.35986, 0.35090, 0.35444, 5000 Philips TL85, Ultralume 50
F11 0.38052, 0.37713, 0.38541, 0.37123, 4000 Philips TL84, Ultralume 40
F12 0.43695, 0.40441, 0.44256, 0.39717, 3000 Philips TL83, Ultralume 30
(*) CCT= correlated color temperature
standard is the 1931 definition
illuinant_in = 'A','B','C', 'D50','D55','D65','D75' etc.
for 1964 version use
illuinant_in = 'A_64','B_64','C_64', 'D50_64','D55_64','D65_64','D75_64' etc.
verification performed by checking
http://brucelindbloom.com/index.html?Eqn_ChromAdapt.html
"""
ind1931 = np.arange(0, 2)
ind1964 = np.arange(2, 4)
illuminants = [ 'A', 'B', 'C', 'D50', 'D55','D65','D75','E','F1', \
'F2','F3','F4','F5','F6','F7','F8','F9','F10','F11','F11']
xy = np.array([[0.44757, 0.40745, 0.45117, 0.40594], \
[0.34842, 0.35161, 0.34980, 0.35270], \
[0.31006, 0.31616, 0.31039, 0.31905], \
[0.34567, 0.35850, 0.34773, 0.35952], \
[0.33242, 0.34743, 0.33411, 0.34877], \
[0.31271, 0.32902, 0.31382, 0.33100], \
[0.29902, 0.31485, 0.29968, 0.31740], \
[1/3 , 1/3 , 1/3 , 1/3 ], \
[0.31310, 0.33727, 0.31811, 0.33559], \
[0.37208, 0.37529, 0.37925, 0.36733], \
[0.40910, 0.39430, 0.41761, 0.38324], \
[0.44018, 0.40329, 0.44920, 0.39074], \
[0.31379, 0.34531, 0.31975, 0.34246], \
[0.37790, 0.38835, 0.38660, 0.37847], \
[0.31292, 0.32933, 0.31569, 0.32960], \
[0.34588, 0.35875, 0.34902, 0.35939], \
[0.37417, 0.37281, 0.37829, 0.37045], \
[0.34609, 0.35986, 0.35090, 0.35444], \
[0.38052, 0.37713, 0.38541, 0.37123], \
[0.43695, 0.40441, 0.44256, 0.39717]])
for i in range(len(illuminants)):
if illuminants[i] == illuminant_in:
index_row = i
index_cols = ind1931
if len(illuminant_in) > 3:
if illuminant_in[-3:] == '_64':
index_cols = ind1964
data = xy[index_row, index_cols]
X, Y, Z = xyY_to_XYZ(data[0], data[1], 1)
return X, Y, Z
def xyY_to_XYZ(x, y, Y):
"""© Ivo Ihrke 2011
Universitaet des Saarlandes / MPI Informatik
"""
Xo = Y * x / y
Yo = Y
Zo = Y * ( 1 - x - y ) / y
return Xo, Yo, Zo
def chromatic_adaptation_xyz(from_illum, to_illum, method='Bradford'):
"""
(c) Ivo Ihrke 2011
Universitaet des Saarlandes / MPI Informatik
computes chromatic adaptation matrix for XYZ space
chromatic adaptation
input: illuminant names (from, to) and method
method = 'XYZScaling', 'Bradford', 'vonKries'
default 'Bradford'
implementation and choice of default according to
http://brucelindbloom.com/index.html?Eqn_ChromAdapt.html
the conversion matrices given on this webpage seem to
use the 'XYZScaling' method which is mentioned as the
worst choice.
"""
fX, fY, fZ = illuminant_xyz( from_illum )
tX, tY, tZ = illuminant_xyz( to_illum )
# setup Ma (cone response domain transform) according to <method>
mselected = 0
if method == 'XYZScaling':
Ma = np.eye(3)
mselected = 1
if method == 'Bradford':
Ma = np.array([ [0.8951000, 0.2664000, -0.1614000], \
[-0.7502000, 1.7135000, 0.0367000], \
[0.0389000, -0.0685000, 1.0296000] ])
mselected = 1
if method == 'vonKries':
Ma = np.array([ [0.4002400, 0.7076000, -0.0808100], \
[-0.2263000, 1.1653200, 0.0457000], \
[0.0000000, 0.0000000, 0.91822000]])
mselected = 1
if not mselected:
# print('chromatic_adaptation_xyz: unknown transform - returning unit matrix');
M = np.eye(3)
else:
# compute transform matrix
# rho, gamma, beta
from_rgb = Ma @ np.array((fX, fY, fZ))
to_rgb = Ma @ np.array((tX, tY, tZ))
M = np.linalg.inv(Ma) @ np.diag(to_rgb / from_rgb) @ Ma
return M
def XYZ_to_RGB_linear(rgb_space='sRGB'):
""" Data from
http://brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html
Name Gamma Reference White Red Primary Green Primary Blue Primary Volume (deltaE^3) Lab Gamut Efficiency Coding Efficiency
x y Y x y Y x y Y
Lab Gamut - D50 - - - - - - - - - 2,381,085 97.0 35.1
Adobe RGB (1998) 2.2 D65 0.6400 0.3300 0.297361 0.2100 0.7100 0.627355 0.1500 0.0600 0.075285 1,208,631 50.6 100.0
Apple RGB 1.8 D65 0.6250 0.3400 0.244634 0.2800 0.5950 0.672034 0.1550 0.0700 0.083332 798,403 33.5 100.0
Best RGB 2.2 D50 0.7347 0.2653 0.228457 0.2150 0.7750 0.737352 0.1300 0.0350 0.034191 2,050,725 77.6 96.5
Beta RGB 2.2 D50 0.6888 0.3112 0.303273 0.1986 0.7551 0.663786 0.1265 0.0352 0.032941 1,717,450 69.3 99.0
Bruce RGB 2.2 D65 0.6400 0.3300 0.240995 0.2800 0.6500 0.683554 0.1500 0.0600 0.075452 988,939 41.5 100.0
CIE RGB 2.2 E 0.7350 0.2650 0.176204 0.2740 0.7170 0.812985 0.1670 0.0090 0.010811 1,725,261 64.3 96.1
ColorMatch RGB 1.8 D50 0.6300 0.3400 0.274884 0.2950 0.6050 0.658132 0.1500 0.0750 0.066985 836,975 35.2 100.0
Don RGB 4 2.2 D50 0.6960 0.3000 0.278350 0.2150 0.7650 0.687970 0.1300 0.0350 0.033680 1,802,358 72.1 98.8
ECI RGB v2 L* D50 0.6700 0.3300 0.320250 0.2100 0.7100 0.602071 0.1400 0.0800 0.077679 1,331,362 55.3 99.7
Ekta Space PS5 2.2 D50 0.6950 0.3050 0.260629 0.2600 0.7000 0.734946 0.1100 0.0050 0.004425 1,623,899 65.7 99.5
NTSC RGB 2.2 C 0.6700 0.3300 0.298839 0.2100 0.7100 0.586811 0.1400 0.0800 0.114350 1,300,252 54.2 99.9
PAL/SECAM RGB 2.2 D65 0.6400 0.3300 0.222021 0.2900 0.6000 0.706645 0.1500 0.0600 0.071334 849,831 35.7 100.0
ProPhoto RGB 1.8 D50 0.7347 0.2653 0.288040 0.1596 0.8404 0.711874 0.0366 0.0001 0.000086 2,879,568 91.2 87.3
SMPTE-C RGB 2.2 D65 0.6300 0.3400 0.212395 0.3100 0.5950 0.701049 0.1550 0.0700 0.086556 758,857 31.9 100.0
sRGB #2.2 D65 0.6400 0.3300 0.212656 0.3000 0.6000 0.715158 0.1500 0.0600 0.072186 832,870 35.0 100.0
Wide Gamut RGB 2.2 D50 0.7350 0.2650 0.258187 0.1150 0.8260 0.724938 0.1570 0.0180 0.016875 2,164,221 77.6 91.9
#2.2 - actual transform is more complex (see http://brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html )
deltaE^3 is the volume of the gamut in Lab space
Coding efficiency is the amount of the gamut inside the horseshoe diagram
implementation verified with Bruce Lindbloom's data
"""
color_spaces = [ 'Adobe RGB (1998)', 'Apple RGB', 'Best RGB', 'Beta RGB', 'Bruce RGB','CIE RGB', \
'ColorMatch RGB','Don RGB 4','ECI RGB v2','Ekta Space PS5','NTSC RGB','PAL/SECAM RGB','ProPhoto RGB','SMPTE-C RGB','sRGB','Wide Gamut RGB']
reference_whites = [ 'D65', 'D65','D50','D50','D65','E','D50','D50','D50','D50','C','D65','D50','D65','D65','D50' ]
xyY_red = np.array([[ 0.6400, 0.3300, 0.297361], \
[0.6250, 0.3400, 0.244634], \
[0.7347, 0.2653, 0.228457], \
[0.6888, 0.3112, 0.303273], \
[0.6400, 0.3300, 0.240995], \
[0.7350, 0.2650, 0.176204], \
[0.6300, 0.3400, 0.274884], \
[0.6960, 0.3000, 0.278350], \
[0.6700, 0.3300, 0.320250], \
[0.6950, 0.3050, 0.260629], \
[0.6700, 0.3300, 0.298839], \
[0.6400, 0.3300, 0.222021], \
[0.7347, 0.2653, 0.288040], \
[0.6300, 0.3400, 0.212395], \
[0.6400, 0.3300, 0.212656], \
[0.7350, 0.2650, 0.258187]])
xyY_green = np.array([[0.2100, 0.7100, 0.627355], \
[0.2800, 0.5950, 0.672034], \
[0.2150, 0.7750, 0.737352], \
[0.1986, 0.7551, 0.663786], \
[0.2800, 0.6500, 0.683554], \
[0.2740, 0.7170, 0.812985], \
[0.2950, 0.6050, 0.658132], \
[0.2150, 0.7650, 0.687970], \
[0.2100, 0.7100, 0.602071], \
[0.2600, 0.7000, 0.734946], \
[0.2100, 0.7100, 0.586811], \
[0.2900, 0.6000, 0.706645], \
[0.1596, 0.8404, 0.711874], \
[0.3100, 0.5950, 0.701049], \
[0.3000, 0.6000, 0.715158], \
[0.1150, 0.8260, 0.724938]])
xyY_blue = np.array([[0.1500, 0.0600, 0.075285], \
[0.1550, 0.0700, 0.083332], \
[0.1300, 0.0350, 0.034191], \
[0.1265, 0.0352, 0.032941], \
[0.1500, 0.0600, 0.075452], \
[0.1670, 0.0090, 0.010811], \
[0.1500, 0.0750, 0.066985], \
[0.1300, 0.0350, 0.033680], \
[0.1400, 0.0800, 0.077679], \
[0.1100, 0.0050, 0.004425], \
[0.1400, 0.0800, 0.114350], \
[0.1500, 0.0600, 0.071334], \
[0.0366, 0.0001, 0.000086], \
[0.1550, 0.0700, 0.086556], \
[0.1500, 0.0600, 0.072186], \
[0.1570, 0.0180, 0.016875 ]])
for i in range(len(color_spaces)):
if color_spaces[i] == rgb_space:
index_row = i
Xr, Yr, Zr = xyY_to_XYZ( xyY_red[index_row, 0], xyY_red[index_row, 1], xyY_red[index_row, 2] )
Xg, Yg, Zg = xyY_to_XYZ( xyY_green[index_row, 0], xyY_green[index_row, 1], xyY_green[index_row, 2] )
Xb, Yb, Zb = xyY_to_XYZ( xyY_blue[index_row, 0], xyY_blue[index_row, 1], xyY_blue[index_row, 2] )
illuminant = reference_whites[index_row]
Xw, Yw, Zw = illuminant_xyz(illuminant)
S = np.linalg.inv( np.array([[Xr, Xg, Xb], [Yr, Yg, Yb], [Zr, Zg, Zb]]) ) @ np.array([Xw, Yw, Zw])
# this matrix is RGB to XYZ
M = np.array([ [S[0]*Xr, S[1]*Xg, S[2]*Xb], [S[0]*Yr, S[1]*Yg, S[2]*Yb], [S[0]*Zr, S[1]*Zg, S[2]*Zb]])
M = np.linalg.inv(M)
return M, illuminant
def apply_color_matrix(I, matrix):
""" Applies a 3x3 color matrix to a 3-channel image I
(c) 2011 Ivo Ihrke
Universitaet des Saarlandes
MPI Informatik
"""
vec = np.reshape(I, (I.shape[0]*I.shape[1], 3))
out_vec = matrix @ vec.T
out = np.reshape(out_vec.T, (I.shape[0], I.shape[1], 3))
return out
def lRGB2XYZ(lRGB):
""" linear RGB to XYZ
"""
invM = XYZ_to_RGB_linear('sRGB')[0]
M = np.linalg.inv(invM)
# linear rgb
R = lRGB[:,:,0]
G = lRGB[:,:,1]
B = lRGB[:,:,2]
# through matrix to XYZ
X = M[0,0] * R + M[0,1] * G + M[0,2] * B
Y = M[1,0] * R + M[1,1] * G + M[1,2] * B
Z = M[2,0] * R + M[2,1] * G + M[2,2] * B
XYZ = np.dstack((X, Y, Z))
return XYZ
def XYZ2lRGB(XYZ):
""" XYZ to linear RGB
"""
invM = XYZ_to_RGB_linear('sRGB')[0]
# XYZ
X = XYZ[:,:,0]
Y = XYZ[:,:,1]
Z = XYZ[:,:,2]
# through matrix to lRGB
R = invM[0,0] * X + invM[0,1] * Y + invM[0,2] * Z
G = invM[1,0] * X + invM[1,1] * Y + invM[1,2] * Z
B = invM[2,0] * X + invM[2,1] * Y + invM[2,2] * Z
RGB = np.dstack((R, G, B))
return RGB
def writeHDR(name, data):
#flip from rgb to bgr for cv2
cv2.imwrite(name, data[:, :, ::-1].astype(np.float32))
def readHDR(name):
raw_in = cv2.imread(name, flags=cv2.IMREAD_ANYDEPTH)
#flip from bgr to rgb
return raw_in[:, :, ::-1]```