Question about colour.sd_to_XYZ definition.

[xkennyz]

Hi all!

First time writing here.
I just wanted to know if it's possible to make it so that the colour.sd_to_XYZ only integrates over the inputted wavelengths. For example, if I have a spectral distribution ranging from 500 to 700nm, the function will extrapolate the first and last values linearly to the left and right, respectively. Then it will perform the calculation.
I've tried setting the other wavelengths to zero (e.g., 495nm and below to 0, 705nm and above to 0), but the function doesn't quite give me what I expect it to. I believe it's because there's some sort of averaging that the function does to integrate (between 495 and 500, 700 and 705).
Basically, I want to somehow find a way to limit the bounds of integration to the range of my spectral distribution rather than it being 360~830 indefinitely.
I'm a little new to this stuff, so I'm sorry if I missed something basic.

Thanks!

[xkennyz]

Never mind! I modified the source code to do what I wanted.
Thanks for the great package!

[tjdcs]

Use sd_to_XYZ(sd.values, shape=sd.shape) and also sd_to_XYZ(sd.values, shape=sd.shape, method="Integration") should do what you want. However, the default behavior is is based on CIE recommendations and other standard practices. In my own work I use the above function calls because they also use a faster code path, but they are considered less accurate.

Your example is relatively bizarre, I'm curious why you don't have any data below 500nm? In this case I don't think there is any right answer as so much of the visual system's sensitivity is below this wavelength. It's hard to say any colorimetric conversions are meaningful.

[xkennyz]

Haha sorry, I was just giving an example of random wavelengths. In hindsight, it does look suspect.
The real data I'm dealing with goes from 480 to 620nm. Not entirely sure why I didn't just say that.

As for putting in the shapes, I've found that no matter what shape I put in, the spectral distribution gets extrapolated out to 360 and 780.
When I looked into the source code, I found that what happens is the spectral distribution itself gets reshaped into the shape of cmfs, RATHER THAN reshaping cmfs to the spectral distribution, which explains why nothing I did was able to get the function to work--it was working, just not the way I expected it to.

[tjdcs]

Yes. That is the recommended procedure by the CIE. However, the code path I posted above will adhere to your preference.

[xkennyz]

Sadly, as I mentioned above, it does not replicate the result I want.

data.txt
Using the data I have as an example (attached), when I run your suggestions, it comes out to be:
XYZ with sd_to_XYZ is: [ 753.45524368 2268.18722396 207.25291574]

The values I'm looking for are:
XYZ is: [ 743.2332, 2264.1323, 206.3468 ]

Using the code I modified, I get:
XYZ is: [ 743.51760662 2264.25109939 206.3403274 ]

Which is close enough to what I wanted.

[KelSolaar]

It is because shape is only used when integrating a n-dimensional array and not when using the colour.continouous.AbstractContinuousSignal sub-classes, e.g., colour.SpectralDistribution. The reason why is it like that is because integrating n-dimensional arrays was implemented much later and it did not occur to me that shape could/should also drive the integration of the initial code path but I think it should, it totally makes sense!

[KelSolaar]

(Note that the docs are explicit about when shape is being used though)

[xkennyz]

Huh. Go figure. Thanks, that explains it!

[tjdcs]

I don't get the same results as you with the data.txt you provided, but here is a notebook to explore the difference more. As @KelSolaar said there is a switch to force the reshaping of the cmfs rather than reshaping the input data based on passing in an n-array instead of an sd here:

colour/colour/colorimetry/tristimulus_values.py

Lines 677 to 715 in 7fb9fcc

	else:
	attest(
	shape is not None,
	"A spectral shape must be explicitly passed with a spectral data array!",
	)
	shape = cast(SpectralShape, shape)
	R = as_float_array(sd)
	shape_R = R.shape
	wl_c_r = R.shape[-1]
	wl_c = len(shape.wavelengths)
	attest(
	wl_c_r == wl_c,
	f"Spectral data array with {wl_c_r} wavelengths is not compatible "
	f"with spectral shape with {wl_c} wavelengths!",
	)
	if cmfs.shape != shape:
	runtime_warning(f'Aligning "{cmfs.name}" cmfs shape to "{shape}".')
	cmfs = reshape_msds(cmfs, shape, copy=False)
	if illuminant.shape != shape:
	runtime_warning(f'Aligning "{illuminant.name}" illuminant shape to "{shape}".')
	illuminant = reshape_sd(illuminant, shape, copy=False)
	XYZ_b = cmfs.values
	S = illuminant.values
	R = np.reshape(R, (-1, wl_c_r))
	d_w = cmfs.shape.interval
	with sdiv_mode():
	k = cast(Real, optional(k, sdiv(100, (np.sum(XYZ_b[..., 1] * S) * d_w))))
	XYZ = k * np.dot(R * S, XYZ_b) * d_w
	XYZ = from_range_100(np.reshape(XYZ, [*list(shape_R[:-1]), 3]))

Unless I'm still misunderstanding @xkennyz variables. I tried debugging with the example data but ended up with totally different XYZ values.

import numpy as np
from colour.colorimetry.spectrum import SpectralDistribution, SpectralShape
from colour.colorimetry.tristimulus_values import sd_to_XYZ

data = [
    0.001447728,
    0.031316602,
    0.281427964,
    1.164928271,
    2.827705839,
    4.683560365,
    6.304900832,
    7.418446126,
    7.92553908,
    7.948622573,
    7.571148614,
    6.952271141,
    6.218380217,
    5.365852514,
    4.528914745,
    3.753913463,
    3.028762779,
    2.392047624,
    1.861467159,
    1.419677822,
    1.066408985,
    0.789787672,
    0.576711539,
    0.415404445,
    0.296450351,
    0.208676513,
    0.144236724,
    0.100071995,
    0.068515426,
]

shape = SpectralShape(480, 620, 5)
sd = SpectralDistribution(data, domain=shape)

# This code path reshapes the cmfs to the input sd.shape
#
# The switch is on line 662 of tristimulus_values.py based on the type of the
# first argument

XYZ1 = sd_to_XYZ(sd.values, shape=sd.shape, method="Integration")

# This code path reshapes sd to match the default CMFs.

XYZ2 = sd_to_XYZ(sd)

print("XYZ1: " + np.array2string(XYZ1, formatter={"float_kind": lambda x: "%.4f" % x}))
print("XYZ2: " + np.array2string(XYZ2, formatter={"float_kind": lambda x: "%.4f" % x}))

> python test.py
XYZ1: [112.7542 343.0918 31.2653]
XYZ2: [103.2273 310.7786 28.3986]

I use this code path pretty often because it is faster during a lot of the optimization code I write, if it's not behaving as described then I'm really confused and pretty concerned about a lot of my results. Ha.

[tjdcs]

note* instead of passing in sd.values you could pass in the data variable directly but there is a bug right now where it fails to calculate a cache key correctly (need to make data an np.ndarray). I'll see if theres a nice place to insert a fix shortly.

XYZ1 = sd_to_XYZ(np.asarray(data), shape=sd.shape, method="Integration")