intdwtdoc.cpp

图像小波变换(包括提升小波变换)及SPIHT压缩解压缩程序

      L3=1<<(level-1);
	  float var7,var8,var9;
   	  var7=0;var8=0;var9=0;
      sum1=0;sum2=0;sum3=0;mean1=0;mean2=0;mean3=0;
	
	  for(i=0;i<(height/L3);i++)
	     for(j=(width/L3);j<2*(width/L3);j++)
		 {
	    	sum1=sum1+image[i*width+j];
		    sum2=sum2+image[(i+iheight/L3)*width+j-width/L3];
		    sum3=sum3+image[(i+height/L3)*width+j];
		 }
	  
	    mean1=sum1/((height/L3)*(width/L3));
	    mean2=sum2/((height/L3)*(width/L3));
	    mean3=sum3/((height/L3)*(width/L3));
	
	  for(i=0;i<(height/L3);i++)
	    for(j=(width/L3);j<2*(width/L3);j++)
		{
		  var7=var7+(image[i*width+j]-mean1)*(image[i*width+j]-mean1);
		  var8=var8+(image[(i+iheight/L3)*iwidth+j-iwidth/L3]-mean2)*(image[(i+iheight/L3)*iwidth+j-iwidth/L3]-mean2);
		  var9=var9+(image[(i+iheight/L3)*iwidth+j]-mean3)*(image[(i+iheight/L3)*iwidth+j]-mean3);
		}

    	var7=var7/((height/L3)*(width/L3));
    	var8=var8/((height/L3)*(width/L3));
    	var9=var9/((height/L3)*(width/L3));
	

	 //计算分辨率为0的各子带的均值和方差
	  L4=1<<level;
	  float var10,var11,var12;
   	  var10=0;var11=0;var12=0;
      sum1=0;sum2=0;sum3=0;mean1=0;mean2=0;mean3=0;
	
	  for(i=0;i<(height/L4);i++)
	     for(j=(width/L4);j<2*(width/L4);j++)
		 {
	    	sum1=sum1+image[i*width+j];
		    sum2=sum2+image[(i+iheight/L4)*width+j-width/L4];
		    sum3=sum3+image[(i+height/L4)*width+j];
		 }
	  
	    mean1=sum1/((height/L4)*(width/L4));
	    mean2=sum2/((height/L4)*(width/L4));
	    mean3=sum3/((height/L4)*(width/L4));
	
	  for(i=0;i<(height/L4);i++)
	    for(j=(width/L4);j<2*(width/L4);j++)
		{
		  var10=var10+(image[i*width+j]-mean1)*(image[i*width+j]-mean1);
		  var11=var11+(image[(i+iheight/L4)*iwidth+j-iwidth/L4]-mean2)*(image[(i+iheight/L4)*iwidth+j-iwidth/L4]-mean2);
		  var12=var12+(image[(i+iheight/L4)*iwidth+j]-mean3)*(image[(i+iheight/L4)*iwidth+j]-mean3);
		}

    	var10=var10/((height/L4)*(width/L4));
    	var11=var11/((height/L4)*(width/L4));
    	var12=var12/((height/L4)*(width/L4));

        //估算噪声的方差
	    subHH1=new float[(height/L1)*(width/L1)]; 
	    for(i=0;i<height/L1;i++)
           for(j=0;j<width/L1;j++)
		   {
		     subHH1[i*(width/L1)+j]=fabs(image[(i+height/L1)*width+j+width/L1]);
             sum=sum+subHH1[i*(width/L1)+j];
		   }
				  
	       
	      sum=sum/((height/L1)*(width/L1));
 
  
   
          var=(sum/0.8)*(sum/0.8);
    
          //计算参数
          float b1,b2,b3,b4;
	      b1=sqrt(log((width/L1)/L1));
	      b2=sqrt(log((width/L2)/L2));
	      b3=sqrt(log((width/L3)/L3));
          b4=sqrt(log((width/L4)/L4));

  	
          //计算各子带阈值
          float T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12;
		  
		  //分辨率为3的三个子带阈值
          T1=(b1*var)/sqrt(var1);
	      T2=(b1*var)/sqrt(var2);
	      T3=(b1*var)/sqrt(var3);
          
		  //分辨率为2的三个子带阈值
	      T4=(b2*var)/sqrt(var4);
          T5=(b2*var)/sqrt(var5);
	      T6=(b2*var)/sqrt(var6);

          //分辨率为1的三个子带阈值
	      T7=(b3*var)/sqrt(var7);
	      T8=(b3*var)/sqrt(var8);
	      T9=(b3*var)/sqrt(var9);
		  
		  //分辨率为0的三个子带阈值
	      T10=(b4*var)/sqrt(var10);
	      T11=(b4*var)/sqrt(var11);
	      T12=(b4*var)/sqrt(var12);

    //利用改进阈值进行去噪

	//分辨率为3的三个子带进行去噪
    for(i=0;i<(height/L1);i++)
	  for(j=(width/L1);j<2*(width/L1);j++)
	  {
		  if(fabs(image[i*width+j])<T1)
		  image[i*width+j]=0;
		  if(fabs(image[i*width+j])>=T1)
		  image[i*width+j]=image[i*width+j]*(1-pow(T1/(fabs(image[i*width+j])),q));
		  
		  if(fabs(image[(i+height/L1)*width+j-width/L1])<T2)
		  image[(i+height/L1)*width+j-width/L1]=0;
		  if(fabs(image[(i+height/L1)*width+j-width/L1])>=T2)
	      image[(i+height/L1)*width+j-width/L1]=image[(i+height/L1)*width+j-width/L1]*(1-pow(T2/(fabs(image[(i+height/L1)*width+j-width/L1])),q));
		  		 
		  if(fabs(image[(i+iheight/L1)*iwidth+j])<T3)
		  image[(i+iheight/L1)*iwidth+j]=0;
		  if(fabs(image[(i+iheight/L1)*iwidth+j])>=T3)
	      image[(i+iheight/L1)*iwidth+j]=image[(i+iheight/L1)*iwidth+j]*(1-pow(T3/(fabs(image[(i+iheight/L1)*iwidth+j])),q));
	      
	  }
	
	//分辨率为2的三个子带进行去噪
	for(i=0;i<(height/L2);i++)
	  for(j=(width/L2);j<2*(width/L2);j++)
	  {
		  if(fabs(image[i*width+j])<T4)
		  image[i*width+j]=0;
		  if(fabs(image[i*width+j])>=T4)
		  image[i*width+j]=image[i*width+j]*(1-pow(T4/(fabs(image[i*width+j])),1/q));
		  
		  if(fabs(image[(i+height/L2)*width+j-width/L2])<T5)
		  image[(i+height/L2)*width+j-width/L2]=0;
		  if(fabs(image[(i+height/L2)*width+j-width/L2])>=T5)
	      image[(i+height/L2)*width+j-width/L2]=image[(i+height/L2)*width+j-width/L2]*(1-pow(T5/(fabs(image[(i+height/L2)*width+j-width/L2])),q));
		  		 
		  if(fabs(image[(i+iheight/L2)*iwidth+j])<T6)
		  image[(i+iheight/L2)*iwidth+j]=0;
		  if(fabs(image[(i+iheight/L2)*iwidth+j])>=T6)
	      image[(i+iheight/L2)*iwidth+j]=image[(i+iheight/L2)*iwidth+j]*(1-pow(T6/(fabs(image[(i+iheight/L2)*iwidth+j])),q));
	      
		  
	  }

     //分辨率为1的三个子带进行去噪
	 for(i=0;i<(height/L3);i++)
	    for(j=(width/L3);j<2*(width/L3);j++)
	  {
		 if(fabs(image[i*width+j])<T7)
		  image[i*width+j]=0;
		  if(fabs(image[i*width+j])>=T7)
		  image[i*width+j]=image[i*width+j]*(1-pow(T7/(fabs(image[i*width+j])),q));
		  
		  if(fabs(image[(i+height/L3)*width+j-width/L3])<T8)
		  image[(i+height/L3)*width+j-width/L3]=0;
		  if(fabs(image[(i+height/L3)*width+j-width/L3])>=T8)
	      image[(i+height/L3)*width+j-width/L3]=image[(i+height/L3)*width+j-width/L3]*(1-pow(T8/(fabs(image[(i+height/L3)*width+j-width/L3])),q));
		  		 
		  if(fabs(image[(i+iheight/L3)*iwidth+j])<T9)
		  image[(i+iheight/L3)*iwidth+j]=0;
		  if(fabs(image[(i+iheight/L3)*iwidth+j])>=T9)
	      image[(i+iheight/L3)*iwidth+j]=image[(i+iheight/L3)*iwidth+j]*(1-pow(T9/(fabs(image[(i+iheight/L3)*iwidth+j])),q));
	   	      
	  }
      
	  //分辨率为0的三个子代带进行去噪
	  for(i=0;i<(height/L4);i++)
	     for(j=(width/L4);j<2*(width/L4);j++)
	  {
		 if(fabs(image[i*width+j])<T10)
		  image[i*width+j]=0;
		  if(fabs(image[i*width+j])>=T10)
		  image[i*width+j]=image[i*width+j]*(1-pow(T10/(fabs(image[i*width+j])),q));
		  
		  if(fabs(image[(i+height/L4)*width+j-width/L4])<T11)
		  image[(i+height/L4)*width+j-width/L4]=0;
		  if(fabs(image[(i+height/L4)*width+j-width/L4])>=T11)
	      image[(i+height/L4)*width+j-width/L4]=image[(i+height/L4)*width+j-width/L4]*(1-pow(T11/(fabs(image[(i+height/L4)*width+j-width/L4])),q));
		  		 
		  if(fabs(image[(i+iheight/L4)*iwidth+j])<T12)
		  image[(i+iheight/L4)*iwidth+j]=0;
		  if(fabs(image[(i+iheight/L4)*iwidth+j])>=T12)
	      image[(i+iheight/L4)*iwidth+j]=image[(i+iheight/L4)*iwidth+j]*(1-pow(T12/(fabs(image[(i+iheight/L4)*iwidth+j])),q));
	   	      
	  }

    

}

