Quick math questions

snibgo · April 10, 2021, 8:26pm

Some links I have found useful:

CIE 1931 color space - Wikipedia Wikipedia: CIE 1931 color space

White point - Wikipedia Wikipedia: White point

Standard illuminant - Wikipedia Wikipedia: Standard illuminant

http://www.brucelindbloom.com/ Bruce Lindbloom. The maths behind colour science.

https://ninedegreesbelow.com/ Nine Degrees Below Photograpy, Elle Stone.
Articles and tutorials on ICC profile color management and free/libre image editing. Includes a large collection of profiles.

A Standard Default Color Space for the Internet - sRGB A Standard Default Color Space for the Internet - sRGB, 1996. Proposed standard.

http://color.org/chardata/rgb/sRGB.pdf How to interpret the sRGB color space (specified in IEC 61966-2-1) for ICC profiles

http://www2.units.it/ipl/students_area/imm2/files/Colore1/sRGB.pdf IEC/4WD 61966-2-1: Colour Measurement and Management in Multimedia Systems and Equipment - Part 2-1: Default RGB Colour Space - sRGB

In quest of well behaved working spaces Nine Degrees Below Photograpy: In quest of well behaved working spaces, Elle Stone.

afre · April 10, 2021, 8:53pm

I am aware of some of these links. I guess the best way is to learn by doing. I will start an exercises thread if need be (when I have more energy).

In the meantime, (if any of you want to participate,)

That is the transform matrix (I think) from xyz to rec2020 with D50 illuminant. As an exercise, I want to change the illuminant to E.

age · April 10, 2021, 9:49pm

Yes.
For me the first thing to try as exercise is to find the srgb (d65) to xyz and the xyz to srgb (d65) matrix given the xy chromaticity coordinates

Here the pratical exercise
Linear Transformation Example for sRGB Space

The second exercise would be find the rec.2020 (d65) to xyz matrix given the xy chromaticity coordinates.

The next step would be learn how to do the Bradford chromatic adaptation transform

snibgo · April 10, 2021, 10:12pm

For Chromatic Adaptation, a useful source is Welcome to Bruce Lindbloom's Web Site

Given the (X,Y,Z) of the reference whites of the source and destination, and cone response matrix MA and its inverse MA^-1, we calculate the linear transformation matrix M.

Lindbloom gives values of MA and MA^-1 for three methods: XYZ scaling, Bradford and Von Kries. Other methods include Sharp, CMCCAT2000, CAT02, BS, BS-PC and Fairchild.

For example, to convert from D65 to E, using the Bradford method, Lindbloom gives M:

 1.0502616  0.0270757 -0.0232523
 0.0390650  0.9729502 -0.0092579
-0.0024047  0.0026446  0.9180873

Consider the numbers on the diagonal. Approx 1.05 will slightly increase the value in the red channel; 0.97 will slightly decrease green, and 0.918 will decrease blue. The overall transformation will make the image slightly redder.

According to Standard illuminant - Wikipedia, illuminant E is a theoretical illuminant roughly similar to D55. We can see from File:Planckian-locus.png - Wikipedia that D55 is slightly redder than D65. So this confirms that the matrix M shifts colours in the right direction.

The matrix xyz_rec2020[3][3] can be multiplied by M to give a conversion from XYZ to REC2020 with standard illuminant E.

afre · April 10, 2021, 10:27pm

In this case, it is from D50 to E.

Thanks for pointing out what happens to the colour. It reinforces the learning process. Indeed, E and D55 are similar. I hope to be able to hop between arbitrary illuminants.

snibgo · April 10, 2021, 10:51pm

Sorry, I misread your post. For D50 to E, Lindbloom gives M as:

 1.0025535  0.0036238  0.0359837
 0.0096914  0.9819125  0.0105947
 0.0089181 -0.0160789  1.2208770

So this shifts lightly to the blue.

Lindbloom has calculated M for our convenience. For programatic calculation, we don’t want to store all those values for M. We do need to store the white points of whatever illuminants we care about, eg from Standard illuminant - Wikipedia , and we need the cone response matrices for whatever methods we want to use.

snibgo · April 10, 2021, 11:02pm

Whenever I see a 3x3 colour transformation matrix, I ask myself: Are the diagonals approx 1.0, and the other values approx 0.0? If so, then it will cause a subtle colour change. Then, which diagonals are above 1.0 and which are below 1.0? This gives the colour shift.

