endecode.cpp

设计并实现了两种分层多描述视频编码器.通过对小波域的运动估计算法进行了分析和研究

	dedata_file.close_file();
	image.coe(1./factor1);
    //运动补偿
	image.Desend_coeff(dimension, DePcoeff);
    Deall_phase_motion_Estimation(dimension, DePcoeff2);
	
	for(l=0;l<144;l++)
		for(m=0;m<176;m++)
			DePcoeff[l][m]+=DePcoeff2[l][m];

	image.send_coeff(dimension, DePcoeff);

//    if(mdc==1) image.send_coeff1(dimension, DePcoeff);
/**/    
    


 
	if (!i) level1 = image.pyramid_levels();
    image.recover(mean);	




    total_time.stop();

    printf("\n  Rate =%6.3f bits/pixel --> ", rate_mult * dedata_file.bytes_used());
	  
	avbit+=rate_mult * dedata_file.bytes_used();
     
	mse = image.compare(num,pic_f_name);
    if (mse > 0.0) {
        psnr = dBW(peak_2 / mse);
        printf("mean squared-error =%8.3f =%7.2f dB PSNR\n", mse, psnr); 
	}
	if(i==0)
		avpsnr1+=psnr;
	else
		avpsnr2+=psnr;
    total_time.display("\n  Total execution time (I/O included) =");


    if(num==1)
	{
		for(l=0;l<144;l++)
				for(m=0;m<176;m++)
				{
					Wcoeff[l][m]=floor(image.readdress(l,m)+0.5);
					if(Wcoeff[l][m]<0) Wcoeff[l][m]=0;
					if(Wcoeff[l][m]>255) Wcoeff[l][m]=255;
				}
		if(i==0)
		{
			
			for(l=0;l<144;l++)
				fwrite(Wcoeff[l],sizeof(unsigned char),176,output1);
		}
		else
		{
			for(l=0;l<144;l++)
				fwrite(Wcoeff[l],sizeof(unsigned char),176,output2);
		}
	}/**/
  }
//  image.recall(level1);
  image.lvset(level1);
//  image.coeprintf();


  image.Desend_coeff1(dimension, DePcoeff);
  image.coeP(1./factor1,DePcoeff); 

  	for(l=0;l<144;l++)
		for(m=0;m<176;m++)
			DePcoeff[l][m]+=DePcoeff2[l][m];
  	image.send_coeff(dimension, DePcoeff);/**/
//  image.coecmp();

  image.recover(mean,recoeff);
    if(num==1)
	{
		for(l=0;l<144;l++)
				for(m=0;m<176;m++)
				{
					Wcoeff[l][m]=floor(image.readdress(l,m)+0.5);
					if(Wcoeff[l][m]<0) Wcoeff[l][m]=0;
					if(Wcoeff[l][m]>255) Wcoeff[l][m]=255;
				}		
		for(l=0;l<144;l++)
			fwrite(Wcoeff[l],sizeof(unsigned char),176,output0);

	}
  
  mse = image.compare(num,pic_f_name);
  if (mse > 0.0) {
	psnr = dBW(peak_2 / mse);
	avpsnr0+=psnr;
	printf("mean squared-error =%8.3f =%7.2f dB PSNR\n", mse, psnr); }


}

//----------------------------------------------------------------------

