convert BGR to Lab without OpenCV

Question

With little experience in color spaces, I used the following code to convert BGR image (array of unsigned characters where each channel ranges from 0 to 255) to lab color space:

double F(double input) // function f(...), which is used for defining L, a and b changes within [4/29,1]
{
    if (input > 0.008856)
        return (pow(input, 0.333333333)); // maximum 1
    else
        return ((841/108)*input + 4/29);  //841/108 = 29*29/36*16
}

// RGB to XYZ
void RGBtoXYZ(uchar R, uchar G, uchar B, double &X, double &Y, double &Z)
{
    // RGB Working Space: sRGB 
    // Reference White: D65
    X = 0.412453*R + 0.357580*G + 0.189423*B; // maximum value = 0.959456 * 255 = 244.66128
    Y = 0.212671*R + 0.715160*G + 0.072169*B; // maximum value = 1 * 255 = 255
    Z = 0.019334*R + 0.119193*G + 0.950227*B; //  maximum value = 1.088754 * 255 = 277.63227
}

// XYZ to CIELab
void XYZtoLab(double X, double Y, double Z, double &L, double &a, double &b)
{
    const double Xo = 244.66128; // reference white
    const double Yo = 255.0;
    const double Zo = 277.63227;
    L = 116 * F(Y / Yo) - 16; // maximum L = 100
    a = 500 * (F(X / Xo) - F(Y / Yo)); // maximum 
    b = 200 * (F(Y / Yo) - F(Z / Zo));
}

// RGB to CIELab
void RGBtoLab(double R, double G, double B, double &L, double &a, double &b)
{
    double X, Y, Z;
    RGBtoXYZ(R, G, B, X, Y, Z);
    XYZtoLab(X, Y, Z, L, a, b);
}

I have re-converted the resulting lab image to BGR (using cvtcolor) to display it using OpenCV using the following code (I don't want to do the conversion using OpenCV, I have just used it to display the results. Basically I wanted to understand how color space conversion works):

// Lchannel, Achannel, Bchannel are arrays of type double
cv::Mat temp64bitL(height, width, CV_64FC1, Lchannel);
cv::Mat temp32bitL;
temp64bitL.convertTo(temp32bitL, CV_32F);

cv::Mat temp64bitA(height, width, CV_64FC1, Achannel);
cv::Mat temp32bitA;
temp64bitA.convertTo(temp32bitA, CV_32F);
cv::Mat temp64bitB(height, width, CV_64FC1, Bchannel);
cv::Mat temp32bitB;
temp64bitB.convertTo(temp32bitB, CV_32F);
    cv::Mat chan[3] = {
    temp32bitL, temp32bitA, temp32bitB
};
cv::Mat merged;
cv::merge(chan, 3, merged);
cv::Mat BGRImage;
cv::cvtColor(merged, BGRImage, CV_Lab2BGR, 3);

However, the computed image is different from the original image. is that due to a problem in the code?

Answer 1

Your code has a bug in double F(double input). It does not work as intended because of the integer division you have. You might be willing to change the function to read something like below. Note the double castings to make the divisions work in the floating-point domain, and the use of cbrt instead of pow.

#include <cmath>

double F(double input) // function f(...), which is used for defining L, a and b
                       // changes within [4/29,1]
{
  if (input > 0.008856)
    return std::cbrt(input); // maximum 1 --- prefer cbrt to pow for cubic root
  else
    return ((double(841) / 108) * input +
            double(4) / 29); // 841/108 = 29*29/36*16
}

Then, another problem could be the reference values you are using for the XYZ space. We have the below reference values, coming from D65 / CIE-1931:

double Xo = 95.047;
double Yo = 100;
double Zo = 108.883;

Then, our RGBtoXYZ conversion was working like this:

template <class float_t> struct Convert<XYZ<float_t>> {
  template <class real_t> static XYZ<float_t> from(const RGB<real_t> &rhs) {
    // Assume RGB has the type invariance satisfied, i.e., channels \in [0,255]
    float_t var_R = float_t(rhs.comp1()) / 255;
    float_t var_G = float_t(rhs.comp2()) / 255;
    float_t var_B = float_t(rhs.comp3()) / 255;

    var_R = (var_R > 0.04045) ? std::pow((var_R + 0.055) / 1.055, 2.4)
                              : var_R / 12.92;
    var_G = (var_G > 0.04045) ? std::pow((var_G + 0.055) / 1.055, 2.4)
                              : var_G / 12.92;
    var_B = (var_B > 0.04045) ? std::pow((var_B + 0.055) / 1.055, 2.4)
                              : var_B / 12.92;

    var_R *= 100;
    var_G *= 100;
    var_B *= 100;

    return XYZ<float_t>{var_R * float_t(0.4124) + var_G * float_t(0.3576) +
                            var_B * float_t(0.1805),
                        var_R * float_t(0.2126) + var_G * float_t(0.7152) +
                            var_B * float_t(0.0722),
                        var_R * float_t(0.0193) + var_G * float_t(0.1192) +
                            var_B * float_t(0.9505)};
  }
};

where RGB was assumed to have its channels inside the valid range, as stated in the comment. Then, the XYZtoLAB function we have is the same except for the cbrt and reference value changes.

EDIT. Above numbers are obtained from EasyRGB's Math page. You can find the conversion from sRGB to XYZ and XYZ to Lab on the page, with a table of XYZ reference values. What we used was the set for "Daylight, sRGB, Adobe RGB."