And we can readily see the visual effect of the transformation (Windows BAT syntax):

rem From D50 to E.

set SMAT=^
1.0025535,0.0036238,0.0359837,^
0.0096914,0.9819125,0.0105947,^
0.0089181,-0.0160789,1.2208770

magick ^
  toes.png ^
  -color-matrix %SMAT% ^
  x.jpg

The toes look cold. Brrr.

Thanatomanic · April 11, 2021, 5:57am

@afre Following the maths on Bruce’s site I arrive at this (fairly straightforward) method to compute the example matrix you picked out. Note that for some odd reason, the variable name in RT is not what I expected. As I understand Bruce’s site, this matrix is used to convert an RGB value in Rec.2020 to a XYZ value and not the other way around. Maybe I’m wrong…

Anyway, here’s the maths. Forgive me for using slightly different notation… force of habit.

Pick your origin color space, here Rec.2020, and note its primaries and white point in xy chromaticity coordinates:

(x_R,y_R) =(0.708,0.292)\\ (x_G,y_G) = (0.170,0.797)\\ (x_B,y_B) = (0.131,0.046)\\ (x_W,y_W) = (0.3127,0.3290)

These values are taken from Wikipedia. The white point is D65 and therefore needs adaptation to D50. Adaptation is a simple matrix multiplication of (X,Y,Z) values with the respective transformation matrix \mathbf{M_A} (see Bruce’s table). For going from D65 to D50 we have,

\mathbf{M_A}=\begin{bmatrix} 1.0478112 & 0.0228866 & -0.0501270\\ 0.0295424 &0.9904844 & -0.0170491 \\ -0.0092345 & 0.0150436 & 0.7521316 \end{bmatrix}

We can now calculate Bradford-adapted (X,Y,Z) values from the chromaticity coordinates:

\begin{bmatrix} X_R'\\ Y_R'\\ Z_R' \end{bmatrix} = \mathbf{M_A} \begin{bmatrix} x_R / y_R\\ 1\\ (1 - x_R - y_R) / y_R \end{bmatrix}

Do the same for the green, blue and white coordinates.
Finally, the conversion matrix \mathbf{M} to go from RGB to XYZ is given here:

\begin{bmatrix} S_R\\ S_G\\ S_B \end{bmatrix}= \begin{bmatrix} X_R' & Y_R' & Z_R'\\ X_G' & Y_G' & Z_G'\\ X_B' & Y_B' & Z_B' \end{bmatrix}^{-1} \begin{bmatrix} X_W'\\ Y_W'\\ Z_W' \end{bmatrix}\\ \mathbf{M} = \begin{bmatrix} S_R X_R' & S_G X_G' & S_B X_B'\\ S_R Y_R' & S_G Y_G' & S_B Y_B'\\ S_R Z_R' & S_G Z_G' & S_B Z_B' \end{bmatrix}

When plugging in the above numbers, I end up with the following matrix, with all values rounded to 7 decimal places:

\mathbf{M}_\textrm{Rec. 2020 to XYZ} = \begin{bmatrix} 0.6734241 & 0.1656411 & 0.1251286\\ 0.2790177 & 0.6753402 & 0.0456377\\ -0.0019300 & 0.0299784 & 0.7973330 \end{bmatrix}

Which is identical to the values in RT’s source code.

age · April 11, 2021, 1:05pm

So this is how I have interpreted the Bruce Lindbloom page, RGB/XYZ Matrices
The difference is that he starts with the XYZ value for white point while I start from the xy coordinates for white point too.

rgbxyz-brucelindbloom.py (1.3 KB)

I think it should be easy to translate in g’mic language.

With numpy I use np.linalg.inv for a matrix inversion, how to do this task with g’mic ?

age · April 11, 2021, 3:13pm

Here finally my rec.2020 to xyz d50

rgb d65 to xyz d65
[[6.36958048e-01 1.44616904e-01 1.68880975e-01]
[2.62700212e-01 6.77998072e-01 5.93017165e-02]
[4.99410657e-17 2.80726930e-02 1.06098506e+00]]

xyz d65 to rgb d65
[[ 1.71665119 -0.35567078 -0.25336628]
[-0.66668435 1.61648124 0.01576855]
[ 0.01763986 -0.04277061 0.94210312]]

rgb d65 to xyz d50
[[ 0.67012458 0.16035917 0.09262633]
[ 0.29474229 0.67845094 0.01987506]
[-0.00896833 0.0437666 0.79752177]]

