postprocessing.c

JM 11.0 KTA 2.1 Source Code

      }
    }//j
  }//i
}



void Filter_training(imgpel **imgY_filt, imgpel **imgY_rec, int fl, int **varImg,int QPindex)
{
  int i, j, ii, jj, iiLimit, jjLimit, k, l, ind,QPind=0;
  int vectorInd, yLocal, ELocal[MAX_SQR_FILT_LENGTH];
  int sqrFiltLength=(fl+1)*(fl+1);
  static int first_training=1;

  int *pPattern=NULL;
  if (fl ==4)
  {
    pPattern=patternFullPel9x9;
  }
  else  if (fl ==3)
  {
    pPattern=patternFullPel7x7;
  }
  else  if (fl ==2)
  {
    pPattern=patternFullPel5x5;
  }



  if (first_training==1)
  {
    if(img->height == 144)
    {
      FILTER_NUMBER = 4;
      FilterSetPerQP = 0;
    }
    else if(img->height == 288)
    {
      FILTER_NUMBER = 8;
    }
    else if(img->height == 720)
    {
      FILTER_NUMBER = 16;
    }
    else
    {
      printf("Currently only QCIF, CIF and 720P are supported!\n");
      exit(1);
    }

    for (QPind=0;QPind<FILTERS_SET_NUMBER;QPind++)
    {

      memset(EGlobal[QPind],0, MAX_FILTER_NUMBER*sqrFiltLength*sqrFiltLength*sizeof(double));
      memset(yGlobal[QPind],0, MAX_FILTER_NUMBER*sqrFiltLength*sizeof(double));

      QP_set[QPind]=0;

    }
  }

  first_training=0;
  for (i = 0; i < input->img_height; ++i){
    for (j = 0; j < input->img_width; ++j){
      k=0;  memset(ELocal, 0, sqrFiltLength*sizeof(int));
      for (ii=i-fl; ii<=i+fl; ii++){
        iiLimit=max(0, min(input->img_height-1,ii));
        for (jj=j-fl; jj<=j+fl; jj++){  
          jjLimit=max(0, min(input->img_width-1,jj));

          vectorInd=pPattern[k++];
          ELocal[vectorInd]+=imgY_rec[iiLimit][jjLimit];
        }
      }
      yLocal=imgY_orgPP[i][j];

      ind=min(varImg[i][j]*FILTER_NUMBER/16, FILTER_NUMBER-1);			

      for (k=0; k<sqrFiltLength; k++)
      {
        for (l=0; l<sqrFiltLength; l++)
        {

          EGlobal[QPindex][ind][k][l]+=(double)(ELocal[k]*ELocal[l]);
        }
        yGlobal[QPindex][ind][k]+=(double)(ELocal[k]*yLocal);
      }
    }
  }
  QP_set[QPindex]=1;
}



int calculate_filtersPP(int fl)
{
  int ind,QPind=0;
  int sqrFiltLength=(fl+1)*(fl+1);
  for (QPind=0;QPind<FILTERS_SET_NUMBER;QPind++)
  {
    if (QP_set[QPind]==1)
    {
      for (ind=0; ind<FILTER_NUMBER; ind++)
      {
        gnsSolveByChol_Post(EGlobal[QPind][ind], yGlobal[QPind][ind], filterCoeffPP[QPind][ind], sqrFiltLength);
      }
    }
  }
  return 1;
}


double find_square_error(imgpel **imgY_filt, imgpel **imgY_orgPP)
{
  int i, j;
  double	seFilt=0,temp;
  for (i = 0; i < input->img_height; ++i)
  {
    for (j = 0; j < input->img_width; ++j)
    {
      temp=(double)(imgY_filt[i][j])-(double)(imgY_orgPP[i][j]);
      seFilt+=temp*temp;
    }
  }
  return seFilt;
}


