Converts image to CIE 1931 Chromaticity diagram

picture_to_cie_diagram.py

"""
Copyright (C) 2024 Rachel030219

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
"""
import os
import sys
import math
from abc import abstractmethod

import numpy
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import PathPatch
from matplotlib.path import Path
import time
import colour
import imageio.v3 as iio
from colour.plotting import *
import exiftool

# here is the svg path used to draw the horseshoe, resized & flipped to fit 1.0x1.0 mpl axes, orig is 512x512
# src: https://commons.wikimedia.org/wiki/File:CIExy1931.svg
# svg_path = 'M 0.74023,0.26172 C 0.53262,0.46934 0.39414,0.60195 0.30527,0.69082 0.22871,0.76738 0.12168,0.83418 0.08203,0.83418 c -0.05605,-0.00000 -0.07715,-0.08809 -0.07715,-0.19238 0.00000,-0.11582 0.03594,-0.44414 0.12891,-0.59844 0.00000,-0.00000 0.02012,-0.03047 0.04238,-0.03906 L 0.74023,0.26172 z'
horseshoe_path = Path(np.array([[0.74023, 0.26172],
       [0.53262, 0.46934],
       [0.39414, 0.60195],
       [0.30527, 0.69082],
       [0.22871, 0.76738],
       [0.12168, 0.83418],
       [0.08203, 0.83418],
       [0.02598, 0.83418],
       [0.00488, 0.74609],
       [0.00488, 0.6418 ],
       [0.00488, 0.52598],
       [0.04082, 0.19766],
       [0.13379, 0.04336],
       [0.13379, 0.04336],
       [0.15391, 0.01289],
       [0.17617, 0.0043 ],
       [0.74023, 0.26172],
       [0.74023, 0.26172]]), np.array([1,  4,  4,  4,  4,  4,  4,  4,  4,  4,  4,  4,  4,  4,  4,  4,  2, 79],
                                      dtype=np.uint8))

# Classes used to parse ICC profiles.
# Referring to ICC.1:2022 standard, section 10.6 *curveType* and 10.18 *parametricCurveType*
class ConversionFunction:
    """
    Base class for conversion functions.
    It should be noted that in ICC profiles all values are encoded as big-endian.
    Since all conversion methods require a conversion from bytes to int, they should all follow byteorder 'big'.
    """

    @abstractmethod
    def convert_trc(self, pixels: np.ndarray) -> np.ndarray:
        pass


def parse_s15Fixed16Number(s15Fixed16Number: bytes) -> float:
    """
    Parse a s15Fixed16Number, which is a 4-byte signed fixed-point number, with 16 fractional bits.
    The first 16 bits are the integer part, and the last 16 bits are the fractional part.
    i.e. 0x00010000 represents 1, 0x80000000 represents -32768, 0x7FFFFFFF represents 32767 + (65535/65536), etc.
    :param s15Fixed16Number: a 4-byte signed fixed-point number, with 16 fractional bits.
    :return: a float number.
    """
    # first calculate the integer part
    integer_bytes = s15Fixed16Number[0:2]
    integer_part = int.from_bytes(integer_bytes, byteorder='big', signed=True)
    # then calculate the fractional part
    fractional_bytes = s15Fixed16Number[2:4]
    fractional_part = int.from_bytes(fractional_bytes, byteorder='big', signed=False) / np.iinfo(np.uint16).max
    if integer_part < 0:
        fractional_part = -fractional_part
    return integer_part + fractional_part