void transform1(Pel_Type** recoeff)
{//小波变换
  int i = NumbTap, j = 0, i1, j1, i2, j2;
  double temp_line[176+2*NumbTap];
  double * t, * in_line;
  int nx, ny, mx = 144, my = 176;
//  FILE *fileout;

  in_line= temp_line + i;
//第1次分解，根据小波系数下采样取奇偶、偶奇样点分解为4个图像
  nx = mx;  mx >>= 1;  ny = my;  my >>= 1;   
  // transformation of columns
  for (j = 0; j < ny; j++) {
      for (i = 0; i < nx; i++) in_line[i] = recoeff[i][j];
      Reflection(in_line, in_line + nx - 1);

      for (i = 0, t = in_line; i < mx; i++) {
        w1[0][0][i][j] = Filter_L(T_LowPass, t++); //低频子带
        w1[1][0][i][j] = Filter_L(T_LowPass, t); //下移一个象素的低频子带
        w1[0][0][i+mx][j] = Filter_H(T_HighPass, t++); //高频子带
        w1[1][0][i+mx][j] = Filter_H(T_HighPass, t); //下移一个象素的高频子带 
	  } 
  }
  // transformation of rows
  for (i = 0; i < nx; i++) {
      for (j = 0; j < ny; j++) in_line[j] = w1[0][0][i][j];
      Reflection(in_line, in_line + ny - 1);
      for (j = 0, t = in_line; j < my; j++) {
        w1[0][0][i][j] = Filter_L(T_LowPass, t++); //低频子带
        w1[0][1][i][j] = Filter_L(T_LowPass, t); //右移一个象素的低频子带
        w1[0][0][i][j+my] = Filter_H(T_HighPass, t++); //高频子带
        w1[0][1][i][j+my] = Filter_H(T_HighPass, t); //右移一个象素的高频子带 
	  } 
      for (j = 0; j < ny; j++) in_line[j] = w1[1][0][i][j];
      Reflection(in_line, in_line + ny - 1);
      for (j = 0, t = in_line; j < my; j++) {
        w1[1][0][i][j] = Filter_L(T_LowPass, t++); //下移一个象素的低频子带
        w1[1][1][i][j] = Filter_L(T_LowPass, t); //对角线移位低频子带
        w1[1][0][i][j+my] = Filter_H(T_HighPass, t++); //下移一个象素的高频子带
        w1[1][1][i][j+my] = Filter_H(T_HighPass, t); //对角线移位高频子带 
	  } 
  }

  for(i=0, i2=0; i<mx; i++, i2++) {
	  for (j=my, j2=0; j<ny; j++) {
		  w[0][0][i2][j2++]=w1[0][0][i][j];
		  w[0][0][i2][j2++]=w1[0][1][i][j];
	  }
	  ++i2; j2=0;
	  for (j=my; j<ny; j++) {
		  w[0][0][i2][j2++]=w1[1][0][i][j];
		  w[0][0][i2][j2++]=w1[1][1][i][j];
	  }
  }
  for(i=mx, i2=0; i<nx; i++, i2++) {
	  for (j=0, j2=0; j<my; j++) {
		  w[0][1][i2][j2++]=w1[0][0][i][j];
		  w[0][1][i2][j2++]=w1[0][1][i][j];
	  }
	  ++i2; j2=0;
	  for (j=0;j<my;j++) {
		  w[0][1][i2][j2++]=w1[1][0][i][j];
		  w[0][1][i2][j2++]=w1[1][1][i][j];
	  }
  }
  for(i=mx, i2=0; i<nx; i++, i2++) {
	  for (j=my, j2=0; j<ny; j++) {
		  w[0][2][i2][j2++]=w1[0][0][i][j];
		  w[0][2][i2][j2++]=w1[0][1][i][j];
	  }
	  ++i2; j2=0;
	  for (j=my;j<ny;j++) {
		  w[0][2][i2][j2++]=w1[1][0][i][j];
		  w[0][2][i2][j2++]=w1[1][1][i][j];
	  }
  }


//第2次分解
  nx = mx;  mx >>= 1;  ny = my;  my >>= 1;
  for (i1 = 0; i1 < 2; i1++) { 
	 for (j1 = 0; j1 < 2; j1++) {
  // transformation of columns
		for (j = 0; j < ny; j++) {
			for (i = 0; i < nx; i++) in_line[i] = w1[i1][j1][i][j];
			Reflection(in_line, in_line + nx - 1);

			for (i = 0, t = in_line; i < mx; i++) {
				w2[i1][j1][i][j] = Filter_L(T_LowPass, t++); //低频子带
				w2[i1+2][j1][i][j] = Filter_L(T_LowPass, t); //下移2个象素的低频子带
				w2[i1][j1][i+mx][j] = Filter_H(T_HighPass, t++); //高频子带
				w2[i1+2][j1][i+mx][j] = Filter_H(T_HighPass, t); //下移2个象素的高频子带 
			}
		} 
  // transformation of rows
		for (i = 0; i < nx; i++) {
			for (j = 0; j < ny; j++) in_line[j] = w2[i1][j1][i][j];
			Reflection(in_line, in_line + ny - 1);
			for (j = 0, t = in_line; j < my; j++) {
				w2[i1][j1][i][j] = Filter_L(T_LowPass, t++); //低频子带
				w2[i1][j1+2][i][j] = Filter_L(T_LowPass, t); //右移2个象素的低频子带
				w2[i1][j1][i][j+my] = Filter_H(T_HighPass, t++); //高频子带
				w2[i1][j1+2][i][j+my] = Filter_H(T_HighPass, t); //右移2个象素的高频子带 
			} 
			for (j = 0; j < ny; j++) in_line[j] = w2[i1+2][j1][i][j];
			Reflection(in_line, in_line + ny - 1);
			for (j = 0, t = in_line; j < my; j++) {
				w2[i1+2][j1][i][j] = Filter_L(T_LowPass, t++); //低频子带
				w2[i1+2][j1+2][i][j] = Filter_L(T_LowPass, t); //对角线移位2个象素的低频子带
				w2[i1+2][j1][i][j+my] = Filter_H(T_HighPass, t++); //高频子带
				w2[i1+2][j1+2][i][j+my] = Filter_H(T_HighPass, t); //对角线移位2个象素的高频子带 
			} 
		}
	} 
  }

  for(i=0,i2=0;i<mx;i++) 
	  for (i1=0;i1<4;i1++,i2++) 
		for (j=my,j2=0;j<ny;j++) 
		  for (j1=0; j1<4; j1++) w[1][0][i2][j2++]=w2[i1][j1][i][j];
  
  for(i=mx,i2=0;i<nx;i++)
	  for (i1=0; i1<4 ;i1++,i2++) 
		for (j=0,j2=0;j<my;j++) 
		  for (j1=0; j1<4; j1++) w[1][1][i2][j2++]=w2[i1][j1][i][j];
  
  for(i=mx,i2=0;i<nx;i++) 
	  for (i1=0; i1<4 ;i1++,i2++) 
		for (j=my,j2=0;j<ny;j++) 
		  for (j1=0; j1<4; j1++) w[1][2][i2][j2++]=w2[i1][j1][i][j];



//第3次分解
  nx = mx;  mx >>= 1;  ny = my;  my >>= 1;    
  // transformation of columns
  for (i1 = 0; i1 < 4; i1++) { 
	 for (j1 = 0; j1 < 4; j1++) {
  // transformation of columns
		for (j = 0; j < ny; j++) {
			for (i = 0; i < nx; i++) in_line[i] = w2[i1][j1][i][j];
			Reflection(in_line, in_line + nx - 1);

			for (i = 0, t = in_line; i < mx; i++) {
				w3[i1][j1][i][j] = Filter_L(T_LowPass, t++); //低频子带
				w3[i1+4][j1][i][j] = Filter_L(T_LowPass, t); //下移4个象素的低频子带
				w3[i1][j1][i+mx][j] = Filter_H(T_HighPass, t++); //高频子带
				w3[i1+4][j1][i+mx][j] = Filter_H(T_HighPass, t); //下移4个象素的高频子带 
			} 
		}
 // transformation of rows
		for (i = 0; i < nx; i++) {
			for (j = 0; j < ny; j++) in_line[j] = w3[i1][j1][i][j];
			Reflection(in_line, in_line + ny - 1);
			for (j = 0, t = in_line; j < my; j++) {
				w3[i1][j1][i][j] = Filter_L(T_LowPass, t++); //低频子带
				w3[i1][j1+4][i][j] = Filter_L(T_LowPass, t); //右移4个象素的低频子带
				w3[i1][j1][i][j+my] = Filter_H(T_HighPass, t++); //高频子带
				w3[i1][j1+4][i][j+my] = Filter_H(T_HighPass, t); //右移4个象素的高频子带 
			} 
			for (j = 0; j < ny; j++) in_line[j] = w3[i1+4][j1][i][j];
			Reflection(in_line, in_line + ny - 1);
			for (j = 0, t = in_line; j < my; j++) {
				w3[i1+4][j1][i][j] = Filter_L(T_LowPass, t++); //低频子带
				w3[i1+4][j1+4][i][j] = Filter_L(T_LowPass, t); //对角线移位4个象素的低频子带
				w3[i1+4][j1][i][j+my] = Filter_H(T_HighPass, t++); //高频子带
				w3[i1+4][j1+4][i][j+my] = Filter_H(T_HighPass, t); //对角线移位4个象素的高频子带 
			} 
		}
	} 
  } 
  for(i=0,i2=0;i<mx;i++) 
	  for (i1=0; i1<8; i1++, i2++) 
		for (j=0,j2=0;j<my;j++) 
		  for (j1=0; j1<8; j1++) w1[0][0][i2][j2++]=w3[i1][j1][i][j];

  for(i=0,i2=0;i<mx;i++) 
	  for (i1=0; i1<8; i1++, i2++) 
		for (j=my,j2=0;j<ny;j++) 
		  for (j1=0; j1<8; j1++) w1[0][1][i2][j2++]=w3[i1][j1][i][j];
  
  for(i=mx,i2=0;i<nx;i++) 
	  for (i1=0; i1<8; i1++, i2++) 
		for (j=0,j2=0;j<my;j++) 
		  for (j1=0; j1<8; j1++) w1[1][0][i2][j2++]=w3[i1][j1][i][j];
  
  for(i=mx,i2=0;i<nx;i++) 
	  for (i1=0; i1<8; i1++, i2++) 
		for (j=my,j2=0;j<ny;j++) 
		  for (j1=0; j1<8; j1++) w1[1][1][i2][j2++]=w3[i1][j1][i][j];
 
}

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

