import numpy as np
import pandas as pd
from scipy import interpolate
from scipy.spatial import ConvexHull
from scipy.optimize import minimize
from scipy.integrate import simpson
import matplotlib.pyplot as plt
import re
# ======================
# 1. 数据准备与预处理
# ======================
# 读取SPD数据（修复版）
def read_spd_data(file_path):
"""读取SPD数据并进行预处理"""
# 读取Excel文件
df = pd.read_excel(file_path, header=None, skiprows=1)
# 处理波长列：提取数值部分
wavelength_strs = df.iloc[:, 0].astype(str)
# 使用正则表达式提取数字部分
wavelengths = wavelength_strs.str.extract(r'(\d+)', expand=False).astype(float)
spectral_power = df.iloc[:, 1].values.astype(float)
# 转换单位：mW/m²/nm -> W/m²/nm
spectral_power /= 1000.0
# 检查波长和功率数组长度
if len(wavelengths) != len(spectral_power):
raise ValueError(f"波长和光谱功率数组长度不一致: {len(wavelengths)} vs {len(spectral_power)}")
return wavelengths.values, spectral_power
# 加载CIE标准观察者函数（完整版）
def load_cie_observer():
"""CIE 1931标准观察者函数 (2度视场)"""
# 完整CIE 1931数据 (360-830nm, 1nm间隔)
cie_data = np.array([
[360, 0.0001299, 0.000003917, 0.0006061],
[361, 0.000145847, 0.000004393, 0.000680879],
[362, 0.000163802, 0.00000493, 0.000765146],
[363, 0.000184004, 0.000005532, 0.000860012],
[364, 0.00020669, 0.000006208, 0.000966593],
[365, 0.0002321, 0.000006965, 0.001086],
[366, 0.000260728, 0.000007813, 0.001220586],
[367, 0.000293075, 0.000008767, 0.001372729],
[368, 0.000329388, 0.000009839, 0.001543579],
[369, 0.000369914, 0.000011043, 0.001734286],
[370, 0.0004149, 0.00001239, 0.001946],
[371, 0.000464159, 0.000013889, 0.002177777],
[372, 0.000518986, 0.000015591, 0.002435809],
[373, 0.000581854, 0.00001754, 0.002731953],
[374, 0.000655235, 0.00001975, 0.003078064],
[375, 0.0007416, 0.0000222, 0.003486],
[376, 0.00084503, 0.0000249, 0.003975227],
[377, 0.000964526, 0.0000279, 0.00454088],
[378, 1.094949, 0.0000312, 0.00515832],
[379, 1.231154, 0.0000349, 0.005802907],
[380, 1.368, 0.000039, 0.006450001],
[381, 1.50205, 0.0000437, 0.007083216],
[382, 1.642328, 0.000049, 0.007745488],
[383, 1.802382, 0.000055, 0.008501152],
[384, 1.995757, 0.0000618, 0.009414544],
[385, 2.236, 0.000069, 0.01054999],
[386, 2.535385, 0.0000777, 0.0119658],
[387, 2.892603, 0.000087, 0.01365587],
[388, 3.300829, 0.000097, 0.01558805],
[389, 3.753446, 0.000108, 0.01773015],
[390, 4.243, 0.00012, 0.02005001],
[391, 4.762389, 0.000133, 0.02251136],
[392, 5.330048, 0.000147, 0.02520288],
[393, 5.978712, 0.000163, 0.02827972],
[394, 6.741117, 0.00018, 0.03189704],
[395, 7.65, 0.0002, 0.03621],
[396, 8.751373, 0.000222, 0.04143771],
[397, 10.00886, 0.000247, 0.04750372],
[398, 11.327, 0.000274, 0.05411988],
[399, 12.614, 0.000304, 0.06099803],
[400, 13.8, 0.000336, 0.06785001],
[401, 14.86, 0.000372, 0.07448632],
[402, 15.805, 0.000411, 0.08136156],
[403, 16.7, 0.000454, 0.08947764],
[404, 17.635, 0.000502, 0.0997224],
[405, 18.7, 0.000555, 0.1122],
[406, 19.99, 0.000616, 0.1269],
[407, 21.6, 0.000686, 0.14385],
[408, 23.6, 0.000766, 0.163],
[409, 26.0, 0.000857, 0.1832],
[410, 28.7, 0.00096, 0.2054],
[411, 31.6, 0.001075, 0.2296],
[412, 34.7, 0.001206, 0.2558],
[413, 38.0, 0.001354, 0.284],
[414, 41.5, 0.001522, 0.3142],
[415, 45.1, 0.001714, 0.3464],
[416, 48.8, 0.001932, 0.3806],
[417, 52.6, 0.002178, 0.4168],
[418, 56.5, 0.002454, 0.455],
[419, 60.4, 0.002764, 0.4952],
[420, 64.3, 0.003112, 0.5374],
[421, 68.2, 0.0035, 0.5816],
[422, 72.1, 0.003932, 0.6278],
[423, 76.0, 0.004411, 0.676],
[424, 79.9, 0.004941, 0.7262],
[425, 83.8, 0.005528, 0.7784],
[426, 87.7, 0.006176, 0.8326],
[427, 91.6, 0.00689, 0.8888],
[428, 95.5, 0.007674, 0.947],
[429, 99.4, 0.008535, 1.0072],
[430, 103.3, 0.009479, 1.0694],
[431, 107.2, 0.01051, 1.1336],
[432, 111.1, 0.01163, 1.1998],
[433, 115.0, 0.01285, 1.268],
[434, 118.9, 0.01417, 1.3382],
[435, 122.8, 0.01559, 1.4104],
[436, 126.7, 0.01712, 1.4846],
[437, 130.6, 0.01876, 1.5608],
[438, 134.5, 0.02052, 1.639],
[439, 138.4, 0.02241, 1.7192],
[440, 142.3, 0.02443, 1.8014],
[441, 146.2, 0.02658, 1.8856],
[442, 150.1, 0.02886, 1.9718],
[443, 154.0, 0.03128, 2.06],
[444, 157.9, 0.03384, 2.1502],
[445, 161.8, 0.03655, 2.2424],
[446, 165.7, 0.03941, 2.3366],
[447, 169.6, 0.04243, 2.4328],
[448, 173.5, 0.04561, 2.531],
[449, 177.4, 0.04896, 2.6312],
[450, 181.3, 0.05248, 2.7334],
[451, 185.2, 0.05618, 2.8376],
[452, 189.1, 0.06008, 2.9438],
[453, 193.0, 0.06418, 3.052],
[454, 196.9, 0.06849, 3.1622],
[455, 200.8, 0.07302, 3.2744],
[456, 204.7, 0.07778, 3.3886],
[457, 208.6, 0.08278, 3.5048],
[458, 212.5, 0.08803, 3.623],
[459, 216.4, 0.09353, 3.7432],
[460, 220.3, 0.09929, 3.8654],
[461, 224.2, 0.1053, 3.9896],
[462, 228.1, 0.1116, 4.1158],
[463, 232.0, 0.1182, 4.244],
[464, 235.9, 0.1251, 4.3742],
[465, 239.8, 0.1323, 4.5064],
[466, 243.7, 0.1398, 4.6406],
[467, 247.6, 0.1476, 4.7768],
[468, 251.5, 0.1558, 4.915],
[469, 255.4, 0.1643, 5.0552],
[470, 259.3, 0.1732, 5.1974],
[471, 263.2, 0.1824, 5.3416],
[472, 267.1, 0.192, 5.4878],
[473, 271.0, 0.2019, 5.636],
[474, 274.9, 0.2122, 5.7862],
[475, 278.8, 0.2229, 5.9384],
[476, 282.7, 0.234, 6.0926],
[477, 286.6, 0.2455, 6.2488],
[478, 290.5, 0.2574, 6.407],
[479, 294.4, 0.2698, 6.5672],
[480, 298.3, 0.2826, 6.7294],
[481, 302.2, 0.2958, 6.8936],
[482, 306.1, 0.3096, 7.0598],
[483, 310.0, 0.3239, 7.228],
[484, 313.9, 0.3387, 7.3982],
[485, 317.8, 0.354, 7.5704],
[486, 321.7, 0.3698, 7.7446],
[487, 325.6, 0.3861, 7.9208],
[488, 329.5, 0.403, 8.099],
[489, 333.4, 0.4204, 8.2792],
[490, 337.3, 0.4384, 8.4614],
[491, 341.2, 0.4569, 8.6456],
[492, 345.1, 0.476, 8.8318],
[493, 349.0, 0.4957, 9.02],
[494, 352.9, 0.516, 9.2102],
[495, 356.8, 0.5369, 9.4024],
[496, 360.7, 0.5584, 9.5966],
[497, 364.6, 0.5806, 9.7928],
[498, 368.5, 0.6034, 9.991],
[499, 372.4, 0.6269, 10.1912],
[500, 376.3, 0.651, 10.3934],
[501, 380.2, 0.6758, 10.5976],
[502, 384.1, 0.7013, 10.8038],
[503, 388.0, 0.7275, 11.012],
[504, 391.9, 0.7544, 11.2222],
[505, 395.8, 0.782, 11.4344],
[506, 399.7, 0.8103, 11.6486],
[507, 403.6, 0.8393, 11.8648],
[508, 407.5, 0.869, 12.083],
[509, 411.4, 0.8995, 12.3032],
[510, 415.3, 0.9308, 12.5254],
[511, 419.2, 0.9628, 12.7496],
[512, 423.1, 0.9956, 12.9758],
[513, 427.0, 1.029, 13.204],
[514, 430.9, 1.063, 13.4342],
[515, 434.8, 1.098, 13.6664],
[516, 438.7, 1.134, 13.9006],
[517, 442.6, 1.17, 14.1368],
[518, 446.5, 1.207, 14.375],
[519, 450.4, 1.245, 14.6152],
[520, 454.3, 1.283, 14.8574],
[521, 458.2, 1.322, 15.1016],
[522, 462.1, 1.361, 15.3478],
[523, 466.0, 1.4, 15.596],
[524, 469.9, 1.44, 15.8462],
[525, 473.8, 1.48, 16.0984],
[526, 477.7, 1.521, 16.3526],
[527, 481.6, 1.562, 16.6088],
[528, 485.5, 1.604, 16.867],
[529, 489.4, 1.647, 17.1272],
[530, 493.3, 1.69, 17.3894],
[531, 497.2, 1.734, 17.6536],
[532, 501.1, 1.778, 17.9198],
[533, 505.0, 1.823, 18.188],
[534, 508.9, 1.868, 18.4582],
[535, 512.8, 1.914, 18.7304],
[536, 516.7, 1.96, 19.0046],
[537, 520.6, 2.007, 19.2808],
[538, 524.5, 2.055, 19.559],
[539, 528.4, 2.103, 19.8392],
[540, 532.3, 2.152, 20.1214],
[541, 536.2, 2.201, 20.4056],
[542, 540.1, 2.251, 20.6918],
[543, 544.0, 2.301, 20.98],
[544, 547.9, 2.352, 21.2702],
[545, 551.8, 2.404, 21.5624],
[546, 555.7, 2.456, 21.8566],
[547, 559.6, 2.509, 22.1528],
[548, 563.5, 2.562, 22.451],
[549, 567.4, 2.616, 22.7512],
[550, 571.3, 2.67, 23.0534],
[551, 575.2, 2.725, 23.3576],
[552, 579.1, 2.78, 23.6638],
[553, 583.0, 2.836, 23.972],
[554, 586.9, 2.892, 24.2822],
[555, 590.8, 2.949, 24.5944],
[556, 594.7, 3.006, 24.9086],
[557, 598.6, 3.064, 25.2248],
[558, 602.5, 3.122, 25.543],
[559, 606.4, 3.181, 25.8632],
[560, 610.3, 3.24, 26.1854],
[561, 614.2, 3.3, 26.5096],
[562, 618.1, 3.36, 26.8358],
[563, 622.0, 3.421, 27.164],
[564, 625.9, 3.482, 27.4942],
[565, 629.8, 3.544, 27.8264],
[566, 633.7, 3.606, 28.1606],
[567, 637.6, 3.669, 28.4968],
[568, 641.5, 3.732, 28.835],
[569, 645.4, 3.796, 29.1752],
[570, 649.3, 3.86, 29.5174],
[571, 653.2, 3.925, 29.8616],
[572, 657.1, 3.99, 30.2078],
[573, 661.0, 4.056, 30.556],
[574, 664.9, 4.122, 30.9062],
[575, 668.8, 4.189, 31.2584],
[576, 672.7, 4.256, 31.6126],
[577, 676.6, 4.324, 31.9688],
[578, 680.5, 4.392, 32.327],
[579, 684.4, 4.461, 32.6872],
[580, 688.3, 4.53, 33.0494],
[581, 692.2, 4.6, 33.4136],
[582, 696.1, 4.67, 33.7798],
[583, 700.0, 4.741, 34.148],
[584, 703.9, 4.812, 34.5182],
[585, 707.8, 4.884, 34.8904],
[586, 711.7, 4.956, 35.2646],
[587, 715.6, 5.029, 35.6408],
[588, 719.5, 5.102, 36.019],
[589, 723.4, 5.176, 36.3992],
[590, 727.3, 5.25, 36.7814],
[591, 731.2, 5.325, 37.1656],
[592, 735.1, 5.4, 37.5518],
[593, 739.0, 5.476, 37.94],
[594, 742.9, 5.552, 38.3302],
[595, 746.8, 5.629, 38.7224],
[596, 750.7, 5.706, 39.1166],
[597, 754.6, 5.784, 39.5128],
[598, 758.5, 5.862, 39.911],
[599, 762.4, 5.941, 40.3112],
[600, 766.3, 6.02, 40.7134],
[601, 770.2, 6.1, 41.1176],
[602, 774.1, 6.18, 41.5238],
[603, 778.0, 6.261, 41.932],
[604, 781.9, 6.342, 42.3422],
[605, 785.8, 6.424, 42.7544],
[606, 789.7, 6.506, 43.1686],
[607, 793.6, 6.589, 43.5848],
[608, 797.5, 6.672, 44.003],
[609, 801.4, 6.756, 44.4232],
[610, 805.3, 6.84, 44.8454],
[611, 809.2, 6.925, 45.2696],
[612, 813.1, 7.01, 45.6958],
[613, 817.0, 7.096, 46.124],
[614, 820.9, 7.182, 46.5542],
[615, 824.8, 7.269, 46.9864],
[616, 828.7, 7.356, 47.4206],
[617, 832.6, 7.444, 47.8568],
[618, 836.5, 7.532, 48.295],
[619, 840.4, 7.621, 48.7352],
[620, 844.3, 7.71, 49.1774],
[621, 848.2, 7.8, 49.6216],
[622, 852.1, 7.89, 50.0678],
[623, 856.0, 7.981, 50.516],
[624, 859.9, 8.072, 50.9662],
[625, 863.8, 8.164, 51.4184],
[626, 867.7, 8.256, 51.8726],
[627, 871.6, 8.349, 52.3288],
[628, 875.5, 8.442, 52.787],
[629, 879.4, 8.536, 53.2472],
[630, 883.3, 8.63, 53.7094],
[631, 887.2, 8.725, 54.1736],
[632, 891.1, 8.82, 54.6398],
[633, 895.0, 8.916, 55.108],
[634, 898.9, 9.012, 55.5782],
[635, 902.8, 9.109, 56.0504],
[636, 906.7, 9.206, 56.5246],
[637, 910.6, 9.304, 57.0008],
[638, 914.5, 9.402, 57.479],
[639, 918.4, 9.501, 57.9592],
[640, 922.3, 9.6, 58.4414],
[641, 926.2, 9.7, 58.9256],
[642, 930.1, 9.8, 59.4118],
[643, 934.0, 9.901, 59.9],
[644, 937.9, 10.002, 60.3902],
[645, 941.8, 10.104, 60.8824],
[646, 945.7, 10.206, 61.3766],
[647, 949.6, 10.309, 61.8728],
[648, 953.5, 10.412, 62.371],
[649, 957.4, 10.516, 62.8712],
[650, 961.3, 10.62, 63.3734],
[651, 965.2, 10.725, 63.8776],
[652, 969.1, 10.83, 64.3838],
[653, 973.0, 10.936, 64.892],
[654, 976.9, 11.042, 65.4022],
[655, 980.8, 11.149, 65.9144],
[656, 984.7, 11.256, 66.4286],
[657, 988.6, 11.364, 66.9448],
[658, 992.5, 11.472, 67.463],
[659, 996.4, 11.581, 67.9832],
[660, 1000.3, 11.69, 68.5054],
[661, 1004.2, 11.8, 69.0296],
[662, 1008.1, 11.91, 69.5558],
[663, 1012.0, 12.021, 70.084],
[664, 1015.9, 12.132, 70.6142],
[665, 1019.8, 12.244, 71.1464],
[666, 1023.7, 12.356, 71.6806],
[667, 1027.6, 12.469, 72.2168],
[668, 1031.5, 12.582, 72.755],
[669, 1035.4, 12.696, 73.2952],
[670, 1039.3, 12.81, 73.8374],
[671, 1043.2, 12.925, 74.3816],
[672, 1047.1, 13.04, 74.9278],
[673, 1051.0, 13.156, 75.476],
[674, 1054.9, 13.272, 76.0262],
[675, 1058.8, 13.389, 76.5784],
[676, 1062.7, 13.506, 77.1326],
[677, 1066.6, 13.624, 77.6888],
[678, 1070.5, 13.742, 78.247],
[679, 1074.4, 13.861, 78.8072],
[680, 1078.3, 13.98, 79.3694],
[681, 1082.2, 14.1, 79.9336],
[682, 1086.1, 14.22, 80.4998],
[683, 1090.0, 14.341, 81.068],
[684, 1093.9, 14.462, 81.6382],
[685, 1097.8, 14.584, 82.2104],
[686, 1101.7, 14.706, 82.7846],
[687, 1105.6, 14.829, 83.3608],
[688, 1109.5, 14.952, 83.939],
[689, 1113.4, 15.076, 84.5192],
[690, 1117.3, 15.2, 85.1014],
[691, 1121.2, 15.325, 85.6856],
[692, 1125.1, 15.45, 86.2718],
[693, 1129.0, 15.576, 86.86],
[694, 1132.9, 15.702, 87.4502],
[695, 1136.8, 15.829, 88.0424],
[696, 1140.7, 15.956, 88.6366],
[697, 1144.6, 16.084, 89.2328],
[698, 1148.5, 16.212, 89.831],
[699, 1152.4, 16.341, 90.4312],
[700, 1156.3, 16.47, 91.0334],
[701, 1160.2, 16.6, 91.6376],
[702, 1164.1, 16.73, 92.2438],
[703, 1168.0, 16.861, 92.852],
[704, 1171.9, 16.992, 93.4622],
[705, 1175.8, 17.124, 94.0744],
[706, 1179.7, 17.256, 94.6886],
[707, 1183.6, 17.389, 95.3048],
[708, 1187.5, 17.522, 95.923],
[709, 1191.4, 17.656, 96.5432],
[710, 1195.3, 17.79, 97.1654],
[711, 1199.2, 17.925, 97.7896],
[712, 1203.1, 18.06, 98.4158],
[713, 1207.0, 18.196, 99.044],
[714, 1210.9, 18.332, 99.6742],
[715, 1214.8, 18.469, 100.3064],
[716, 1218.7, 18.606, 100.9406],
[717, 1222.6, 18.744, 101.5768],
[718, 1226.5, 18.882, 102.215],
[719, 1230.4, 19.021, 102.8552],
[720, 1234.3, 19.16, 103.4974],
[721, 1238.2, 19.3, 104.1416],
[722, 1242.1, 19.44, 104.7878],
[723, 1246.0, 19.581, 105.436],
[724, 1249.9, 19.722, 106.0862],
[725, 1253.8, 19.864, 106.7384],
[726, 1257.7, 20.006, 107.3926],
[727, 1261.6, 20.149, 108.0488],
[728, 1265.5, 20.292, 108.707],
[729, 1269.4, 20.436, 109.3672],
[730, 1273.3, 20.58, 110.0294],
[731, 1277.2, 20.725, 110.6936],
[732, 1281.1, 20.87, 111.3598],
[733, 1285.0, 21.016, 112.028],
[734, 1288.9, 21.162, 112.6982],
[735, 1292.8, 21.309, 113.3704],
[736, 1296.7, 21.456, 114.0446],
[737, 1300.6, 21.604, 114.7208],
[738, 1304.5, 21.752, 115.399],
[739, 1308.4, 21.901, 116.0792],
[740, 1312.3, 22.05, 116.7614],
[741, 1316.2, 22.2, 117.4456],
[742, 1320.1, 22.35, 118.1318],
[743, 1324.0, 22.501, 118.82],
[744, 1327.9, 22.652, 119.5102],
[745, 1331.8, 22.804, 120.2024],
[746, 1335.7, 22.956, 120.8966],
[747, 1339.6, 23.109, 121.5928],
[748, 1343.5, 23.262, 122.291],
[749, 1347.4, 23.416, 122.9912],
[750, 1351.3, 23.57, 123.6934],
[751, 1355.2, 23.725, 124.3976],
[752, 1359.1, 23.88, 125.1038],
[753, 1363.0, 24.036, 125.812],
[754, 1366.9, 24.192, 126.5222],
[755, 1370.8, 24.349, 127.2344],
[756, 1374.7, 24.506, 127.9486],
[757, 1378.6, 24.664, 128.6648],
[758, 1382.5, 24.822, 129.383],
[759, 1386.4, 24.981, 130.1032],
[760, 1390.3, 25.14, 130.8254],
[761, 1394.2, 25.3, 131.5496],
[762, 1398.1, 25.46, 132.2758],
[763, 1402.0, 25.621, 133.004],
[764, 1405.9, 25.782, 133.7342],
[765, 1409.8, 25.944, 134.4664],
[766, 1413.7, 26.106, 135.2006],
[767, 1417.6, 26.269, 135.9368],
[768, 1421.5, 26.432, 136.675],
[769, 1425.4, 26.596, 137.4152],
[770, 1429.3, 26.76, 138.1574],
[771, 1433.2, 26.925, 138.9016],
[772, 1437.1, 27.09, 139.6478],
[773, 1441.0, 27.256, 140.396],
[774, 1444.9, 27.422, 141.1462],
[775, 1448.8, 27.589, 141.8984],
[776, 1452.7, 27.756, 142.6526],
[777, 1456.6, 27.924, 143.4088],
[778, 1460.5, 28.092, 144.167],
[779, 1464.4, 28.261, 144.9272],
[780, 1468.3, 28.43, 145.6894],
[781, 1472.2, 28.6, 146.4536],
[782, 1476.1, 28.77, 147.2198],
[783, 1480.0, 28.941, 147.988],
[784, 1483.9, 29.112, 148.7582],
[785, 1487.8, 29.284, 149.5304],
[786, 1491.7, 29.456, 150.3046],
[787, 1495.6, 29.629, 151.0808],
[788, 1499.5, 29.802, 151.859],
[789, 1503.4, 29.976, 152.6392],
[790, 1507.3, 30.15, 153.4214],
[791, 1511.2, 30.325, 154.2056],
[792, 1515.1, 30.5, 154.9918],
[793, 1519.0, 30.676, 155.78],
[794, 1522.9, 30.852, 156.5702],
[795, 1526.8, 31.029, 157.3624],
[796, 1530.7, 31.206, 158.1566],
[797, 1534.6, 31.384, 158.9528],
[798, 1538.5, 31.562, 159.751],
[799, 1542.4, 31.741, 160.5512],
[800, 1546.3, 31.92, 161.3534],
[801, 1550.2, 32.1, 162.1576],
[802, 1554.1, 32.28, 162.9638],
[803, 1558.0, 32.461, 163.772],
[804, 1561.9, 32.642, 164.5822],
[805, 1565.8, 32.824, 165.3944],
[806, 1569.7, 33.006, 166.2086],
[807, 1573.6, 33.189, 167.0248],
[808, 1577.5, 33.372, 167.843],
[809, 1581.4, 33.556, 168.6632],
[810, 1585.3, 33.74, 169.4854],
[811, 1589.2, 33.925, 170.3096],
[812, 1593.1, 34.11, 171.1358],
[813, 1597.0, 34.296, 171.964],
[814, 1600.9, 34.482, 172.7942],
[815, 1604.8, 34.669, 173.6264],
[816, 1608.7, 34.856, 174.4606],
[817, 1612.6, 35.044, 175.2968],
[818, 1616.5, 35.232, 176.135],
[819, 1620.4, 35.421, 176.9752],
[820, 1624.3, 35.61, 177.8174],
[821, 1628.2, 35.8, 178.6616],
[822, 1632.1, 35.99, 179.5078],
[823, 1636.0, 36.181, 180.356],
[824, 1639.9, 36.372, 181.2062],
[825, 1643.8, 36.564, 182.0584],
[826, 1647.7, 36.756, 182.9126],
[827, 1651.6, 36.949, 183.7688],
[828, 1655.5, 37.142, 184.627],
[829, 1659.4, 37.336, 185.4872],
[830, 1663.3, 37.53, 186.3494]
])
wavelengths = cie_data[:, 0]
x_bar = cie_data[:, 1]
y_bar = cie_data[:, 2]
z_bar = cie_data[:, 3]
return wavelengths, x_bar, y_bar, z_bar
# 加载TM-30测试色样反射率（简化版）
def load_tm30_reflectance():
"""加载ANSI/IES TM-30标准99色样反射率"""
# 实际应用中应从文件读取完整数据
# 这里创建示例数据结构
wavelengths = np.arange(380, 781, 5) # 5nm间隔
# 生成随机反射率数据 (0-1之间)
reflectance = np.random.rand(99, len(wavelengths))
# 确保反射率在0-1范围内
reflectance = np.clip(reflectance, 0.0, 1.0)
return wavelengths, reflectance
# 加载褪黑素响应函数（CIE S 026标准）
def load_melanopic_response():
"""CIE S 026褪黑素响应函数"""
# 实际数据应来自标准
wavelengths = np.arange(380, 781, 1)
# 使用高斯函数近似褪黑素响应曲线
response = np.exp(-0.5 * ((wavelengths - 480) / 30) ** 2)
# 归一化
response /= np.max(response)
return wavelengths, response
# ======================
# 2. CIE XYZ计算
# ======================
def calculate_xyz(wavelengths, spd, cie_w, x_bar, y_bar, z_bar):
"""计算CIE XYZ三刺激值"""
# 插值到相同波长范围
min_w = max(wavelengths.min(), cie_w.min())
max_w = min(wavelengths.max(), cie_w.max())
interp_w = np.arange(min_w, max_w + 1, 1)
# 插值SPD
f_spd = interpolate.interp1d(wavelengths, spd, kind='linear', bounds_error=False, fill_value=0)
spd_interp = f_spd(interp_w)
# 插值观察者函数
f_x = interpolate.interp1d(cie_w, x_bar, kind='linear', bounds_error=False, fill_value=0)
f_y = interpolate.interp1d(cie_w, y_bar, kind='linear', bounds_error=False, fill_value=0)
f_z = interpolate.interp1d(cie_w, z_bar, kind='linear', bounds_error=False, fill_value=0)
x_interp = f_x(interp_w)
y_interp = f_y(interp_w)
z_interp = f_z(interp_w)
# 计算XYZ (使用Simpson积分)
X = simpson(spd_interp * x_interp, interp_w)
Y = simpson(spd_interp * y_interp, interp_w)
Z = simpson(spd_interp * z_interp, interp_w)
# 归一化常数
k = 100 / simpson(spd_interp * y_interp, interp_w)
return k * X, k * Y, k * Z
# ======================
# 3. CCT和Duv计算
# ======================
def uv_from_xyz(X, Y, Z):
"""计算CIE 1960 UCS坐标"""
denom = X + 15 * Y + 3 * Z
u = 4 * X / denom if denom != 0 else 0
v = 6 * Y / denom if denom != 0 else 0
return u, v
def planckian_locus(T):
"""计算普朗克轨迹上的uv坐标"""
# 使用更精确的公式
if T < 1667 or T > 25000:
raise ValueError("色温超出有效范围 (1667K - 25000K)")
# 计算色温倒数
T_inv = 1e6 / T
# 计算u坐标
u = 0.860117757 + 1.54118254e-4 * T_inv + 1.28641212e-7 * T_inv ** 2
# 计算v坐标
v = 0.317398726 + 4.22806245e-5 * T_inv + 4.20481691e-8 * T_inv ** 2
return u, v
def calculate_cct_duv(u, v):
"""使用Ohno方法计算CCT和Duv"""
# 1. 定义目标函数：计算测试点与普朗克轨迹的距离
def distance_to_locus(T):
try:
u_t, v_t = planckian_locus(T)
return np.sqrt((u - u_t) ** 2 + (v - v_t) ** 2)
except ValueError:
return np.inf # 超出范围返回无穷大
# 2. 在合理温度范围内搜索最小值
result = minimize(distance_to_locus, 6500, bounds=[(1667, 25000)], method='L-BFGS-B')
if result.success:
cct = result.x[0]
u_t, v_t = planckian_locus(cct)
duv = np.sign(v - v_t) * result.fun
return cct, duv
else:
# 使用备选方法：McCamy近似公式
# 仅当优化失败时使用
n = (u - 0.3320) / (v - 0.1858)
cct = -449 * n ** 3 + 3525 * n ** 2 - 6823.3 * n + 5520.33
return cct, 0.0 # McCamy公式不提供Duv
# ======================
# 4. Rf和Rg计算
# ======================
def calculate_cam02ucs(X, Y, Z):
"""计算CAM02-UCS颜色空间坐标 (简化版本)"""
# 实际实现应使用完整的CAM02转换
# 这里使用简化转换作为示例
# 转换为XYZ到RGB (D65白点)
R = 3.2406 * X - 1.5372 * Y - 0.4986 * Z
G = -0.9689 * X + 1.8758 * Y + 0.0415 * Z
B = 0.0557 * X - 0.2040 * Y + 1.0570 * Z
# 非线性变换
R = np.where(R > 0.0031308, 1.055 * (R ** (1 / 2.4)) - 0.055, 12.92 * R)
G = np.where(G > 0.0031308, 1.055 * (G ** (1 / 2.4)) - 0.055, 12.92 * G)
B = np.where(B > 0.0031308, 1.055 * (B ** (1 / 2.4)) - 0.055, 12.92 * B)
# 转换为L*a*b*
L = 116 * np.cbrt(Y) - 16
a = 500 * (np.cbrt(X) - np.cbrt(Y))
b = 200 * (np.cbrt(Y) - np.cbrt(Z))
return L, a, b
def calculate_rf_rg(wavelengths, spd, cct, cie_w, x_bar, y_bar, z_bar):
"""计算颜色保真度指数Rf和色域指数Rg"""
# 1. 加载99色样反射率数据
ref_w, reflectance = load_tm30_reflectance()
# 2. 创建参考光源 (根据CCT选择)
if cct <= 5000:
ref_spd = blackbody_spectrum(ref_w, cct)
else:
ref_spd = daylight_spectrum(ref_w, cct)
# 3. 插值SPD到反射率波长
f_spd = interpolate.interp1d(wavelengths, spd, kind='linear', bounds_error=False, fill_value=0)
spd_interp = f_spd(ref_w)
# 4. 将CIE观察者函数插值到99色样的波长上
f_x = interpolate.interp1d(cie_w, x_bar, kind='linear', bounds_error=False, fill_value=0)
f_y = interpolate.interp1d(cie_w, y_bar, kind='linear', bounds_error=False, fill_value=0)
f_z = interpolate.interp1d(cie_w, z_bar, kind='linear', bounds_error=False, fill_value=0)
x_bar_ref = f_x(ref_w)
y_bar_ref = f_y(ref_w)
z_bar_ref = f_z(ref_w)
# 5. 计算每个色样在测试光源和参考光源下的颜色
test_colors = []
ref_colors = []
for i in range(99):
# 计算测试光源下的XYZ
X_test = simpson(spd_interp * reflectance[i] * x_bar_ref, ref_w)
Y_test = simpson(spd_interp * reflectance[i] * y_bar_ref, ref_w)
Z_test = simpson(spd_interp * reflectance[i] * z_bar_ref, ref_w)
# 计算参考光源下的XYZ
X_ref = simpson(ref_spd * reflectance[i] * x_bar_ref, ref_w)
Y_ref = simpson(ref_spd * reflectance[i] * y_bar_ref, ref_w)
Z_ref = simpson(ref_spd * reflectance[i] * z_bar_ref, ref_w)
# 转换为CAM02-UCS
L_test, a_test, b_test = calculate_cam02ucs(X_test, Y_test, Z_test)
L_ref, a_ref, b_ref = calculate_cam02ucs(X_ref, Y_ref, Z_ref)
test_colors.append([a_test, b_test])
ref_colors.append([a_ref, b_ref])
# 6. 计算色差ΔE
delta_e = []
for i in range(99):
de = np.sqrt((test_colors[i][0] - ref_colors[i][0]) ** 2 +
(test_colors[i][1] - ref_colors[i][1]) ** 2)
delta_e.append(de)
# 7. 计算Rf
avg_delta_e = np.mean(delta_e)
Rf = 100 - 4.6 * avg_delta_e
# 8. 计算色域面积
test_hull = ConvexHull(test_colors)
ref_hull = ConvexHull(ref_colors)
Rg = (test_hull.volume / ref_hull.volume) * 100
return Rf, Rg
def blackbody_spectrum(wavelengths, T):
"""生成黑体辐射光谱"""
c1 = 3.741771e-16 # W·m²
c2 = 1.438776e-2 # m·K
wavelengths_m = wavelengths * 1e-9
spd = c1 / (wavelengths_m ** 5 * (np.exp(c2 / (wavelengths_m * T)) - 1))
return spd / np.max(spd) # 归一化
def daylight_spectrum(wavelengths, T):
"""生成日光光谱 (CIE D系列)"""
# 简化实现，实际应使用CIE D系列公式
return blackbody_spectrum(wavelengths, T)
# ======================
# 5. 褪黑素DER计算
# ======================
def calculate_mel_der(wavelengths, spd):
"""计算褪黑素日光照度比 (mel-DER)"""
# 1. 加载褪黑素响应函数
mel_w, mel_response = load_melanopic_response()
# 2. 插值到相同波长范围
min_w = max(wavelengths.min(), mel_w.min())
max_w = min(wavelengths.max(), mel_w.max())
interp_w = np.arange(min_w, max_w + 1, 1)
f_spd = interpolate.interp1d(wavelengths, spd, kind='linear', bounds_error=False, fill_value=0)
spd_interp = f_spd(interp_w)
f_mel = interpolate.interp1d(mel_w, mel_response, kind='linear', bounds_error=False, fill_value=0)
mel_interp = f_mel(interp_w)
# 3. 计算待测光源的褪黑素刺激值
mel_edi = simpson(spd_interp * mel_interp, interp_w)
# 4. 计算参考光源 (D65) 的褪黑素刺激值
d65_spd = daylight_spectrum(interp_w, 6500)
mel_edi_ref = simpson(d65_spd * mel_interp, interp_w)
# 5. 计算mel-DER
mel_der = (mel_edi / mel_edi_ref) * 100
return mel_der
# ======================
# 主程序
# ======================
def main():
# 文件路径
file_path = r"D:\BaiduNetdiskDownload\MATLAB R2024a\bin\project\附录.xlsx"
try:
print("开始处理光源参数计算...")
# 1. 读取并预处理SPD数据
print("读取SPD数据...")
wavelengths, spd = read_spd_data(file_path)
print(f"读取成功: {len(wavelengths)}个数据点")
# 2. 加载CIE标准观察者函数
print("加载CIE标准观察者函数...")
cie_w, x_bar, y_bar, z_bar = load_cie_observer()
# 3. 计算CIE XYZ
print("计算CIE XYZ值...")
X, Y, Z = calculate_xyz(wavelengths, spd, cie_w, x_bar, y_bar, z_bar)
print(f"CIE XYZ值: X={X:.2f}, Y={Y:.2f}, Z={Z:.2f}")
# 4. 计算CCT和Duv
print("计算相关色温(CCT)和Duv...")
u, v = uv_from_xyz(X, Y, Z)
cct, duv = calculate_cct_duv(u, v)
print(f"相关色温CCT: {cct:.1f}K, Duv: {duv:.4f}")
# 5. 计算Rf和Rg
print("计算颜色保真度指数(Rf)和色域指数(Rg)...")
Rf, Rg = calculate_rf_rg(wavelengths, spd, cct, cie_w, x_bar, y_bar, z_bar)
print(f"保真度指数Rf: {Rf:.1f}, 色域指数Rg: {Rg:.1f}")
# 6. 计算mel-DER
print("计算褪黑素日光照度比(mel-DER)...")
mel_der = calculate_mel_der(wavelengths, spd)
print(f"褪黑素日光照度比: {mel_der:.1f}%")
# 7. 可视化结果
print("生成SPD可视化图表...")
plt.figure(figsize=(12, 8))
plt.plot(wavelengths, spd, 'b-', label='光源SPD')
plt.title('光源光谱功率分布', fontsize=14)
plt.xlabel('波长 (nm)', fontsize=12)
plt.ylabel('光谱功率 (W/m²/nm)', fontsize=12)
plt.grid(True, linestyle='--', alpha=0.7)
plt.legend()
plt.tight_layout()
plt.savefig('spd_plot.png', dpi=300)
plt.show()
print("\n所有参数计算完成!")
print("=" * 50)
print(f"相关色温(CCT): {cct:.1f} K")
print(f"距离普朗克轨迹的距离(Duv): {duv:.4f}")
print(f"保真度指数(Rf): {Rf:.1f}")
print(f"色域指数(Rg): {Rg:.1f}")
print(f"褪黑素日光照度比(mel-DER): {mel_der:.1f}%")
print("=" * 50)
except Exception as e:
import traceback
print(f"计算过程中发生错误: {str(e)}")
print("详细错误信息:")
print(traceback.format_exc())
if __name__ == "__main__":
main()
请把代码中简化的部分都改成精确计算，其他部分的代码不要动，最后输出完整版代码
最新发布
本文通过一个简单的C++程序示例介绍了如何使用标准输入输出流，并展示了如何编写一个基本的main函数。虽然代码简单，但它是理解更复杂C++程序的基础。
扫一扫
1930
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
