Replicate Photoshop sRGB to LAB Conversion

Question

The task I want to achieve is to replicate Photoshop RGB to LAB conversion.
For simplicity, I will describe what I did to extract only the L Channel.

Extracting Photoshop's L Channel

Here is RGB Image which includes all RGB colors (Please click and download):

In order to extract Photoshop's LAB what I did is the following:

1. Loaded the image into Photoshop.
2. Set Mode to LAB.
3. Selected the L Channel in the Channel Panel.
4. Set Mode to Grayscale.
5. Set mode to RGB.
6. Saved as PNG.

This is the L Channel of Photoshop (This is exactly what seen on screen when L Channel is selected in LAB Mode):

sRGB to LAB Conversion

My main reference is Bruce Lindbloom great site.
Also known is that Photoshop is using D50 White Point in its LAB Mode (See also Wikipedia's LAB Color Space Page).

Assuming the RGB image is in sRGB format the conversion is given by:

sRGB -> XYZ (White Point D65) -> XYZ (White Point D50) -> LAB

Assuming data is in Float within the [0, 1] range the stages are given by:

1. Transform sRGB into XYZ.
   The conversion Matrix is given by RGB -> XYZ Matrix (See sRGB D65).
2. Converting from XYZ D65 to XYZ D50
   The conversion is done using Chromatic Adaptation Matrix. Since the previous step and this are Matrix Multiplication they can be combined into one Matrix which goes from sRGB -> XYZ D50 (See the bottom of RGB to XYZ Matrix). Note that Photoshop uses Bradford Adaptation Method.
3. Convert from XYZ D50 to LAB
   The conversion is done using the XYZ to LAB Steps.

MATLAB Code

Since, for start, I'm only after the L Channel things are a bit simpler. The images are loaded into MATLAB and converted into Float [0, 1] range.

This is the code:

%% Setting Enviorment Parameters

INPUT_IMAGE_RGB             = 'RgbColors.png';
INPUT_IMAGE_L_PHOTOSHOP     = 'RgbColorsL.png';


%% Loading Data

mImageRgb   = im2double(imread(INPUT_IMAGE_RGB));
mImageLPhotoshop     = im2double(imread(INPUT_IMAGE_L_PHOTOSHOP));
mImageLPhotoshop     = mImageLPhotoshop(:, :, 1); %<! All channels are identical


%% Convert to L Channel

mImageLMatlab = ConvertRgbToL(mImageRgb, 1);


%% Display Results
figure();
imshow(mImageLPhotoshop);
title('L Channel - Photoshop');

figure();
imshow(mImageLMatlab);
title('L Channel - MATLAB');

Where the function ConvertRgbToL() is given by:

function [ mLChannel ] = ConvertRgbToL( mRgbImage, sRgbMode )

OFF = 0;
ON  = 1;

RED_CHANNEL_IDX     = 1;
GREEN_CHANNEL_IDX   = 2;
BLUE_CHANNEL_IDX    = 3;

RGB_TO_Y_MAT = [0.2225045, 0.7168786, 0.0606169]; %<! D50

Y_CHANNEL_THR = 0.008856;

% sRGB Compensation
if(sRgbMode == ON)
    vLinIdx = mRgbImage < 0.04045;

    mRgbImage(vLinIdx)  = mRgbImage(vLinIdx) ./ 12.92;
    mRgbImage(~vLinIdx) = ((mRgbImage(~vLinIdx) + 0.055) ./ 1.055) .^ 2.4;
end

% RGB to XYZ (D50)
mY = (RGB_TO_Y_MAT(1) .* mRgbImage(:, :, RED_CHANNEL_IDX)) + (RGB_TO_Y_MAT(2) .* mRgbImage(:, :, GREEN_CHANNEL_IDX)) + (RGB_TO_Y_MAT(3) .* mRgbImage(:, :, BLUE_CHANNEL_IDX));

vYThrIdx = mY > Y_CHANNEL_THR;

mY3 = mY .^ (1 / 3);

mLChannel = ((vYThrIdx .* (116 * mY3 - 16.0)) + ((~vYThrIdx) .* (903.3 * mY))) ./ 100;


end

As one could see the results are different.
Photoshop is much darker for most colors.

Anyone knows how to replicate Photoshop's LAB conversion?
Anyone can spot issue in this code?

Thank You.

Answer 1

Latest answer (we know that it is wrong now, waiting for a proper answer)

Photoshop is a very old and messy software. There's no clear documentation as to why this or that happens to the pixel values when you are performing conversions from a mode to another.

Your problem happens because when you are converting the selected L* channel to Greyscale in Adobe Photoshop, there's a change in gamma. Natively, the conversion uses a gamma of 1.74 for single channel to greyscale conversion. Don't ask me why, I would guess this is related to old laser printers (?).

Anyway, this is the best way I found to do it:

Open your file, turn it to LAB mode, select the L channel only

Then go to:

Edit > Convert to profile

You will select "custom gamma" and enter the value 2.0 (don't ask me why 2.0 works better, I have no idea what's in the mind of Adobe's software makers...) This operation will turn your picture into a greyscale one with only one channel

Then you can convert it to RGB mode.

If you compare the result with your result, you will see differences up to 4 dot something % - all located in the darkest areas.

I suspect this is because the gamma curve application does not appy to LAB mode in the dark values (Cf. as you know, all XYZ values below 0.008856 are linear in LAB)

CONCLUSION:

As far as I know, there is no proper implemented way in Adobe Photoshop to extract the L channel from LAB mode to grey mode!

Previous answer

this is the result I get with my own method:

It seems to be exactly the same result as the Adobe Photoshop one.

I am not sure what went wrong on your side since the steps that you are describing are exactly the same ones that I followed and that I would have advised you to follow. I don't have Matlab so I used python:

import cv2, Syn

# your file
fn = "EASA2.png"

#reading the file
im = cv2.imread(fn,-1)

#openCV works in BGR, i'm switching to RGB
im = im[:,:,::-1]

#conversion to XYZ
XYZ = Syn.sRGB2XYZ(im)

#white points D65 and D50
WP_D65 = Syn.Yxy2XYZ((100,0.31271, 0.32902))
WP_D50 = Syn.Yxy2XYZ((100,0.34567, 0.35850))

#bradford
XYZ2 = Syn.bradford_adaptation(XYZ, WP_D65, WP_D50) 

#conversion to L*a*b*
LAB = Syn.XYZ2Lab(XYZ2, WP_D50)

#picking the L channel only
L = LAB[:,:,0] /100. * 255.

#image output
cv2.imwrite("result.png", L)

the Syn library is my own stuff, here are the functions (sorry for the mess):

def sRGB2XYZ(sRGB):

    sRGB = np.array(sRGB)
    aShape = np.array([1,1,1]).shape
    anotherShape = np.array([[1,1,1],[1,1,1]]).shape
    origShape = sRGB.shape

    if sRGB.shape == aShape:
        sRGB = np.reshape(sRGB, (1,1,3))

    elif len(sRGB.shape) == len(anotherShape):
        h,d = sRGB.shape
        sRGB = np.reshape(sRGB, (1,h,d))

    w,h,d = sRGB.shape

    sRGB = np.reshape(sRGB, (w*h,d)).astype("float") / 255.

    m1 = sRGB[:,0] > 0.04045
    m1b = sRGB[:,0] <= 0.04045
    m2 = sRGB[:,1] > 0.04045
    m2b = sRGB[:,1] <= 0.04045
    m3 = sRGB[:,2] > 0.04045
    m3b = sRGB[:,2] <= 0.04045

    sRGB[:,0][m1] = ((sRGB[:,0][m1] + 0.055 ) / 1.055 ) ** 2.4
    sRGB[:,0][m1b] = sRGB[:,0][m1b] / 12.92

    sRGB[:,1][m2] = ((sRGB[:,1][m2] + 0.055 ) / 1.055 ) ** 2.4
    sRGB[:,1][m2b] = sRGB[:,1][m2b] / 12.92

    sRGB[:,2][m3] = ((sRGB[:,2][m3] + 0.055 ) / 1.055 ) ** 2.4
    sRGB[:,2][m3b] = sRGB[:,2][m3b] / 12.92

    sRGB *= 100. 

    X = sRGB[:,0] * 0.4124 + sRGB[:,1] * 0.3576 + sRGB[:,2] * 0.1805
    Y = sRGB[:,0] * 0.2126 + sRGB[:,1] * 0.7152 + sRGB[:,2] * 0.0722
    Z = sRGB[:,0] * 0.0193 + sRGB[:,1] * 0.1192 + sRGB[:,2] * 0.9505

    XYZ = np.zeros_like(sRGB)

    XYZ[:,0] = X
    XYZ[:,1] = Y
    XYZ[:,2] = Z

    XYZ = np.reshape(XYZ, origShape)

    return XYZ

def Yxy2XYZ(Yxy):

    Yxy = np.array(Yxy)
    aShape = np.array([1,1,1]).shape
    anotherShape = np.array([[1,1,1],[1,1,1]]).shape
    origShape = Yxy.shape

    if Yxy.shape == aShape:
        Yxy = np.reshape(Yxy, (1,1,3))

    elif len(Yxy.shape) == len(anotherShape):
        h,d = Yxy.shape
        Yxy = np.reshape(Yxy, (1,h,d))

    w,h,d = Yxy.shape

    Yxy = np.reshape(Yxy, (w*h,d)).astype("float")

    XYZ = np.zeros_like(Yxy)

    XYZ[:,0] = Yxy[:,1] * ( Yxy[:,0] / Yxy[:,2] )
    XYZ[:,1] = Yxy[:,0]
    XYZ[:,2] = ( 1 - Yxy[:,1] - Yxy[:,2] ) * ( Yxy[:,0] / Yxy[:,2] )

    return np.reshape(XYZ, origShape)

def bradford_adaptation(XYZ, Neutral_source, Neutral_destination):
    """should be checked if it works properly, but it seems OK"""

    XYZ = np.array(XYZ)
    ashape = np.array([1,1,1]).shape
    siVal = False

    if XYZ.shape == ashape:


        XYZ = np.reshape(XYZ, (1,1,3))
        siVal = True


    bradford = np.array(((0.8951000, 0.2664000, -0.1614000),
                          (-0.750200, 1.7135000,  0.0367000),
                          (0.0389000, -0.068500,  1.0296000)))

    inv_bradford = np.array(((0.9869929, -0.1470543, 0.1599627),
                              (0.4323053,  0.5183603, 0.0492912),
                              (-.0085287,  0.0400428, 0.9684867)))

    Xs,Ys,Zs = Neutral_source
    s = np.array(((Xs),
                   (Ys),
                   (Zs)))

    Xd,Yd,Zd = Neutral_destination
    d = np.array(((Xd),
                   (Yd),
                   (Zd)))


    source = np.dot(bradford, s)
    Us,Vs,Ws = source[0], source[1], source[2]

    destination = np.dot(bradford, d)
    Ud,Vd,Wd = destination[0], destination[1], destination[2]

    transformation = np.array(((Ud/Us, 0, 0),
                                (0, Vd/Vs, 0),
                                (0, 0, Wd/Ws)))

    M = np.mat(inv_bradford)*np.mat(transformation)*np.mat(bradford)

    w,h,d = XYZ.shape
    result = np.dot(M,np.rot90(np.reshape(XYZ, (w*h,d)),-1))
    result = np.rot90(result, 1)
    result = np.reshape(np.array(result), (w,h,d))

    if siVal == False:
        return result
    else:
        return result[0,0]

def XYZ2Lab(XYZ, neutral):
    """transforms XYZ to CIE Lab
    Neutral should be normalized to Y = 100"""

    XYZ = np.array(XYZ)
    aShape = np.array([1,1,1]).shape
    anotherShape = np.array([[1,1,1],[1,1,1]]).shape
    origShape = XYZ.shape

    if XYZ.shape == aShape:
        XYZ = np.reshape(XYZ, (1,1,3))

    elif len(XYZ.shape) == len(anotherShape):
        h,d = XYZ.shape
        XYZ = np.reshape(XYZ, (1,h,d))

    N_x, N_y, N_z = neutral
    w,h,d = XYZ.shape

    XYZ = np.reshape(XYZ, (w*h,d)).astype("float")

    XYZ[:,0] = XYZ[:,0]/N_x
    XYZ[:,1] = XYZ[:,1]/N_y
    XYZ[:,2] = XYZ[:,2]/N_z

    m1 = XYZ[:,0] > 0.008856
    m1b = XYZ[:,0] <= 0.008856
    m2 = XYZ[:,1] > 0.008856 
    m2b = XYZ[:,1] <= 0.008856
    m3 = XYZ[:,2] > 0.008856
    m3b = XYZ[:,2] <= 0.008856

    XYZ[:,0][m1] = XYZ[:,0][XYZ[:,0] > 0.008856] ** (1/3.0)
    XYZ[:,0][m1b] = ( 7.787 * XYZ[:,0][m1b] ) + ( 16 / 116.0 )

    XYZ[:,1][m2] = XYZ[:,1][XYZ[:,1] > 0.008856] ** (1/3.0)
    XYZ[:,1][m2b] = ( 7.787 * XYZ[:,1][m2b] ) + ( 16 / 116.0 )

    XYZ[:,2][m3] = XYZ[:,2][XYZ[:,2] > 0.008856] ** (1/3.0)
    XYZ[:,2][m3b] = ( 7.787 * XYZ[:,2][m3b] ) + ( 16 / 116.0 )

    Lab = np.zeros_like(XYZ)

    Lab[:,0] = (116. * XYZ[:,1] ) - 16.
    Lab[:,1] = 500. * ( XYZ[:,0] - XYZ[:,1] )
    Lab[:,2] = 200. * ( XYZ[:,1] - XYZ[:,2] )

    return np.reshape(Lab, origShape)