CFA Strength and Color Discrimination

[Jack Hogan]

I have often wondered what people mean when they say that CFA strength has weakened in recent times, resulting in poorer color discrimination than in the old days.

How are 'strength' and 'discrimination' defined? Recently there was a suggestion that higher strength could be represented by narrower Spectral Sensitivity Functions in the CFA and that a ballpark proxy for discrimination could be provided by SMI, the Sensitivity Metamerism Index (just another way of showing the average deltaE off a number of color patches).

Below I show some preliminary findings* based on a D50 illuminant, a Colorchecker 24 target and simulated gaussian CFA SSFs of varying standard deviation. In a nutshell, absent any mistakes on my part* it appears that with my simulated setup SMI peaks when the standard deviation allows a fair amount of overlap - but not too much. This makes intuitive sense to me.



Below are the resulting dE values and CFA SSFs for standard deviations of 10, 20, 40, and 60 nm. Peak SMI with my setup occurred at a standard deviation of 39nm.


Standard Deviation of 10nm --> SMI 53.7


Standard Deviation of 20nm --> SMI 66.1


Standard Deviation of 40nm --> SMI 86.3


Standard Deviation of 60nm --> SMI 66.4

Here is the matrix from white balanced 'raw' data to XYZ for the highest SMI, at a standard deviation of 39nm:


Correlated Color Temperature is 4919K

In case anyone is wondering, below is the CFA obtained by averaging a number of Nikon camera SSFs, taken from RIT *.



Happy to discuss these findings and/or try different combinations should anyone be interested. And especially interested in obtaining feedback on the procedure I followed below.

Jack

*The low SMI obtained from the Nikon average CFA is somewhat surprising to me, especially given that I measured my own camera, a Nikon D610, and got an SMI of 85 (see the very bottom here ). Makes me wonder whether I made any procedural mistakes, please check it below:

1) D50 illuminant and CC24 reflectivity generated by matlab plug-in Optprop in energy units
2) multiply together illuminant and reflectivity for each CC24 patch
3) multiply the results of 2) by CIE 2006 2 degree LMS color matching functions to get XYZ values for each CC24 patch
4) convert results in 3) to Lab, these will be the reference values for the patches

5) multiply the result of 2) above by wavelength in order to convert energy to photons
6) generate gaussian CFA SSFs centered on 470, 530 and 590 nm with the desired standard deviation
7) multiply the result of 5) above by CFA SSFs to get 'raw' RGB values
8) white balance raw data based on patch D4
9) find the matrix that minimizes deltaE to reference Lab values in 4)
10) calculate SMI = 100-5.5*average deltaE of color patches in CC24

 

[JimKasson]

Good work, Jack. I think the next step might be to explore shapes of the curves.

Jim

 

[Iliah Borg]

There may be other reasons for SMI difference between your camera and calculations made from RIT data, too; but I never saw spectral response curves for sensors published on the web to align well with my own measurements. Part of it maybe that they don't read raw data correctly.

 

[Jack Hogan]

There may be other reasons for SMI difference between your camera and calculations made from RIT data, too; but I never saw spectral response curves for sensors published on the web to align well with my own measurements. Part of it maybe that they don't read raw data correctly.

Yeah, for one they seldom specify the units... On the other hand I've been known to make huge mistakes all by myself

 

[Iliah Borg]

There may be other reasons for SMI difference between your camera and calculations made from RIT data, too; but I never saw spectral response curves for sensors published on the web to align well with my own measurements. Part of it maybe that they don't read raw data correctly.

Yeah, for one they seldom specify the units...

I came across arbitrary and log scales, as well as biased data.

On the other hand I've been known to make huge mistakes all by myself

Who isn't. Goes under "other reasons" I started with.

--
http://www.libraw.org/

 

[J A C S]

You can try to mimic the human response by shifting the R,G,B Gaussians so that the peaks are where the human ones are (somewhere around 570, 540, 440) and their widths are close to the human ones. One could even do some curve fitting with those 6 parameters.

 

[JimKasson]

You can try to mimic the human response by shifting the R,G,B Gaussians so that the peaks are where the human ones are (somewhere around 570, 540, 440) and their widths are close to the human ones. One could even do some curve fitting with those 6 parameters.