void CINTDWTDoc::OnSaveimage() 
{
	// TODO: Add your command handler code here
	
	int i,j;
	CString str;
   	BMPFile=lpDst1;

	//建立保存路径
	CString FolderName="E:\\";
    //保存文件的名字
	CString FileName="Spihtimage";

   	//输出图象
	for(i=0;i<iheight;i++)
		for(j=0;j<iwidth;j++)
		{
			BMPFile.SetXY8(j,i,data2[i*iwidth+j]);
			str.Format("%s\\%s.bmp",FolderName,FileName);
		    BMPFile.Save(str);
		}
	
	
}

void CINTDWTDoc::OnOpenoriginalimage() 
{
	// TODO: Add your command handler code here
	CFileDialog fileopenbox(TRUE,NULL,"*.bmp",OFN_HIDEREADONLY,"files(*.bmp)|*.bmp|",NULL);	
    fileopenbox.m_ofn.lpstrTitle="小波变换";
	if(fileopenbox.DoModal()!=IDOK) return;
	CString strPathname=fileopenbox.GetPathName();
	
	if(IsBmp(strPathname))
	{
	
		lpSrc.Read(strPathname);
		
	}
    else
	{
		AfxMessageBox("不是bmp文件！");
	}
	
	int i,j;
    CSize m_size;
	m_size.cx = lpSrc.Width();
	m_size.cy = lpSrc.Height();
	iwidth=m_size.cx ;
	iheight=m_size.cy;
		
    data1=new int[iwidth * iheight];
	data2=new double[iwidth * iheight];
	frame=new double[m_size.cx * m_size.cy];
	
	for( i = 0; i <iheight; i++ )
		for( j = 0; j < iwidth; j++ )
		{
		   data1[i*iwidth+j]=(unsigned char)lpSrc.GetXY8(j,i)-128;
		   data2[i*iwidth+j]=(unsigned char)lpSrc.GetXY8(j,i)-128;
		   frame[i*iwidth+j]=(unsigned char)lpSrc.GetXY8(j,i);
		}
	
	UpdateAllViews(FALSE);
	return ;
}



