About RGB Colourspace Models Performance

Introduction

The purpose of this document is to compare RGB colourspaces operations performance against spectral reference.

In a RGB colourspace, the basis vectors are defined by the primaries. Mathematical operations are dependent on the chosen basis vectors. As an example, a triplet of tristimulus values, converted to three different RGB colourspaces value, subsequently exponentiated, and converted back to CIE XYZ colourpace will produce different tristimulus values triplets, one for each of the three RGB colourspaces. In this context, the result of the multiplication, division and power operations is dependent on the RGB colourspace primaries, producing three different final tristimulus values while addition and subtraction operations would produce three equal tristimulus values.

The rgb_colourspace_primaries_dependency definition below demonstrates this behaviour: Two colours are selected from the colour rendition chart, they are then used as operands for different operations happening respectively in sRGB and Rec. 2020 colourspaces. The two resulting colours are then converted back to CIE XYZ colourspace and compared.

Note: A companion spreadsheet is available at the following url: http://goo.gl/akrH5m

In [1]:
import operator
from pprint import pprint

import colour


def rgb_colourspace_primaries_dependency(operator=operator.add):
    name, data, illuminant = colour.COLOURCHECKERS['ColorChecker 2005']

    # Aliases for *sRGB* colourspace.
    sRGB_w = colour.sRGB_COLOURSPACE.whitepoint
    sRGB_XYZ_to_RGB = colour.sRGB_COLOURSPACE.XYZ_to_RGB_matrix
    sRGB_RGB_to_XYZ = colour.sRGB_COLOURSPACE.RGB_to_XYZ_matrix

    # Aliases for *Rec. 2020* colourspace.
    REC_2020_w = colour.REC_2020_COLOURSPACE.whitepoint
    REC_2020_XYZ_to_RGB = colour.REC_2020_COLOURSPACE.XYZ_to_RGB_matrix
    REC_2020_RGB_to_XYZ = colour.REC_2020_COLOURSPACE.RGB_to_XYZ_matrix

    # Selecting various colours from the colour rendition chart:
    # *Dark Skin*
    _, _, x, y, Y = data[0]
    XYZ_r1 = colour.xyY_to_XYZ((x, y, Y))
    # *Green*
    _, _, x, y, Y = data[13]
    XYZ_r2 = colour.xyY_to_XYZ((x, y, Y))

    # Defining the *sRGB* *Dark Skin* colour.
    sRGB_r1 = colour.XYZ_to_RGB(XYZ_r1,
                                illuminant,
                                sRGB_w,
                                sRGB_XYZ_to_RGB)

    # Defining the *sRGB* *Green* colour.
    sRGB_r2 = colour.XYZ_to_RGB(XYZ_r2,
                                illuminant,
                                sRGB_w,
                                sRGB_XYZ_to_RGB)

    # Performing *sRGB* *Dark Skin* and *Green* colours operation.
    sRGB_m = operator(sRGB_r1, sRGB_r2)

    # Converting resulting colour to *CIE* *XYZ* colourspace.
    XYZ_sRGB_m1 = colour.RGB_to_XYZ(sRGB_m,
                                    sRGB_w,
                                    illuminant,
                                    sRGB_RGB_to_XYZ)

    # Defining the *Rec. 2020* *Dark Skin* colour.
    REC_2020_r1 = colour.XYZ_to_RGB(XYZ_r1,
                                    illuminant,
                                    REC_2020_w,
                                    REC_2020_XYZ_to_RGB)

    # Defining the *Rec. 2020* *Green* colour.
    REC_2020_r2 = colour.XYZ_to_RGB(XYZ_r2,
                                    illuminant,
                                    REC_2020_w,
                                    REC_2020_XYZ_to_RGB)

    # Performing the *Rec. 2020* *Dark Skin* and *Green* colours operation.
    REC_2020_m = operator(REC_2020_r1, REC_2020_r2)

    # Converting resulting colour to *CIE* *XYZ* colourspace.
    XYZ_REC_2020_m1 = colour.RGB_to_XYZ(REC_2020_m,
                                        REC_2020_w,
                                        illuminant,
                                        REC_2020_RGB_to_XYZ)

    return (XYZ_sRGB_m1, XYZ_REC_2020_m1)