Without messing with the widths:





Thanks to Jack for the code.

Jim

--
http://blog.kasson.com

 

[Iliah Borg]

You can try to mimic the human response by shifting the R,G,B Gaussians so that the peaks are where the human ones are (somewhere around 570, 540, 440) and their widths are close to the human ones.

You may find SMI is about 90 for "synthetic humans". As long as it it a matter of a matrix, simple shifting does not result in much of a change.

One could even do some curve fitting with those 6 parameters.

 

[J A C S]

You can try to mimic the human response by shifting the R,G,B Gaussians so that the peaks are where the human ones are (somewhere around 570, 540, 440) and their widths are close to the human ones.

You may find SMI is about 90 for "synthetic humans".

I am not sure what that means. Isn't LMS ideal by definition?

As long as it it a matter of a matrix, simple shifting does not result in much of a change.

It should. After all, if you shift too much, you can have R and G equal and you do not want that.

 

[JimKasson]

You can try to mimic the human response by shifting the R,G,B Gaussians so that the peaks are where the human ones are (somewhere around 570, 540, 440) and their widths are close to the human ones.

You may find SMI is about 90 for "synthetic humans".

That appears to be the case with Jack's program. Here's what I got when I set the CFA to the reference XYZ:







--

the last word the last word - Photography meets digital computer technology. Photography wins -- most of the time.

Photography meets digital computer technology. Photography wins -- most of the time.

blog.kasson.com

 

[J A C S]

You can try to mimic the human response by shifting the R,G,B Gaussians so that the peaks are where the human ones are (somewhere around 570, 540, 440) and their widths are close to the human ones.

You may find SMI is about 90 for "synthetic humans".

That appears to be the case with Jack's program. Here's what I got when I set the CFA to the reference XYZ:

Shouldn't you get 100% by definition?

 

[JimKasson]

You can try to mimic the human response by shifting the R,G,B Gaussians so that the peaks are where the human ones are (somewhere around 570, 540, 440) and their widths are close to the human ones.

You may find SMI is about 90 for "synthetic humans".

That appears to be the case with Jack's program. Here's what I got when I set the CFA to the reference XYZ:

Shouldn't you get 100% by definition?

I would think so, but I don't know what Jack is doing in any great detail.

Jack, what say you?

Forward Matrix to XYZ D50/2 =

1.0348 -0.1331 0.0641
0.0595 0.8915 0.0490
0.0794 -0.0988 0.8262

Shouldn't that be all zeros except for ones on the diagonals?

Here's how I got the CFA:

if XYZFlag
CFA = csvread('lin2012xyz2e_1_7sf.csv');
CFA,1) = [];
end

Jim

 

[Iliah Borg]

You can try to mimic the human response by shifting the R,G,B Gaussians so that the peaks are where the human ones are (somewhere around 570, 540, 440) and their widths are close to the human ones.

You may find SMI is about 90 for "synthetic humans".

That appears to be the case with Jack's program. Here's what I got when I set the CFA to the reference XYZ:

Sometimes it is consistency, especially when a non-CIE1931 observer is used, while the reference values are XYZ or Lab calculated with CIE1931, this happens.

Can be for a number of other reasons tho.

--
http://www.libraw.org/

 

[Jack Hogan]

You can try to mimic the human response by shifting the R,G,B Gaussians so that the peaks are where the human ones are (somewhere around 570, 540, 440) and their widths are close to the human ones.

You may find SMI is about 90 for "synthetic humans".

That appears to be the case with Jack's program. Here's what I got when I set the CFA to the reference XYZ:

Shouldn't you get 100% by definition?

I would think so, but I don't know what Jack is doing in any great detail.

Jack, what say you?

1) Reference XYZ is obtained by multiplying the D50 illuminant in units of energy by the reflectivity of the cc24 patches, then multiplying that by the linear energy CMFs.

2) The way the code is set up the D50 illuminant is in units of energy; it then is multiplied by wavelength to get to units proportional to photons; the photons are reflected based on patch reflectivity as-is, filtered by the CFAs as-is and counted/saved by the sensor in the raw data (we could perform the conversion to photons last and it would make no difference). A matrix is then found to approximate 1). Does this sound right?