They are different from the rawtherapee matrices but I’m confident that my code is correct because I could match the d65 matrices here
https://colour.readthedocs.io/en/feature-read_the_docs/colour.models.dataset.rec_2020.html

Rawtherapee has matrices more close to this online calculator
http://www.russellcottrell.com/photo/matrixCalculator.htm

Thanatomanic · April 11, 2021, 4:19pm

I believe one of us has mixed up the order of operations. In your final step you have rec2020toxyzd50=np.matmul(M,d65tod50) , so you seem to do \mathbf{M} \cdot \mathbf{M_A} while I effectively do the reverse, \mathbf{M_A} \cdot \mathbf{M}. Since matrix multiplication doesn’t commute, we get different results.

age · April 11, 2021, 6:35pm

Ah!
Thanks @Thanatomanic
Fixed
rgbxyz-brucelindbloom-2020.py (1.5 KB)

Reptorian · May 12, 2021, 2:37am

I have a challenging question (there’s probably a simple solution.)

Using pixel coordinate space: I would like to be able to generate a quadrilateral gradient that is basically a “loft of two lines”.

To expand, I’d like a solution to this:

Note that the each points do have the same colors, but as you approach the center of colinear lines, the colors are very much different.

EDIT: Added separate channels.

How can I generate arbitrary quads without that problem in pixel-coordinate space? Every point in colinear lines should have the same color basically.

It should look more like this, but with shifted center point of quads with C0 continuity:

EDIT: I think this will help - math - Relative position of a point within a quadrilateral - Stack Overflow

Thanatomanic · May 12, 2021, 5:44am

@Reptorian Indeed, to figure out the values at arbitrary points you’ll have to solve a quadratic equation to go from global (arbitrary quadrilateral) to local (square) coordinates. This might help too: Inigo Quilez :: fractals, computer graphics, mathematics, shaders, demoscene and more

Then just basically do this for all four quadrants? I’ll see if I can cook up a working example…

Edit: this works, but I’m not satisfied with my code yet to share it

Reptorian · May 12, 2021, 4:24pm

@Thanatomatic, I just had finished testing the result. It works, almost. Some problem in the mismatching face, it isn’t continuous like the other one which is not acceptable, and two, it doesn’t match. I will try with the other approach to see if it gets anywhere. Continuous output is more important on the mismatching face. In the case of no solution to this problem, I think this still can be a new filter for gmic.

Test Code

$ sp cat +rep_powertweak 55%,60%
rep_powertweak:

f "begin(
  const ww=w-1;
  const hh=h-1;
  const cx=ww/2;
  const cy=hh/2;
  const pc=.5;
  const px=.5+$1*.5;
  const py=1-(.5+$2*.5);
  const pc_px=ww*px;
  const pc_py=hh*py;
  const inc_x1=(pc_py-cy)/pc_px;
  const inc_x2=(pc_py-cy)/(ww-pc_px);
  const inc_y1=(pc_px-cx)/pc_py;
  const inc_y2=(pc_px-cx)/(hh-pc_py);
  xcross(a,b)=a[0]*b[1]-a[1]*b[0];
  invBilinear(p,a,b,c,d)=(
   
   e = b-a;
   f = d-a;
   g = a-b+c-d;
   h = p-a;
   
   k2 = xcross( g, f );
   k1 = xcross( e, f ) + xcross( h, g );
   k0 = xcross( h, e );
   
   if(abs(k2)<0.001,
    [(h[0]*k1+f[0]*k0)/(e[0]*k1-g[0]*k0),-k0/k1];
   ,
    m = k1*k1 - 4.0*k0*k2;
    if(m<0.0,[-1,-1],
     m = sqrt( m );
     ik2 = 0.5/k2;
     v = (-k1 - m)*ik2;
     u = (h[0] - f[0]*v)/(e[0] + g[0]*v);
     v = (-k1 + m)*ik2;
     u = (h[0] - f[0]*v)/(e[0] + g[0]*v);
     [u,v];
    );
   );
  );
);
xp=x/w;
yp=y/h;

side_y=x<pc_px?(y+inc_x1*(pc_px-x))>pc_py:(y+inc_x2*(x-pc_px))>pc_py;
side_x=y<pc_py?(x+inc_y1*(pc_py-y))>pc_px:(x+inc_y2*(y-pc_py))>pc_px;
side=side_y*2+side_x;
 
