HDR color-space transformations result in negative RGB values (Yxy to XYZ to sRGB)

Question

I'm currently adding HDR to an old engine and stumbled on a color-space transformation problem.

1. I'm defining my lights in the Yxy color space
2. Then I'm converting Yxy to XYZ
3. XYZ to sRGB transformation.
4. Using RGB values > 1.0 when rendering and normalizing the result with tone mapping in the end.

I'm working with rather huge numbers since the main light source is the sun having up to 150k Lux Illuminance.

YxyToXYZ function

osg::Vec3 YxyToXYZ( const osg::Vec3& Yxy )
{
    if ( Yxy[2] > 0.0f )
    {
        return osg::Vec3( Yxy[0] * Yxy[1] / Yxy[2] , Yxy[0] , Yxy[0] * ( 1.0f - Yxy[1] - Yxy[2])  / Yxy[2] );
    }
    else
    {
        return osg::Vec3( 0.0f , 0.0f , 0.0f );   
    }
}

XYZtosRGB

osg::Vec3 XYZToSpectralRGB( const osg::Vec3& XYZ )
{
    // Wikipedia matrix
    osg::Vec3 rgb;
    rgb[0] = 3.240479  * XYZ[0] - 1.537150 * XYZ[1] - 0.498535 * XYZ[2];
    rgb[1] = -0.969256 * XYZ[0] + 1.875992 * XYZ[1] + 0.041556 * XYZ[2];
    rgb[2] = 0.055648 * XYZ[0]  - 0.204043 * XYZ[1] + 1.057311 * XYZ[2];

    std::cout << "newrgb rgb r:" << rgb[0] << " g:" << rgb[1] << " b:" << rgb[2] << std::endl;

    // The matrix in pbrt book p. 235 gives unexpected results. We expect that if we have
    // x = y = 0.33333 we get a white pixel but that matrix gives red. Hence we use a different
    // matrix that is often used by 3D people
    rgb[0] = 2.5651  * XYZ[0] -1.1665 * XYZ[1] -0.3986 * XYZ[2];
    rgb[1] = -1.0217 * XYZ[0] + 1.9777 * XYZ[1] + 0.0439 * XYZ[2];
    rgb[2] = 0.0753 * XYZ[0]  -0.2543 * XYZ[1] + 1.1892 * XYZ[2];

    std::cout << "oldrgb rgb r:" << rgb[0] << " g:" << rgb[1] << " b:" << rgb[2] << std::endl;

    return rgb;
}

Test samples:

Yxy Y:1 x:1 y:1
XYZ X:1 Y:1 Z:-1
newrgb rgb r:2.20186 g:0.86518 b:-1.20571
oldrgb rgb r:1.7972  g:0.9121  b:-1.3682

Yxy Y:25 x:0.26 y:0.28
XYZ X:23.2143 Y:25 Z:41.0714
newrgb rgb r:16.3211 g:26.106  b:39.616
oldrgb rgb r:14.0134 g:27.5275 b:44.2327

Yxy Y:3100 x:0.27 y:0.29
XYZ X:2886.21 Y:3100 Z:4703.45
newrgb rgb r:2242.7  g:3213.56 b:4501.09
oldrgb rgb r:1912.47 g:3388.51 b:5022.34

Yxy Y:6e+06 x:0.33 y:0.33
XYZ X:6e+06 Y:6e+06 Z:6.18182e+06
newrgb rgb r:7.13812e+06 g:5.69731e+06 b:5.64573e+06
oldrgb rgb r:5.92753e+06 g:6.00738e+06 b:6.27742e+06

Question:

1. I suppose negative values should be just clamped away? Or do i have a error in my calculations?
2. The 2 matrices produce similar but different values (primaries very close to sRGB). I would like to replace the old matrix (which I have no idea where it's coming from) with the Wiki one. Does anyone know where the old matrix is coming from and which would be the correct one?
3. I found a partial answer @ Yxy to RGB conversion which sounds like from a former colleague :)... but doesn't solve my problems.

Many Thanks

Answer 1

Your newrgb rgb computations are correct, I get the same output using Colour:

import colour

xyY = (1.0, 1.0, 1.0)
XYZ = colour.xyY_to_XYZ(xyY)
print(colour.XYZ_to_sRGB(XYZ, apply_encoding_cctf=False))

xyY = (0.26, 0.28, 25.0)
XYZ = colour.xyY_to_XYZ(xyY)
print(colour.XYZ_to_sRGB(XYZ, apply_encoding_cctf=False))

xyY = (0.27, 0.29, 3100.0)
XYZ = colour.xyY_to_XYZ(xyY)
print(colour.XYZ_to_sRGB(XYZ, apply_encoding_cctf=False))

xyY = (0.33, 0.33, 6e+06)
XYZ = colour.xyY_to_XYZ(xyY)
print(colour.XYZ_to_sRGB(XYZ, apply_encoding_cctf=False))

# [ 2.2020461   0.86530782 -1.20530687]
# [ 16.31921754  26.10605109  39.60507773]
# [ 2242.49706509  3213.58480355  4499.85141917]
# [ 7138073.69325106  5697605.86069197  5644291.15301836]

You get negative values in the first conversion because your xy chromaticity coordinates are outside the spectral locus, thus they represent imaginary colours.

The Wikipedia matrix is correct and is the one from IEC 61966-2-1:1999 which is the official standard for sRGB colourspace.