void CINTDWTDoc::OnIntdwt53() 
{
	// TODO: Add your command handler code here
	int i,j;

	//建立对话框
	CLEVEL Dialog;
    Dialog.m_level = 0;
	if (Dialog.DoModal()!=IDOK)
    return;
	level=Dialog.m_level;
	
	// 删除对话框
	delete Dialog;	
	    
	lpDst1=lpSrc;
   	image1=new int[iwidth*iheight];
	
	//整数小波变换(5\3小波)	
	INTDWT53.Dwt2D(data1, iheight,iwidth,level);

	for(i=0;i<iheight;i++)
	  for(j=0;j<iwidth;j++)
	    image1[i*iwidth+j]=data1[i*iwidth+j];
	
	//规一化小波分解后的图像
	uint8(image1,iheight,iwidth);
   	
     for(i=0;i<iheight;i++)
        for(j=0;j<iwidth;j++)
        {
 	       lpDst1.SetXY8(j,i,image1[i*iwidth+j]);
		}
	
	//释放内存		
	delete []image1;
		
	UpdateAllViews(NULL);
	SetModifiedFlag();
}

void CINTDWTDoc::OnIntidwt53() 
{
	// TODO: Add your command handler code here
	int i,j;
	
	lpDst2=lpSrc;
    
	//整数小波反变换(5\3小波)
	INTDWT53.IDwt2D(data1, iheight,iwidth,level);

	//规一化重建图像
	uint8(data1,iheight,iwidth);

    for(i=0;i<iheight;i++)
        for(j=0;j<iwidth;j++)
        {
 	       lpDst2.SetXY8(j,i,data1[i*iwidth+j]);
		}
	
    //清除内存
	delete []data1;
   	
	UpdateAllViews(NULL);
	SetModifiedFlag();
	
}

