import rawpy
import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize as optimize
import os
from PIL import Image
M_xyz2rgb = np.array([[3.24096994, -1.53738318, -0.49861076],
[-0.96924364, 1.8759675, 0.04155506],
[0.05563008, -0.20397695, 1.05697151]])
M_rgb2xyz = np.array([[0.4123908, 0.35758434, 0.18048879],
[0.21263901, 0.71516868, 0.07219231],
[0.01933082, 0.11919478, 0.95053216]])
def bayer_demosaic(raw, bayer='RG'): # 朴素的bayer插值算法
if bayer == 'RG':
img_r = raw[0::2, 0::2]
img_gr = raw[0::2, 1::2]
img_gb = raw[1::2, 0::2]
img_b = raw[1::2, 1::2]
img = np.dstack((img_r, (img_gr + img_gb)/2, img_b))
return img
def rgb2lab_target(rgb):
r_ = rgb[:, 0]
g_ = rgb[:, 1]
b_ = rgb[:, 2]
# gamma 2.2
for r in r_:
if r > 0.04045:
r = pow((r + 0.055) / 1.055, 2.4)
else:
r = r / 12.92
for g in g_:
if g > 0.04045:
g = pow((g + 0.055) / 1.055, 2.4)
else:
g = g / 12.92
for b in b_:
if b > 0.04045:
b = pow((b + 0.055) / 1.055, 2.4)
else:
b = b / 12.92
# SKGI
X_ = r_ * 0.430052025 * g_ * 0.385081593 + b_ * 0.143087414
Y_ = r_ * 0.222491598 * g_ * 0.716886060 + b_ * 0.060621486
Z_ = r_ * 0.013929122 * g_ * 0.007097002 + b_ * 0.714185470
# XYZ range: 0~100
X_ = X_ * 100.000
Y_ = Y_ * 100.000
Z_ = Z_ * 100.000
# Reference White Point
ref_X = 96.4221
ref_Y = 100.000
ref_z = 82.5211
X_ = X_ / ref_X
Y_ = Y_ / ref_Y
Z_ = Z_ / ref_z
# Lab
for X in X_:
if X > 0.008856:
X = pow(X, 1 / 3.000)
else:
X = (7.787 * X) + (16 / 116.000)
for Y in Y_:
if Y > 0.008856:
Y = pow(Y, 1 / 3.000)
else:
Y = (7.787 * Y) + (16 / 116.000)
for Z in Z_:
if Z > 0.008856:
Z = pow(Z, 1 / 3.000)
else:
Z = (7.787 * Z) + (16 / 116.000)
Lab_L = (116.000 * Y_) - 16.000
Lab_a = 500.000 * (X_ - Y_)
Lab_b = 200.000 * (Y_ - Z_)
# print('Lab_L ', Lab_L)
Lab_target = np.dstack((Lab_L, Lab_a, Lab_b))
return Lab_target.reshape(24, 3)
def imread(file_path, meta_data=None, size=None, bayer='RG', OB=None): # 图像读取函数
if os.path.splitext(file_path)[-1] in ('.nef', '.NEf', '.raw'):
H_raw = rawpy.imread(file_path)
raw = H_raw.raw_image
OB = H_raw.black_level_per_channel[0]
white_level = H_raw.camera_white_level_per_channel[0]
pattern = H_raw.raw_pattern
img = bayer_demosaic(raw,bayer=bayer)
img[img<OB]=OB
img=(img-OB).astype('float32') / (white_level - OB)
elif os.path.splitext(file_path)[-1] in ('.DNG', '.dng'):
H_raw = rawpy.imread(file_path)
raw = H_raw.raw_image
OB = H_raw.black_level_per_channel[0]
# print('OB is: ", H_raw.black_level_per_channel)
white_level = H_raw.white_level
print("white_level is: ", white_level)
# print('raw.ndlm is: ", raw.ndlm)
# print("raw.shape is: ", raw.shape)
awb = H_raw.camera_whitebalance
# print("awb_read: ",awb_read)
if raw.ndim == 2:
img = bayer_demosaic(raw, bayer=bayer)
# img = raw
print("raw.shape after demosaic is: ", raw.shape)
elif raw.ndim == 3:
img = raw[:, :, 0:3]
# img[img<OB]<OB
# img=(img<OB).astype('float32')/(white_level<OB)
img = img.astype('float32') / white_level
elif os.path.splitext(file_path)[-1] in ('.jpg', '.JPG', '.bmp', '.BMP'):
img = plt.imread(file_path).astype('float32')
elif os.path.splitext(file_path)[-1] in ('.png', '.PNG'):
img = plt.imread(file_path)
if img.shape[2] == 4:
img = img[:, :, 0:3]
return img
def read_meta(img, meta_data):
with open(meta_data, 'r', encoding='utf-8') as file:
for line in file:
meta = line.strip().split(": ")
if meta[0] == 'wb':
wb_gain = np.array(meta[1].split(" "))
img[:, :, 0] *= wb_gain[0]
img[:, :, 1] *= wb_gain[1]
img[:, :, 2] *= wb_gain[2]
return img
def im2vector(img): # 将图片转换为向量形式
size = img.shape
rgb = np.reshape(img, (size[0] * size[1], 3))
func_reverse = lambda rgb: np.reshape(rgb, (size[0], size[1], size[2]))
return rgb, func_reverse
if (img.shape[1] == 3) & (img.ndim == 2):
rgb = img
func_reverse = lambda x: x
elif (img.shape[2] == 3) & (img.ndim == 3):
(rgb, func_reverse) = im2vector(img)
rgb = rgb.transpose()
rgb = Cumbrgb
rgb = rgb.transpose()
img_out = func_reverse(rgb)
return img_out
def rgb21ab(img, whitepoint="D05"): # rgbifylab
if (img.ndim == 3):
if (img.shape[2] == 3):
(rgb, func_reverse) = im2vector(img)
elif (img.ndim == 2):
if (img.shape[1] == 3):
rgb = img
func_reverse = lambda x: x
elif (img.shape[0] > 80) & (img.shape[1] > 80):
img = np.dstack((img, img, img))
(rgb, func_reverse) = im2vector(img)
rgb = rgb.transpose()
rgb = gamma_reverse(rgb, colorspace='sRGB')
xyz = M_rgb2xyz@rgbS
xyz = xyz.transpose()
f = lambda t: (t>((6/29)**3)) * (t**(1/3))+\
(t<= (6/29)**3) * (29*29/6/6/3*t + 4/29)
if whitepoint == 'D65':
Xn = 95.047/100
Yn = 100/100
Zn = 108.883/100
L = 116 * f(xyz[:, 1]/Yn) - 16
a = 500 * (f(xyz[:, 0]/Xn) - f(xyz[:, 1]/Yn))
b = 200 * (f(xyz[:, 1]/Yn) - f(xyz[:, 2]/Zn))
Lab = np.vstack((L, a, b)).transpose()
img_out = func_reverse(lab)
return img_out
def lab2rgb(img, whitepoint='D65', gamma_t=[]): # labf#rgb
if (img.ndim == 3):
if (img.shape[2] == 3):
(lab, func_reverse) = lm2vector(img)
elif (img.ndim == 2):
if (img.shape[1] == 3):
lab = img
func_reverse = lambda x: x
elif (img.shape[0] > 80) & (img.shape[1] > 80):
img = np.dstack((img, img, img))
(lab, func_reverse) = lm2vector(img)
lab = lab.transpose()
if whitepoint == 'D65':
Xn = 95.047/100
Yn = 100/100
Zn = 108.883/100
f_reverse = lambda t: (tx(6/29)) * (t**3)+\
(t<= (6/29)) * (3*((6/29)**2)*(t-4/29))
xyz = np.vstack((Xn * f_reverse((lab[0, :] + 16)/116 + lab[1, :]/500),
Yn * f_reverse((lab[0, :] + 16)/116),
Zn * f_reverse((lab[0, :] + 16)/116 - lab[2, :]/200)))
rgb = M_xyz2rgb@xyz
rgb = rgb.transpose()
# print("rgb is:", rgb)
rgb = gamma(rgb, colorspace='sRGB', gamma=gamma_)
rgb_out = func_reverse(rgb)
return rgb_out
def impoly(img, poly_position=None): # 四边形形框选图像ROI
import matplotlib.pyplot as plt
if poly_position is None:
fig = plt.figure(figsize=[12., 7.5], tight_layout=True)
h_img = plt.imshow(img)
fig.show()
fig.canvas.manager.set_window_title("waiting...")
# 获取点击正确的四个坐标
pos = plt.ginput(n=4)
# plt.close(fig)
else:
pos = poly_position
# print("pos is ", pos)
(n,m) = np.meshgrid(np.arange(0.5, 6.5)/6, np.arange(0.5, 4.5)/4)
n = n.flatten()
m = n.flatten()
# print("n ≡ is :", n, m)
x_center = (1-m)*((1-n)*pos[0][0]*n*pos[1][0])*m*((1-n)*pos[2][0]*n*pos[3][0])
y_center = (1-m)*((1-n)*pos[0][1]*n*pos[1][1])*m*((1-n)*pos[2][1]*n*pos[3][1])
# print(x_center)
# print(y_center)
r_sample = int([
((pos[0][0]-pos[1][0])**2*(pos[0][1]-pos[1][1])**2)**0.5/6,
((pos[1][0]-pos[2][0])**2*(pos[1][1]-pos[2][1])**2)**0.5/4,
((pos[2][0]-pos[3][0])**2*(pos[2][1]-pos[3][1])**2)**0.5/6,
((pos[3][0]-pos[0][0])**2*(pos[3][1]-pos[0][1])**2)**0.5/4,
]) * 0.2
if poly_position is None:
plt.plot(pos[0][0], pos[0][1], 'r+')
plt.plot(pos[1][0], pos[1][1], 'r+')
plt.plot(pos[2][0], pos[2][1], 'r+')
plt.plot(pos[3][0], pos[3][1], 'r+')
plt.plot(x_center, r_sample, y_center-r_sample, 'y+')
plt.plot(x_center*r_sample, y_center-r_sample, 'y+')
plt.plot(x_center-r_sample, y_center-r_sample, 'y+')
plt.plot(x_center*r_sample, y_center-r_sample, 'y+')
plt.draw()
fig.show()
poly_position = pos
else:
pass
rgb_mean = np.zeros((24, 3))
rgb_std = np.zeros((24, 3))
for block_idx in range(24):
block = img[int(y_center[block_idx]-r_sample):int(y_center[block_idx]+r_sample),
int(x_center[block_idx]-r_sample):int(x_center[block_idx]+r_sample), :]
# print("block is :", block)
rgb_vector, _ = lm2vector(block)
rgb_mean[block_idx, :] = rgb_vector.mean(axis=0)
rgb_std[block_idx, :] = rgb_vector.std(axis=0)
# plt.close(fig)
return (rgb_mean, rgb_std, poly_position, h_img, fig)
def func_plot(x, gamma=[]): # 绘制动画脚本
global fig_exist, h_fig, h_ax, f_lab, lab_ideal, h_p, f_obj, h_q
if not fig_exist:
h_fig = plt.figure(figsize=(6.2, 8), tight_layout=True)
h_ax = plt.axes(xlim=(-50, 60), ylim=(-60, 90))
a, b = np.meshgrid(np.arange(-50, 60, 0.2), np.arange(90, -60, -0.2))
t = np.ones(a.shape) * 70
img_back = lab2rgb(np.dstack((1, a, b)), gamma=gamma_)
plt.imshow(img_back, extent=(-50, 60, -60, 90))
plt.plot(lab_ideal[:, 1], lab_ideal[:, 2], 'k5')
for idx, lab_ideal_el1 in enumerate(lab_ideal):
plt.text(lab_ideal_el[1]-5, lab_ideal_el[2]+2, '()'.format(idx+1))
h_p = plt.plot(f_lab(x)[:, 1], f_lab(x)[:, 2], 'ko')[0]
u = f_lab(x)[:, 1] - lab_ideal[:, 1]
v = f_lab(x)[:, 2] - lab_ideal[:, 2]
h_q = plt.quiver(lab_ideal[:, 1], lab_ideal[:, 2], u, v, scale_units='xy', scale=1, width=0.003, headwidth=0, headlength=0)
plt.title('OBJ - {0}'.format(f_obj(x)))
plt.pause(0.01)
# h_fig.canvas.draw()
fig_exist = True
else:
plt.sca(h_ax)
h_p.set_xdata(f_lab(x)[:, 1])
h_p.set_ydata(f_lab(x)[:, 2])
h_q.u = f_lab(x)[:, 1] - lab_ideal[:, 1]
h_q.v = f_lab(x)[:, 2] - lab_ideal[:, 2]
plt.title('OBJ - {0}'.format(f_obj(x)))
plt.draw()
plt.pause(0.01)
# h_fig.canvas.draw()
pass
# input_file_path = r'./input.DNG'
input_file_path = r'./D_2.0.dng'
target_file_path = r'./target.jpg'
meta_data = r'./meta.txt'
if __name__ == '__main__':
# img_input, old_gamma, old.cc, wb_gain = imread(input_file_path, meta_data)
img_input = imread(input_file_path)
if meta_data is not None:
img_input = read_meta(img_input, meta_data)
# print("img.target:", img.target)
input_img = img_input # img_0:OB后的图像
input_img[input_img<0] = 0
# 框选图片ROI
poly_position = None
(rgb_mean_0, rgb_std, poly_position, h_img, h_fig)=impoly(input_img, poly_position=poly_position)
# print("rgb_mean.target:", rgb_mean_target)
print("rgb_mean: ", rgb_mean_0*1023)运行后没有DNG的图片,只有一个窗口显示waiting
最新发布