colour-science / colour

Colour Science for Python
https://www.colour-science.org
BSD 3-Clause "New" or "Revised" License
2.08k stars 259 forks source link

About interpolation versus zero filling in "colour.spectral_to_XYZ" definition. #232

Closed eyeforcolor closed 8 years ago

eyeforcolor commented 8 years ago

First of all, a big thank you for creating this excellent python library for color research! Something like this was sorely missing and this is so well written!

While I was testing it, I think I may have stumbled across a bug in the spectral_to_XYZ() conversion. Currently, inside this conversion, if the shape of the SPD does not match the shape of color matching functions, the "missing" SPD wavelength data is being zero filled. In my opinion, this should be replaced by interpolation. I have done some testing using the color checker SPD data included in the library.

Values below are CIELAB calculated using spectral_to_XYZ then XYZ_to_Lab functions using D65 illuminant for "BabelColor Average" and "ColorChecker N Ohta" SPDs. (I have trimmed it to the first 5 patches only).

Using clone.zeros() in the spectral_to_XYZ() method

Babel CIELAB Ohta CIELAB DEab
25.17 1.16 -4.37 36.47 1.80 -5.63 11.39
37.04 1.59 -5.28 52.18 2.74 -6.75 15.26
73.31 3.15 -6.65 95.47 4.28 -10.20 22.46
62.00 6.29 59.48 82.36 7.72 73.84 24.96
61.16 2.54 -7.27 80.96 4.03 -9.35 19.97

Change to clone.align() in the spectral_to_XYZ() method

Babel CIELAB Ohta CIELAB DEab
35.87 1.49 -5.53 36.47 1.80 -5.63 0.68
50.82 2.04 -6.68 52.18 2.74 -6.75 1.53
96.53 4.03 -8.44 95.47 4.28 -10.20 2.07
82.28 7.97 74.91 82.36 7.72 73.83 1.11
81.21 3.25 -9.21 80.96 4.03 -9.36 0.83

As you will notice, there is a big difference between the CIELAB values of Babel ColorChecker data versus the Ohta ColorChecker data, when zero filling the missing wavelengths, something I did not expect. I see that the Ohta CIELAB values are calculated from 5 nm data whereas the Babel values are calculated from 10 nm data. So I changed the zero filling to interpolate and they are much closer (I am aware that these are from different measurements so they cannot be the same).

Apologies for the long explanation. I would be happy to create a PR for this fix, let me know!

Thanks! Rohit Color Scientist

KelSolaar commented 8 years ago

Hi!

Thanks for the kind words and digging around the spectral interpolation jazz! :+1:

This is a currently known behaviour, and one we are not happy with! We are also working on it (more on this later).

In the current implementation of colour.spectral_to_XYZ definition we don't account for the differential distance dx when doing the operations. That's the ending d\lambda in those equations.

As a result the Luminance computed is different which should explain discrepancies you find in computing the CIE Lab values. That being said the colour should be almost the same (not accounting for the precision loss because of the different bin size and provided you normalise the values):

import colour
from colour.plotting import *

spd = colour.COLOURCHECKERS_SPDS['cc_ohta']['dark skin']
spd_c = spd.clone().interpolate(colour.SpectralShape(steps=10))
XYZ_a = colour.spectral_to_XYZ(spd) / 100
XYZ_b = colour.spectral_to_XYZ(spd_c) / 100

print('CIE XYZ:\n\t{0}\n\t{1}'.format(XYZ_a, XYZ_b))
print('CIE xy:\n\t{0}\n\t{1}'.format(colour.XYZ_to_xy(XYZ_a), 
                                     colour.XYZ_to_xy(XYZ_b)))
print('CIE Lab:\n\t{0}\n\t{1}'.format(colour.XYZ_to_Lab(XYZ_a), 
                                      colour.XYZ_to_Lab(XYZ_b)))
print('CIE Lab (Y Nornalised):\n\t{0}\n\t{1}'.format(colour.XYZ_to_Lab(XYZ_a / XYZ_a[1]), 
                                                     colour.XYZ_to_Lab(XYZ_b / XYZ_b[1])))

spds_CIE_1931_chromaticity_diagram_plot((spd, spd_c))
CIE XYZ:
    [ 0.02386821  0.0199886   0.01118792]
    [ 0.01193872  0.01000525  0.00557721]
