Color profile conversion - image not correct

I have jpg image with embeded ICC profile. I have extracted ICC profile from the file and now I need to convert decompressed image to sRGB.

My formulas are based on http://www.brucelindbloom.com/index.html?Math.html

I use this algorithm:

1. Convert image from RGB -> XYZ
2. Do chromatic adaptation on XYZ
3. Convert XYZ -> sRGB

For the tasks, I have this code:

Image2d<float> ColorSpace::ConvertRgbToXyz(const Image2d<uint8_t> & input,
    const IccProfile & inputIcc) {
    size_t len = input.GetPixelsCount();
    std::vector<float> data;
    data.resize(len * 3);

    for (size_t i = 0; i < len; i++) {
        const uint8_t * rgb = input.GetPixelStart(i);
        float r = (rgb[0] / 255.0f);
        float g = (rgb[1] / 255.0f);
        float b = (rgb[2] / 255.0f);

        r = (r <= 0) ? 0 : pow(r, inputIcc.rGamma);
        g = (g <= 0) ? 0 : pow(g, inputIcc.gGamma);
        b = (b <= 0) ? 0 : pow(b, inputIcc.bGamma);

        //sRGB
        float x = inputIcc.rXYZ[0] * r + inputIcc.gXYZ[0] * g + inputIcc.bXYZ[0] * b;
        float y = inputIcc.rXYZ[1] * r + inputIcc.gXYZ[1] * g + inputIcc.bXYZ[1] * b;
        float z = inputIcc.rXYZ[2] * r + inputIcc.gXYZ[2] * g + inputIcc.bXYZ[2] * b;

        data[i * 3 + 0] = x;
        data[i * 3 + 1] = y;
        data[i * 3 + 2] = z;

    }

    return Image2d<float>(input.GetWidth(),
        input.GetHeight(),
        std::move(data),
        PixelFormat::XYZ);
}

Chromatic conversion from white point [Xws, Yws, Zws] to target white point [Xwd, Ywd, Zwd] (used white for D65 => Xwd = 0.95047f, Ywd = 1.0f, Zwd = 1.08883f)

void ColorSpace::ChromaticAdaptation(Image2d<float> & input,
    float Xws, float Yws, float Zws,
    float Xwd, float Ywd, float Zwd,
    ChromaticMethod method) {
    size_t len = input.GetPixelsCount();

    Matrix3x3 m = Matrix3x3(
        (Xwd / Xws), 0, 0,
        0, (Ywd / Yws), 0,
        0, 0, (Zwd / Zws)
    );

    if (method == ChromaticMethod::Bradford) {
        const Matrix3x3 mA = Matrix3x3(
            0.8951000, 0.2664000, -0.1614000,
            -0.7502000, 1.7135000, 0.0367000,
            0.0389000, -0.0685000, 1.0296000
        );

        const Matrix3x3 mAInv = Matrix3x3(
            0.9869929, -0.1470543, 0.1599627,
            0.4323053, 0.5183603, 0.0492912,
            -0.0085287, 0.0400428, 0.9684867
        );

        float xs = mA.M[0][0] * Xws + mA.M[0][1] * Yws + mA.M[0][2] * Zws;
        float ys = mA.M[1][0] * Xws + mA.M[1][1] * Yws + mA.M[1][2] * Zws;
        float zs = mA.M[2][0] * Xws + mA.M[2][1] * Yws + mA.M[2][2] * Zws;

        float xd = mA.M[0][0] * Xwd + mA.M[0][1] * Ywd + mA.M[0][2] * Zwd;
        float yd = mA.M[1][0] * Xwd + mA.M[1][1] * Ywd + mA.M[1][2] * Zwd;
        float zd = mA.M[2][0] * Xwd + mA.M[2][1] * Ywd + mA.M[2][2] * Zwd;

        m = Matrix3x3(
            (xd / xs), 0, 0,
            0, (yd / ys), 0,
            0, 0, (zd / zs)
        );

        Matrix3x3 res = mAInv;
        res *= m;
        res *= mA;

        m = res;
    }


    for (size_t i = 0; i < len; i++) {
        float * xyz = input.GetPixelStart(i);
        float x = xyz[0];
        float y = xyz[1];
        float z = xyz[2];

        float xd = m.M[0][0] * x + m.M[0][1] * y + m.M[0][2] * z;
        float yd = m.M[1][0] * x + m.M[1][1] * y + m.M[1][2] * z;
        float zd = m.M[2][0] * x + m.M[2][1] * y + m.M[2][2] * z;

        xyz[0] = xd;
        xyz[1] = yd;
        xyz[2] = zd;

    }
}