for operator in (operator.add,
                 operator.sub,
                 operator.mul,
                 operator.div,
                 operator.pow):
    print('{0}:'.format(operator.__name__))
    pprint(rgb_colourspace_primaries_dependency(operator))
    print('\n')
add:
(array([ 0.26503479,  0.3326    ,  0.12989551]),
 array([ 0.26503479,  0.3326    ,  0.12989551]))


sub:
(array([-0.03466529, -0.131     , -0.02810806]),
 array([-0.03466529, -0.131     , -0.02810806]))


mul:
(array([ 0.01419744,  0.02001609,  0.00525311]),
 array([ 0.01744617,  0.02218764,  0.00496058]))


div:
(array([ 1.60211763,  0.94783487,  0.64260732]),
 array([ 0.81894892,  0.5243635 ,  0.54464893]))


pow:
(array([ 0.69403138,  0.58944147,  0.64315778]),
 array([ 0.69452189,  0.58949601,  0.63399593]))


Methodology

We have decided to test multiplication operation because it is affected by RGB colourspace primaries dependencies and especially because it is easily possible to design synthetic spectral power distributions that can be applied on existing spectral data, like the one from a colour rendition chart and also on the same data converted to RGB upon conversion to tristimulus values.

In early versions of this document, the synthetic spectral power distributions were named filters as they were designed to mimic existing filter sets from real life. However we discovered that referencing existing real life filters was introducing an important bias in the results. We were basically finding the best colourspace for a given set of filters thus we decided to change the filters design so that we have a very large amount of different synthetic spectral power distributions.

The methodology implementation code is available in the Process section and is illustrated in the following figure:

In [2]:
from IPython.core.display import Image

Image(filename="resources/images/About Colourspaces & Colour Rendition Charts - Overall.png")
Out[2]:

The spectral test data is coming from N. Ohta Colour Rendition Chart Classic and X-Rite Colour Rendition Chart SG colour rendition charts.

The filters design will be detailed in the Filters Design section but let's assume that we have a large number of random spectral power distributions providing coverage of the spectral locus.

Picking a test sample, a random filter and a RGB colourspace model, we follow two paths:

  • The first path is the spectral path where the test sample is weighted by the random filter and then converted to tristimulus values under the RGB colourspace model illuminant. The tristimulus values are then converted to CIE Lab and labeled as the Reference Lab matrix.
  • The second path is the RGB path where the test sample and the random filter are converted to tristimulus values under the RGB colourspace model illuminant, then respectively converted to RGB. The test sample RGB triplet is then weighted by the random filter RGB triplet. The resulting RGB triplet is also converted to CIE Lab and labeled as the Test Lab matrix.

$\Delta_{E_{ab}}$ CIE 2000 colour difference is then computed with the Reference Lab and Test Lab matrices.

This process is repeated for each sample, each random filter and each RGB colourspace. The resulting $\Delta_{E_{ab}}$ values are averaged together to get an overall value for each RGB colourspace. We also take that opportunity to compute standard deviation $\sigma$ and min / max $\Delta_{E_{ab}}$ values for each RGB colourspace.

Filters Design

In order to improve upon earlier versions, the filters are now quasi-randomly generated. We have defined 3 base distributions as follows:

  • Normal distribution: A basic normal distribution that can be centered on an arbitrary wavelength.
  • Bowl distribution: A distribution with the shape of a bowl and that can be centered on an arbitrary wavelength. Its purpose is to provide coverage for the spectral locus line of purples that cannot be reached using a normal distribution.
  • Random Spline distribution: This distribution is used to generate random splines distribution with very low points count.