/*
* gnsCholeskyDec_Post:
*
* Parameters:
*     inpMatr          I: 6x6 symmetric input matrix
*     outMatr          O: 6x6 upper triangular output matrix
*
* Function:
*     Performs Cholesky decomposition of inpMatr, i.e. calculates
*     upper triangular outMatr such that outMatr'*outMatr = inpMatr
* Returns:
*     1 if inpMatr is positive definite, i.e. factorization has been done
*     0 if inpMatr is not positive definite, i.e. factorization has failed
*
*/

int gnsCholeskyDec_Post(double inpMatr[MAX_SQR_FILT_LENGTH][MAX_SQR_FILT_LENGTH] , double outMatr[MAX_SQR_FILT_LENGTH][MAX_SQR_FILT_LENGTH], int noEq)
{ 
  int 
    i, j, k;     /* Looping Variables */
  double 
    scale;       /* scaling factor for each row */
  double 
    invDiag[MAX_SQR_FILT_LENGTH];  /* Vector of the inverse of diagonal entries of outMatr */


  /*
  *  Cholesky decomposition starts
  */

  for(i = 0; i < noEq; i++)
    for(j = i; j < noEq; j++)
    {
      /* Compute the scaling factor */
      scale=inpMatr[i][j];
      if ( i > 0) for( k = i - 1 ; k >= 0 ; k--)
        scale -= outMatr[k][j] * outMatr[k][i];

      /* Compute i'th row of outMatr */
      if(i==j){
        if(scale <= 0 ){  /* If inpMatr is singular */
          return(0);
        }
        else              /* Normal operation */
          invDiag[i] =  1.0/(outMatr[i][i]=sqrt(scale));
      }
      else{
        outMatr[i][j] = scale*invDiag[i]; /* Upper triangular part          */
        outMatr[j][i] = 0.0;              /* Lower triangular part set to 0 */
      }                    
    }

    return(1); /* Signal that Cholesky factorization is successfully performed */
}

/*
* gnsTransposeBacksubstitution_post:
*
* Parameters:
*     U          I: a 6x6 upper triangular matrix
*     rhs        I: a vector of size 6
*     x          O: vector that satisfies U'*x = rhs
*
* Function:
*     Solves U'*x = rhs for x
*
* Returns:
*     Nothing
*
* History
*     8.4.97         Initial version.
*/

void gnsTransposeBacksubstitution_post(double U[MAX_SQR_FILT_LENGTH][MAX_SQR_FILT_LENGTH], double rhs[], 
                                       double x[], int order)
{
  int 
    i,j;              /* Looping variables */
  double 
    sum;              /* Holds backsubstitution from already handled rows */

  /* Backsubstitution starts */
  x[0] = rhs[0]/U[0][0];               /* First row of U'                   */
  for (i = 1; i < order; i++){         /* For the rows 1..order-1           */

    for (j = 0, sum = 0.0; j < i; j++) /* Backsubst already solved unknowns */
      sum += x[j]*U[j][i];

    x[i]=(rhs[i] - sum)/U[i][i];       /* i'th component of solution vect.  */
  }
}


/*
* gnsTransposeBacksubstitution_post:
*
* Parameters:
*     R          I: a 6x6 upper triangular matrix
*     z          I: a vector of size 6
*     A          O: vector that satisfies R*A = z
*
* Function:
*     Solves by backsubstitution the system given by the matrix
*     equation: R * A = z, with R being upper-triangular.
*
* Returns:
*     Nothing
*
*/

void gnsBacksubstitution_Post(double R[MAX_SQR_FILT_LENGTH][MAX_SQR_FILT_LENGTH], double z[MAX_SQR_FILT_LENGTH], int R_size, double A[MAX_SQR_FILT_LENGTH])
{
  int
    i, j;

  double
    sum;

  R_size--;

  A[R_size] = z[R_size] / R[R_size][R_size];

  for (i = R_size-1; i >= 0; i--) {

    for (j = i+1, sum = 0.0; j <= R_size; j++)
      sum += R[i][j] * A[j];

    A[i] = (z[i] - sum) / R[i][i];

  }
}