void CINTDWTDoc::uint8(double *Imagebuffer, int height, int width)
{
	int i,j,max,min;
	max=min=0;
	
	//搜索最大值和最小值
	for(i=0;i<height;i++)
	  for(j=0;j<width;j++)
	  {
		if (*(Imagebuffer+i*width+j)>max)
	    	  max=*(Imagebuffer+i*width+j);
		if (*(Imagebuffer+i*width+j)<min)
		      min=*(Imagebuffer+i*width+j);
	  }

    //进行规一化
	  for(i=0;i<iheight;i++)
	    for(j=0;j<iwidth;j++)
         *(Imagebuffer+i*iwidth+j)=(*(Imagebuffer+i*iwidth+j)*255)/(max-min)+128;

	//像素幅值控制
	for(i=0;i<height;i++)
    	for(j=0;j<width;j++)
		{
           if(*(Imagebuffer+i*width+j)>255)
		    *(Imagebuffer+i*width+j)=255;
		   if(*(Imagebuffer+i*width+j)<0)
			*(Imagebuffer+i*width+j)=0;	
		}

}

void CINTDWTDoc::uint8(int *Imagebuffer, int height, int width)
{
	int i,j,max,min;
	
	//搜索最大值和最小值
	max=min=0;
	for(i=0;i<height;i++)
	  for(j=0;j<width;j++)
	  {
		if (*(Imagebuffer+i*width+j)>max)
	    	  max=*(Imagebuffer+i*width+j);
		if (*(Imagebuffer+i*width+j)<min)
		      min=*(Imagebuffer+i*width+j);
	  }
    
	//进行规一化
    for(i=0;i<iheight;i++)
	    for(j=0;j<iwidth;j++)
         *(Imagebuffer+i*iwidth+j)=(*(Imagebuffer+i*iwidth+j)*255)/(max-min)+128;
	
	  
	//像素幅值控制
	for(i=0;i<height;i++)
    	for(j=0;j<width;j++)
		{
           if(*(Imagebuffer+i*width+j)>255)
		    *(Imagebuffer+i*width+j)=255;
		   if(*(Imagebuffer+i*width+j)<0)
			*(Imagebuffer+i*width+j)=0;	
		}

}


void CINTDWTDoc::OnSpihtcoding() 
{
	// TODO: Add your command handler code here
    int i,j,bpp;
	Level=3; //设置小波变换的层数
	bpp=0.5;  //每个象素的比特数，通过该参数控制压缩后文件大小
	image1=new int[iheight*iwidth];

			
	//浮点型小波变换
    FLOATDWT.Dwt2D(data2, iheight,iwidth,Level);

	for(i=0;i<iheight;i++)
	  for(j=0;j<iwidth;j++)
	      image1[i*iwidth+j]=int(data2[i*iwidth+j]);
	 	 
	  CSPIHTCoder SPIHTCoder;
	  SPIHTCoder.SetBitRate(bpp);  //设置spiht编码时的码率
	  SPIHTCoder.Encode(image1,iwidth,iheight,Level);  //进行spiht编码
     
	  delete []data2;
	  delete []image1;

	  UpdateAllViews(NULL);
      SetModifiedFlag();
     
	
}



void CINTDWTDoc::OnSpihtdecoding() 
{
	// TODO: Add your command handler code here
	int m,n,MSE;
	float PSNR;
	CString str;
	double *image;
    lpDst1=lpSrc;
	
	CSPIHTCoder SPIHTCoder;
	SPIHTCoder.Decode();     //spiht解码

    iheight=SPIHTCoder.header.height; 
	iwidth=SPIHTCoder.header.width;
    image=new double[iwidth*iheight];

	for(m=0;m<iheight;m++)
	  for(n=0;n<iwidth;n++)
	      image[m*iwidth+n]=SPIHTCoder.MR->m[m][n];

	 //小波反变换
	 FLOATDWT.IDwt2D(image, iheight,iwidth,Level);  

     //将图像归一化
	 uint8(image, iheight, iwidth);

	  
	 for(m=0;m<iheight;m++)
        for(n=0;n<iwidth;n++)
        {
 	       lpDst1.SetXY8(n,m,int(image[m*iwidth+n]));
		}

      //计算PSNR
	   MSE=0;
	   for(m=0;m<iheight;m++)
	      for(n=0;n<iwidth;n++)
		      MSE=MSE+(image[m*iwidth+n]-frame[m*iwidth+n])*(image[m*iwidth+n]-frame[m*iwidth+n]);
		      MSE=MSE/(iwidth*iheight);
		      PSNR=10*log10(255*255/MSE);  
			  str.Format("PSNR为%f",PSNR);
	          AfxMessageBox(str);
	 
   //清除内存
	delete []image;
   	
	UpdateAllViews(NULL);
	SetModifiedFlag();
	
}