Quick math questions

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

x
The toes look cold. Brrr.

@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… :man_shrugging:

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.

1 Like

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 :thinking: ?

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

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.

2 Likes

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

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

@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
image

1 Like

@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

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.

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.

1 Like

@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.

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?

Yes. This is the G’MIC syntax.

ah got it, I just assumed it was C code :slight_smile:

1 Like

FYI, it is an expression separator, where the last assigned value is returned.

@Thanatomanic I found that this is wrong:

C:\Windows\System32>gmic echo {n=16;base=4;logb(n)=log(n)/log(base);n">"=base?(t=floor(logb(n));n+(n+1)*t+(base*(1-base*t))/(base-1)):n;}
[gmic]-0./ Start G'MIC interpreter.
[gmic]-0./ 40.666666666666664
[gmic]-0./ End G'MIC interpreter.

The correct is suppose to be 30.

EDIT: Oh, it was supposed to be in power function. My bad.

I went back to see if I can make the gradient in gmic, but this shows up as math processing error.

Edit: Now I see it in quote. Sorry about that.

I think I’m realizing the answer involve using newson-raphson method. This way I can directly map a Archimedes spiral ramp in 2D array.

Does anyone has any idea how the coefficient are found: How To Get A Spline's Length Using Gauss–Legendre Quadrature? | Noah Zuo's Blog

Derivative and something else, but missing something.

Read Gauss–Legendre quadrature - Wikipedia. See references.