/*
* gnsSolveByChol_Post:
*
* Parameters:
*     LHS              I: Left Hand Side Matrix(6x6) of least squares eqn.
*     rhs              I: Right Hand Side vector(6x1) of least squares eqn.
*     x                O: vector that satisfies LHS*x = rhs
*
* Function:
*     Solves LHS*x = rhs for x. Uses Cholesky decomposition + two
*     backsubstitutions for the solution. In case LHS is singular, 
*     regularizes LHS by adding a small positive value to the
*     diagonals of LHS
*
* Returns:
*     Nothing
*
*/

void gnsSolveByChol_Post(double LHS[][MAX_SQR_FILT_LENGTH], double rhs[], double x[], int noEq)
{ 
  double
    aux[MAX_SQR_FILT_LENGTH],     /* Auxiliary vector */
    U[MAX_SQR_FILT_LENGTH][MAX_SQR_FILT_LENGTH];    /* Upper triangular Cholesky factor of LHS */
  int
    i;          /* Looping variable */

  /* The equation to be solved is LHSx = rhs */

  /* Compute upper triangular U such that U'*U = LHS */
  if(gnsCholeskyDec_Post(LHS, U, noEq)){ /* If Cholesky decomposition has been successful */

    /* Now, the equation is  U'*U*x = rhs, where U is upper triangular
    * Solve U'*aux = rhs for aux
    */
    gnsTransposeBacksubstitution_post(U, rhs, aux, noEq);         

    /* The equation is now U*x = aux, solve it for x (new motion coefficients) */
    gnsBacksubstitution_Post(U, aux, noEq, x);   

  }
  else{ /* LHS was singular */ 

    /* Regularize LHS */
    for(i=0; i<noEq; i++)
      LHS[i][i] += REG;

    /* Compute upper triangular U such that U'*U = regularized LHS */
    gnsCholeskyDec_Post(LHS, U, noEq);

    /* Solve  U'*aux = rhs for aux */  
    gnsTransposeBacksubstitution_post(U, rhs, aux, noEq);   

    /* Solve U*x = aux for x */
    gnsBacksubstitution_Post(U, aux, noEq, x);
  }  
}



void generateSqrtTable()
{
  int i,j;
  sqrtTable[0]=0;
  for (i=1;i<16;++i)
  {
    int start=i*i;
    for (j=start;j<start+i+1;++j)
    {
      sqrtTable[j]=i;
    }
    for (j=start+i+1;j<(i+1)*(i+1);++j)
    {
      sqrtTable[j]=i+1;
    }
  }
}



//Global functions

int postProcessing(Picture *frame_pic)
{
  static imgpel   **imgY_filt;
  static int **imgY_var;
  imgpel   **refY;
  int memory_size = 0;


  int apply_filter=0;
  static int first=0;
  int img_width, img_height;  

  static int frameNo=0;
  //static int FrameBCountinGOP=1;
  int DisplayOrder_frameNo=enc_picture->frameNumberinFile;
  //Bitstream *tempStream=NULL;
  currentQP = frame_pic->slices[0]->qp;

  img_width  = input->img_width;
  img_height = input->img_height;
  refY       = enc_picture->imgY;

  imgY_orgPP   = imgY_org_frm;

  if (first==0)
  {
    memory_size += get_mem2Dpel(&imgY_filt, img->height, img->width);

    QP_map=calloc((input->no_frames-1)*(input->successive_Bframe + 1) + 1, sizeof(int));
    SeqNo_map=calloc((input->no_frames-1)*(input->successive_Bframe + 1) + 1, sizeof(int));

    memory_size += get_mem2Dint(&imgY_var, img->height, img->width);
    generateSqrtTable();
    first=1;
  }

  calcVar(imgY_var, refY, FILTER_LENGTH>>1);

  if (img->type !=B_SLICE || FilterSetPerQP==0)
  {
    QP_index=0;
  }
  else if (img->nal_reference_idc==0)
  {
    QP_index=1;
  }
  else
  {
    QP_index=2;
  }

  apply_filter=0;
  Filter_training(imgY_filt,refY,FILTER_LENGTH>>1,imgY_var,QP_index);
  QP_map[frameNo] =QP_index;
  SeqNo_map[DisplayOrder_frameNo] = frameNo;


  frameNo++;


  return (apply_filter);
}





