Why does ColorConvert give a different answer than other websites when converting from RGB to XYZ?

Hopefully someone else can shed more light on this, someone with more expertise on colorspaces, but I was able to find out a little bit more. I found Bruce Lindbloom's website which lists many different transformation matrices for going from RGB to XYZ. The matrix used in the above-linked wiki page corresponds to converting between sRGB and XYZ, using a reference white of D65. I know not what the reference white value means, but it will apparently be important here.

I tried the hard way to find out what matrix Mathematica uses to convert between RGB and XYZ, just by taking 3 random colors and converting them and then solving for the matrix,

varmat = 
  Array[ToExpression[
     "m" <> IntegerString[#1] <> IntegerString[#2]] &, {3, 3}];
eqn = {x, y, z} == varmat.(If[# > .04045,
        ((# + 0.055)/1.055)^2.4,
        #/12.92] & /@ {r, g, b});
eqns = Table[
   rgb = RandomReal[1, 3];
   xyz = List @@ ColorConvert[RGBColor @@ rgb, XYZColor];
   eqn /. (#1 -> #2 & @@@ Transpose@{{r, g, b}, rgb}) /. (#1 -> #2 & @@@
       Transpose@{{x, y, z}, xyz})
   , {3}];
varmat /. First@Solve[eqns, Flatten@varmat] // MatrixForm

And this matches exactly the transformation listed on Lindbloom's site for sRGB with a reference white of D50. I had not seen anything listed on the help pages for RGBColor or XYZColor about a reference white. But then I found the following on the page for ColorConvert

ColorConvert automatically performs chromatic (white point) adaptation. D50 white point is assumed for "XYZ", "LAB", "LUV", and "LCH" and D65 for "RGB", "CMYK", "HSB", and "Grayscale".

So that solves this problem, I post it here in case anyone in the future is trying to convert colors and notices the discrepancy.

I do wonder why Mathematica uses a different reference white than seems to be the norm for sRGB.

Edit: I was also able to find that Mathematica uses the D50 Reference White when converting to CIELAB from XYZ, which is again different than many other sites that have conversion utilities.

Just for reference purposes, here is a manual re-implementation of ColorConvert[color, "RGB" -> "XYZ"]:

(* whitepoints *)
d65 = {0.95047, 1., 1.08883}; d50 = {0.96422, 1., 0.82521};

(* cone response domain *)
bradford = {{0.8951, 0.2664, -0.1614},
            {-0.7502, 1.7135, 0.0367},
            {0.0389, -0.0685, 1.0296}};

(* inverse companding function *)
InversesRGBGamma = Function[x, Piecewise[{{x/12.92, x <= 0.04045}},
                                         ((x + 0.055)/1.055)^2.4], Listable];

(* chromaticity coordinates of primaries (RGB) *)
srgb = {{0.64, 0.33}, {0.3, 0.6}, {0.15, 0.06}};

(* chromatic adaptation matrix, D65 to D50 *)
d65tod50 = LinearSolve[bradford, DiagonalMatrix[(bradford.d50)/(bradford.d65)].bradford];

(* RGB -> XYZ matrix *)
r2x = With[{xyz = Transpose[{#1/#2, 1., (1 - #1 - #2)/#2} & @@@ srgb]}, 
           xyz.DiagonalMatrix[LinearSolve[xyz, d65]]];

rgb2xyz[col_RGBColor] := d65tod50.r2x.InversesRGBGamma[List @@ col]

Example:

rgb2xyz /@ ColorData[61, "ColorList"]
   {{0.199667, 0.106547, 0.00721723}, {0.073949, 0.056383, 0.147498},
    {0.823078, 0.923995, 0.234938}, {0.226765, 0.233149, 0.365262},
    {0.099347, 0.136449, 0.027346}, {0.521978, 0.511432, 0.0633633},
    {0.0405846, 0.033663, 0.113543}, {0.3308, 0.169661, 0.0191982},
    {0.671893, 0.718828, 0.33747}}

List @@@ ColorConvert[ColorData[61, "ColorList"], XYZColor]
   {{0.199667, 0.106547, 0.00721724}, {0.0739489, 0.0563829, 0.147498},
    {0.823078, 0.923995, 0.234938}, {0.226765, 0.233149, 0.365262},
    {0.099347, 0.136449, 0.027346}, {0.521978, 0.511432, 0.0633633},
    {0.0405845, 0.033663, 0.113543}, {0.3308, 0.16966, 0.0191983},
    {0.671893, 0.718828, 0.33747}}

There are discrepancies, but they seem slight.

---

For completeness, here is the manual "XYZ" -> "RGB" conversion:

sRGBGamma = Function[x, With[{z = Abs[x]},
                             Sign[x] Piecewise[{{12.92 z, z <= 0.0031308}},
                                               1.055 z^(1/2.4) - 0.055]],
                     Listable];

d50tod65 = LinearSolve[bradford, DiagonalMatrix[(bradford.d65)/(bradford.d50)].bradford];

xyz2rgb[col_XYZColor] :=
   Clip[sRGBGamma[LinearSolve[r2x, d50tod65.(List @@ col)]], {0., 1.}]