If so I would guess that when using CMFs as CFAs there might be one conversion missing (or alternatively one too many): the linear CMFs refer to energy, not photons. Therefore to use them as CFAs I would think we need to compensate by dividing CMFs by wavelength.

In such a case the full trip would then be: D50 in units of energy; times wavelength to get to units of photons; times cc24 reflectance; times CMFs/wavelength = raw data. Find the matrix and we should have a perfect match to 1). Well, minus imperfections...

Makes sense or am I off to left field?

Forward Matrix to XYZ D50/2 =

1.0348 -0.1331 0.0641
0.0595 0.8915 0.0490
0.0794 -0.0988 0.8262

Shouldn't that be all zeros except for ones on the diagonals?

Here's how I got the CFA:

if XYZFlag
CFA = csvread('lin2012xyz2e_1_7sf.csv');
CFA,1) = [];
end

 

[JimKasson]

1) Reference XYZ is obtained by multiplying the D50 illuminant in units of energy by the reflectivity of the cc24 patches, then multiplying that by the linear energy CMFs.

2) The way the code is set up the D50 illuminant is in units of energy; it then is multiplied by wavelength to get to units proportional to photons; the photons are reflected based on patch reflectivity as-is, filtered by the CFAs as-is and counted/saved by the sensor in the raw data (we could perform the conversion to photons last and it would make no difference). A matrix is then found to approximate 1). Does this sound right?

If so I would guess that when using CMFs as CFAs there might be one conversion missing (or alternatively one too many): the linear CMFs refer to energy, not photons. Therefore to use them as CFAs I would think we need to compensate by dividing CMFs by wavelength.

In such a case the full trip would then be: D50 in units of energy; times wavelength to get to units of photons; times cc24 reflectance; times CMFs/wavelength = raw data. Find the matrix and we should have a perfect match to 1). Well, minus imperfections...

Makes sense or am I off to left field?

I tried weighting the XYZ CMFs as you suggested, but that didn't fix it. I PM'd you about that.

So then I tried changing the seed to

seed = diag([1 1 1]);

with the unweighted XYZ:



No joy.

Then I applied your suggested weighting:







Bingo!

The good news: this is a check on your methodology, and it looks good.

The bad news: the goal metric is multimodal. That means that we'll need stronger stuff than Nelder-Mead. I could try the genetic algorithm but that will make my code a lot less portable , since it will then need an exotic Matlab toolkit.

I await your advice.

Jim

--
http://blog.kasson.com

 

[Iliah Borg]

You can try to mimic the human response by shifting the R,G,B Gaussians so that the peaks are where the human ones are (somewhere around 570, 540, 440) and their widths are close to the human ones.

You may find SMI is about 90 for "synthetic humans".

That appears to be the case with Jack's program. Here's what I got when I set the CFA to the reference XYZ:

Shouldn't you get 100% by definition?

I would think so, but I don't know what Jack is doing in any great detail.

Jack, what say you?

1) Reference XYZ is obtained by multiplying the D50 illuminant in units of energy by the reflectivity of the cc24 patches, then multiplying that by the linear energy CMFs.

2) The way the code is set up the D50 illuminant is in units of energy; it then is multiplied by wavelength to get to units proportional to photons; the photons are reflected based on patch reflectivity as-is, filtered by the CFAs as-is and counted/saved by the sensor in the raw data (we could perform the conversion to photons last and it would make no difference). A matrix is then found to approximate 1). Does this sound right?

If so I would guess that when using CMFs as CFAs there might be one conversion missing (or alternatively one too many): the linear CMFs refer to energy, not photons. Therefore to use them as CFAs I would think we need to compensate by dividing CMFs by wavelength.

I would say your number 1 also doesn't mention normalization (observer y by illuminant SPD for each wavelength, summed).

 

[Iliah Borg]

1) Reference XYZ is obtained by multiplying the D50 illuminant in units of energy by the reflectivity of the cc24 patches, then multiplying that by the linear energy CMFs.

2) The way the code is set up the D50 illuminant is in units of energy; it then is multiplied by wavelength to get to units proportional to photons; the photons are reflected based on patch reflectivity as-is, filtered by the CFAs as-is and counted/saved by the sensor in the raw data (we could perform the conversion to photons last and it would make no difference). A matrix is then found to approximate 1). Does this sound right?

