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);
}