image.c

JM 11.0 KTA 2.1 Source Code

  // cropping for chroma
  if (p->frame_cropping_flag)
  {
    crop_left   = p->frame_cropping_rect_left_offset;
    crop_right  = p->frame_cropping_rect_right_offset;
    crop_top    = ( 2 - p->frame_mbs_only_flag ) *  p->frame_cropping_rect_top_offset;
    crop_bottom = ( 2 - p->frame_mbs_only_flag ) *   p->frame_cropping_rect_bottom_offset;
  }
  else
  {
    crop_left = crop_right = crop_top = crop_bottom = 0;
  }

  if ((p->chroma_format_idc==YUV400) && input->write_uv)
  {
    size_x_cr = p->size_x/2;
    size_y_cr = p->size_y/2;
  }
  else
  {
    size_x_cr = p->size_x_cr - crop_left - crop_right;
    size_y_cr = p->size_y_cr - crop_top  - crop_bottom;
  }

  framesize_in_bytes = (((int64)size_y*size_x) + ((int64)size_y_cr*size_x_cr)*2) * symbol_size_in_bytes;

  if (psnrPOC==0 && img->psnr_number)
    img->idr_psnr_number = img->number*img->ref_poc_gap/(input->poc_scale);

  img->psnr_number=max(img->psnr_number,img->idr_psnr_number+psnrPOC);

  frame_no = img->idr_psnr_number+psnrPOC;

  // KS: this buffer should actually be allocated only once, but this is still much faster than the previous version
  buf = malloc ( size_y * size_x * symbol_size_in_bytes );

  if (NULL == buf)
  {
    no_mem_exit("find_snr: buf");
  }

  status = lseek (p_ref, framesize_in_bytes * frame_no, SEEK_SET);
  if (status == -1)
  {
    fprintf(stderr, "Error in seeking frame number: %d\n", frame_no);
    free (buf);
    return;
  }

  if(rgb_output)
    lseek (p_ref, framesize_in_bytes/3, SEEK_CUR);

  read(p_ref, buf, size_y * size_x * symbol_size_in_bytes);
  buf2img(imgY_ref, buf, size_x, size_y, symbol_size_in_bytes);

  if (p->chroma_format_idc != YUV400)
  {
    for (uv=0; uv < 2; uv++)
    {
      if(rgb_output && uv==1)
        lseek (p_ref, -framesize_in_bytes, SEEK_CUR);

      read(p_ref, buf, size_y_cr * size_x_cr*symbol_size_in_bytes);
      buf2img(imgUV_ref[uv], buf, size_x_cr, size_y_cr, symbol_size_in_bytes);
    }
  }

  if(rgb_output) 
    lseek (p_ref, framesize_in_bytes*2/3, SEEK_CUR);

  free (buf);

  img->quad[0]=0;
  diff_y=0;
  for (j=0; j < size_y; ++j)
  {
    for (i=0; i < size_x; ++i)
    {
#ifdef  INTERNAL_BIT_DEPTH_INCREASE
      diff_y += SQR_DEPTH(p->imgY[j][i], imgY_ref[j][i], input->output_bitdepth, img->BitDepthIncrease);
#else
      diff_y += img->quad[abs(p->imgY[j][i]-imgY_ref[j][i])];
#endif
    }
  }

  // Chroma
  diff_u=0;
  diff_v=0;

  if (p->chroma_format_idc != YUV400)
  {
    for (j=0; j < size_y_cr; ++j)
    {
      for (i=0; i < size_x_cr; ++i)
      {
#ifdef  INTERNAL_BIT_DEPTH_INCREASE
        diff_u += SQR_DEPTH(p->imgUV[0][j][i], imgUV_ref[0][j][i], input->output_bitdepth, img->BitDepthIncreaseChroma);
        diff_v += SQR_DEPTH(p->imgUV[1][j][i], imgUV_ref[1][j][i], input->output_bitdepth, img->BitDepthIncreaseChroma);
#else
        diff_u += img->quad[abs(imgUV_ref[0][j][i]-p->imgUV[0][j][i])];
        diff_v += img->quad[abs(imgUV_ref[1][j][i]-p->imgUV[1][j][i])];
#endif
      }
    }
  }

#if ZEROSNR
  if (diff_y == 0)
    diff_y = 1;
  if (diff_u == 0)
    diff_u = 1;
  if (diff_v == 0)
    diff_v = 1; 