double Apply_postprocessing(imgpel **refY,imgpel **imgY_orgPP,Bitstream *bitstream)
{ 
  int j;
  double seFilt,seOrig;
  static imgpel   **imgY_filt;
  static int first1=0;
  static int **imgY_var;
  int memory_size = 0;
  int apply_filter=0;
  static int frameNo1=0;
  filterLengthPP = FILTER_LENGTH>>1;
  if (first1==0)
  {
    memory_size += get_mem2Dpel(&imgY_filt, input->img_height, input->img_width);
    memory_size += get_mem2Dint(&imgY_var, input->img_height, input->img_width);
    Apply_filter_map=calloc((input->no_frames-1)*(input->successive_Bframe + 1) + 1, sizeof(int));
    first1=1;
  }


  calcVar(imgY_var, refY, filterLengthPP);
  filter_frame(imgY_filt,refY,filterLengthPP,imgY_var,QP_map[SeqNo_map[frameNo1]]);
  seOrig =find_square_error(refY,imgY_orgPP);
  seFilt =find_square_error(imgY_filt,imgY_orgPP);


  if (seOrig>seFilt)
  {
    apply_filter=1;
    for (j = 0; j < input->img_height; j++)
      memcpy(refY[j],imgY_filt[j],sizeof(imgpel)*input->img_width);
  }


  Apply_filter_map[SeqNo_map[frameNo1]] = apply_filter;


  frameNo1++;

  if (apply_filter)
    return 10*log10(255.0*255.0*input->img_width*input->img_height/seFilt);
  else
    return 10*log10(255.0*255.0*input->img_width*input->img_height/seOrig);
}




void finish_bitstream(Bitstream *bitstream)
{

  int i;
  for (i=0;i<((input->no_frames-1)*(input->successive_Bframe + 1) + 1);i++)
  {
    u_1("apply post filter",Apply_filter_map[i],bitstream);
  }

  // try to put 0 at the end of slice
  if (bitstream->bits_to_go != 8 )
  {
    while(bitstream->bits_to_go != 8)
    {
      u_1("QL: paddle", 1, bitstream); 
    }
  }
  if (QP_map)
    free(QP_map);
  if (SeqNo_map)
    free(SeqNo_map);
  if (Apply_filter_map)
    free(Apply_filter_map);
}

int writePostfilterSEIUnit(Bitstream* currStream)
{

  NALU_t *n ;
  int outbits, NALUlen;
  byte rbsp[MAXRBSPSIZE];
  clear_sei_message(NORMAL_SEI);
  write_sei_message(NORMAL_SEI, currStream->streamBuffer, currStream->byte_pos, SEI_POST_FILTER_HINTS);
  finalize_sei_message(NORMAL_SEI);
  memcpy(rbsp,sei_message[NORMAL_SEI].data,sei_message[NORMAL_SEI].payloadSize);


  n= AllocNALU(MAXRBSPSIZE);

  NALUlen = RBSPtoNALU (rbsp, n, sei_message[NORMAL_SEI].payloadSize, NALU_TYPE_SEI, NALU_PRIORITY_HIGHEST, 0, 0);
  n->startcodeprefix_len = 3;
  n->forbidden_bit = 0;

  start_post_sequence();
  outbits = WriteNALU (n);
  stats->bit_ctr += outbits;
  terminate_post_sequence();
  FreeNALU(n);

  return (outbits);
}


void direct_output_postfilter(StorablePicture *p, int p_out)
{
  // we have a frame (or complementary field pair)
  // so output it directly
  write_picture (p, p_out, FRAME);
}

#endif