cvtColor() does not correctly convert LAB images to RGB (Patch #1511)

The code in the Lab2RGB_f struct (in modules/imgproc/src/color.cpp) does not calculate the inverse of the RGB-LAB mapping equations described in the documentation. The following code snippet demonstrates the result of generating an BGR image, converting it to LAB, converting it back to BGR, and then writing the results to disk.

 1const int sz=100;
 2const int szm=sz-1;
 3float pi2=2*3.1415f;
 4cv::Mat orig(sz,sz,CV_32FC3);
 5
 6// Generate a pretty test image
 7for(int i=0; i<sz; i++) {
 8  for(int j=0; j<sz; j++) {
 9    float b=(1+cos((szm-i)*(szm-j)*pi2/(10*float(szm))))/2;
10    float g=(1+cos((szm-i)*j*pi2/(10*float(szm))))/2;
11    float r=(1+sin(i*j*pi2/(10*float(szm))))/2;
12
13    // The following lines aren't necessary, but just to prove that
14    // the BGR values all lie in [0,1]...
15    if(b<0) b=0; else if(b>1) b=1;
16    if(g<0) g=0; else if(g>1) g=1;
17    if(r<0) r=0; else if(r>1) r=1;
18    orig.at<cv::Vec3f>(i,j)=cv::Vec3f(b,g,r);
19  }
20}
21
22// Convert test image to LAB
23cv::Mat lab;
24cv::cvtColor(orig,lab,CV_BGR2Lab);
25// Convert LAB image back to BGR
26cv::Mat recons;
27cv::cvtColor(lab,recons,CV_Lab2BGR);
28
29// Write original test image
30cv::Mat orig8u;
31orig.convertTo(orig8u,CV_8U,255);
32cv::imwrite("orig.tiff",orig8u);
33cv::Mat recons8u;
34// Write image converted back to RGB from LAB space
35recons.convertTo(recons8u,CV_8U,255);
36cv::imwrite("recons.tiff",recons8u);
37

I have written my own code that correctly inverts the equations in the documentation, but I have not spent enough time studying the OpenCV code to work out how to incorporate it. Nonetheless, here it is (it assumes that the template arguments are float or double):

 1template<class S,class T>
 2void lab2rgb(
 3    const S li,
 4    const S ai,
 5    const S bi,
 6    T& ro,
 7    T& go,
 8    T& bo)
 9{
10  const S lThresh=0.008856*903.3;
11
12  S y, fy;
13  if(li<=lThresh) {
14    y=li/903.3;
15    fy=7.787*y+16.0/116.0;
16  } else {
17    fy=(li+16.0)/116.0;
18    y=fy*fy*fy;
19  }
20  S fxz[]={S(ai/500.0+fy),S(bi/-200.0+fy)};
21  S xzUnnormr2;
22  const S fThresh=7.787*0.008856+16.0/116.0;
23
24  for(int i=0; i<2; i++) {
25    S f=fxz[i];
26    if(f<=fThresh) {
27      xzUnnorm[i]=(f-16.0/116.0)/7.787;
28    } else {
29      xzUnnorm[i]=f*f*f;
30    }
31  }
32
33  S x=xzUnnormr0*0.950456, z=xzUnnormr1*1.088754;
34  ro=3.240479*x-1.53715*y-0.498535*z;
35  go=-0.969256*x+1.875991*y+0.041556*z;
36  bo=0.055648*x-0.204043*y+1.057311*z;
37}