#endif

  // Collecting SNR statistics
  if (diff_y != 0)
    snr->snr_y=(float)(10*log10(max_pix_value_sqd*(double)((double)(size_x)*(size_y) / diff_y)));        // luma snr for current frame
  else
    snr->snr_y=0.0;
  if (diff_u != 0)
    snr->snr_u=(float)(10*log10(max_pix_value_sqd_uv*(double)((double)(size_x_cr)*(size_y_cr) / (diff_u))));    //  chroma snr for current frame
  else
    snr->snr_u=0.0;
  if (diff_v != 0)
    snr->snr_v=(float)(10*log10(max_pix_value_sqd_uv*(double)((double)(size_x_cr)*(size_y_cr) / (diff_v))));    //  chroma snr for current frame
  else
    snr->snr_v=0;

  if (img->number == 0) // first
  {
    snr->snr_ya=snr->snr_y1=snr->snr_y;                                                        // keep luma snr for first frame
    snr->snr_ua=snr->snr_u1=snr->snr_u;                                                        // keep chroma snr for first frame
    snr->snr_va=snr->snr_v1=snr->snr_v;                                                        // keep chroma snr for first frame

  }
  else
  {
    snr->snr_ya=(float)(snr->snr_ya*(snr->frame_ctr)+snr->snr_y)/(snr->frame_ctr+1); // average snr chroma for all frames
    snr->snr_ua=(float)(snr->snr_ua*(snr->frame_ctr)+snr->snr_u)/(snr->frame_ctr+1); // average snr luma for all frames
    snr->snr_va=(float)(snr->snr_va*(snr->frame_ctr)+snr->snr_v)/(snr->frame_ctr+1); // average snr luma for all frames
  } 

  // picture error concealment
  if(p->concealed_pic)
  {
    fprintf(stdout,"%04d(P)  %8d %5d %5d %7.4f %7.4f %7.4f  %s %5d\n", 
      frame_no, p->frame_poc, p->pic_num, p->qp, 
      snr->snr_y, snr->snr_u, snr->snr_v, yuv_types[p->chroma_format_idc], 0);      

  }
}

#ifdef EIGHTH_PEL
void get_fullpel_block(int ref_frame, StorablePicture **list, int x_pres, int y_pres, int max_x, int max_y, int block[BLOCK_SIZE][BLOCK_SIZE])
{
  int i, j;

  for (j = 0; j < BLOCK_SIZE; j++)
    for (i = 0; i < BLOCK_SIZE; i++)
      block[i][j] = list[ref_frame]->imgY[max(0,min(y_pres+j,max_y))][max(0,min(x_pres+i,max_x))];
}

void interp_block_X(int ref_frame, StorablePicture **list, int pres_x, int pres_y, const int *coefx, int max_x, int max_y, int block[BLOCK_SIZE][BLOCK_SIZE])
{
  int i, j, x;
  int result;
  int rounding = 128;
  int shifting = 8;

  for( j = 0; j < BLOCK_SIZE; j++ )
    for( i = 0; i < BLOCK_SIZE; i++ )
    {
      for(result = 0, x = -3; x < 5; x++ )
        result += list[ref_frame]->imgY[Clip3(0, max_y, pres_y+j)][Clip3(0, max_x, pres_x+i+x)]*coefx[x+3];
      block[i][j] = Clip3(0, img->max_imgpel_value, ( result + rounding ) >> shifting );
    }
}

void interp_block_Y(int ref_frame, StorablePicture **list, int pres_x, int pres_y, const int *coefy, int max_x, int max_y, int block[BLOCK_SIZE][BLOCK_SIZE])
{
  int i, j, y;
  int result;
  int rounding = 128;
  int shifting = 8;

  for (j = 0; j < BLOCK_SIZE; j++)
    for (i = 0; i < BLOCK_SIZE; i++)
    {
      for(result = 0, y = -3; y < 5; y++)
        result += list[ref_frame]->imgY[Clip3(0, max_y, pres_y+j+y)][Clip3(0, max_x, pres_x+i)]*coefy[y+3];
      block[i][j] = Clip3(0, img->max_imgpel_value, ( result + rounding ) >> shifting );
    }
}

