Quantify colour differences when converting between RGB and YCbCR, as well as SDR and HDR

colourspace_differences.py

#!/usr/bin/env python3

# Script output:
#
#     ΔE (CAM16-UCS) [50, 95, 99, 99.9] percentiles
#
#     sRGB (8-bit full) → sRGB (8-bit full)
#     (For reference only)
#     Quantisation error:   [ 0.11  0.23  0.32  0.44]
#     ±1 1ch corruption:    [ 0.16  0.43  0.53  0.7 ]
#     ±1 3ch corruption:    [ 0.34  0.66  0.88  1.15]
#
#     sRGB (8-bit full) → BT.709 YCbCr (below) → sRGB (8-bit full)
#     (e.g., SDR screen capture played back on SDR display)
#      8-bit limited:       [ 0.24  0.46  0.62  0.86]
#      8-bit full:          [ 0.2   0.42  0.53  0.7 ]
#     10-bit limited:       [ 0.  0.  0.  0.]
#
#     sRGB (8-bit full) → BT.2100 PQ YCbCr (below) → sRGB (8-bit full)
#     (e.g., SDR screen capture in HDR format played back on SDR display)
#     10-bit limited:       [ 0.01  0.06  0.1   0.17]
#     12-bit limited:       [ 0.    0.01  0.03  0.05]
#     14-bit limited:       [ 0.    0.    0.    0.02]
#
#     sRGB (8-bit full) → BT.2100 PQ (below)
#     (e.g., SDR content on HDR display)
#     RGB 10-bit full:      [ 0.09  0.2   0.27  0.4 ]
#     RGB 12-bit full:      [ 0.02  0.05  0.07  0.09]
#     RGB 14-bit full:      [ 0.01  0.01  0.02  0.02]
#     YCbCr 10-bit limited: [ 0.16  0.32  0.43  0.61]
#     YCbCr 12-bit limited: [ 0.04  0.08  0.11  0.16]
#     YCbCr 14-bit limited: [ 0.01  0.02  0.03  0.04]
#
#     HDR10 RGB (full) → BT.2100 PQ YCbCr (below) → HDR10 RGB (full)
#     (e.g., HDR screen capture played back on HDR display)
#     10-bit limited:       [ 0.15  0.38  0.52  0.89]
#     12-bit limited:       [ 0.  0.  0.  0.]
#     14-bit limited:       [ 0.  0.  0.  0.]
#
# Observations:
#  - For lossless recording of sRGB in YCbCr, use 10 bits
#  - For lossless recording of HDR10 in YCbCr, use 12 bits
#  - The conversion loss in (sRGB → 8-bit YCbCr) and (HDR10 → 10-bit YCbCr) is
#    very similar, and both are undetectable by a human (<1 ΔE)
#  - The conversion loss in sRGB → HDR10 RGB (SDR → HDR) is about the same as
#    the quantisation error in sRGB itself

# pip install colour-science
import colour, numpy, random

ITERATIONS = 10000
PERCENTILES = [50, 95, 99, 99.9]
SDR_MAX_NITS = 100

print(f"ΔE (CAM16-UCS) {PERCENTILES} percentiles")

def sRGB_deltaE(srgb1, srgb2):
    xyz1 = colour.sRGB_to_XYZ(srgb1 / 255)
    xyz2 = colour.sRGB_to_XYZ(srgb2 / 255)
    jab1 = colour.XYZ_to_CAM16UCS(xyz1)
    jab2 = colour.XYZ_to_CAM16UCS(xyz2)
    return colour.delta_E(jab1, jab2, "CAM16-UCS")

def print_deltaEs(title, deltaEs):
    title += ":"
    print(f"{title:<21s} {numpy.round(numpy.percentile(deltaEs, PERCENTILES), 2)}")

# https://extremelearning.com.au/unreasonable-effectiveness-of-quasirandom-sequences/
phi3 = 1.2207440846057596
alpha = 1 / numpy.array([phi3, phi3**2, phi3**3])
def quasirandom_RGB(n):
    return alpha * n % 1

########################################

print("\nsRGB (8-bit full) → sRGB (8-bit full)")
print("(For reference only)")

# sRGB quantisation error
deltaEs = []
for n in range(ITERATIONS):
    rgb1 = numpy.clip(quasirandom_RGB(n) * 256, 0, 255)
    rgb2 = numpy.round(rgb1)
    deltaEs.append(sRGB_deltaE(rgb1, rgb2))
print_deltaEs("Quantisation error", deltaEs)

# sRGB → ±1 1ch corruption
deltaEs = []
for n in range(ITERATIONS):
    rgb1 = numpy.floor(quasirandom_RGB(n) * 256)
    rgb2 = numpy.copy(rgb1)
    ch = random.randint(0, 2)
    rgb2[ch] = numpy.clip(rgb2[ch] + random.randint(-1, 1), 0, 255)
    deltaEs.append(sRGB_deltaE(rgb1, rgb2))