def parse_trc(trc_param: bytes) -> ConversionFunction:
    """
    Parse the TRC parameter.
    :param trc_param: TRC parameter, a 1D array, representing the sampled one-dimensional function.
    :return: a conversion function.
    """
    if not trc_param:
        return GammaConversion(gamma=2.2)
    trc_type = trc_param[0:4].decode('utf-8')
    if trc_type == 'curv':
        """
        For curv types, according to the standard,
        When n is equal to 0, an identity response is assumed.
        When n is equal to 1, then the curve value shall be interpreted as a gamma value, encoded as a u8Fixed8Number.
         Gamma shall be interpreted as the exponent in the equation y = xγ and not as an inverse.
        When n greater than 1, the curve values (a sampled one-dimensional function) are defined as follows:
         The first entry is located at 0,0,
         the last entry at 1,0,
         intermediate entries are uniformly spaced using an increment of 1,0/(n−1).
         Entries are encoded as uInt16 (values represented by the entries in the range 0,0 to 1,0 are from 0 to 65 535).
         Function values between the entries shall be obtained through linear interpolation.
        """
        n = int.from_bytes(trc_param[8:12], byteorder='big', signed=False)
        if n == 0:
            return GammaConversion(gamma=2.2)
        elif n == 1:
            return GammaConversion(gamma=(trc_param[12] + trc_param[13] / 256))
        elif n > 1:
            matrix = np.zeros(n)
            for i in range(n):
                matrix[i] = (int.from_bytes(trc_param[12 + i * 2: 12 + (i + 1) * 2], byteorder='big', signed=False)
                             / np.iinfo(np.uint16).max)
            return CurveConversion(matrix)
        else:
            raise ValueError(f'n should be greater than 0, got n: {n}')
    elif trc_type == 'para':
        curve_type = int.from_bytes(trc_param[8:10], byteorder='big', signed=False)
        length = math.ceil(len(trc_param[10:]) / 4)
        params = np.zeros(length)
        for i in range(length):
            params[i] = parse_s15Fixed16Number(trc_param[12 + i * 4: 12 + (i + 1) * 4])
        return ParametricConversion(curve_type, params)
    else:
        raise ValueError(f'TRC type {trc_type} not supported')


