vec3 convert_ciexyz_to_linear_rgb(vec3 xyz) {
    // http://www.brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html sRGB M^-1 matrix
    //    mat3 M = mat3(
    //        3.2404542, -1.5371385, -0.4985314,
    //        -0.9692660, 1.8760108, 0.0415560,
    //        0.0556434, -0.2040259, 1.0572252
    //    );
    // TODO: figure out which is correct, GLSL uses column-major for matrices

    mat3 M = mat3(
    3.2404542, -0.9692660, 0.0556434,
    -1.5371385, 1.8760108, -0.2040259,
    -0.4985314, 0.0415560, 1.0572252
    );
    return M * xyz;
}


vec3 convert_cielab_to_xyz(vec3 cielab, vec3 ref_white) {
    float Xn = ref_white.x;
    float Yn = ref_white.y;
    float Zn = ref_white.z;

    float L = cielab.x;
    float a = cielab.y;
    float b = cielab.z;

    float epsilon = 216./24289.;
    float kappa = 24389./27.;

    float fy = (L + 16.) / 116.;
    float fz = fy - b / 200.;
    float fx = a / 500. + fy;

    float xr;
    float yr;
    float zr;

    if (pow(fx, 3.) > epsilon) {
        xr = pow(fx, 3.);
    } else {
        xr = (116. * fx - 16.) / kappa;
    }

    if (L > kappa * epsilon) {
        yr = pow((L + 16.) / 116., 3.);
    } else {
        yr = L / kappa;
    }

    if (pow(fz, 3.) > epsilon) {
        zr = pow(fz, 3.);
    } else {
        zr = (116. * fz - 16.) / kappa;
    }

    float X = xr * Xn;
    float Y  = yr * Yn;
    float Z = zr * Zn;
    return vec3(X, Y, Z);
}


vec3 convert_cielab_to_rgb(vec3 cielab) {
    // using D65
    float Xn = .950489;
    float Yn = 1.;
    float Zn = 1.0888;
    vec3 xyz = convert_cielab_to_xyz(cielab, vec3(Xn, Yn, Zn));
    vec3 linear_rgb = convert_ciexyz_to_linear_rgb(xyz);
    return linear_rgb;
}

vec3 convert_cieluv_to_rgb(vec3 cieluv) {
    return cieluv;// TODO
}

vec3 convert_oklab_to_rgb(vec3 oklab) {
    return oklab;// TODO
}