From the base distributions the filters are designed:

  • The first step is to generate the initial filters using the Normal and Bowl distributions, the wavelengths have been cherry picked to provide as much as possible an evenly spaced coverage of the spectral locus boundaries as seen in the Initial Spectral Power Distributions section.
  • We then start designing the final filters: the initial filters are randomly offset along the spectrum and averaged with a Random Spline distribution using ratio 0.85 / 0.15 in order to provide variation.
  • The resulting set of filters is combined with a pool of pure Random Spline distributions with a 0.5 / 0.5 ratio, we label that pool of filters Random Filters.
  • The Random Filters are then saturated to regain spectral locus coverage. The ratio is 1/3 default Random Filters, 1/3 square saturated Random Filters, 1/3 cubic saturated Random Filters.
  • The Random Filters are checked for uniqueness and all-most-equal ones are discarded.

The Random Filters can be seen in the Random Filters - Design section and their design process is illustrated in the below figure:

In [3]:
Image(filename="resources/images/About Colourspaces & Colour Rendition Charts - Filters Design.png")
Out[3]:

Implementation

Colour Rendition Charts

In [4]:
%matplotlib inline

import logging
reload(logging)

import matplotlib.pyplot
import numpy as np
import os
import pylab
import re
from IPython.display import HTML
from mpl_toolkits.mplot3d import Axes3D, art3d
from copy import deepcopy

from colour.characterisation.dataset.colour_checkers.spds import (
    COLOURCHECKER_INDEXES_TO_NAMES_MAPPING)
from colour.plotting import *

# Logging configuration.
LOGGER = logging.getLogger('colour-science')
LOGGER.setLevel(logging.DEBUG)

LOGGING_FILE = 'about_colourspaces_colour_rendition_charts.log'

os.path.exists(LOGGING_FILE) and os.remove(LOGGING_FILE)

HANDLER = logging.FileHandler(LOGGING_FILE)
HANDLER.setFormatter(logging.Formatter('%(asctime)s - %(levelname)-8s: %(message)s'))

LOGGER.addHandler(HANDLER) 

LOGGER.info('*' * 79)
LOGGER.info('About Colourspaces & Colour Rendition Charts')
LOGGER.info('*' * 79)

pylab.rcParams['figure.figsize'] = 28, 14

CMFS = colour.CMFS['CIE 1931 2 Degree Standard Observer']
SHAPE = CMFS.shape


def read_xrite_spds(path):
    xrite_reflectances_file = open(path)
    xrite_reflectances_data = xrite_reflectances_file.read().strip().split(
        '\n')

    xrite_spds = []
    is_spectral_data_format, is_spectral_data = False, False
    for line in xrite_reflectances_data:
        line = line.strip()

        if line == 'END_DATA_FORMAT':
            is_spectral_data_format = False

        if line == 'END_DATA':
            is_spectral_data = False

        if is_spectral_data_format:
            wavelengths = [float(x) for x in re.findall('nm(\d+)', line)]
            index = len(wavelengths)

        if is_spectral_data:
            tokens = line.split()
            values = [float(x) for x in tokens[-index:]]
            xrite_spds.append(colour.SpectralPowerDistribution(tokens[1],
                                                               dict(zip(
                                                                   wavelengths,
                                                                   values))))

        if line == 'BEGIN_DATA_FORMAT':
            is_spectral_data_format = True

        if line == 'BEGIN_DATA':
            is_spectral_data = True
    return xrite_spds


COLOUR_CHECKER = 'SG'

LOGGER.info('Colour checker: \'{0}\''.format(COLOUR_CHECKER))

if COLOUR_CHECKER == 'Classic':
    # Preparing N. Ohta Colour Rendition Chart Classic (24 samples).
    COLOUR_CHECKER = deepcopy(
        colour.COLOURCHECKERS_SPDS['ColorChecker N Ohta'])
    COLOUR_CHECKER_MAPPING = COLOURCHECKER_INDEXES_TO_NAMES_MAPPING