class ICCProfile:
    """
    Class to parse ICC profile.
    """
    _trc_param = {}
    _matrix = None
    _profile_description = None
    _cat = None
    _pcs_illuminant_XYZ = None
    _pcs_illuminant_conversion_matrix = None
    _matrix_missing = False

    def __init__(self, icc_path: str):
        self.icc_path = icc_path

    def get_trc(self) -> dict:
        """
        Get the TRC function for a given channel.
        Reason for the three return values is that it is unable to guarantee that all channel are of the same type.
        Therefore, the caller should apply conversion on each channel separately.
        :return: three conversion functions, corresponding to R, G, and B channels.
        """
        if not self._trc_param:
            if self.icc_path is not None and os.path.exists(self.icc_path):
                with exiftool.ExifTool() as et:
                    trc_red_param = et.execute(*['-icc_profile:RedTRC', '-b', self.icc_path], raw_bytes=True)
                    trc_green_param = et.execute(*['-icc_profile:GreenTRC', '-b', self.icc_path], raw_bytes=True)
                    trc_blue_param = et.execute(*['-icc_profile:BlueTRC', '-b', self.icc_path], raw_bytes=True)
                    self._trc_param = {'R': parse_trc(trc_red_param), 'G': parse_trc(trc_green_param),
                                       'B': parse_trc(trc_blue_param)}
            else:
                self._trc_param = {'R': GammaConversion(gamma=2.2), 'G': GammaConversion(gamma=2.2),
                                   'B': GammaConversion(gamma=2.2)}
        return self._trc_param

    def get_matrix(self) -> np.ndarray:
        """
        Get the matrix for R, G and B channels. Matrices read from ICC profiles should be regarded as columns.
        This method returns an array like:
            [[rX, gX, bX],
             [rY, gY, bY],
             [rZ, gZ, bZ]]
        Then, numpy's matrix dot product can be used to convert from RGB to XYZ: 'XYZ = matrix dot RGB', where matrix is
        the matrix returned by this method, and RGB is a 2D array with shape of (n, 3).
        In dot product, (A·B)^T = B^T·A^T, therefore, the matrix returned by this method could be transposed before
        using it in dot product, to avoid transposing a very large RGB array.
        Furthermore, after performing such conversion, the result array is in PCSXYZ color space with illuminant of D50,
        to get D65 results another adaptation transformation is needed.
        :return: a 2D array with shape of (3, 3), representing the matrix.
        """
        if self._matrix is None:
            if self.icc_path is not None and os.path.exists(self.icc_path):
                with exiftool.ExifTool() as et:
                    red_string = et.execute(*['-icc_profile:RedMatrixColumn', '-b', self.icc_path])
                    green_string = et.execute(*['-icc_profile:GreenMatrixColumn', '-b', self.icc_path])
                    blue_string = et.execute(*['-icc_profile:BlueMatrixColumn', '-b', self.icc_path])
                    if red_string and green_string and blue_string:
                        red_column = list(
                            map(float, et.execute(*['-icc_profile:RedMatrixColumn', '-b', self.icc_path]).split(' ')))
                        green_column = list(
                            map(float, et.execute(*['-icc_profile:GreenMatrixColumn', '-b', self.icc_path]).split(' ')))
                        blue_column = list(
                            map(float, et.execute(*['-icc_profile:BlueMatrixColumn', '-b', self.icc_path]).split(' ')))
                        self._matrix = np.vstack([red_column, green_column, blue_column]).transpose()
                        self._matrix_missing = False
                        print(f'using {self._matrix} for RGB to XYZ')
                    else:
                        self._matrix = colour.models.rgb.RGB_COLOURSPACE_sRGB.matrix_RGB_to_XYZ
                        self._matrix_missing = True
                        print(f'profile matrix not found, using sRGB matrix instead')
            else:
                self._matrix = colour.models.rgb.RGB_COLOURSPACE_sRGB.matrix_RGB_to_XYZ
                self._matrix_missing = True
        return self._matrix

    def get_chromatic_adaptation_transformation(self) -> np.ndarray | None:
        """
        Get the matrix for chromatic adaptation transformation (CAT). This is used to convert from PCS white point (D65)
        to the actual white point, but *NOT USED* since the actual white point is not consistent.
        :return: a 2D array with shape of (3, 3), representing the matrix.
        """
        if self._cat is None:
            if self.icc_path is not None and os.path.exists(self.icc_path):
                with exiftool.ExifTool() as et:
                    cat_string = et.execute(*['-icc_profile:ChromaticAdaptation', '-b', self.icc_path])
                    if cat_string:
                        cat = list(map(float, cat_string.split(' ')))
                        self._cat = np.vstack([cat[0:3], cat[3:6], cat[6:9]])
                        print(f"using {self._cat} for chromatic adaptation")
                    else:
                        self._cat = None
                        print(f"chromatic adaptation matrix missing, not performing transformation")
            else:
                self._cat = None
                print(f"icc profile missing, not performing transformation")
        return self._cat

    def get_pcs_illuminant(self) -> np.ndarray:
        """
        Get the white point of PCS.
        :return: a 1D XYZ array.
        """
        if self._pcs_illuminant_XYZ is None:
            if self.icc_path is not None and os.path.exists(self.icc_path):
                with exiftool.ExifTool() as et:
                    pcs_illuminant = et.execute(*['-icc_profile:ConnectionSpaceIlluminant', '-b', self.icc_path])
                    if pcs_illuminant:
                        self._pcs_illuminant_XYZ = list(map(float, pcs_illuminant.split(' ')))
                        print(f"extracted PCS illuminant {self._pcs_illuminant_XYZ}")
                    else:
                        # use D65 when no matrix is defined, otherwise use D50
                        if self._matrix_missing:
                            self._pcs_illuminant_XYZ = colour.xyY_to_XYZ(np.array([0.31271, 0.32902, 1]))
                            print(f"using default illuminant D65")
                        else:
                            self._pcs_illuminant_XYZ = colour.xyY_to_XYZ(np.array([0.34567, 0.35850, 1]))
                            print(f"matrix found but no PCS illuminant is present, using D50")
            else:
                self._pcs_illuminant_XYZ = colour.xyY_to_XYZ(np.array([0.31271, 0.32902, 1]))
        return self._pcs_illuminant_XYZ

    def get_conversion_matrix_to_D65(self) -> np.ndarray:
        """
        Get the matrix for chromatic adaptation transformation (CAT) from PCS white point to D65. This method first
        extracts the white point from the ICC profile, then performs
        :return:
        """
        if self._pcs_illuminant_conversion_matrix is None:
            self._pcs_illuminant_conversion_matrix = colour.adaptation.vonkries.matrix_chromatic_adaptation_VonKries(
                    self.get_pcs_illuminant(),  # use PCS white point
                    colour.xyY_to_XYZ(np.array([0.31271, 0.32902, 1])),  # D65 in 2-degree observer
                    transform='Bradford')
        print(f'using {self._pcs_illuminant_conversion_matrix} for chromatic adaptation to D65')
        return self._pcs_illuminant_conversion_matrix

    def get_profile_description(self) -> str:
        """
        Get the profile description.
        :return: a string, the profile description.
        """
        if self._profile_description is None:
            if self.icc_path is not None and os.path.exists(self.icc_path):
                with exiftool.ExifToolHelper() as et:
                    metadata = et.get_metadata(self.icc_path)
                    try:
                        self._profile_description = metadata[0]["ICC_Profile:ProfileDescription"]
                    except KeyError:
                        self._profile_description = 'sRGB'
            else:
                return 'sRGB'
        return self._profile_description