print_deltaEs("±1 1ch corruption", deltaEs)

# sRGB → ±1 3ch corruption
deltaEs = []
for n in range(ITERATIONS):
    rgb1 = numpy.floor(quasirandom_RGB(n) * 256)
    rgb2 = numpy.clip(rgb1 + numpy.floor(numpy.random.random_sample(3) * 3) - 1, 0, 255)
    deltaEs.append(sRGB_deltaE(rgb1, rgb2))
print_deltaEs("±1 3ch corruption", deltaEs)

########################################

print("\nsRGB (8-bit full) → BT.709 YCbCr (below) → sRGB (8-bit full)")
print("(e.g., SDR screen capture played back on SDR display)")

def sRGB_to_BT709_to_sRGB_deltaEs(bits, *, limited):
    deltaEs = []
    for n in range(ITERATIONS):
        rgb1 = numpy.floor(quasirandom_RGB(n) * 256)
        yuv = colour.RGB_to_YCbCr(
            rgb1, K=colour.WEIGHTS_YCBCR["ITU-R BT.709"],
            in_bits=8, in_int=True, in_legal=False,
            out_bits=bits, out_int=True, out_legal=limited,
        )
        rgb2 = colour.YCbCr_to_RGB(
            yuv, K=colour.WEIGHTS_YCBCR["ITU-R BT.709"],
            in_bits=bits, in_int=True, in_legal=limited,
            out_bits=8, out_int=True, out_legal=False,
        )
        deltaEs.append(sRGB_deltaE(rgb1, rgb2))
    return deltaEs

# sRGB → BT.709 (yuv444p) → sRGB
print_deltaEs(" 8-bit limited", sRGB_to_BT709_to_sRGB_deltaEs(8, limited=True))

# sRGB → BT.709 (yuvj444p) → sRGB
print_deltaEs(" 8-bit full", sRGB_to_BT709_to_sRGB_deltaEs(8, limited=False))

# sRGB → BT.709 (yuv444p10le) → sRGB
print_deltaEs("10-bit limited", sRGB_to_BT709_to_sRGB_deltaEs(10, limited=True))

########################################

print("\nsRGB (8-bit full) → BT.2100 PQ YCbCr (below) → sRGB (8-bit full)")
print("(e.g., SDR screen capture in HDR format played back on SDR display)")

def sRGB_to_BT2100PQ_to_sRGB_deltaEs(bits):
    deltaEs = []
    for n in range(ITERATIONS):
        rgb1 = numpy.floor(quasirandom_RGB(n) * 256) / 255
        rgb2 = colour.RGB_to_RGB(rgb1, "sRGB", "ITU-R BT.2020", apply_cctf_decoding=True, apply_cctf_encoding=True)
        rgb2 = colour.eotf_inverse(SDR_MAX_NITS * colour.eotf(rgb2, "ITU-R BT.1886"), "ITU-R BT.2100 PQ")

        yuv2 = colour.RGB_to_YCbCr(
            rgb2, K=colour.WEIGHTS_YCBCR["ITU-R BT.2020"],
            in_int=False, in_legal=False,
            out_bits=bits, out_int=True, out_legal=True,
        )
        rgb2 = colour.YCbCr_to_RGB(
            yuv2, K=colour.WEIGHTS_YCBCR["ITU-R BT.2020"],
            in_bits=bits, in_int=True, in_legal=True,
            out_int=False, out_legal=False,
        )

        rgb2 = colour.eotf_inverse(colour.eotf(rgb2, "ITU-R BT.2100 PQ") / SDR_MAX_NITS, "ITU-R BT.1886")
        rgb2 = colour.RGB_to_RGB(rgb2, "ITU-R BT.2020", "sRGB", apply_cctf_decoding=True, apply_cctf_encoding=True)
        rgb2 = numpy.clip(numpy.round(rgb2 * 255), 0, 255) / 255

        deltaEs.append(sRGB_deltaE(rgb1, rgb2))
    return deltaEs

# sRGB → BT.2100 PQ YCbCr (yuv444p10le) → sRGB
print_deltaEs("10-bit limited", sRGB_to_BT2100PQ_to_sRGB_deltaEs(10))

# sRGB → BT.2100 PQ YCbCr (yuv444p12le) → sRGB
print_deltaEs("12-bit limited", sRGB_to_BT2100PQ_to_sRGB_deltaEs(12))

# sRGB → BT.2100 PQ YCbCr (yuv444p14le) → sRGB
print_deltaEs("14-bit limited", sRGB_to_BT2100PQ_to_sRGB_deltaEs(14))

########################################

print("\nsRGB (8-bit full) → BT.2100 PQ (below)")
print("(e.g., SDR content on HDR display)")

def BT2100PQ_to_XYZ(rgb):
    rgb = colour.eotf_inverse(colour.eotf(rgb, "ITU-R BT.2100 PQ") / SDR_MAX_NITS, "ITU-R BT.1886")
    return colour.RGB_to_XYZ(rgb, "ITU-R BT.2020", apply_cctf_decoding=True)

