YUV subsampling chroma channels in python

When halving chroma channels width and/or height, what's the correct way to subsamble? To take one chroma pixel for each 2x2 luma pixel, when sampling from a full-resolution chroma source, which chroma pi开发者_如何学Goxel do we take - top left? Average of all 4? Doesn't matter?

This is my current code

if field == 'top':
    i = 0
elif field == 'bottom':
    i = 1

U = fromstring(Udata, dtype='uint8', count=w/2*h).reshape(h,w/2)

# halve chroma height (it's already half-width from UYVY source) by line skipping
U = U[i::2]

# scale chroma by a factor of 0.5 (2x2 pixels in -> 1 pixel out)
U = (U[0::2, 0::2]>>2) + (U[0::2, 1::2]>>2) + (U[1::2, 0::2]>>2) + (U[1::2, 1::2]>>2) + \
  (((U[0::2, 0::2]%4)  + (U[0::2, 1::2]%4)  + (U[1::2, 0::2]%4)  + (U[1::2, 1::2]%4)) >> 2)

---

Ideally you would interpolate the data to produce the best possible image, however I think that your question is best answered with a quick overview of WHERE the luma and chroma samples are located (imagine a camera with one grid for luma and the other for chroma). For subsampling there are multiple standards for selecting the location of chroma relative to luma samples.

most often the chroma is either co-located with the "top left" luma pixel,

XO  X   XO  X

X   X   X   X
 
XO  X   XO  X
 
X   X   X   X

or is located in the center of the square,

X   X   X   X
  O      O
X   X   X   X

X   X   X   X
  O      O
X   X   X   X

or is located on the left side of the square between the two left luma samples.

X   X   X   X
O       O
X   X   X   X

X   X   X   X
O       O
X   X   X   X

notice that in all cases the chroma value co-located with at least some luma samples will be a combination of several chroma values, this is generally done with a FIR filter, which is why another filter is recommended for the reverse operation (however you can get OK results with averaging)

---

The correct way to convert from YUV 4:2:0 to YUV 4:2:2 is to use a 6-tap FIR filter.

The src for this comes from the mpeg2 reference implementation that you can find here http://www.mpeg.org/MPEG/video/mssg-free-mpeg-software.html

Go for the "mpegv12.zip" file

Implementation in C:

/* vertical 1:2 interpolation filter */
static void conv420to422(unsigned char* src, unsigned char* dst)
{
  int w, h, i, j, j2;
  int jm6, jm5, jm4, jm3, jm2, jm1, jp1, jp2, jp3, jp4, jp5, jp6, jp7;

  w = 352>>1;
  h = 288>>1;

  printf("hello \n");

  if (1)
  {
    /* intra frame */
    for (i=0; i<w; i++)
    {
      for (j=0; j<h; j++)
      {
    //printf("%d,%d\n", i, j);
        j2 = j<<1;
        jm3 = (j<3) ? 0 : j-3;
        jm2 = (j<2) ? 0 : j-2;
        jm1 = (j<1) ? 0 : j-1;
        jp1 = (j<h-1) ? j+1 : h-1;
        jp2 = (j<h-2) ? j+2 : h-1;
        jp3 = (j<h-3) ? j+3 : h-1;

        /* FIR filter coefficients (*256): 5 -21 70 228 -37 11 */
        /* New FIR filter coefficients (*256): 3 -16 67 227 -32 7 */
        dst[w*j2] =     Clip[(int)(  3*src[w*jm3]
                             -16*src[w*jm2]
                             +67*src[w*jm1]
                            +227*src[w*j]
                             -32*src[w*jp1]
                             +7*src[w*jp2]+128)>>8];

        dst[w*(j2+1)] = Clip[(int)(  3*src[w*jp3]
                             -16*src[w*jp2]
                             +67*src[w*jp1]
                            +227*src[w*j]
                             -32*src[w*jm1]
                             +7*src[w*jm2]+128)>>8];
      }
      src++;
      dst++;
    }
  }
  else
  {
    /* intra field */
    for (i=0; i<w; i++)
    {
      for (j=0; j<h; j+=2)
      {
        j2 = j<<1;

        /* top field */
        jm6 = (j<6) ? 0 : j-6;
        jm4 = (j<4) ? 0 : j-4;
        jm2 = (j<2) ? 0 : j-2;
        jp2 = (j<h-2) ? j+2 : h-2;
        jp4 = (j<h-4) ? j+4 : h-2;
        jp6 = (j<h-6) ? j+6 : h-2;

        /* Polyphase FIR filter coefficients (*256): 2 -10 35 242 -18 5 */
        /* New polyphase FIR filter coefficients (*256): 1 -7 30 248 -21 5 */
        dst[w*j2] = Clip[(int)(  1*src[w*jm6]
                         -7*src[w*jm4]
                         +30*src[w*jm2]
                        +248*src[w*j]
                         -21*src[w*jp2]
                          +5*src[w*jp4]+128)>>8];

        /* Polyphase FIR filter coefficients (*256): 11 -38 192 113 -30 8 */
        /* New polyphase FIR filter coefficients (*256):7 -35 194 110 -24 4 */
        dst[w*(j2+2)] = Clip[(int)( 7*src[w*jm4]
                             -35*src[w*jm2]
                            +194*src[w*j]
                            +110*src[w*jp2]
                             -24*src[w*jp4]
                              +4*src[w*jp6]+128)>>8];

        /* bottom field */
        jm5 = (j<5) ? 1 : j-5;
        jm3 = (j<3) ? 1 : j-3;
        jm1 = (j<1) ? 1 : j-1;
        jp1 = (j<h-1) ? j+1 : h-1;
        jp3 = (j<h-3) ? j+3 : h-1;
        jp5 = (j<h-5) ? j+5 : h-1;
        jp7 = (j<h-7) ? j+7 : h-1;

        /* Polyphase FIR filter coefficients (*256): 11 -38 192 113 -30 8 */
        /* New polyphase FIR filter coefficients (*256):7 -35 194 110 -24 4 */
        dst[w*(j2+1)] = Clip[(int)( 7*src[w*jp5]
                             -35*src[w*jp3]
                            +194*src[w*jp1]
                            +110*src[w*jm1]
                             -24*src[w*jm3]
                              +4*src[w*jm5]+128)>>8];

        dst[w*(j2+3)] = Clip[(int)(  1*src[w*jp7]
                             -7*src[w*jp5]
                             +30*src[w*jp3]
                            +248*src[w*jp1]
                             -21*src[w*jm1]
                              +5*src[w*jm3]+128)>>8];
      }
      src++;
      dst++;
    }
  }
}

and in Python:

def conv420to422(src, dst):
    """420 to 422 - vertical 1:2 interpolation filter """

    width = 352                  # 352
    height = 288                 # 288
    w = width >> 1               # 176
    h = height >> 1              # 144

    n = 0
    k = 0
    for i in range(0, w):
        for j in range(0, h):
            j2 = j<<1
            jm3 = 0 if (j < 3) else j - 3
            jm2 = 0 if (j < 2) else j - 2
            jm1 = 0 if (j < 1) else j - 1
            jp1 = j + 1 if (j < h - 1) else h - 1
            jp2 = j + 2 if (j < h - 2) else h - 1
            jp3 = j + 3 if (j < h - 3) else h - 1

            a = (3*src[n+w*jm3]-16*src[n+w*jm2]+67*src[n+w*jm1]+227*src[n+w*j]-32*src[n+w*jp1]+7*src[n+w*jp2]+128)>>8
            dst[k+w*j2] = clip(a)
            b = (3*src[n+w*jp3]-16*src[n+w*jp2]+67*src[n+w*jp1]+227*src[n+w*j]-32*src[n+w*jm1]+7*src[n+w*jm2]+128)>>8
            dst[k+w*(j2+1)] = clip(b)
        n += 1
        k += 1
    return dst