class CurveConversion(ConversionFunction):
    """
    Class used to carry out curve type conversion.
    """

    def __init__(self, matrix: np.array):
        """
        :param matrix: a 1D array, representing the sampled one-dimensional function.
        The matrix contains n values, splitting the range 0,0 to 1,0 into n−1 equal sized intervals.
        Each value represents a point of exactly a fraction
        Values between those fraction points should be calculated via linear interpolation.
        """
        self.matrix = matrix

    def convert_trc(self, pixels: np.ndarray) -> np.ndarray:
        """
        Convert the pixels using the curve function.
        :param pixels: a 1D array, the channel of pixels to be converted, whose dtype is float32.
        :return: a 1D array, converted and clipped to [0, 1].
        The standard says that the domain and range of the curve function is [0, 1], but it does not specify
        what to do when the output is out of range. Therefore, we clip it to [0, 1].
        """
        n = len(self.matrix)
        # np.interp does not check for out of range values, so we need to do it manually
        print(f'using matrix {self.matrix} for TRC curve conversion')
        if np.min(pixels) < 0 or np.max(pixels) > 1:
            raise ValueError(f'pixels should be in range [0, 1], got {np.min(pixels)} and {np.max(pixels)}')
        result_array = np.interp(pixels, np.linspace(0, 1, n), self.matrix)
        return np.clip(result_array, 0, 1)


class GammaConversion(ConversionFunction):
    """
    CLass used for gamma correction.
    Typically, for an image with gamma value of 2.2, the conversion function is y = x^2.2.
    """

    def __init__(self, gamma: float = 2.2):
        self.gamma = gamma

    def convert_trc(self, pixels: np.ndarray) -> np.ndarray:
        print(f'using gamma {self.gamma} for TRC conversion')
        if self.gamma != 1:
            return pixels ** self.gamma
        else:
            return pixels


class ParametricConversion(ConversionFunction):
    """
    Class used for parametric curve conversion.
    """

    def __init__(self, curve_type: int, params: np.ndarray):
        self.curve_type = curve_type
        self.params = params

    def convert_trc(self, pixels: np.ndarray) -> np.ndarray:
        """
        Convert the pixels using the parametric curve function. There are five types, referring to
        Table 68 in ICC:1-2022 standard.
        :param pixels: a 1D array, the channel of pixels to be converted, whose dtype is float32.
        :return: a 1D array, converted and clipped to [0, 1]. The standard says that:
        > Any function value outside the range shall be clipped to the range of the function.
        Therefore, we do the clipping before returning.
        """
        print(f'using type {self.curve_type} and params {self.params} for TRC parametric curve conversion')
        if self.curve_type == 0:
            return GammaConversion(gamma=self.params[0].astype(float)).convert_trc(pixels)
        elif self.curve_type == 1:
            g = self.params[0]
            a = self.params[1]
            b = self.params[2]
            result_array = np.piecewise(pixels, [pixels < -b / a, pixels >= -b / a],
                                        [lambda x: 0, lambda x: (a * x + b) ** g])
        elif self.curve_type == 2:
            g = self.params[0]
            a = self.params[1]
            b = self.params[2]
            c = self.params[3]
            result_array = np.piecewise(pixels, [pixels < -b / a, pixels >= -b / a],
                                        [lambda x: c, lambda x: ((a * x + b) ** g + c)])
        elif self.curve_type == 3:
            g = self.params[0]
            a = self.params[1]
            b = self.params[2]
            c = self.params[3]
            d = self.params[4]
            result_array = np.piecewise(pixels, [pixels < d, pixels >= d],
                                        [lambda x: c * x, lambda x: (a * x + b) ** g])
        elif self.curve_type == 4:
            g = self.params[0]
            a = self.params[1]
            b = self.params[2]
            c = self.params[3]
            d = self.params[4]
            e = self.params[5]
            f = self.params[6]
            result_array = np.piecewise(pixels, [pixels < d, pixels >= d],
                                        [lambda x: (c * x + f), lambda x: ((a * x + b) ** g + e)])
        else:
            raise ValueError(f'curve type {self.curve_type} not supported')
        return np.clip(result_array, 0, 1)


