1 /*
2  *  Copyright (C) 2005-2018 Team Kodi
3  *  This file is part of Kodi - https://kodi.tv
4  *
5  *  SPDX-License-Identifier: LGPL-2.1-or-later
6  *  See LICENSES/README.md for more information.
7  */
8 
9 #include "ConversionMatrix.h"
10 
11 //------------------------------------------------------------------------------
12 // constants for primaries and transfers functions and color models
13 //------------------------------------------------------------------------------
14 
15 // source: https://www.khronos.org/registry/DataFormat/specs/1.2/dataformat.1.2.html#PRIMARY_CONVERSION
16 
17 struct ConvYCbCr
18 {
19   float Kr, Kb;
20 };
21 
22 struct Primaries
23 {
24   float primaries[3][2];
25   float whitepoint[2];
26 };
27 
28 const ConvYCbCr BT709YCbCr = {0.2126, 0.0722};
29 const ConvYCbCr BT601YCbCr = {0.299, 0.114};
30 const ConvYCbCr BT2020YCbCr = {0.2627, 0.0593};
31 const ConvYCbCr ST240YCbCr = {0.212, 0.087};
32 
33 const Primaries PrimariesBT709 = {{{0.640, 0.330}, {0.300, 0.600}, {0.150, 0.060}},
34   {0.3127, 0.3290} };
35 const Primaries PrimariesBT610_525 = {{{0.640, 0.340}, {0.310, 0.595}, {0.155, 0.070}},
36   {0.3127, 0.3290} };
37 const Primaries PrimariesBT610_625 = {{{0.640, 0.330}, {0.290, 0.600}, {0.150, 0.060}},
38   {0.3127, 0.3290} };
39 const Primaries PrimariesBT2020 = {{{0.708, 0.292}, {0.170, 0.797}, {0.131, 0.046}},
40   {0.3127, 0.3290} };
41 
42 //------------------------------------------------------------------------------
43 // Matrix helpers
44 //------------------------------------------------------------------------------
45 // source: http://timjones.io/blog/archive/2014/10/20/the-matrix-inverted
46 
47 template <unsigned Order>
48 float CalculateDeterminant(float (&src)[Order][Order]);
49 
50 template <unsigned Order>
GetSubmatrix(float (& dest)[Order-1][Order-1],float (& src)[Order][Order],unsigned row,unsigned col)51 void GetSubmatrix(float (&dest)[Order-1][Order-1], float (&src)[Order][Order], unsigned row, unsigned col)
52 {
53   unsigned colCount = 0;
54   unsigned rowCount = 0;
55 
56   for (unsigned i = 0; i < Order; i++)
57   {
58     if (i != row)
59     {
60       colCount = 0;
61       for (unsigned j = 0; j < Order; j++)
62       {
63         if (j != col)
64         {
65           dest[rowCount][colCount] = src[i][j];
66           colCount++;
67         }
68       }
69       rowCount++;
70     }
71   }
72 }
73 
74 template <unsigned Order>
CalculateMinor(float (& src)[Order][Order],unsigned row,unsigned col)75 float CalculateMinor(float (&src)[Order][Order], unsigned row, unsigned col)
76 {
77   float sub[Order-1][Order-1];
78   GetSubmatrix<Order>(sub, src, row, col);
79   return CalculateDeterminant<Order - 1>(sub);
80 }
81 
82 template <unsigned Order>
CalculateDeterminant(float (& src)[Order][Order])83 float CalculateDeterminant(float (&src)[Order][Order])
84 {
85   float det = 0.0f;
86 
87   for (unsigned i = 0; i < Order; i++)
88   {
89     // Get minor of element (0, i)
90     float minor = CalculateMinor<Order>(src, 0, i);
91 
92     // If this is an odd-numbered row, negate the value.
93     float factor = (i % 2 == 1) ? -1.0f : 1.0f;
94 
95     det += factor * src[0][i] * minor;
96   }
97   return det;
98 }
99 
100 template <>
CalculateDeterminant(float (& src)[2][2])101 float CalculateDeterminant<2>(float (&src)[2][2])
102 {
103   return src[0][0] * src[1][1] - src[0][1] * src[1][0];
104 }
105 
106 //------------------------------------------------------------------------------
107 // Matrix classes
108 //------------------------------------------------------------------------------
109 
110 template <unsigned Order>
CMatrix(float (& src)[Order][Order])111 CMatrix<Order>::CMatrix(float (&src)[Order][Order])
112 {
113   Copy(m_mat, src);
114 }
115 
116 template <unsigned Order>
CMatrix(float (& src)[Order-1][Order-1])117 CMatrix<Order>::CMatrix(float (&src)[Order-1][Order-1])
118 {
119   *this = src;
120 }
121 
122 template <unsigned Order>
operator =(const CMatrix & src)123 CMatrix<Order>& CMatrix<Order>::operator=(const CMatrix& src)
124 {
125   Copy(m_mat, src.m_mat);
126   return *this;
127 }
128 
129 template <unsigned Order>
operator =(const float (& src)[Order-1][Order-1])130 CMatrix<Order>& CMatrix<Order>::operator=(const float (&src)[Order-1][Order-1])
131 {
132   for (unsigned i=0; i<Order-1; ++i)
133     for (unsigned j=0; j<Order-1; ++j)
134       m_mat[i][j] = src[i][j];
135 
136   for (unsigned i=0; i<Order; ++i)
137     m_mat[i][Order-1] = 0;
138 
139   for (unsigned i=0; i<Order; ++i)
140     m_mat[Order-1][i] = 0;
141 
142   return *this;
143 }
144 
145 template <unsigned Order>
146 float (&CMatrix<Order>::Get())[Order][Order]
__anon951207f00102null147 {
148   return m_mat;
149 }
150 
151 template <unsigned Order>
Invert()152 CMatrix<Order> CMatrix<Order>::Invert()
153 {
154   CMatrix<Order> ret;
155   Invert(ret.m_mat, m_mat);
156   return ret;
157 }
158 
159 template <unsigned Order>
operator *(const CMatrix & other)160 CMatrix<Order> CMatrix<Order>::operator*(const CMatrix& other)
161 {
162   return *this * other.m_mat;
163 }
164 
165 template <unsigned Order>
operator *=(const CMatrix & other)166 CMatrix<Order> CMatrix<Order>::operator*=(const CMatrix& other)
167 {
168   CMatrix<Order> tmp = *this * other.m_mat;
169   *this = tmp;
170   return *this;
171 }
172 
173 template <unsigned Order>
operator *(const float (& other)[Order][Order])174 CMatrix<Order> CMatrix<Order>::operator*(const float (&other)[Order][Order])
175 {
176   CMatrix<Order> ret;
177   for (unsigned i=0; i<Order; ++i)
178     for (unsigned j=0; j<Order; ++j)
179       for (unsigned k=0; k<Order; ++k)
180         ret.m_mat[i][j] += m_mat[i][k] * other[k][j];
181 
182   return ret;
183 }
184 
185 template <unsigned Order>
Copy(float (& dst)[Order][Order],const float (& src)[Order][Order])186 void CMatrix<Order>::Copy(float (&dst)[Order][Order], const float (&src)[Order][Order])
187 {
188   for (unsigned i=0; i<Order; ++i)
189     for (unsigned j=0; j<Order; ++j)
190       dst[i][j] = src[i][j];
191 }
192 
193 template <unsigned Order>
Invert(float (& dst)[Order][Order],float (& src)[Order][Order])194 void CMatrix<Order>::Invert(float (&dst)[Order][Order], float (&src)[Order][Order])
195 {
196   // Calculate the inverse of the determinant of src.
197   float det = CalculateDeterminant<Order>(src);
198   float inverseDet = 1.0f / det;
199 
200   for (unsigned j = 0; j < Order; j++)
201   {
202     for (unsigned i = 0; i < Order; i++)
203     {
204       // Get minor of element (j, i) - not (i, j) because
205       // this is where the transpose happens.
206       float minor = CalculateMinor<Order>(src, j, i);
207 
208       // Multiply by (−1)^{i+j}
209       float factor = ((i + j) % 2 == 1) ? -1.0f : 1.0f;
210       float cofactor = minor * factor;
211 
212       dst[i][j] = inverseDet * cofactor;
213     }
214   }
215 }
216 
CGlMatrix(float (& src)[3][3])217 CGlMatrix::CGlMatrix(float (&src)[3][3]) : CMatrix<4>(src)
218 {
219 
220 }
221 
operator *(const float (& other)[4][4])222 CGlMatrix::CMatrix CGlMatrix::operator*(const float (&other)[4][4])
223 {
224   CGlMatrix ret;
225 
226   float (&left)[4][4] = m_mat;
227   const float (&right)[4][4] = other;
228 
229   ret.m_mat[0][0] = left[0][0] * right[0][0] + left[0][1] * right[1][0] + left[0][2] * right[2][0];
230   ret.m_mat[0][1] = left[0][0] * right[0][1] + left[0][1] * right[1][1] + left[0][2] * right[2][1];
231   ret.m_mat[0][2] = left[0][0] * right[0][2] + left[0][1] * right[1][2] + left[0][2] * right[2][2];
232   ret.m_mat[0][3] = left[0][0] * right[0][3] + left[0][1] * right[1][3] + left[0][2] * right[2][3] + left[0][3];
233   ret.m_mat[1][0] = left[1][0] * right[0][0] + left[1][1] * right[1][0] + left[1][2] * right[2][0];
234   ret.m_mat[1][1] = left[1][0] * right[0][1] + left[1][1] * right[1][1] + left[1][2] * right[2][1];
235   ret.m_mat[1][2] = left[1][0] * right[0][2] + left[1][1] * right[1][2] + left[1][2] * right[2][2];
236   ret.m_mat[1][3] = left[1][0] * right[0][3] + left[1][1] * right[1][3] + left[1][2] * right[2][3] + left[1][3];
237   ret.m_mat[2][0] = left[2][0] * right[0][0] + left[2][1] * right[1][0] + left[2][2] * right[2][0];
238   ret.m_mat[2][1] = left[2][0] * right[0][1] + left[2][1] * right[1][1] + left[2][2] * right[2][1];
239   ret.m_mat[2][2] = left[2][0] * right[0][2] + left[2][1] * right[1][2] + left[2][2] * right[2][2];
240   ret.m_mat[2][3] = left[2][0] * right[0][3] + left[2][1] * right[1][3] + left[2][2] * right[2][3] + left[2][3];
241 
242   return ret;
243 }
244 
CScale(float x,float y,float z)245 CScale::CScale(float x, float y, float z)
246 {
247   m_mat[0][0] = x;
248   m_mat[1][1] = y;
249   m_mat[2][2] = z;
250   m_mat[3][3] = 1;
251 }
252 
CTranslate(float x,float y,float z)253 CTranslate::CTranslate(float x, float y, float z)
254 {
255   m_mat[0][0] = 1;
256   m_mat[1][1] = 1;
257   m_mat[2][2] = 1;
258   m_mat[3][3] = 1;
259   m_mat[0][3] = x;
260   m_mat[1][3] = y;
261   m_mat[2][3] = z;
262 }
263 
264 //------------------------------------------------------------------------------
265 // Conversion classes
266 //------------------------------------------------------------------------------
267 
ConversionToRGB(float Kr,float Kb)268 ConversionToRGB::ConversionToRGB(float Kr, float Kb)
269 {
270   float Kg = 1-Kr-Kb;
271   a11 = Kr;
272   a12 = Kg;
273   a13 = Kb;
274   CbDen = 2*(1-Kb);
275   CrDen = 2*(1-Kr);
276 
277   m_mat[0][0] = a11;       m_mat[0][1] = a12;       m_mat[0][2] = a13;
278   m_mat[1][0] = -Kr/CbDen; m_mat[1][1] = -Kg/CbDen; m_mat[1][2] = 0.5;
279   m_mat[2][0] = 0.5;       m_mat[2][1] = -Kg/CrDen; m_mat[2][2] = -Kb/CrDen;
280 
281   CMatrix<3> inv(m_mat);
282   Copy(m_mat, inv.Invert().Get());
283 };
284 
operator =(const float (& src)[3][3])285 ConversionToRGB& ConversionToRGB::operator=(const float (&src)[3][3])
286 {
287   Copy(m_mat, src);
288   return *this;
289 }
290 
PrimaryToXYZ(const float (& primaries)[3][2],const float (& whitepoint)[2])291 PrimaryToXYZ::PrimaryToXYZ(const float (&primaries)[3][2], const float (&whitepoint)[2])
292 {
293   float By = CalcBy(primaries, whitepoint);
294   float Gy = CalcGy(primaries, whitepoint, By);
295   float Ry = CalcRy(By, Gy);
296 
297   m_mat[0][0] = Ry*primaries[0][0]/primaries[0][1];
298   m_mat[0][1] = Gy*primaries[1][0]/primaries[1][1];
299   m_mat[0][2] = By*primaries[2][0]/primaries[2][1];
300   m_mat[1][0] = Ry;
301   m_mat[1][1] = Gy;
302   m_mat[1][2] = By;
303   m_mat[2][0] = Ry/primaries[0][1] * (1- primaries[0][0] - primaries[0][1]);
304   m_mat[2][1] = Gy/primaries[1][1] * (1- primaries[1][0] - primaries[1][1]);
305   m_mat[2][2] = By/primaries[2][1] * (1- primaries[2][0] - primaries[2][1]);
306 }
307 
CalcBy(const float p[3][2],const float w[2])308 float PrimaryToXYZ::CalcBy(const float p[3][2], const float w[2])
309 {
310   float val = ((1-w[0])/w[1] - (1-p[0][0])/p[0][1]) * (p[1][0]/p[1][1] - p[0][0]/p[0][1]) -
311   (w[0]/w[1] - p[0][0]/p[0][1]) * ((1-p[1][0])/p[1][1] - (1-p[0][0])/p[0][1]);
312 
313   val /= ((1-p[2][0])/p[2][1] - (1-p[0][0])/p[0][1]) * (p[1][0]/p[1][1] - p[0][0]/p[0][1]) -
314   (p[2][0]/p[2][1] - p[0][0]/p[0][1]) * ((1-p[1][0])/p[1][1] - (1-p[0][0])/p[0][1]);
315 
316   return val;
317 }
318 
CalcGy(const float p[3][2],const float w[2],const float By)319 float PrimaryToXYZ::CalcGy(const float p[3][2], const float w[2], const float By)
320 {
321   float val = w[0]/w[1] - p[0][0]/p[0][1] - By * (p[2][0]/p[2][1] - p[0][0]/p[0][1]);
322   val /= p[1][0]/p[1][1] - p[0][0]/p[0][1];
323 
324   return val;
325 }
326 
CalcRy(const float By,const float Gy)327 float PrimaryToXYZ::CalcRy(const float By, const float Gy)
328 {
329   return 1.0 - Gy - By;
330 }
331 
PrimaryToRGB(float (& primaries)[3][2],float (& whitepoint)[2])332 PrimaryToRGB::PrimaryToRGB(float (&primaries)[3][2], float (&whitepoint)[2]) : PrimaryToXYZ(primaries, whitepoint)
333 {
334   CMatrix<3> inv(m_mat);
335   Copy(m_mat, inv.Invert().Get());
336 }
337 
338 //------------------------------------------------------------------------------
339 
SetColParams(AVColorSpace colSpace,int bits,bool limited,int textuteBits)340 void CConvertMatrix::SetColParams(AVColorSpace colSpace, int bits, bool limited, int textuteBits)
341 {
342   if (m_colSpace == colSpace &&
343       m_srcBits == bits &&
344       m_limitedSrc == limited &&
345       m_srcTextureBits == textuteBits &&
346       m_pMat)
347     return;
348 
349   m_colSpace = colSpace;
350   m_srcBits = bits;
351   m_limitedSrc = limited;
352   m_srcTextureBits = textuteBits;
353 
354   GenMat();
355 }
356 
SetColPrimaries(AVColorPrimaries dst,AVColorPrimaries src)357 void CConvertMatrix::SetColPrimaries(AVColorPrimaries dst, AVColorPrimaries src)
358 {
359   m_colPrimariesDst = dst;
360   m_colPrimariesSrc = src;
361 
362   if (m_colPrimariesDst != m_colPrimariesSrc)
363   {
364     Primaries primToRGB;
365     Primaries primToXYZ;
366     switch (m_colPrimariesSrc)
367     {
368       case AVCOL_PRI_BT709:
369         primToXYZ = PrimariesBT709;
370         m_gammaSrc = 2.2;
371         break;
372       case AVCOL_PRI_BT470BG:
373         primToXYZ = PrimariesBT610_625;
374         m_gammaSrc = 2.2;
375         break;
376       case AVCOL_PRI_SMPTE170M:
377       case AVCOL_PRI_SMPTE240M:
378         primToXYZ = PrimariesBT610_525;
379         m_gammaSrc = 2.2;
380         break;
381       case AVCOL_PRI_BT2020:
382         primToXYZ = PrimariesBT2020;
383         m_gammaSrc = 2.4;
384         break;
385       default:
386         primToXYZ = PrimariesBT709;
387         m_gammaSrc = 2.2;
388         break;
389     }
390     switch (m_colPrimariesDst)
391     {
392       case AVCOL_PRI_BT709:
393         primToRGB = PrimariesBT709;
394         m_gammaDst = 2.2;
395         break;
396       case AVCOL_PRI_BT470BG:
397         primToRGB = PrimariesBT610_625;
398         m_gammaDst = 2.2;
399         break;
400       case AVCOL_PRI_SMPTE170M:
401       case AVCOL_PRI_SMPTE240M:
402         primToRGB = PrimariesBT610_525;
403         m_gammaDst = 2.2;
404         break;
405       case AVCOL_PRI_BT2020:
406         primToRGB = PrimariesBT2020;
407         m_gammaDst = 2.4;
408         break;
409       default:
410         primToRGB = PrimariesBT709;
411         m_gammaDst = 2.2;
412         break;
413     }
414     PrimaryToXYZ toXYZ(primToXYZ.primaries, primToXYZ.whitepoint);
415     PrimaryToRGB toRGB(primToRGB.primaries, primToRGB.whitepoint);
416 
417     CMatrix<3> tmp = toRGB*toXYZ;
418     m_pMatPrim.reset(new CMatrix<3>(tmp));
419   }
420   else
421   {
422     m_pMatPrim.reset();
423   }
424 }
425 
SetParams(float contrast,float black,bool limited)426 void CConvertMatrix::SetParams(float contrast, float black, bool limited)
427 {
428   m_contrast = contrast;
429   m_black = black;
430   m_limitedDst = limited;
431 }
432 
GenMat()433 void CConvertMatrix::GenMat()
434 {
435   ConvYCbCr convYCbCr;
436   switch (m_colSpace)
437   {
438     case AVCOL_SPC_BT709:
439       convYCbCr = BT709YCbCr;
440       break;
441     case AVCOL_SPC_BT470BG:
442     case AVCOL_SPC_SMPTE170M:
443       convYCbCr = BT601YCbCr;
444       break;
445     case AVCOL_SPC_SMPTE240M:
446       convYCbCr = ST240YCbCr;
447       break;
448     case AVCOL_SPC_BT2020_NCL:
449     case AVCOL_SPC_BT2020_CL:
450       convYCbCr = BT2020YCbCr;
451       break;
452     default:
453       convYCbCr = BT709YCbCr;
454       break;
455   }
456 
457   ConversionToRGB mConvRGB(convYCbCr.Kr, convYCbCr.Kb);
458 
459   m_pMat.reset(new CGlMatrix(mConvRGB.Get()));
460 
461   CTranslate trans(0, -0.5, -0.5);
462   *m_pMat *= trans;
463 
464   if (m_limitedSrc)
465   {
466     if (m_srcBits >= 12)
467     {
468       CScale scale(4080.0f / (3760 - 256), 4080.0f / (3840 - 256), 4080.0f / (3840 - 256));
469       CTranslate trans(- 256.0f / 4080, - 256.0f / 4080, - 256.0f / 4080);
470       *m_pMat *= scale;
471       *m_pMat *= trans;
472     }
473     else if (m_srcBits == 10)
474     {
475       CScale scale(1020.0f / (940 - 64), 1020.0f / (960 - 64), 1020.0f / (960 - 64));
476       CTranslate trans(- 64.0f / 1020, - 64.0f / 1020, - 64.0f / 1020);
477       *m_pMat *= scale;
478       *m_pMat *= trans;
479     }
480     else
481     {
482       CScale scale(255.0f / (235 - 16), 255.0f / (240 - 16), 255.0f / (240 - 16));
483       CTranslate trans(- 16.0f / 255, - 16.0f / 255, - 16.0f / 255);
484       *m_pMat *= scale;
485       *m_pMat *= trans;
486     }
487   }
488 
489   if (m_srcTextureBits > 8)
490   {
491     float val = 65535.0f / ((1 << m_srcTextureBits) - 1);
492     CScale scale(val, val, val);
493     *m_pMat *= scale;
494   }
495 }
496 
GetYuvMat(float (& mat)[4][4])497 void CConvertMatrix::GetYuvMat(float (&mat)[4][4])
498 {
499   if (!m_pMat)
500     return;
501 
502   CScale contrast(m_contrast, m_contrast, m_contrast);
503   CTranslate black(m_black, m_black, m_black);
504 
505   CGlMatrix ret = contrast;
506   ret *= black;
507   if (m_limitedDst)
508   {
509     float valScale = (235 - 16) / 255.0f;
510     float valTrans = 16.0f / 255;
511     CScale scale(valScale, valScale, valScale);
512     CTranslate trans(valTrans, valTrans, valTrans);
513     ret *= trans;
514     ret *= scale;
515   }
516 
517   ret *= m_pMat->Get();
518   float (&src)[4][4] = ret.Get();
519 
520   for (int i=0; i<4; ++i)
521     for (int j=0; j<4; ++j)
522       mat[i][j] = src[j][i];
523 
524   mat[0][3] = 0.0f;
525   mat[1][3] = 0.0f;
526   mat[2][3] = 0.0f;
527   mat[3][3] = 1.0f;
528 }
529 
GetPrimMat(float (& mat)[3][3])530 bool CConvertMatrix::GetPrimMat(float (&mat)[3][3])
531 {
532   if (!m_pMatPrim)
533     return false;
534 
535   float (&src)[3][3] = m_pMatPrim->Get();
536 
537   for (int i=0; i<3; ++i)
538     for (int j=0; j<3; ++j)
539       mat[i][j] = src[j][i];
540 
541   return true;
542 }
543 
GetGammaSrc()544 float CConvertMatrix::GetGammaSrc()
545 {
546   return m_gammaSrc;
547 }
548 
GetGammaDst()549 float CConvertMatrix::GetGammaDst()
550 {
551   return m_gammaDst;
552 }
553 
GetRGBYuvCoefs(AVColorSpace colspace,float (& coefs)[3])554 bool CConvertMatrix::GetRGBYuvCoefs(AVColorSpace colspace, float (&coefs)[3])
555 {
556   switch (colspace)
557   {
558     case AVCOL_SPC_BT709:
559       coefs[0] = BT709YCbCr.Kr;
560       coefs[1] = 1 - BT709YCbCr.Kr - BT709YCbCr.Kb;
561       coefs[2] = BT709YCbCr.Kb;
562       break;
563     case AVCOL_SPC_BT470BG:
564     case AVCOL_SPC_SMPTE170M:
565       coefs[0] = BT601YCbCr.Kr;
566       coefs[1] = 1 - BT601YCbCr.Kr - BT601YCbCr.Kb;
567       coefs[2] = BT601YCbCr.Kb;
568       break;
569     case AVCOL_SPC_SMPTE240M:
570       coefs[0] = ST240YCbCr.Kr;
571       coefs[1] = 1 - ST240YCbCr.Kr - ST240YCbCr.Kb;
572       coefs[2] = ST240YCbCr.Kb;
573       break;
574     case AVCOL_SPC_BT2020_NCL:
575     case AVCOL_SPC_BT2020_CL:
576       coefs[0] = BT2020YCbCr.Kr;
577       coefs[1] = 1 - BT2020YCbCr.Kr - BT2020YCbCr.Kb;
578       coefs[2] = BT2020YCbCr.Kb;
579       break;
580     default:
581       return false;
582   }
583   return true;
584 }
585