Given a hyperspectral image and the sensitivity data of photoreceptors, I need to obtain the responses of the long, medium, and short photoreceptors for each pixel at every wavelength. Then, I have to convert from LMS to XYZ and then to RGB.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from spectral import open_image

img = open_image('food_chili.hdr').load()
metadata = img.metadata

if 'wavelength' in metadata:
    wavelength = np.array(metadata['wavelength'], dtype=float)
else:
    raise ValueError("not found!")

sensitivity_data = pd.read_csv("linss2_10e_5_8dp.csv")
sensitivity_data.columns = sensitivity_data.columns.str.strip()

wavelengths_LMS = sensitivity_data.iloc[:, 0].values
L_sensitivity = sensitivity_data.iloc[:, 1].values
M_sensitivity = sensitivity_data.iloc[:, 2].values
S_sensitivity = np.nan_to_num(sensitivity_data.iloc[:, 3].values)


L_interpolated = np.interp(wavelength, wavelengths_LMS, L_sensitivity)
M_interpolated = np.interp(wavelength, wavelengths_LMS, M_sensitivity)
S_interpolated = np.interp(wavelength, wavelengths_LMS, S_sensitivity)

#normalization LMS
L_interpolated /= np.trapz(L_interpolated, wavelength)
M_interpolated /= np.trapz(M_interpolated, wavelength)
S_interpolated /= np.trapz(S_interpolated, wavelength)


min_wavelength, max_wavelength = 400, 700
visible_bands = np.where((wavelength >= min_wavelength) & (wavelength <= max_wavelength))[0]
visible_wavelength = wavelength[visible_bands]
L_interpolated = L_interpolated[visible_bands
M_interpolated = M_interpolated[visible_bands]
S_interpolated = S_interpolated[visible_bands]

# response LMS
response = np.zeros((img.shape[0], img.shape[1], 3)) 
for i in range(img.shape[0]):
    for j in range(img.shape[1]):
        pixel = img[i, j, :].squeeze()[bande_visibili] 
        response[i, j, 0] = np.sum(pixel * L_interpolated)  
        response[i, j, 1] = np.sum(pixel * M_interpolated) 
        response[i, j, 2] = np.sum(pixel * S_interpolated) 

print("LMS min and max:")
print("L min:", np.min(response[:, :, 0]), "L max:", np.max(response[:, :, 0]))
print("M min:", np.min(response[:, :, 1]), "M max:", np.max(response[:, :, 1]))
print("S min:", np.min(response[:, :, 2]), "S max:", np.max(response[:, :, 2]))

matrix_LMS_to_XYZ = np.array([
    [0.4002, 0.7075, -0.0808],
    [-0.2263, 1.1653, 0.0457],
    [0.0000, 0.0000, 0.9182]
])
response_XYZ = np.dot(response, matrix_LMS_to_XYZ.T)

print("xyz min and max:")
print("X min:", np.min(response_XYZ[:, :, 0]), "X max:", np.max(response_XYZ[:, :, 0]))
print("Y min:", np.min(response_XYZ[:, :, 1]), "Y max:", np.max(response_XYZ[:, :, 1]))
print("Z min:", np.min(response_XYZ[:, :, 2]), "Z max:", np.max(response_XYZ[:, :, 2]))


matrix_XYZ_to_RGB = np.array([
    [3.2406, -1.5372, -0.4986],
    [-0.9689, 1.8758, 0.0415],
    [0.0557, -0.2040, 1.0570]
])
rgb_image = np.dot(response, matrix_LMS_to_RGB.T)

rgb_image[rgb_image < 0] = 0
rgb_image /= np.max(rgb_image) 
rgb_image = (rgb_image * 255).astype(np.uint8)
plt.imshow(rgb_image)
plt.title('Immagine RGB a partire da LMS')
plt.axis('off')
plt.show()
plt.imsave('rgb_image_from_lms.png', rgb_image)

The problem is that the resulting image has unbalanced colors; for example, green tends to shift towards red.