def plot_xy_coordinates_with_color(xy_array, output_png_path):
    start_time = time.time()

    # convert xy_array to xyY, then XYZ and RGB ('sRGB' color space for plotting)
    xy2_array = colour.xy_to_xyY(xy_array, 0.6)
    xyz_array = colour.xyY_to_XYZ(xy2_array)
    rgb_array = colour.XYZ_to_RGB(xyz_array, colour.RGB_COLOURSPACES['sRGB'])
    rgb_array_clipped = np.clip(rgb_array, 0, 1)

    # plot
    plt.style.use('dark_background')
    fig, ax = plt.subplots(figsize=(8, 9))
    ax.set_title('CIE 1931 Chromaticity Diagram')
    # draw horseshoe
    horseshoe_patch = PathPatch(horseshoe_path, facecolor='none', edgecolor='#DDDDDD', linewidth=0.5)
    ax.add_patch(horseshoe_patch)
    # draw xy coordinates
    ax.scatter(xy2_array[:, 0], xy2_array[:, 1], color=rgb_array_clipped, s=0.1, edgecolors=None, linewidths=0)
    # setup axes
    ax.set_xlim(0, 0.8)
    ax.set_ylim(0, 0.9)
    ax.xaxis.set_ticks(np.arange(0, 0.9, 0.1))
    ax.yaxis.set_ticks(np.arange(0, 1.0, 0.1))
    # draw white point and annotation
    ax.scatter(0.3127, 0.3290, color='#DDDDDD', s=5, edgecolors=None, linewidths=0)
    ax.scatter(0.3127, 0.3290, color=(0, 0, 0, 0), s=20, edgecolors='#DDDDDD', linewidths=0.5)
    ax.text(0.325, 0.319, 'D65', color='#DDDDDD', fontsize=6, ha='center', va='center')
    ax.scatter(0.3457, 0.3585, color='#DDDDDD', s=5, edgecolors=None, linewidths=0)
    ax.scatter(0.3457, 0.3585, color=(0, 0, 0, 0), s=20, edgecolors='#DDDDDD', linewidths=0.5)
    ax.text(0.353, 0.365, 'D50', color='#DDDDDD', fontsize=6, ha='center', va='center')

    plt.savefig(output_png_path, format='jpg', dpi=500, facecolor='#000000')
    plt.close()
    print('Drawing Chromaticity Diagram spent: {:.2f} seconds'.format(time.time() - start_time))


def image_to_cie_xy(image_path) -> np.ndarray:
    start_time = time.time()
    # read image file
    try:
        img = iio.imread(image_path)
    except FileNotFoundError:
        print("file not found!")
        exit(1)
    # extract image color space using exiftool
    image_data = np.array(img, dtype=np.float32)
    # if image is integer, normalize it to floating [0, 1]
    if np.issubdtype(img.dtype, np.integer):
        print(f"original image is {str(img.dtype):s}, normalizing to float32 [0, 1]")
        image_data = image_data / np.iinfo(img.dtype).max
    else:
        print(f"original image is {str(img.dtype):s}, no need to normalize")

    # if image is RGBA, convert it to RGB by removing alpha channel
    if image_data.shape[2] == 4:
        print("original image is RGBA, removing alpha channel")
        image_data = image_data[:, :, :3]
    # reshape image data to 2D array
    pixels = image_data.reshape(-1, 3)

    # convert RGB (if format is RGB, from corresponding color space, or 'sRGB' by default) to XYZ and then xy
    icc_profile = ICCProfile(image_path)
    pixels_corrected = np.zeros(pixels.shape)
    for i, channel in enumerate(['R', 'G', 'B']):
        pixels_corrected[:, i] = icc_profile.get_trc()[channel].convert_trc(pixels[:, i])

    conversion_matrix = icc_profile.get_matrix()
    pixels_xyz = np.dot(pixels_corrected, conversion_matrix.transpose())
    # perform chromatic adaptation transformation
    pixels_xyz = np.dot(pixels_xyz, icc_profile.get_conversion_matrix_to_D65().transpose())
    # cat_matrix = icc_profile.get_chromatic_adaptation_transformation()
    # if cat_matrix is not None:
    #     pixels_xyz = np.dot(pixels_xyz, np.linalg.inv(cat_matrix).transpose())
    xy_array = colour.XYZ_to_xy(pixels_xyz)
    print('Computing XYZ and xy spent: {:.2f} seconds'.format(time.time() - start_time))

    return xy_array


