Gaussian-Based Spectral Upsampling

Colour Developers

2026-01-08 08:00

Introduction¶

Spectral reflectance upsampling, or reflectance recovery, is the process of converting tristimulus values into spectral reflectance distributions.

This document explores the Gaussian-based spectral upsampling method implemented in Colour, develop branch as of this writing. The method uses clamped Gaussian basis spectra with the hierarchical decomposition algorithm from Smits (1999).

In [1]:

import matplotlib.pyplot as plt
import numpy as np
import timeit

import colour
from colour.characterisation import ColourChecker

colour.plotting.colour_style()

colour.utilities.describe_environment();

===============================================================================
*                                                                             *
*   Interpreter :                                                             *
*       python : 3.13.0 (main, Oct  7 2024, 23:47:22) [Clang 18.1.8 ]         *
*                                                                             *
*   colour-science.org :                                                      *
*       colour : 0.4.7                                                        *
*                                                                             *
*   Runtime :                                                                 *
*       imageio : 2.37.2                                                      *
*       matplotlib : 3.10.8                                                   *
*       networkx : 3.6.1                                                      *
*       numpy : 2.4.0                                                         *
*       scipy : 1.16.3                                                        *
*       tqdm : 4.67.1                                                         *
*       OpenImageIO : 3.1.9.0                                                 *
*                                                                             *
===============================================================================

Smits (1999) Reflectance Recovery¶

A simplistic reflectance recovery approach would create spectra only for red, green, and blue, then combine them to generate other reflectance spectra. However, the sum of useful R, G, and, B spectra does not produce a uniform white:

In [ ]:

MSDS_SMITS1999 = colour.recovery.datasets.smits1999.MSDS_SMITS1999

SD_R_p_G_p_B = (
    MSDS_SMITS1999.signals["red"]
    + MSDS_SMITS1999.signals["green"]
    + MSDS_SMITS1999.signals["blue"]
)
SD_R_p_G_p_B.display_name = "R + G + B"

colour.plotting.plot_multi_sds(
    [MSDS_SMITS1999.signals["white"], SD_R_p_G_p_B],
);

Smits (1999) solves the issue by pre-computing seven independently optimised basis spectra: white, cyan, magenta, yellow, red, green, and blue.

Any colour can be expressed as a weighted sum of white spectrum, plus one secondary spectrum (cyan, magenta, or yellow), and finally one primary spectrum (red, green, or blue), using as much of the wider spectra as possible:

Spectrum RGBToSpectrum(R, G, B):
    Spectrum result = 0
    
    if R <= G and R <= B:
        result += R * S_W
        if G <= B:
            result += (G - R) * S_C
            result += (B - G) * S_B
        else:
            result += (B - R) * S_C
            result += (G - B) * S_G
            
    else if G <= R and G <= B:
        result += G * S_W
        if R <= B:
            result += (R - G) * S_M
            result += (B - R) * S_B
        else:
            result += (B - G) * S_M
            result += (R - B) * S_R
            
    else:  # B <= R and B <= G
        result += B * S_W
        if R <= G:
            result += (R - B) * S_Y
            result += (G - R) * S_G
        else:
            result += (G - B) * S_Y
            result += (R - G) * S_R
    
    return result

Where $R$, $G$, $B$ are scalar RGB values, and $S_W$, $S_C$, $S_M$, $S_Y$, $S_R$, $S_G$, $S_B$ are the pre-computed basis spectra (white, cyan, magenta, yellow, red, green, blue).

Smits (1999) provides spectral distributions for 10 bins between 380nm and 720nm:

bin	white	cyan	magenta	yellow	red	green	blue
1	1.0000	0.9710	1.0000	0.0001	0.1012	0.0000	1.0000
2	1.0000	0.9426	1.0000	0.0000	0.0515	0.0000	1.0000
3	0.9999	1.0007	0.9685	0.1088	0.0000	0.0273	0.8916
4	0.9993	1.0007	0.2229	0.6651	0.0000	0.7937	0.3323
5	0.9992	1.0007	0.0000	1.0000	0.0000	1.0000	0.0000
6	0.9998	1.0007	0.0458	1.0000	0.0000	0.9418	0.0000
7	1.0000	0.1564	0.8369	0.9996	0.8325	0.1719	0.0003
8	1.0000	0.0000	1.0000	0.9586	1.0149	0.0000	0.0369
9	1.0000	0.0000	1.0000	0.9685	1.0149	0.0000	0.0483
10	1.0000	0.0000	0.9959	0.9840	1.0149	0.0025	0.0496