elif COLOUR_CHECKER == 'SG':
    # Preparing X-Rite Colour Rendition Chart SG (140 samples).
    xrite_spds = read_xrite_spds(
        'resources/others/Digital_ColorChecker_SG.txt')
    xrite_spds_names = [xrite_spd.name for xrite_spd in xrite_spds]
    xrite_spds_indexes = range(len(xrite_spds_names))
    COLOUR_CHECKER = dict(zip(xrite_spds_names, xrite_spds))
    COLOUR_CHECKER_MAPPING = dict(zip(xrite_spds_indexes, xrite_spds_names))

[spd.align(SHAPE) for spd in COLOUR_CHECKER.values()]


def get_colour_checker_spds(colour_checker=COLOUR_CHECKER,
                            mapping=COLOUR_CHECKER_MAPPING):
    spds = []
    for index, sample in sorted(mapping.items()):
        spds.append(colour_checker[sample].clone())
    return spds


# TODO: Merge in *colour*.
def spds_CIE_1931_chromaticity_diagram_plot(
        spds,
        cmfs='CIE 1931 2 Degree Standard Observer',
        annotate=True,
        **kwargs):
    CIE_1931_chromaticity_diagram_plot(standalone=False,
                                       **kwargs)

    cmfs, name = get_cmfs(cmfs), cmfs
    cmfs_shape = cmfs.shape

    for spd in spds:
        spd = spd.clone().align(cmfs_shape)
        XYZ = colour.spectral_to_XYZ(spd) / 100
        xy = colour.XYZ_to_xy(XYZ)

        pylab.plot(xy[0], xy[1], 'o', color='white')

        if spd.name is not None and annotate:
            pylab.annotate(spd.name,
                           xy=xy,
                           xytext=(50, 30),
                           textcoords='offset points',
                           arrowprops=dict(arrowstyle='->',
                                           connectionstyle='arc3, rad=0.2'))

    display(standalone=True)


def spds_CIE_1960_UCS_chromaticity_diagram_plot(
        spds,
        cmfs='CIE 1931 2 Degree Standard Observer',
        annotate=True,
        **kwargs):
    CIE_1960_UCS_chromaticity_diagram_plot(standalone=False,
                                           **kwargs)

    cmfs, name = get_cmfs(cmfs), cmfs
    cmfs_shape = cmfs.shape

    for spd in spds:
        spd = spd.clone().align(cmfs_shape)
        XYZ = colour.spectral_to_XYZ(spd) / 100
        uv = colour.UCS_to_uv(colour.XYZ_to_UCS(XYZ))

        pylab.plot(uv[0], uv[1], 'o', color='white')

        if spd.name is not None and annotate:
            pylab.annotate(spd.name,
                           xy=uv,
                           xytext=(50, 30),
                           textcoords='offset points',
                           arrowprops=dict(arrowstyle='->',
                                           connectionstyle='arc3, rad=0.2'))

    display(standalone=True)


def spds_CIE_1976_UCS_chromaticity_diagram_plot(
        spds,
        cmfs='CIE 1931 2 Degree Standard Observer',
        annotate=True,
        **kwargs):
    CIE_1976_UCS_chromaticity_diagram_plot(standalone=False,
                                           **kwargs)

    cmfs, name = get_cmfs(cmfs), cmfs
    cmfs_shape = cmfs.shape

    for spd in spds:
        spd = spd.clone().align(cmfs_shape)
        XYZ = colour.spectral_to_XYZ(spd) / 100
        uv = colour.Luv_to_uv(colour.XYZ_to_Luv(XYZ))

        pylab.plot(uv[0], uv[1], 'o', color='white')

        if spd.name is not None and annotate:
            pylab.annotate(spd.name,
                           xy=uv,
                           xytext=(50, 30),
                           textcoords='offset points',
                           arrowprops=dict(arrowstyle='->',
                                           connectionstyle='arc3, rad=0.2'))

    display(standalone=True)


multi_spd_plot(COLOUR_CHECKER.values(), legend=False, use_spds_colours=True)

spds_CIE_1931_chromaticity_diagram_plot(COLOUR_CHECKER.values(),
                                        annotate=False)