void interp_block_XY(int ref_frame, StorablePicture **list, int pres_x, int pres_y, const int *coefx, const int *coefy, int max_x, int max_y, int block[BLOCK_SIZE][BLOCK_SIZE])
{
  int i, j, n;
  int result;
  int tmp_block[BLOCK_SIZE+8][BLOCK_SIZE];

  int rounding = 32768;
  int shifting = 16;

  // horizontal interpolation
  for (j = -3; j < BLOCK_SIZE+5; j++)
    for (i = 0; i < BLOCK_SIZE; i++)
      for (tmp_block[j+3][i] = 0, n = -3; n < 5; n++)
        tmp_block[j+3][i] += list[ref_frame]->imgY[Clip3(0, max_y, pres_y+j)][Clip3(0,max_x, pres_x+i+n)]*coefx[n+3];

  // vertical interpolation
  for( j = 0; j < BLOCK_SIZE; j++ )
    for( i = 0; i < BLOCK_SIZE; i++ )
    {
      for( result = 0, n = -3; n < 5; n++ )
        result += tmp_block[j+n+3][i] * coefy[n+3];
      block[i][j] = Clip3(0, img->max_imgpel_value, ( result + rounding ) >> shifting );
    }
}

void average_block(int block[BLOCK_SIZE][BLOCK_SIZE], int block2[BLOCK_SIZE][BLOCK_SIZE]) 
{
  int i, j;

  for (j = 0; j < BLOCK_SIZE; j++)
    for (i = 0; i < BLOCK_SIZE; i++)
      block[i][j] = (block[i][j] + block2[i][j]) >> 1;
}

void average_block4(int block[BLOCK_SIZE][BLOCK_SIZE], int block2[BLOCK_SIZE][BLOCK_SIZE], int block3[BLOCK_SIZE][BLOCK_SIZE], int block4[BLOCK_SIZE][BLOCK_SIZE]) 
{
  int i, j;

  for (j = 0; j < BLOCK_SIZE; j++)
    for (i = 0; i < BLOCK_SIZE; i++)
      block[i][j] = (block[i][j] + block2[i][j]+ block3[i][j]+ block4[i][j])/4;
}

void average_block_2(int block[MB_BLOCK_SIZE][MB_BLOCK_SIZE], int block2[MB_BLOCK_SIZE][MB_BLOCK_SIZE]) 
{
  int i, j;

  for (j = 0; j < BLOCK_SIZE; j++)
    for (i = 0; i < BLOCK_SIZE; i++)
      block[i][j] = (block[i][j] + block2[i][j]) >> 1;
}

void average_block4_2(int block[MB_BLOCK_SIZE][MB_BLOCK_SIZE], int block2[MB_BLOCK_SIZE][MB_BLOCK_SIZE], int block3[MB_BLOCK_SIZE][MB_BLOCK_SIZE], int block4[MB_BLOCK_SIZE][MB_BLOCK_SIZE]) 
{
  int i, j;

  for (j = 0; j < BLOCK_SIZE; j++)
    for (i = 0; i < BLOCK_SIZE; i++)
      block[i][j] = (block[i][j] + block2[i][j]+ block3[i][j]+ block4[i][j])/4;
}