Given the interval size (~38nm), they exhibit abrupt transitions:

In [3]:

colour.plotting.plot_multi_sds(
    colour.colorimetry.reshape_msds(colour.recovery.datasets.smits1999.MSDS_SMITS1999,
                                    colour.SPECTRAL_SHAPE_DEFAULT
    ),
    plot_kwargs={
        "use_sd_colours": True,
    },
);

No description has been provided for this image

Gaussian-Based Spectral Upsampling¶

Colour's default spectral shape follows the ASTM E308-15 standard practice (360-780nm at 1nm intervals). At this resolution, increased smoothness in the basis spectra is desirable for better reconstruction performance.

The Gaussian-based method replaces Smits (1999) basis spectra with clamped super-Gaussians that have independently optimised parameters for round-trip colorimetric accuracy with sRGB primaries.

The Gaussian basis parameters are found by minimising a combined error. Given parameters $\boldsymbol{\theta} = (\lambda_R, \lambda_G, \lambda_B, \lambda_C, \lambda_M, \lambda_Y, \sigma_R, \sigma_G, \sigma_B, \sigma_C, \sigma_M, \sigma_Y, n_R, n_G, n_B, n_C, n_M, n_Y)$ representing peak wavelengths, FWHM values, and exponents for all six chromatic basis spectra, the objective function is:

$$\min_{\boldsymbol{\theta}} \left( \sum_{i} \left\| \mathbf{XYZ}_i - \mathbf{XYZ}'_i(\boldsymbol{\theta}) \right\|^2 + \overline{\Delta E^*_{00}} + \sum_{j} (n_j - 2)^2 \right)$$

where:

$\mathbf{XYZ}_i$ are the test colours (RGB primaries, CMY secondaries, and grey)
$\mathbf{XYZ}'_i$ is obtained by converting each test colour to RGB, recovering a spectrum using the Smits (1999) decomposition with the Gaussian basis defined by $\boldsymbol{\theta}$, and integrating back to CIE XYZ tristimulus values
$\overline{\Delta E^*_{00}}$ is the mean CIE DE2000 error on ColorChecker patches
The last term is a smoothness penalty that penalises deviation from the standard Gaussian exponent ($n = 2$)

In [4]:

colour.plotting.plot_multi_sds(
    colour.recovery.gaussian.MSDS_GAUSSIAN_BASIS,
    plot_kwargs={
        "use_sd_colours": True,
    },
);

With the basis spectra obtained, we can then spectrally upsample a ColorChecker Classic and compare the recovered spectra from both methods to the reference ones.

In [5]:

SDS_COLORCHECKER_CLASSIC = list(colour.SDS_COLOURCHECKERS["ISO 17321-1"].values())

RGB_COLORCHECKER_CLASSIC = colour.XYZ_to_sRGB(
    [
        colour.sd_to_XYZ(sd, illuminant=colour.SDS_ILLUMINANTS["D65"]) / 100
        for sd in SDS_COLORCHECKER_CLASSIC
    ],
    apply_cctf_encoding=False,
)

MSDS_SMITS1999_ALIGNED = colour.colorimetry.reshape_msds(
    colour.recovery.datasets.smits1999.MSDS_SMITS1999, colour.SPECTRAL_SHAPE_DEFAULT
)

SDS_COLORCHECKER_CLASSIC_SMITS1999 = colour.recovery.RGB_to_msds_Smits1999(
    RGB_COLORCHECKER_CLASSIC, basis=MSDS_SMITS1999_ALIGNED
)
SDS_COLORCHECKER_CLASSIC_GAUSSIAN = colour.recovery.RGB_to_msds_Gaussian(
    RGB_COLORCHECKER_CLASSIC
)

W = colour.SPECTRAL_SHAPE_DEFAULT.wavelengths

FIGURE = plt.figure(figsize=(18, 12))
GS = FIGURE.add_gridspec(4, 6, hspace=0, wspace=0)
AXES = GS.subplots(sharex=True, sharey=True)

for i, (sd_original, axes) in enumerate(zip(SDS_COLORCHECKER_CLASSIC, AXES.flatten())):
    sd_original_aligned = colour.colorimetry.reshape_sd(
        sd_original, colour.SPECTRAL_SHAPE_DEFAULT
    )

    axes.plot(
        W,
        sd_original_aligned.values,
        label="Original",
        linewidth=1.5,
        color="r",
        alpha=0.5,
    )
    axes.plot(
        W,
        SDS_COLORCHECKER_CLASSIC_SMITS1999[i],
        label="Smits (1999)",
        linewidth=1.5,
        color="g",
        alpha=0.5,
    )
    axes.plot(
        W,
        SDS_COLORCHECKER_CLASSIC_GAUSSIAN[i],
        label="Gaussian",
        linewidth=1.5,
        color="b",
        alpha=0.5,
    )

    axes.text(
        0.5,
        0.9,
        sd_original.name,
        horizontalalignment="center",
        transform=axes.transAxes,
        fontsize=8,
    )
    axes.set_ylim(0, 1)

    if i == 0:
        axes.legend(fontsize=6, loc="upper right")

# Hide tick labels on inner subplots
for axes in AXES.flatten():
    axes.label_outer()

# Remove left ticks on inner columns
for axes in AXES[:, 1:].flatten():
    axes.tick_params(left=False, which="both")

# Hide overlapping y-tick labels at row boundaries (first column only)
for axes in AXES[1:, 0]:
    axes.yaxis.get_major_ticks()[-1].label1.set_visible(False)
for axes in AXES[:-1, 0]:
    axes.yaxis.get_major_ticks()[0].label1.set_visible(False)

FIGURE.supxlabel("Wavelength (nm)")
FIGURE.supylabel("Reflectance")
FIGURE.suptitle("ColorChecker Classic: Reference vs Recovered Spectra", fontsize=14)
plt.tight_layout()
plt.show();

In [6]:

SD_DARK_SKIN = colour.colorimetry.reshape_sd(SDS_COLORCHECKER_CLASSIC[0], colour.SPECTRAL_SHAPE_DEFAULT)

colour.plotting.plot_multi_sds(
    [
        SD_DARK_SKIN,
        colour.SpectralDistribution(
            SDS_COLORCHECKER_CLASSIC_SMITS1999[0], W, name="Smits (1999)"
        ),
        colour.SpectralDistribution(
            SDS_COLORCHECKER_CLASSIC_GAUSSIAN[0], W, name="Gaussian"
        ),
    ],
);

In [7]:

XYZ_COLORCHECKER_CLASSIC_REFERENCE = [
    colour.sd_to_XYZ(sd, illuminant=colour.SDS_ILLUMINANTS["D65"]) / 100
    for sd in SDS_COLORCHECKER_CLASSIC
]

XYZ_COLORCHECKER_CLASSIC_SMITS1999 = [
    colour.sd_to_XYZ(
        colour.SpectralDistribution(sd, W),
        illuminant=colour.SDS_ILLUMINANTS["D65"],
    )
    / 100
    for sd in SDS_COLORCHECKER_CLASSIC_SMITS1999
]


XYZ_COLORCHECKER_CLASSIC_GAUSSIAN = [
    colour.sd_to_XYZ(
        colour.SpectralDistribution(sd, W),
        illuminant=colour.SDS_ILLUMINANTS["D65"],
    )
    / 100
    for sd in SDS_COLORCHECKER_CLASSIC_GAUSSIAN
]

xyY_REFERENCE = colour.XYZ_to_xyY(XYZ_COLORCHECKER_CLASSIC_REFERENCE)
xyY_SMITS1999 = colour.XYZ_to_xyY(XYZ_COLORCHECKER_CLASSIC_SMITS1999)
xyY_GAUSSIAN = colour.XYZ_to_xyY(XYZ_COLORCHECKER_CLASSIC_GAUSSIAN)

CCS_D65 = colour.CCS_ILLUMINANTS["CIE 1931 2 Degree Standard Observer"]["D65"]
NAMES = [sd.name for sd in SDS_COLORCHECKER_CLASSIC]

colour.plotting.plot_multi_colour_checkers(
    [
        ColourChecker(
            "Reference",
            dict(zip(NAMES, xyY_REFERENCE)),
            CCS_D65,
            4,
            6,
        ),
        ColourChecker(
            "Smits (1999)",
            dict(zip(NAMES, xyY_SMITS1999)),
            CCS_D65,
            4,
            6,
        ),
    ]
)

colour.plotting.plot_multi_colour_checkers(
    [
        ColourChecker(
            "Reference",
            dict(zip(NAMES, xyY_REFERENCE)),
            CCS_D65,
            4,
            6,
        ),
        ColourChecker(
            "Gaussian",
            dict(zip(NAMES, xyY_GAUSSIAN)),
            CCS_D65,
            4,
            6,
        ),
    ]
);

In [8]:

Lab_REFERENCE = colour.XYZ_to_Lab(XYZ_COLORCHECKER_CLASSIC_REFERENCE)
Lab_SMITS1999 = colour.XYZ_to_Lab(XYZ_COLORCHECKER_CLASSIC_SMITS1999)
Lab_GAUSSIAN = colour.XYZ_to_Lab(XYZ_COLORCHECKER_CLASSIC_GAUSSIAN)

dE_SMITS1999 = colour.delta_E(Lab_REFERENCE, Lab_SMITS1999, method="CIE 2000")
dE_GAUSSIAN = colour.delta_E(Lab_REFERENCE, Lab_GAUSSIAN, method="CIE 2000")

print("CIE DE2000: Reference vs Smits (1999)")
print(f"  Mean: {np.mean(dE_SMITS1999):.2f}")
print(f"  Max:  {np.max(dE_SMITS1999):.2f}")

print()

print("CIE DE2000: Reference vs Gaussian")
print(f"  Mean: {np.mean(dE_GAUSSIAN):.2f}")
print(f"  Max:  {np.max(dE_GAUSSIAN):.2f}")

CIE DE2000: Reference vs Smits (1999)
  Mean: 1.05
  Max:  3.18

CIE DE2000: Reference vs Gaussian
  Mean: 0.99
  Max:  2.61

The colour.recovery.gaussian.optimise_gaussian_basis_parameters and colour.recovery.gaussian.generate_gaussian_basis definitions allow optimisation of a new set of basis spectra for a different illuminant, colour matching functions, or even colourspace.

It is, for example, possible to optimise basis spectra for the sRGB colourspace under the CIE Standard Illuminant D65 (the default is sRGB primaries under illuminant E):

In [9]:

# Optimise basis for sRGB colourspace (D65 illuminant)
PEAK_WL_D65, FWHM_D65, EXPONENT_D65 = (
    colour.recovery.gaussian.optimise_gaussian_basis_parameters(
        colourspace=colour.models.RGB_COLOURSPACE_sRGB,
    )
)

MSDS_GAUSSIAN_BASIS_D65 = colour.recovery.gaussian.generate_gaussian_basis(
    peak_wavelengths=PEAK_WL_D65,
    fwhm=FWHM_D65,
    exponent=EXPONENT_D65,
)

SDS_COLORCHECKER_CLASSIC_GAUSSIAN_D65 = colour.recovery.RGB_to_msds_Smits1999(
    RGB_COLORCHECKER_CLASSIC, basis=MSDS_GAUSSIAN_BASIS_D65
)

XYZ_COLORCHECKER_CLASSIC_GAUSSIAN_D65 = [
    colour.sd_to_XYZ(
        colour.SpectralDistribution(sd, W),
        illuminant=colour.SDS_ILLUMINANTS["D65"],
    )
    / 100
    for sd in SDS_COLORCHECKER_CLASSIC_GAUSSIAN_D65
]

xyY_GAUSSIAN_D65 = colour.XYZ_to_xyY(XYZ_COLORCHECKER_CLASSIC_GAUSSIAN_D65)

colour.plotting.plot_multi_colour_checkers(
    [
        ColourChecker(
            "Reference",
            dict(zip(NAMES, xyY_REFERENCE)),
            CCS_D65,
            4,
            6,
        ),
        ColourChecker(
            "Gaussian (sRGB/D65)",
            dict(zip(NAMES, xyY_GAUSSIAN_D65)),
            CCS_D65,
            4,
            6,
        ),
    ]
);

Lab_GAUSSIAN_D65 = colour.XYZ_to_Lab(XYZ_COLORCHECKER_CLASSIC_GAUSSIAN_D65)
dE_GAUSSIAN_D65 = colour.delta_E(Lab_REFERENCE, Lab_GAUSSIAN_D65, method="CIE 2000")

print("CIE DE2000: Reference vs Gaussian (sRGB/D65)")
print(f"  Mean: {np.mean(dE_GAUSSIAN_D65):.2f}")
print(f"  Max:  {np.max(dE_GAUSSIAN_D65):.2f}")

CIE DE2000: Reference vs Gaussian (sRGB/D65)
  Mean: 0.77
  Max:  2.24

Performance Benchmarks¶

An advantage of the Gaussian-based method is its support for multi-spectral distribution recovery via RGB_to_msds_Gaussian. Other methods offered by Colour, excluding Smits (1999), only support a single spectral distribution at a time.

In [10]:

# Use ColorChecker XYZ values repeated 10x for benchmarking
XYZ_BENCHMARK = np.tile(XYZ_COLORCHECKER_CLASSIC_REFERENCE, (10, 1))
N_SAMPLES = len(XYZ_BENCHMARK)

# Benchmark msds methods (batch processing)
TIME_GAUSSIAN_MSDS = timeit.timeit(
    lambda: colour.recovery.XYZ_to_msds(XYZ_BENCHMARK, method="Gaussian"),
    number=10,
)

TIME_SMITS1999_MSDS = timeit.timeit(
    lambda: colour.recovery.XYZ_to_msds(XYZ_BENCHMARK, method="Smits 1999"),
    number=10,
)

# Benchmark all sd methods (single colour at a time)
SD_METHODS = ["Gaussian", "Smits 1999", "Mallett 2019", "Meng 2015", "Jakob 2019", "Otsu 2018"]
SD_TIMES = {}

for method in SD_METHODS:
    time_taken = timeit.timeit(
        lambda m=method: [colour.recovery.XYZ_to_sd(xyz, method=m) for xyz in XYZ_BENCHMARK],
        number=1,
    )
    SD_TIMES[method] = time_taken

print("Performance Benchmarks")
print("=" * 22)
print()
print(f"MultiSpectral Distributions Methods - {N_SAMPLES} samples, 10 iterations:")
print(f"  Gaussian:     {TIME_GAUSSIAN_MSDS * 100:.2f} ms")
print(f"  Smits (1999): {TIME_SMITS1999_MSDS * 100:.2f} ms")

print()

print(f"Spectral Distribution Methods - {N_SAMPLES} samples, 1 iteration:")
for method, time_taken in sorted(SD_TIMES.items(), key=lambda x: x[1]):
    print(f"  {method:14s} {time_taken * 1000:.2f} ms")

print()

print("Speedup: `XYZ_to_msds` vs `XYZ_to_sd`")
SPEEDUP_GAUSSIAN = (SD_TIMES["Gaussian"] / N_SAMPLES) / (TIME_GAUSSIAN_MSDS / 10 / N_SAMPLES)
SPEEDUP_SMITS1999 = (SD_TIMES["Smits 1999"] / N_SAMPLES) / (TIME_SMITS1999_MSDS / 10 / N_SAMPLES)
print(f"  Gaussian:     {SPEEDUP_GAUSSIAN:.0f}x faster with `XYZ_to_msds`")
print(f"  Smits (1999): {SPEEDUP_SMITS1999:.0f}x faster with `XYZ_to_msds`")

Performance Benchmarks
======================

MultiSpectral Distributions Methods - 240 samples, 10 iterations:
  Gaussian:     3.17 ms
  Smits (1999): 0.54 ms

Spectral Distribution Methods - 240 samples, 1 iteration:
  Mallett 2019   48.31 ms
  Otsu 2018      97.97 ms
  Gaussian       114.29 ms
  Smits 1999     133.90 ms
  Jakob 2019     3513.36 ms
  Meng 2015      113833.13 ms

Speedup: `XYZ_to_msds` vs `XYZ_to_sd`
  Gaussian:     36x faster with `XYZ_to_msds`
  Smits (1999): 249x faster with `XYZ_to_msds`

Conclusion¶

The Gaussian-based spectral upsampling method provides an effective approach for reflectance recovery with some advantages:

Smooth Spectra: Clamped Gaussian basis spectra produce smoother reflectance spectra compared to Smits (1999) basis spectra
Colourimetric Accuracy: Optimised parameters yield low $\Delta E^*_{00}$ errors on ColorChecker patches
Performance: The RGB_to_msds_Gaussian function enables efficient multi-spectral distribution processing, achieving significant speedups over spectral distribution methods
Flexibility: Basis spectra can be re-optimised for different illuminants, colour matching functions, or colourspaces via optimise_gaussian_basis_parameters

Introduction¶

Smits (1999) Reflectance Recovery¶

Gaussian-Based Spectral Upsampling¶

Performance Benchmarks¶

Conclusion¶

Comments