If so I would guess that when using CMFs as CFAs there might be one conversion missing (or alternatively one too many): the linear CMFs refer to energy, not photons. Therefore to use them as CFAs I would think we need to compensate by dividing CMFs by wavelength.

In such a case the full trip would then be: D50 in units of energy; times wavelength to get to units of photons; times cc24 reflectance; times CMFs/wavelength = raw data. Find the matrix and we should have a perfect match to 1). Well, minus imperfections...

Makes sense or am I off to left field?

I tried weighting the XYZ CMFs as you suggested, but that didn't fix it. I PM'd you about that.

So then I tried changing the seed to

seed = diag([1 1 1]);

with the unweighted XYZ:



No joy.

Then I applied your suggested weighting:



Bingo!

The good news: this is a check on your methodology, and it looks good.

The bad news: the goal metric is multimodal. That means that we'll need stronger stuff than Nelder-Mead. I could try the genetic algorithm but that will make my code a lot less portable , since it will then need an exotic Matlab toolkit.

From the standard:

"Capture the CIE Publication 13[6] eight patches", where [6] is CIE13.3-1995

Patches are ( http://munsell.com/faqs/list-of-colors-by-notation-name/ ) :

7,5R 6/4

5Y 6/4

5GY 6/8

2,5G 6/6

10BG 6/4

5PB 6/8

2,5P 6/8

10P 6/8

Calculate the tentatively optimized linear matrix using Equation (B.6):

A= T*S' * (S*S')^(-1) where ' is transpose, ^(-1) is inversion

T is XYZ for 8 patches, as a matrix

X1 ... Xi ... X8

Y1 ... Yi ... Y8

Z1 ... Zi ... Z8

S is a matrix of sensor outputs for respective patches

Now, a new set of estimated tristimulus values is calculated using:

Tk = A * Sk

--

About LibRaw | LibRaw

www.libraw.org

 

[Erik Kaffehr]

1) Reference XYZ is obtained by multiplying the D50 illuminant in units of energy by the reflectivity of the cc24 patches, then multiplying that by the linear energy CMFs.

2) The way the code is set up the D50 illuminant is in units of energy; it then is multiplied by wavelength to get to units proportional to photons; the photons are reflected based on patch reflectivity as-is, filtered by the CFAs as-is and counted/saved by the sensor in the raw data (we could perform the conversion to photons last and it would make no difference). A matrix is then found to approximate 1). Does this sound right?

If so I would guess that when using CMFs as CFAs there might be one conversion missing (or alternatively one too many): the linear CMFs refer to energy, not photons. Therefore to use them as CFAs I would think we need to compensate by dividing CMFs by wavelength.

In such a case the full trip would then be: D50 in units of energy; times wavelength to get to units of photons; times cc24 reflectance; times CMFs/wavelength = raw data. Find the matrix and we should have a perfect match to 1). Well, minus imperfections...

Makes sense or am I off to left field?

I tried weighting the XYZ CMFs as you suggested, but that didn't fix it. I PM'd you about that.

So then I tried changing the seed to

seed = diag([1 1 1]);

with the unweighted XYZ:



No joy.

Then I applied your suggested weighting:



Bingo!

The good news: this is a check on your methodology, and it looks good.

The bad news: the goal metric is multimodal. That means that we'll need stronger stuff than Nelder-Mead. I could try the genetic algorithm but that will make my code a lot less portable , since it will then need an exotic Matlab toolkit.

I await your advice.

Jim

--
http://blog.kasson.com

Hi Jack and Jim,

Excellent stuff!

What I would ask a bit about is the importance of the left peak on the red curve. How would rendition be affected by zeroing it out?

My layman's conclusion from the work presented:

• There seems to be a strong relation between SMI and colour rendition accuracy using an optimized matrix. No big surprise really, as SMI really measures accuracy.
• It is not obvious that small overlaps are advantageous to colour accuracy, at least not in the SMI sense.
• The shape and distribution of the curves matter. That is also expected.
• It also seems that it is possible to get very high SMI for a CFA design that would approximate the SML curves for human vision.

What is your take on the "Colour Science" presented by Phase One?


