RGB转Lab的那些事(一)

Matlab与OpenCV都有RGB转Lab的具体实现，然而两种版本给出的结果似乎并不一样，那么两者有何区别呢？

首先了解RGB转Lab的理论知识：

一般地，由RGB转Lab都需要先将RGB转为XYZ，再由XYZ转为Lab。通常我们获取的图片往往是sRGB格式的

sRGB (“standard RGB” [41]) was developed (jointly by Hewlett-Packard and Microsoft) with the goal of creating a precisely specified color space for these applications, based on standardized mappings with respect to the col-orimetric CIE XYZ color space.

Several standard image formats, including EXIF (JPEG) and PNG are based on sRGB color data, which makes sRGB the de facto standard for dig-ital still cameras, color  rinters, an d other imaging devices at the consumer level.

Thus, in practice, working with any RGB color data almost always means dealing with sRGB.

sRGB is a nonlinear color space with respect to the XYZ coordinate system, and it is important to carefully distinguish between the linear and nonlinear RGB component values. The nonlinear values (denoted R', G', B') represent the actual color tuples, the data valu es read from an image file or received from a digital camera.

也就是说我们通常获取的图片是非线性颜色空间的sRGB格式，需要将sRGB转为线性颜色空间的RGB

格式（这一步即为gamma矫正过程）

 

注意：Matlab与OpenCV的gamma常数好像都是取得0.04045，而不是0.03928

XYZ转Lab的公式如下：

至此，RGB转Lab的理论说明已经很清楚，其实大致的步骤和网络上的资源也很类似，只是细节方面有一些不一样的地方。Matlab与OpenCV的差别就在于

线性映射矩阵M_RGB与白色参考点，Matlab用的是D50，而OpenCV用的是D65。

首先来看一看Matlab的RGB转Lab的D65、D50实现：

function [L, a, b] = sRGB2LabD50(sR, sG, sB)
% SRGB2LAB convert sRGB to CIELAB(Matlab implemention).
% 
if nargin == 1
    sB = double(sR(:,:,3));
    sG = double(sR(:,:,2));
    sR = double(sR(:,:,1));
end

% normalize for [0 1]
if max(max(sR)) > 1.0 || max(max(sG)) > 1.0 || max(max(sB)) > 1.0
    sR = double(sR) / 255;
    sG = double(sG) / 255;
    sB = double(sB) / 255;
end

[M, N] = size(sR);
s = M * N;
sRL = reshape(sR, 1, s);
sGL = reshape(sG, 1, s);
sBL = reshape(sB, 1, s);

% gamma correction (gamma 2.2)
gammaConst = 0.04045; 
sRL = ( sRL<=gammaConst ) .* sRL  ./ 12.92 +  ...
      ( ( sRL>gammaConst  ) .* (0.055 + sRL) ./ 1.055 )  .^2.4;
sGL = ( sGL<=gammaConst ) .* sGL  ./ 12.92 +  ...
      ( ( sGL>gammaConst  ) .* (0.055 + sGL) ./ 1.055 )  .^2.4;
sBL = ( sBL<=gammaConst ) .* sBL  ./ 12.92 +  ...
      ( ( sBL>gammaConst  ) .* (0.055 + sBL) ./ 1.055 )  .^2.4;
  
RGB = [sRL; sGL; sBL];

% linear map
MAT = [0.436052025 0.385081593 0.143087414
       0.222491598 0.716886060 0.060621486
       0.013929122 0.097097002 0.714185470];  % D50
   
XYZ = MAT * RGB;
XYZ = XYZ';

% normalize for white point
wp = [0.964296 1.0 0.825106];         % D50

XYZ = bsxfun(@rdivide, XYZ, wp);

% cube root normalized xyz
fxyz_n = XYZ .^ (1/3);
% if normalized x, y, or z less than or equal to 216 / 24389 apply function 2  
L = XYZ <= 216 / 24389;
% function 2
k = 841 / 108;
fxyz_n(L) = k * XYZ(L) + 16/116;
clear XYZ L;

LAB = zeros(s,3);
% calculate L*  
LAB(:,1) = 116 * fxyz_n(:,2) - 16;
% calculate a*  
LAB(:,2) = 500 * (fxyz_n(:,1) - fxyz_n(:,2));
% calculate b*  
LAB(:,3) = 200 * (fxyz_n(:,2) - fxyz_n(:,3));

% do the scale and offset and cast to uint8
LAB = round([(255 * (LAB(:,1)/100)) LAB(:,2)+128 LAB(:,3)+128]);
LAB = max(0, min(255, LAB));
LAB = uint8(LAB);

L = reshape(LAB(:,1), M, N);
a = reshape(LAB(:,2), M, N);
b = reshape(LAB(:,3), M, N);

if nargout < 2
    L = cat(3, L, a, b);
end

 

function [L, a, b] = sRGB2LabD65(sR, sG, sB)
% SRGB2LAB convert sRGB to CIELAB(Matlab implemention).
% 
if nargin == 1
    sB = double(sR(:,:,3));
    sG = double(sR(:,:,2));
    sR = double(sR(:,:,1));
end

% normalize for [0 1]
if max(max(sR)) > 1.0 || max(max(sG)) > 1.0 || max(max(sB)) > 1.0
    sR = double(sR) / 255;
    sG = double(sG) / 255;
    sB = double(sB) / 255;
end

[M, N] = size(sR);
s = M * N;
sRL = reshape(sR, 1, s);
sGL = reshape(sG, 1, s);
sBL = reshape(sB, 1, s);

% gamma correction (gamma 2.2)
gammaConst = 0.04045; 
sRL = ( sRL<=gammaConst ) .* sRL  ./ 12.92 +  ...
      ( ( sRL>gammaConst  ) .* (0.055 + sRL) ./ 1.055 )  .^2.4;
sGL = ( sGL<=gammaConst ) .* sGL  ./ 12.92 +  ...
      ( ( sGL>gammaConst  ) .* (0.055 + sGL) ./ 1.055 )  .^2.4;
sBL = ( sBL<=gammaConst ) .* sBL  ./ 12.92 +  ...
      ( ( sBL>gammaConst  ) .* (0.055 + sBL) ./ 1.055 )  .^2.4;
  
RGB = [sRL; sGL; sBL];

MAT = [0.412453 0.357580 0.180423;
       0.212671 0.715160 0.072169;
       0.019334 0.119193 0.950227];   % D65
   
XYZ = MAT * RGB;
XYZ = XYZ';

% normalize for white point
wp = [0.950456 1.0 1.088754];         % D65

XYZ = bsxfun(@rdivide, XYZ, wp);

% cube root normalized xyz
fxyz_n = XYZ .^ (1/3);
% if normalized x, y, or z less than or equal to 216 / 24389 apply function 2  
L = XYZ <= 216 / 24389;
% function 2
k = 841 / 108;
fxyz_n(L) = k * XYZ(L) + 16/116;
clear XYZ L;

LAB = zeros(s,3);
% calculate L*  
LAB(:,1) = 116 * fxyz_n(:,2) - 16;
% calculate a*  
LAB(:,2) = 500 * (fxyz_n(:,1) - fxyz_n(:,2));
% calculate b*  
LAB(:,3) = 200 * (fxyz_n(:,2) - fxyz_n(:,3));

% do the scale and offset and cast to uint8
LAB = round([(255 * (LAB(:,1)/100)) LAB(:,2)+128 LAB(:,3)+128]);
LAB = max(0, min(255, LAB));
LAB = uint8(LAB);

L = reshape(LAB(:,1), M, N);
a = reshape(LAB(:,2), M, N);
b = reshape(LAB(:,3), M, N);

% keyboard
if nargout < 2
    L = cat(3, L, a, b);
end

测试程序：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc; clear all;
close all;

% initial image
N = 100;
I = uint8( zeros(N, N, 3) );
level = 0 : N-1;
for ii = 1: length( level )
    I(ii, :, 1) =  circshift(level, [0 ii]);
    I(ii, :, 2) =  circshift(level, [0 -ii]);
    I(ii, :, 3) =  circshift(level, [0 ii+2]);
end

% rgb2lab (matlab built-in function)
sRGB2LabForm = makecform('srgb2lab');
labImageBuiltIn = applycform(I, sRGB2LabForm);
figure('Name', 'Matlab built-in function');  imshow(labImageBuiltIn);

% rgb to lab function 
labImageForMatlab = sRGB2LabD50(I);
figure('Name', 'One step function'); imshow(labImageBuiltIn);

% difference of two version
diffValue = double(labImageBuiltIn) - double(labImageForMatlab);
h = figure('Name', 'Difference of two version', 'NumberTitle', 'off'); 
subplot(2, 2, 1);  mesh(diffValue(:, :, 1));  title('Difference L'); view([0 0]);
subplot(2, 2, 2);  mesh(diffValue(:, :, 2));  title('Difference A'); view([0 0]);
subplot(2, 2, [3 4]);  mesh(diffValue(:, :, 3));  title('Difference B'); view([0 0]);

 

实验结果：

 

 

                                    (a) matlab自带函数的sRGB转Lab                                                                        (c) sRGB转Lab D50实现

                                                                                                             

                                   (c) sRGB转Lab D65实现                                                                             (d) 对比结果 D50实现

 

 

(5) 对比结果 D65实现

 

从实验结果可以发现，D50的实现与Matlab自带的RGB转Lab差别都在1像素之内(这可以理解为中间计算误差)，而D65的实现与Matlab自带的误差很明显。

参考资料：

Principles of Digital Image Processing: Fundamental Techniques.  Wilhelm Burger, Mark J. Burge