p=[xp,yp];

uv_s0=invBilinear(p,[0,0],[.5,0],[px,py],[0,.5]);
uv_s1=invBilinear(p,[.5,0],[1,0],[1,.5],[px,py]);
uv_s2=invBilinear(p,[0,.5],[px,py],[.5,1],[0,1]);
uv_s3=invBilinear(p,[px,py],[1,.5],[1,1],[.5,1]);

uv_s0=[uv_s0[0]*cx,   uv_s0[1]*cy];
uv_s1=[uv_s1[0]*cx+cx,uv_s1[1]*cy];
uv_s2=[uv_s2[0]*cx,   uv_s2[1]*cy+cy];
uv_s3=[uv_s3[0]*cx+cx,uv_s3[1]*cy+cy];

side==3?I(uv_s3[0],uv_s3[1],0,1,3):
side==2?I(uv_s2[0],uv_s2[1],0,1,3):
side==1?I(uv_s1[0],uv_s1[1],0,1,3):
        I(uv_s0[0],uv_s0[1],0,1,3);
"

EDIT: Another test image

Reptorian · January 20, 2022, 6:18pm

I got a new challenging math question and programming question. I would like to know how to find out the number of digits from numbers concatenated with a base where the counts of number is defined in base 10.

concat_consec_digits_count(n in base 10,b);

concat_consec_digits_count(23,12):
1-2-3-4-5-6-7-8-9-10-11-12-13-14-15-16-17-18-19-20-21-22-23
1-2-3-4-5-6-7-8-9-A -B -10-11-12-13-14-15-16-17-18-19-1A-1B
concat_consec_digits_count(23,12)=11+12*2=35

concat_consec_digits_count(12,5):
1-2-3-4- 5 - 6 - 7 -  8 - 9 -10-11-12
1-2-3-4- 10- 11-12 - 13 -14 -20-21-22
concat_consec_digits_count(12,5)=4*1+8*2=20

In a chatroom, someone found the answer in OEIS - A058183 - OEIS . All I need to do is to replace log_10 with log_base.

Thanatomanic · January 21, 2022, 7:09am

Is this still a challenge if you found a sequence on OEIS? Edit: I didn’t look at the sequence, but tried deriving it myself. I came up with something equivalent, though the formula on OEIS seems somewhat shorter.

After some manual inspection of such base-digit-sequences, I came up with an expression for the length L of the string of concatenated digits in base b from 1 to the largest number with n digits:

T(b,n) = \frac{1-b^n-nb^n+nb^{n+1}}{b-1}

So, for example, there are T(5,2) = 44 characters in the string “12341011121314202122232430313233344041424344” and T(3,3) = 68 characters in the string “12101112202122100101102110111112120121122200201202210211212220221222”.

To find the length S of the string for an arbitrary number N in base b, we first need to know the number of digits d for that number, which is given by:

d = \lfloor\log_b N\rfloor + 1

The length of the string is shortened by (N_\text{max}-N)d digits, where N_\text{max} is the largest number with d digits, N_\text{max} = b^d-1.

We finally have,

S(N,b,d) = T(b,d) - (b^d-1-N)d

Or, to match your definition of concat_consec_digits_count, we have after some simplification and rearrangement:

S(N,b) = N+(N+1)\lfloor\log_bN\rfloor+\frac{b\left(1-b^{\lfloor\log_bN\rfloor}\right)}{b-1}

It’s easy to verify that S(23,12) = 35 and S(12,5)=20 as in your examples.

Reptorian · January 21, 2022, 8:38pm

@Thanatomanic I was able to find a optimized solution for element-wise subtraction with that definition of concat_conset_digits_count. It doesn’t work if N is lower than b. This is my final formula:

concat_consec_digits_count(v)=v>=base?(
    t=floor(logb(v));
    v+(v+1)*t+(base*(1-base^t))/(base-1);
  )
  :v;

So, if v is less than base, then the answer is v, otherwise, it use that formula.

EDIT: Changed (1-base*t) to (1-base^t), and now it works for all cases.

Matt_Maguire · January 21, 2022, 10:31pm

Is it allowed to use a semi-colon (;) like that in the middle of a ? : expression? I’m pretty sure it was forbidden back in the days when I learned C, and we would probably have used a comma (,) instead. But maybe the C standards have evolved since then?

Reptorian · January 21, 2022, 10:32pm

Yes. This is the G’MIC syntax.