import numpy as np
import pandas as pd
from scipy.optimize import minimize_scalar
from scipy.spatial import ConvexHull
# -------------------------- Basic Function Definitions --------------------------
def trapezoidal_integral(f, lambda_nm):
"""
Trapezoidal integration for spectral data
:param f: Integrand array
:param lambda_nm: Wavelength array
:return: Integration result
"""
return np.sum((f[:-1] + f[1:]) / 2 * np.diff(lambda_nm))
def compute_xyz(spd, x_bar, y_bar, z_bar, lambda_nm):
"""
Calculate CIE XYZ tristimulus values
:param spd: Spectral power distribution
:param x_bar, y_bar, z_bar: CIE 1931 standard observer functions
:param lambda_nm: Wavelength array
:return: X, Y, Z tristimulus values
"""
X = trapezoidal_integral(spd * x_bar, lambda_nm)
Y = trapezoidal_integral(spd * y_bar, lambda_nm)
Z = trapezoidal_integral(spd * z_bar, lambda_nm)
return X, Y, Z
def planckian_locus_xy(T):
"""
Calculate chromaticity coordinates (x_b, y_b) on Planckian locus for given color temperature T
:param T: Color temperature (K)
:return: x_b, y_b chromaticity coordinates on Planckian locus
"""
if T < 1000 or T > 25000:
raise ValueError("Color temperature should be between 1000-25000K")
# Calculate x_b (Planckian locus x coordinate)
if T <= 4000:
x_b = -0.2661239e9 / T**3 - 0.2343580e6 / T**2 + 0.8776956e3 / T + 0.179910
else:
x_b = -3.0258469e9 / T**3 + 2.1070379e6 / T**2 + 0.2226347e3 / T + 0.240390
# Calculate y_b (Planckian locus y coordinate)
if T <= 2222:
y_b = -1.1063814 * x_b**3 - 1.34811020 * x_b**2 + 2.18555832 * x_b - 0.20219683
elif T <= 4000:
y_b = -0.9549476 * x_b**3 - 1.37418593 * x_b**2 + 2.09137015 * x_b - 0.16748867
else:
y_b = 3.0817580 * x_b**3 - 5.87338670 * x_b**2 + 3.75112997 * x_b - 0.37001483
return x_b, y_b
def ciede2000_delta_e(lab_sample, lab_standard):
"""
Calculate CIEDE2000 color difference
:param lab_sample: LAB coordinates of test color sample under test light source [L*, a*, b*]
:param lab_standard: LAB coordinates of test color sample under standard light source [L*0, a*0, b*0]
:return: delta_E00 color difference
"""
L, a, b = lab_sample
L0, a0, b0 = lab_standard
# Lightness difference
delta_L = L - L0
# Chroma calculation
C = np.sqrt(a**2 + b**2)
C0 = np.sqrt(a0**2 + b0**2)
C_avg = (C + C0) / 2
# Chroma difference
delta_C = C - C0
# Hue angle calculation
h = np.degrees(np.arctan2(b, a)) if C != 0 else 0
h0 = np.degrees(np.arctan2(b0, a0)) if C0 != 0 else 0
# Hue difference
delta_h = h - h0
if C * C0 == 0:
delta_H = 0
else:
if np.abs(delta_h) <= 180:
delta_H = delta_h
else:
delta_H = delta_h - 360 if delta_h > 0 else delta_h + 360
# Hue difference in polar coordinates
delta_H = 2 * np.sqrt(C * C0) * np.sin(np.radians(delta_H) / 2)
# Average hue angle
h_avg = (h + h0) / 2 if np.abs(h - h0) <= 180 else (h + h0 + 360) / 2 if h + h0 < 360 else (h + h0 - 360) / 2
# Calculate T for hue rotation term
T = 1 - 0.17 * np.cos(np.radians(h_avg - 30)) + 0.24 * np.cos(np.radians(2 * h_avg)) + \
0.32 * np.cos(np.radians(3 * h_avg + 6)) - 0.20 * np.cos(np.radians(4 * h_avg - 63))
# Rotation term
delta_theta = 30 * np.exp(-(((h_avg - 275) / 25) ** 2))
R_C = 2 * np.sqrt((C_avg ** 7) / (C_avg ** 7 + 25 ** 7))
R_T = -R_C * np.sin(np.radians(2 * delta_theta))
# Calculate parametric factors
L_avg = (L + L0) / 2
S_L = 1 + (0.015 * (L_avg - 50) ** 2) / np.sqrt(20 + (L_avg - 50) ** 2)
S_C = 1 + 0.045 * C_avg
S_H = 1 + 0.015 * C_avg
# Calculate delta_E00
delta_E00 = np.sqrt(
(delta_L / (1 * S_L)) ** 2 +
(delta_C / (1 * S_C)) ** 2 +
(delta_H / (1 * S_H)) ** 2 +
R_T * (delta_C / S_C) * (delta_H / S_H)
)
return delta_E00
def xy_to_uv_prime(x, y):
"""
Convert xy chromaticity coordinates to CIE u'v' chromaticity coordinates
:param x, y: xy chromaticity coordinates
:return: u', v' chromaticity coordinates
"""
denominator = -2 * x + 12 * y + 3 # Corrected denominator formula
u_prime = 4 * x / denominator
v_prime = 9 * y / denominator
return u_prime, v_prime
def shoelace_area(points):
"""
Calculate polygon area using shoelace formula
:param points: Polygon vertex coordinate array (n, 2)
:return: Polygon area
"""
points = np.vstack([points, points[0]]) # Close the polygon
return 0.5 * np.abs(np.sum(points[:-1, 0] * points[1:, 1] - points[1:, 0] * points[:-1, 1]))
# -------------------------- Main Calculation Program --------------------------
def main():
"""Main program: Read data and calculate all spectral characteristic parameters"""
# -------------------------- Data Reading --------------------------
# Read Excel data (assuming all data is in different sheets of the same Excel file)
# In practical applications, modify the file path and sheet names according to data storage location
excel_file = "Problem 1.xlsx"
try:
# Read LED spectral power distribution (SPD)
led_spd_df = pd.read_excel(excel_file, sheet_name="LED_SPD")
lambda_nm = led_spd_df["Wavelength"].values # Wavelength array (380-780nm)
spd_led = led_spd_df["SPD"].values # LED spectral power distribution
# Read CIE 1931 standard observer functions
cie_xyz_df = pd.read_excel(excel_file, sheet_name="CIE_1931_XYZ")
x_bar = cie_xyz_df["x_bar"].values # xbar function
y_bar = cie_xyz_df["y_bar"].values # ybar function
z_bar = cie_xyz_df["z_bar"].values # zbar function
# Read 24 CIE test color sample reflectances
color_samples_df = pd.read_excel(excel_file, sheet_name="CIE_24_Colors")
rho = color_samples_df.iloc[:, 1:25].values # Reflectance array (401x24)
# Read standard light source D65 SPD
d65_spd_df = pd.read_excel(excel_file, sheet_name="D65_SPD")
spd_d65 = d65_spd_df["SPD"].values # D65 spectral power distribution
# Read melanopic efficiency function
melanopic_df = pd.read_excel(excel_file, sheet_name="Melanopic_S")
s_mel = melanopic_df["S"].values # melanopic efficiency function
except FileNotFoundError:
print(f"Error: Excel file '{excel_file}' not found.")
return
except KeyError as e:
print(f"Error: Missing column in Excel file - {e}")
return
# Validate data arrays have the same length
data_arrays = [lambda_nm, spd_led, x_bar, y_bar, z_bar, spd_d65, s_mel]
if len(set(len(arr) for arr in data_arrays)) != 1:
print("Error: All spectral data arrays must have the same length")
return
# -------------------------- Color Characteristic Parameters Calculation (CCT, Duv) --------------------------
# Calculate LED's XYZ tristimulus values and chromaticity coordinates
X_led, Y_led, Z_led = compute_xyz(spd_led, x_bar, y_bar, z_bar, lambda_nm)
sum_xyz_led = X_led + Y_led + Z_led
if sum_xyz_led == 0:
print("Error: XYZ sum is zero, cannot calculate chromaticity coordinates")
return
x_led = X_led / sum_xyz_led
y_led = Y_led / sum_xyz_led
# Convert to u', v' coordinates for more accurate CCT calculation
u_led, v_led = xy_to_uv_prime(x_led, y_led)
# Function to calculate distance to Planckian locus in u', v' coordinates
def distance_to_locus(T):
x_b, y_b = planckian_locus_xy(T)
u_b, v_b = xy_to_uv_prime(x_b, y_b)
return (u_led - u_b)**2 + (v_led - v_b)** 2
# Find CCT by minimizing distance to Planckian locus
try:
result = minimize_scalar(distance_to_locus, bounds=(1000, 25000), method='bounded')
cct = result.x # Correlated color temperature
# Calculate Duv (distance from Planckian locus)
x_b, y_b = planckian_locus_xy(cct)
u_b, v_b = xy_to_uv_prime(x_b, y_b)
duv = (v_led - v_b) / 0.0014 # 0.0014 is the scaling factor for Duv
except ValueError as e:
print(f"Error calculating CCT: {e}")
cct, duv = None, None
# -------------------------- Color Rendering Parameters Calculation (Rf, Rg) --------------------------
n_samples = 24 # Number of CIE test color samples
xyz_led_samples = np.zeros((n_samples, 3)) # XYZ values of color samples under LED
xyz_d65_samples = np.zeros((n_samples, 3)) # XYZ values of color samples under D65
# Calculate XYZ tristimulus values for all color samples under LED and D65
for i in range(n_samples):
rho_i = rho[:, i] # Reflectance of the i-th color sample
# XYZ under LED
spd_led_i = spd_led * rho_i
X, Y, Z = compute_xyz(spd_led_i, x_bar, y_bar, z_bar, lambda_nm)
xyz_led_samples[i] = (X, Y, Z)
# XYZ under D65
spd_d65_i = spd_d65 * rho_i
X, Y, Z = compute_xyz(spd_d65_i, x_bar, y_bar, z_bar, lambda_nm)
xyz_d65_samples[i] = (X, Y, Z)
# Calculate XYZ tristimulus values of D65 (white reference)
X_d65, Y_d65, Z_d65 = compute_xyz(spd_d65, x_bar, y_bar, z_bar, lambda_nm)
sum_xyz_d65 = X_d65 + Y_d65 + Z_d65
# Calculate LAB coordinates for all color samples
def xyz_to_lab(xyz, xyz_ref):
"""XYZ to LAB color space conversion"""
x, y, z = xyz
x_ref, y_ref, z_ref = xyz_ref
# Avoid division by zero
if x_ref == 0 or y_ref == 0 or z_ref == 0:
return [0, 0, 0]
x /= x_ref
y /= y_ref
z /= z_ref
# Non-linear transformation function
def f(t):
return t**(1/3) if t > (6/29)**3 else (t * (29/6)**2)/3 + 4/29
L = 116 * f(y) - 16
a = 500 * (f(x) - f(y))
b = 200 * (f(y) - f(z))
return np.array([L, a, b])
# Calculate LAB values using respective white references
lab_led = np.array([xyz_to_lab(xyz, (X_led, Y_led, Z_led)) for xyz in xyz_led_samples])
lab_d65 = np.array([xyz_to_lab(xyz, (X_d65, Y_d65, Z_d65)) for xyz in xyz_d65_samples])
# Calculate fidelity index Rf
rf = None
try:
delta_e00_list = [ciede2000_delta_e(lab_led[i], lab_d65[i]) for i in range(n_samples)]
avg_delta_e00 = np.mean(delta_e00_list)
rf = 100 - 4.6 * avg_delta_e00 # ANSI/IES TM-30-15 formula
rf = max(0, min(100, rf)) # Ensure Rf is within 0-100 range
except Exception as e:
print(f"Error calculating Rf: {e}")
# Calculate gamut index Rg
rg = None
try:
def get_uv_prime_points(xyz_samples):
"""Convert XYZ color samples to u'v' coordinate array"""
uv_points = []
for xyz in xyz_samples:
X, Y, Z = xyz
sum_xyz = X + Y + Z
if sum_xyz == 0:
uv_points.append([0, 0])
continue
x = X / sum_xyz
y = Y / sum_xyz
u, v = xy_to_uv_prime(x, y)
uv_points.append([u, v])
return np.array(uv_points)
# Calculate gamut areas using convex hulls for LED and D65
uv_led = get_uv_prime_points(xyz_led_samples)
uv_d65 = get_uv_prime_points(xyz_d65_samples)
# Compute convex hulls to get the gamut boundaries
hull_led = ConvexHull(uv_led)
hull_d65 = ConvexHull(uv_d65)
# Calculate areas using shoelace formula on convex hull vertices
area_led = shoelace_area(uv_led[hull_led.vertices])
area_d65 = shoelace_area(uv_d65[hull_d65.vertices])
rg = (area_led / area_d65) * 100 if area_d65 != 0 else 0
except Exception as e:
print(f"Error calculating Rg: {e}")
# -------------------------- Circadian Rhythm Parameter Calculation (mel-DER) --------------------------
mel_der = None
try:
# Calculate melanopic irradiance for LED
e_mel_led = trapezoidal_integral(spd_led * s_mel, lambda_nm)
# Calculate melanopic irradiance for D65 for comparison
e_mel_d65 = trapezoidal_integral(spd_d65 * s_mel, lambda_nm)
# melanopic Daylight Equivalent Ratio
mel_der = e_mel_led / e_mel_d65 if e_mel_d65 != 0 else 0
except Exception as e:
print(f"Error calculating mel-DER: {e}")
# -------------------------- Result Output --------------------------
print("LED Light Source Spectral Characteristic Parameters Calculation Results:")
print(f"1. Correlated Color Temperature (CCT): {cct:.1f} K" if cct is not None else "1. Correlated Color Temperature (CCT): Error")
print(f"2. Distance from Planckian Locus (Duv): {duv:.4f}" if duv is not None else "2. Distance from Planckian Locus (Duv): Error")
print(f"3. Fidelity Index (Rf): {rf:.1f}" if rf is not None else "3. Fidelity Index (Rf): Error")
print(f"4. Gamut Index (Rg): {rg:.1f} %" if rg is not None else "4. Gamut Index (Rg): Error")
print(f"5. Melatonin Daylight Illuminance Ratio (mel-DER): {mel_der:.3f}" if mel_der is not None else "5. Melatonin Daylight Illuminance Ratio (mel-DER): Error")
# -------------------------- Program Entry --------------------------
if __name__ == "__main__":
main() 能把他转化为MATLAB的代码并加上注释
最新发布
本文提供了一种高效的算法来输出指定范围内的第k个素数，并优化了时间复杂度至50ms内完成计算，通过筛法优化避免重复计算，适用于大规模数据处理。
扫一扫
331
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
