typedef struct {

  double m[5][4];

} ClMat;

typedef void (*ClmInitFunc)(ClMat *m);

void clmZero(ClMat *m) {

  for(int r = 0; r < 5; ++r)

    for(int c = 0; c < 4; ++c)

      m->m[r][c] = 0;

}

void clmIdentity(ClMat *m) {

  for(int r = 0; r < 5; ++r)

    for(int c = 0; c < 4; ++c)

      m->m[r][c] = r == c;

}

void clmCopy(ClMat *dst, ClMat *src)

  { if (dst != src) memcpy(dst, src, sizeof(*dst)); }

void clmMulColor(ClMat *m, double val) {

  for(int r = 0; r < 5; ++r)

    for(int c = 0; c < 4; ++c)

      m->m[r][c] = r == c ? (r < 3 ? val : 1) : 0;

}

void clmMulPix(double *dst, const ClMat *m, const double *pix) {

  if (dst == pix) {

    double p[4];

    memcpy(p, pix, sizeof(p));

    clmMulPix(dst, m, p);

    return;

  }

  dst[0] = m->m[0][0]*pix[0] + m->m[1][0]*pix[1] + m->m[2][0]*pix[2] + m->m[3][0]*pix[3] + m->m[4][0];

  dst[1] = m->m[0][1]*pix[0] + m->m[1][1]*pix[1] + m->m[2][1]*pix[2] + m->m[3][1]*pix[3] + m->m[4][1];

  dst[2] = m->m[0][2]*pix[0] + m->m[1][2]*pix[1] + m->m[2][2]*pix[2] + m->m[3][2]*pix[3] + m->m[4][2];

  dst[3] = m->m[0][3]*pix[0] + m->m[1][3]*pix[1] + m->m[2][3]*pix[2] + m->m[3][3]*pix[3] + m->m[4][3];

}

void clmMul(ClMat *dst, const ClMat *a, const ClMat *b) {

  if (dst == a || dst == b) {

    ClMat m;

    clmCopy(&m, dst);

    clmMul(dst, dst == a ? &m : a, dst == b ? &m : b);

    return;

  }

  clmMulPix(dst->m[0], a, b->m[0]);

  clmMulPix(dst->m[1], a, b->m[1]);

  clmMulPix(dst->m[2], a, b->m[2]);

  clmMulPix(dst->m[3], a, b->m[3]);

  clmMulPix(dst->m[4], a, b->m[4]);

}

void clmGray(ClMat *m) {

  for(int r = 0; r < 5; ++r)

    for(int c = 0; c < 4; ++c)

      m->m[r][c] = r < 3 && c < 3 ? 1/3.0 : r == c;

}

void clmToYUV(ClMat *m) {

  m->m[0][0] = 0.299; m->m[0][1] = -0.168736; m->m[0][2] =  0.5;      m->m[0][3] = 0;

  m->m[1][0] = 0.587; m->m[1][1] =  0.331264; m->m[1][2] = -0.418688; m->m[1][3] = 0;

  m->m[2][0] = 0.114; m->m[2][1] =  0.5;      m->m[2][2] = -0.081312; m->m[2][3] = 0;

  m->m[3][0] = 0;     m->m[3][1] =  0;        m->m[3][2] =  0;        m->m[3][3] = 1;

  m->m[4][0] = 0;     m->m[4][1] =  0;        m->m[4][2] =  0;        m->m[4][3] = 0;

}

void clmFromYUV(ClMat *m) {

  m->m[0][0] = 1;     m->m[0][1] =  1;        m->m[0][2] = 1;     m->m[0][3] = 0;

  m->m[1][0] = 0;     m->m[1][1] = -0.344136; m->m[1][2] = 1.772; m->m[1][3] = 0;

  m->m[2][0] = 0.402; m->m[2][1] = -0.714136; m->m[2][2] = 0;     m->m[2][3] = 0;

  m->m[3][0] = 0;     m->m[3][1] =  0;        m->m[3][2] = 0;     m->m[3][3] = 1;

  m->m[4][0] = 0;     m->m[4][1] =  0;        m->m[4][2] = 0;     m->m[4][3] = 0;

}

void clmScreenToCompositeYUV(ClMat *m) {

  // to YUV

  clmToYUV(m);

  // Alpha = Y-Luminance, Y-Luminance = 1

  m->m[0][3] = m->m[0][0];

  m->m[1][3] = m->m[1][0];

  m->m[2][3] = m->m[2][0];

  m->m[3][3] = m->m[3][0];

  m->m[0][0] = 0;

  m->m[1][0] = 0;

  m->m[2][0] = 0;

  m->m[4][0] = 1;

  // YUV -> RGB

  ClMat mm;

  clmFromYUV(&mm);

  clmMul(m, &mm, m);

}

void filterMatrix(Img *img, const ClMat *m) {

  if (!img->data) return;

  for(double *p = img->data, *e = p + 4*img->w*img->h; p < e; p += 4)

    clmMulPix(p, m, p);

}