CIE xy:
    [ 0.43361472  0.3631338 ]
    [ 0.43380111  0.36354738]
CIE Lab:
    [ 15.48126316  10.02488061   6.58787558]
    [ 8.99581624  7.93233569  4.98396404]
CIE Lab (Y Nornalised):
    [ 100.           36.93899272   24.27455226]
    [ 100.           36.81219815   24.51422506]

image

Now as I said, @MichaelMauderer and I are not happy about that and we are working on a complete overhaul of the spectral code. The BIG problem is that we actually can't just say: "Hey! Let's interpolate! \o/". Doing so implies that we know which interpolation type to use, etc... As we worked on that, we realised that we have opened an interesting can of worm and that it will also require a substantial amount of work.

So where we are now? Code wise, this is a prototype for the new spectral code: http://nbviewer.ipython.org/github/colour-science/colour-ramblings/blob/master/spectral_signal.ipynb We are only handling Spectrum object for now and need to discuss how we want to build the Spectra object. I got side-tracked those last months but I will get back to that in January.

The interpolation can of worm analysis draft is here: http://nbviewer.ipython.org/github/colour-science/colour-ramblings/blob/master/analysis_of_spectral_data_interpolation.ipynb

Voila! :)

Ron024 commented 8 years ago

ASTM E308 has a procedure for calculating tristimulius values with band pass corrections for 5 and 10nm sampling. This is the procedure most desktop spectrophotmers implement for their calculations.

KelSolaar commented 8 years ago

I'll take a look and see what they advise. The CIE usually recommends Sprague (1880) interpolation for uniformly spaced measurements and cubic spline for non uniform ones, although from our initial tests those are not necessarily the best interpolation methods.

KelSolaar commented 8 years ago

I will carefully re-read the designations but at first glance in ASTM E308−15 the computations are roughly done as follows:

ASTM E308−15 provides some tables with the tristimulus weighting factors for various illuminant observer combinations but refers to ASTM E2022-11 in order to construct the tables, this involves a Lagrange interpolator. The data is assumed to have been already corrected for spectral bandpass dependence.

Here is a related interesting thread:

References

KelSolaar commented 8 years ago

I have started to work on that a few days ago, the tristimulus weighting factors computation part is essentially done, I get exact same values than ASTM E308−15 for CIE Standard Illuminant A for 10nm and 20nm, I haven't checked other tables yet, however I need to double check but I couldn't get the same table as the one presented in ASTM E2022−11:

https://github.com/colour-science/colour-ramblings/blob/master/astm_e2022.ipynb

eyeforcolor commented 8 years ago

I apologize for the long delay in response (was on vacation).

Seems like you are already on top of this! I knew I was opening a can of worms here. As you mentioned above, there are at least two things that need to be addressed:

  1. Interpolation: I am not debating the trickiness in figuring out the correct interpolation, but I was a little surprised that zero filling was chosen over interpolation (any type)! Frankly, I had never seen this being done except maybe when extrapolating. On a related note, also wanted to post the link to a recent article on CRA from Dr. Luo's team, where in they have done some quantitative analysis of different interpolation and extrapolation methods as it relates to calculation of the tristimulus values.
  2. Integration: I did notice this. In the past I was able to get by using scipy and the excellent trapezoidal integration algorithm the library has to do these calculations (along with their cubic spline interpolation routines) when unequal wavelength data was involved. But, these were all for our own internal calculations and did not have the strict requirements that this excellent library is aiming for.

Look forward to seeing how the library evolves and let me know if and how I can be of help!

KelSolaar commented 8 years ago

No worries :)

Glad you opened the can of worms (It was actually opened internally but it is even better with outside feedback)!

I'm not sure why I choose zero filling originally and in most cases, I'm often aligning stuff to the CMFS prior usage. Thanks for the paper link I'll drop a serious eye onto it as we have done some work around that.

We have a Gitter room by the way, feel free to drop by if you have anything to talk about that doesn't warrant an issue!

KelSolaar commented 8 years ago

As I'm implementing the new code into the API (https://github.com/colour-science/colour/compare/feature/ASTM_E308_15), I remember why I originally choose the zero filling method, it was essentially for speed reasons. I have some performance issues (especially with CQS, CRI and CCT computations) and I'll likely need to add more caching mechanism.

KelSolaar commented 8 years ago

The new ASTM code has been merged into develop.