void get_eighthpel_block(int ref_frame, StorablePicture **list, int x_pos, int y_pos, struct img_par *img, int block[BLOCK_SIZE][BLOCK_SIZE])
{
  static const int COEF[3][8] =
  {
    { -3, 12, -37, 229,  71, -21,  6, -1 },
    { -3, 12, -39, 158, 158, -39, 12, -3 },
    { -1,  6, -21,  71, 229, -37, 12, -3 }
  };
  int rounding  = 32768;
  int rounding2 = 128;
  int shifting  = 16;
  int shifting2 = 8;

  int block2[BLOCK_SIZE][BLOCK_SIZE];
  int dx = x_pos & 7, x = 0;
  int dy = y_pos & 7, y = 0;
  int pres_x = x_pos >> 3;
  int pres_y = y_pos >> 3;

  int tmp[BLOCK_SIZE][BLOCK_SIZE+8];
  int result;
  int i, j;

  int max_x = dec_picture->size_x - 1;
  int max_y = dec_picture->size_y - 1;

  if (dec_picture->mb_field[img->current_mb_nr])
    max_y = dec_picture->size_y/2 - 1;

  // choose filter depending on the sub-pel position
  if( dx == 0 && dy == 0 )
  { // full-pel position, just return full-pel block
    get_fullpel_block(ref_frame, list, pres_x, pres_y, max_x, max_y, block);
  }
  else if( dy == 0 )
  { // horizontal interpolation only
    if (dx == 1) // get a block just left to searched position
      get_fullpel_block(ref_frame, list, pres_x, pres_y, max_x, max_y, block);
    else
      interp_block_X(ref_frame, list, pres_x, pres_y, COEF[(dx>>1)-1], max_x, max_y, block);

    if( dx&1 )
    { // 1/8-pel position, needs 2nd block right to searched position for bilinear interpolation
      if( dx == 7 )
        get_fullpel_block(ref_frame, list, pres_x+1, pres_y, max_x, max_y, block2);
      else
        interp_block_X(ref_frame, list, pres_x, pres_y, COEF[dx>>1], max_x, max_y, block2);

      average_block(block, block2);
    }
  }
  else if( dx == 0 )
  { // vertical interpolation only
    if (dy == 1) // get a block just above searched position
      get_fullpel_block(ref_frame, list, pres_x, pres_y, max_x, max_y, block);
    else
      interp_block_Y(ref_frame, list, pres_x, pres_y, COEF[(dy>>1)-1], max_x, max_y, block);

    if( dy&1 )
    { // 1/8-pel position, needs 2nd block below searched position for bilinear interpolation
      if (dy == 7)
        get_fullpel_block(ref_frame, list, pres_x, pres_y+1, max_x, max_y, block2);
      else
        interp_block_Y(ref_frame, list, pres_x, pres_y, COEF[(dy>>1)], max_x, max_y, block2);

      average_block(block, block2);
    }
  }
  else if( !(dx&1) )
  {  // horizontal 1/4-pel
    for( j = -3; j < BLOCK_SIZE + 5; j++ )
      for( i = 0; i < BLOCK_SIZE; i++ )
      {
        for( tmp[i][j+3] = 0, x = -3; x < 5; x++ )
          tmp[i][j+3] += list[ref_frame]->imgY[Clip3(0, max_y, pres_y+j )][Clip3(0, max_x, pres_x+i+x )]*COEF[(dx>>1)-1][x+3];
      }

      if( dy == 1 )
      {
        for( j = 0; j < BLOCK_SIZE; j++ )
          for( i = 0; i < BLOCK_SIZE; i++ )
            block[i][j] = Clip3(0, img->max_imgpel_value, ( tmp[i][j+3] + rounding2 ) >> shifting2 );
      }
      else
      {
        for( j = 0; j < BLOCK_SIZE; j++ )
          for( i = 0; i < BLOCK_SIZE; i++ )
          {
            for( result = 0, y = -3; y < 5; y++ )
              result += tmp[i][j+y+3] * COEF[(dy>>1)-1][y+3];
            block[i][j] = Clip3(0, img->max_imgpel_value, ( result + rounding ) >> shifting );
          }
      }

      if( dy&1 )
      { // 1/8-pel position, lower 1/4-pel position for averaging
        if( dy == 7 )
          for( j = 0; j < BLOCK_SIZE; j++ )
            for( i = 0; i < BLOCK_SIZE; i++ )
              block2[i][j] = Clip3(0, img->max_imgpel_value, ( tmp[i][j+4] + rounding2 ) >> shifting2 );
        else
        {
          for( j = 0; j < BLOCK_SIZE; j++ )
            for( i = 0; i < BLOCK_SIZE; i++ )
            {
              for(result = 0, y = -3; y < 5; y++)
                result += tmp[i][j+y+3] * COEF[dy>>1][y+3];
              block2[i][j] = Clip3(0, img->max_imgpel_value, ( result + rounding ) >> shifting );
            }
        }
        average_block(block, block2);
      }
  }
  else if( !(dy&1) )
  { // vertical 1/4-pel and horizontal 1/8-pel interpolation
    for( j = 0; j < BLOCK_SIZE; j++ )        // y
      for( i = -3; i < BLOCK_SIZE + 5; i++ ) // x
      { // get vertical 1/4-pel positions
        for( tmp[j][i+3] = 0, y = -3; y < 5; y++ )
          tmp[j][i+3] += list[ref_frame]->imgY[Clip3(0, max_y, pres_y+j+y)][Clip3(0, max_x, pres_x+i)]*COEF[(dy>>1)-1][y+3];

      }

      // left 1/4-pel position for averaging
      if( dx == 1 )