int main(int numb_arg, char * arg[])
{
  int i;
  min_bits=2;
  factor1=0.6962;//1;//1.078;//0.88;//.967;//0.6962;//0.504;//0.437;//1;//1.9;//有final_bit的倍乘因子
//  factor1=2.006;//1.024;//0.53;//去掉final_bit的倍乘因子
  pic_f_name="input.y";//"lena.512"
  dimension.x=144;//512
  dimension.y=176;//512
  pel_bytes=1;
  bit_rate=0.8;//2.25;//0.5;//

  
  output0=fopen("reconstruct0","wb");
  output1=fopen("reconstruct1","wb");
  output2=fopen("reconstruct2","wb");

  Pel_Type ** recoeff;
  NEW_VECTOR(recoeff, 144, Pel_Type *, "error");//512
  for (i = 0; i < 144; i++) NEW_VECTOR(recoeff[i], 176, Pel_Type, "error");//512


  rate_mult = 8.0 / (double(dimension.x) * double(dimension.y));
  if (Max(dimension.x, dimension.y) > 8192) Error("image is too large!");
  byte_budget =65535;//long(0.125 * (bit_rate * dimension.x) *dimension.y) - 1;//


  for (i=0; i<2; i++) {
	if (i>0) {
		cod_f_name="input2.cmp"; 
		mdc=2;
	}
	else {
		cod_f_name="input1.cmp"; 
		mdc=1;
	}
    total_time.start();
    image.read_pic(dimension, pic_f_name, pel_bytes);
    image.transform();
    mean = image.transform_mean();
    image.coe(factor1); 
//	if (mdc==1) image.coewrite(); //调试用，存系数中心重建时读出比较

//	if (i) image.decimate2(1<<(min_bits+3));
//	image.coeprintf();
    Chronometer max_time;
    max_time.start("\n  Starting computation of maximum magnitudes...");
    threshold_bits = Calc_Max_Image();
    max_time.display(" Maximum magnitudes computed in");
  
    data_file.open_file(cod_f_name);
    data_file.code_bits(11, (dimension.x >> 2) - 1);
    data_file.code_bits(11, (dimension.y >> 2) - 1);
    data_file.code_bits(5, threshold_bits);
    data_file.code_bits(1, pel_bytes - 1);
    data_file.code_bits(8, int(mean & 0xFFL));
    data_file.code_bits(8, int((mean >> 8) & 0xFFL));
    if (pel_bytes == 2)
    data_file.code_bits(8, int((mean >> 16) & 0xFFL));

    code_time.start("\n  Starting image compression...");
    Process_Image();
    code_time.display(" Image compressed in");

    data_file.close_file();

    printf("\n  Compressed file size = %ld bytes (%6.3f bits/pixel).\n",
      data_file.bytes_used(), rate_mult * data_file.bytes_used());

    total_time.display("\n  Total execution time (I/O included) =");
    puts(" ");
    total_time.stop();    
  }
//---------------------------reconstruct-----------------------------------------
  reconstruct(output0,recoeff);
  transform1(recoeff);
  
  for(int j=1;j<20;j++)
  {
	 printf("j=%d\n",j);
	 P_Frame(j);
     DeP_Frame(j,recoeff);
	 transform1(recoeff);
  }
  fclose(output2);
  fclose(output1);
  fclose(output0);
  
  avpsnr0=avpsnr0/(j-1);
  avpsnr1=avpsnr1/(j-1);
  avpsnr2=avpsnr2/(j-1);
  avbit=avbit/(j-1);
  printf("**************************对中心描述进行运动估计************************\n");
  printf("P帧个数＝%d\n",j-1);
  printf("第一个描述的平均psnr=%7.2f\n",avpsnr1);
  printf("第二个描述的平均psnr=%7.2f\n",avpsnr2);
  printf("平均bit=%8.3f\n",avbit);
  printf("中心描述的平均psnr=%7.2f\n",avpsnr0);


  return 0;
}


// end of file < CodeTree.C >