def sRGB_to_BT2100PQ_deltaEs(bits, *, ycbcr):
    deltaEs = []
    for n in range(ITERATIONS):
        rgb1 = numpy.floor(quasirandom_RGB(n) * 256) / 255
        rgb2 = colour.RGB_to_RGB(rgb1, "sRGB", "ITU-R BT.2020", apply_cctf_decoding=True, apply_cctf_encoding=True)
        rgb2 = colour.eotf_inverse(SDR_MAX_NITS * colour.eotf(rgb2, "ITU-R BT.1886"), "ITU-R BT.2100 PQ")

        if ycbcr:
            yuv2 = colour.RGB_to_YCbCr(
                rgb2, K=colour.WEIGHTS_YCBCR["ITU-R BT.2020"],
                in_int=False, in_legal=False,
                out_bits=bits, out_int=True, out_legal=True,
            )
            rgb2 = colour.YCbCr_to_RGB(
                yuv2, K=colour.WEIGHTS_YCBCR["ITU-R BT.2020"],
                in_bits=bits, in_int=True, in_legal=True,
                out_int=False, out_legal=False,
            )
        else:
            max = (1 << bits) - 1
            rgb2 = numpy.clip(numpy.round(rgb2 * max), 0, max) / max

        xyz1 = colour.sRGB_to_XYZ(rgb1)
        xyz2 = BT2100PQ_to_XYZ(rgb2)
        jab1 = colour.XYZ_to_CAM16UCS(xyz1)
        jab2 = colour.XYZ_to_CAM16UCS(xyz2)
        deltaEs.append(colour.delta_E(jab1, jab2, "CAM16-UCS"))
    return deltaEs

# sRGB → BT.2100 PQ RGB (gbrp10le)
print_deltaEs("RGB 10-bit full", sRGB_to_BT2100PQ_deltaEs(10, ycbcr=False))

# sRGB → BT.2100 PQ RGB (gbrp12le)
print_deltaEs("RGB 12-bit full", sRGB_to_BT2100PQ_deltaEs(12, ycbcr=False))

# sRGB → BT.2100 PQ RGB (gbrp14le)
print_deltaEs("RGB 14-bit full", sRGB_to_BT2100PQ_deltaEs(14, ycbcr=False))

# sRGB → BT.2100 PQ YCbCr (yuv444p10le)
print_deltaEs("YCbCr 10-bit limited", sRGB_to_BT2100PQ_deltaEs(10, ycbcr=True))

# sRGB → BT.2100 PQ YCbCr (yuv444p12le)
print_deltaEs("YCbCr 12-bit limited", sRGB_to_BT2100PQ_deltaEs(12, ycbcr=True))

# sRGB → BT.2100 PQ YCbCr (yuv444p14le)
print_deltaEs("YCbCr 14-bit limited", sRGB_to_BT2100PQ_deltaEs(14, ycbcr=True))

########################################

print("\nHDR10 RGB (full) → BT.2100 PQ YCbCr (below) → HDR10 RGB (full)")
print("(e.g., HDR screen capture played back on HDR display)")

def HDR10_to_BT2100PQ_to_HDR10_deltaEs(bits):
    deltaEs = []
    for n in range(ITERATIONS):
        rgb1 = numpy.floor(quasirandom_RGB(n) * 1024)

        yuv2 = colour.RGB_to_YCbCr(
            rgb1, K=colour.WEIGHTS_YCBCR["ITU-R BT.2020"],
            in_bits=10, in_int=True, in_legal=False,
            out_bits=bits, out_int=True, out_legal=True,
        )
        rgb2 = colour.YCbCr_to_RGB(
            yuv2, K=colour.WEIGHTS_YCBCR["ITU-R BT.2020"],
            in_bits=bits, in_int=True, in_legal=True,
            out_bits=10, out_int=True, out_legal=False,
        )

        xyz1 = BT2100PQ_to_XYZ(rgb1 / 1023)
        xyz2 = BT2100PQ_to_XYZ(rgb2 / 1023)
        jab1 = colour.XYZ_to_CAM16UCS(xyz1)
        jab2 = colour.XYZ_to_CAM16UCS(xyz2)
        deltaEs.append(colour.delta_E(jab1, jab2, "CAM16-UCS"))
    return deltaEs

# HDR10 → BT.2100 PQ YCbCr (yuv444p10le) → HDR10
print_deltaEs("10-bit limited", HDR10_to_BT2100PQ_to_HDR10_deltaEs(10))

# HDR10 → BT.2100 PQ YCbCr (yuv444p12le) → HDR10
print_deltaEs("12-bit limited", HDR10_to_BT2100PQ_to_HDR10_deltaEs(12))

# HDR10 → BT.2100 PQ YCbCr (yuv444p14le) → HDR10
print_deltaEs("14-bit limited", HDR10_to_BT2100PQ_to_HDR10_deltaEs(14))