And finally

Image2d<uint8_t> ColorSpace::ConvertXyzToSRgb_D65(const Image2d<float> & input) {
    size_t len = input.GetPixelsCount();
    std::vector<uint8_t> data;
    data.resize(len * 3);

    for (size_t i = 0; i < len; i++) {
        const float * xyz = input.GetPixelStart(i);
        float x = xyz[0];
        float y = xyz[1];
        float z = xyz[2];

        //sRGB
        float r =  3.2404542f * x  + -1.5371385f * y + -0.4985314f * z;
        float g = -0.9692660f * x  +  1.8760108f * y +  0.0415560f * z;
        float b =  0.0556434f * x  + -0.2040259f * y +  1.0572252f * z;

        //color companding
        //for sRGB
        r = (r > 0.0031308f) ? 1.055f * std::pow(r, 1 / 2.4f) - 0.055f : 12.92f * r;
        g = (g > 0.0031308f) ? 1.055f * std::pow(g, 1 / 2.4f) - 0.055f : 12.92f * g;
        b = (b > 0.0031308f) ? 1.055f * std::pow(b, 1 / 2.4f) - 0.055f : 12.92f * b;

        data[i * 3 + 0] = static_cast<uint8_t>(255.0f * std::clamp(r, 0.0f, 1.0f));
        data[i * 3 + 1] = static_cast<uint8_t>(255.0f * std::clamp(g, 0.0f, 1.0f));
        data[i * 3 + 2] = static_cast<uint8_t>(255.0f * std::clamp(b, 0.0f, 1.0f));    
    }

    return Image2d<uint8_t>(input.GetWidth(),
        input.GetHeight(),
        std::move(data),
        PixelFormat::RGB);
}

However, the result image is not correct compared with the image opened in image viewer (IrfanView with enabled ICC profile, Firefox or Photoshop).

Left is, how it should look. Right is my output.

Without any color conversion, the image is loaded as:

So my version is closer to the "how it should look", but not quite there.

I am not sure if ConvertRgbToXyz is correct. I take rXYZ / gXYZ / bXYZ values directly from ICC, based on documentation http://www.color.org/specification/ICC1v43_2010-12.pdf - section 9.2.44 for rXYZ -

This tag contains the first column in the matrix, which is used in matrix/TRC transforms.

My ICC profile is this:

Header:
  size         = 25908 bytes
  CMM          = 'argl'
  Version      = 2.2.0
  Device Class = Input
  Color Space  = RGB
  Conn. Space  = XYZ
  Date, Time   = 10 Mar 2017, 13:05:37
  Platform     = Microsoft
  Flags        = Not Embedded Profile, Use anywhere
  Dev. Mnfctr. = 0x0
  Dev. Model   = 0x0
  Dev. Attrbts = Reflective, Glossy, Positive, Color
  Rndrng Intnt = Perceptual
  Illuminant   = 0.964203, 1.000000, 0.824905    [Lab 100.000000, 0.000498, -0.000436]
  Creator      = 'argl'

'desc': CoCa 10.3.2017
'cprt': Copyright_free
'dmnd': proserv 600i
'dmdd': Model

'wtpt': 0.768753, 0.822220, 0.665100    [Lab 92.672746, -4.781618, 1.218722]
'bkpt': 0.003693, 0.004395, 0.003387    [Lab 3.969564, -2.199059, 0.448583]

'rXYZ': 0.956482, 0.326324, 0.010391    [Lab 63.862187, 154.428754, 91.157829]
'gXYZ': 0.070267, 0.995911, -0.271652    [Lab 99.841662, -290.469249, 685.018939]
'bXYZ': 0.136230, -0.105103, 1.321716    [Lab -94.938734, 600.671989, -370.133480]

'rTRC': Curve is gamma of 1.234375
'gTRC': Curve is gamma of 1.296875
'bTRC': Curve is gamma of 1.312500

I have found the problem.

I was using wrong reference white for chromatic adaptation. I was using value from wtpt, but Illuminant should be used instead as an input white.

Output white is still set to D65.