if __name__ == "__main__":
    '''
    The Python script converts image RGB data to CIE xy chromaticity coordinates and RGB values
    using imageio and NumPy (function image_to_cie_xy). 
    It then visualizes these coordinates on the CIE 1931 Chromaticity Diagram 
    with Matplotlib, coloring each point according to 
    its RGB value (function plot_xy_coordinates_with_color). 
    The development focused on optimizing data processing and visualization, 
    utilizing vectorized operations for efficiency, 
    ensuring accurate color representation, and producing visual output.
    '''
    # use first param as path
    image_path = "test.jpg"
    if len(sys.argv) > 1 and sys.argv[1] != "":
        print('Using image file: ' + sys.argv[1])
        image_path = sys.argv[1]
    else:
        print('Using default image file: ' + image_path)
    xy_coordinates = image_to_cie_xy(image_path)

    output_png_path = os.path.splitext(image_path)[0] + '_diagram.jpg'
    plot_xy_coordinates_with_color(xy_coordinates, output_png_path)

Install requirements:

pip3 install numpy matplotlib colour-science imageio PyExifTool

Furthermore, you'll need ExifTool by Phil Harvey installed to extract metadata and most importantly color space from image file.

Usage:

python3 picture_to_cie_diagram.py <image_path>

---

Details

This script could be used as a way to visualize color distribution in a image. Multiple formats and color spaces are supported, credits to colour-science and imageio .

The original script is from Anthony-Hoo, posted here , I corrected & improved the conversion procedures and standardized it using colour-science . I will later write a blog on this script and provide a way to apply this script in reality.

For an early example, here is a test image file:

Once processed by the script, you will get a diagram like:

Note:

There is still much to be improved in current version, and I'm going to release a patch in a few days.

It works, currently, with “standard” color profiles, like Adobe RGB (1998) and Display P3 and so on, however when worked with photos with a custom icc profile (like those taken with stock camera on Xiaomi or Google phones) the color mapping will be slightly different from what it should be. This is, in most cases, not really fatal, but it can cause problems, especially when you need to quantify colors as what this script is supposed to do. I'll later release a blog (Chinese) commenting on what have gone wrong in the original script and update THIS gist, but for now, unless you are sure you will only process TIFFs with STANDARD color spaces, DO NOT TRUST THIS.

Update: Now the script should work correctly after the 6th revision.

Thank you! I'm very grateful to the both of you (@Anthony-Hoo).

My use-case is to visualize the output of colorization deep learning models to check the richness of the predicted colors.
I'm not sure it makes too much sense since it doesn't have any color profiles so I use the standard RGB.

My use-case is to visualize the output of colorization deep learning models to check the richness of the predicted colors. I'm not sure it makes too much sense since it doesn't have any color profiles so I use the standard RGB.

@Dawars Then it should be OK. Anyway, I decided to post improved script before the blog is finished, and this gist have just been updated. Enjoy it ;)

Which version of colour are you using?

The parameters of colour.XYZ_to_RGB() in line 400 have changed and it requires now a matrix and illumination. https://colour.readthedocs.io/en/develop/generated/colour.XYZ_to_RGB.html
I'm not sure if my fix is correct:

illuminant = np.array([0.34570, 0.35850])  # hardcoded
rgb_array = colour.XYZ_to_RGB(xyz_array, illuminant, illuminant, colour.RGB_COLOURSPACES['sRGB'].XYZ_to_RGB_matrix)    

Which version of colour are you using?

The parameters of colour.XYZ_to_RGB() in line 400 have changed and it requires now a matrix and illumination. https://colour.readthedocs.io/en/develop/generated/colour.XYZ_to_RGB.html I'm not sure if my fix is correct:

illuminant = np.array([0.34570, 0.35850])  # hardcoded
rgb_array = colour.XYZ_to_RGB(xyz_array, illuminant, illuminant, colour.RGB_COLOURSPACES['sRGB'].XYZ_to_RGB_matrix)    

@Dawars I'm using colour-science 0.4.4 , which I believe to be the latest version. From your link (and that page of latest/v0.4.4/master branch) it doesn't seem that two illuminant parameters are needed. The only illuminant parameter exists for adapting to different white points of different color spaces, and it should be ok to leave it empty. Anyway, if the method requires two illuminant parameters, then your fix is right, avoiding any unexpected conversion.

Sorry, my bad. I used an older version because I had an older Python environment.🙏