The standard sensor, which I would assume describes the IQ3100 MP seems to have very high response in green and blue channels in near infrared, probably because it lacks an IR filter. But it does have the left peak on the R-channel which I have not seen on any of RIT curves I have plotted.

Would be interesting to see SMI calculated for the Trichromatic sensor. It could of course be that the curves are just "artist's vision" reminiscent of NASA's images of exoplanets and black holes colliding.

Best regards

Erik

--
Erik Kaffehr
Website: http://echophoto.dnsalias.net
Magic uses to disappear in controlled experiments…
Gallery: http://echophoto.smugmug.com
Articles: http://echophoto.dnsalias.net/ekr/index.php/photoarticles

 

[Joofa]

1) Reference XYZ is obtained by multiplying the D50 illuminant in units of energy by the reflectivity of the cc24 patches, then multiplying that by the linear energy CMFs.

In such a case the full trip would then be: D50 in units of energy; times wavelength to get to units of photons; times cc24 reflectance; times CMFs/wavelength = raw data. Find the matrix and we should have a perfect match to 1). Well, minus imperfections...

I tried weighting the XYZ CMFs as you suggested, but that didn't fix it. I PM'd you about that.

So then I tried changing the seed to

seed = diag([1 1 1]);

with the unweighted XYZ:



No joy.

Then I applied your suggested weighting:



Bingo!

The good news: this is a check on your methodology, and it looks good.

I don't understand. You used XYZ CMFs for reference, and it appears then used the same CMFs as test. And, that of course will result into zero error. What is the surprise here???

--
Dj Joofa
http://www.djjoofa.com

 

[Jack Hogan]

1) Reference XYZ is obtained by multiplying the D50 illuminant in units of energy by the reflectivity of the cc24 patches, then multiplying that by the linear energy CMFs.

2) The way the code is set up the D50 illuminant is in units of energy; it then is multiplied by wavelength to get to units proportional to photons; the photons are reflected based on patch reflectivity as-is, filtered by the CFAs as-is and counted/saved by the sensor in the raw data (we could perform the conversion to photons last and it would make no difference). A matrix is then found to approximate 1). Does this sound right?

If so I would guess that when using CMFs as CFAs there might be one conversion missing (or alternatively one too many): the linear CMFs refer to energy, not photons. Therefore to use them as CFAs I would think we need to compensate by dividing CMFs by wavelength.

In such a case the full trip would then be: D50 in units of energy; times wavelength to get to units of photons; times cc24 reflectance; times CMFs/wavelength = raw data. Find the matrix and we should have a perfect match to 1). Well, minus imperfections...

Makes sense or am I off to left field?

I tried weighting the XYZ CMFs as you suggested, but that didn't fix it. I PM'd you about that.

So then I tried changing the seed to

seed = diag([1 1 1]);

with the unweighted XYZ:



No joy.

Then I applied your suggested weighting:



Bingo!

The good news: this is a check on your methodology, and it looks good.

The bad news: the goal metric is multimodal. That means that we'll need stronger stuff than Nelder-Mead. I could try the genetic algorithm but that will make my code a lot less portable , since it will then need an exotic Matlab toolkit.

From the standard:

"Capture the CIE Publication 13[6] eight patches", where [6] is CIE13.3-1995

Patches are ( http://munsell.com/faqs/list-of-colors-by-notation-name/ ) :

7,5R 6/4

5Y 6/4

5GY 6/8

2,5G 6/6

10BG 6/4

5PB 6/8

2,5P 6/8

10P 6/8

Calculate the tentatively optimized linear matrix using Equation (B.6):

A= T*S' * (S*S')^(-1) where ' is transpose, ^(-1) is inversion

T is XYZ for 8 patches, as a matrix

X1 ... Xi ... X8

Y1 ... Yi ... Y8

Z1 ... Zi ... Z8

S is a matrix of sensor outputs for respective patches

Now, a new set of estimated tristimulus values is calculated using:

Tk = A * Sk

--
http://www.libraw.org/

Thanks Iliah, yes good idea when the search doesn't work. I used that approach on the CC24 patches in the past but I discarded it because it produces a suboptimal matrix (it minimizes distance to the XYZ reference as opposed to Lab) . When determining the best compromise matrix I tend to use dE2000, here just dE , and I couldn't figure out how to work dE easily into the normal equation. Any ideas?