intdwtdoc.cpp

该程序实现各种小波提升分解

	  //计算PSNR值
	  for(i=0;i<iwidth;i++)
          for(j=0;j<iheight;j++)
              MSE=MSE+(frame[i*iwidth+j]-data2[i*iwidth+j])*(frame[i*iwidth+j]-data2[i*iwidth+j]);
	          MSE=MSE/(iwidth*iheight);
		      PSNR=10*log10(255*255/MSE);
              str.Format("PSNR为%f",PSNR);
              AfxMessageBox(str);	  
		      
	  UpdateAllViews(NULL);
	  SetModifiedFlag();

}





void CINTDWTDoc::Imagedenoising(int width, int height, int level, float q, double *image)
{
	  int i,j,L1,L2,L3,L4,var;
	  float *subHH1,sum,mean1,mean2,mean3,sum1,sum2,sum3;
	  sum1=0;sum2=0;sum3=0;mean1=0;mean2=0;mean3=0;
	  sum=0; 

	  //计算分辨率为3的各子带的均值和方差
	  L1=1<<(level-3);
	  float var1,var2,var3;
	  var1=0;var2=0;var3=0;
     
      for(i=0;i<(height/L1);i++)
	    for(j=(width/L1);j<2*(width/L1);j++)
	  {
		  sum1=sum1+image[i*width+j];
		  sum2=sum2+image[(i+height/L1)*width+j-(width/L1)];
		  sum3=sum3+image[(i+height/L1)*width+j];
	  }
	  
	  mean1=sum1/((height/L1)*(width/L1));
	  mean2=sum2/((height/L1)*(width/L1));
	  mean3=sum3/((height/L1)*(width/L1));
    
	  for(i=0;i<(height/L1);i++)
	    for(j=(width/L1);j<2*(width/L1);j++)
	  {
		  var1=var1+(image[i*width+j]-mean1)*(image[i*width+j]-mean1);
		  var2=var2+(image[(i+height/L1)*width+j-width/L1]-mean2)*(image[(i+height/L1)*width+j-width/L1]-mean2);
		  var3=var3+(image[(i+height/L1)*width+j]-mean3)*(image[(i+height/L1)*width+j]-mean3);
	  }
	 
	 var1=var1/((height/L1)*(width/L1));
	 var2=var2/((height/L1)*(width/L1));
	 var3=var3/((height/L1)*(width/L1));
	  
     //计算分辨率为2的各子带的均值和方差
	  
	  L2=1<<(level-2);
	  float var4,var5,var6;
	  var4=0;var5=0;var6=0;
      sum1=0;sum2=0;sum3=0;mean1=0;mean2=0;mean3=0;
	
	  for(i=0;i<(height/L2);i++)
	     for(j=(width/L2);j<2*(width/L2);j++)
	  {
		  sum1=sum1+image[i*width+j];
		  sum2=sum2+image[(i+height/L2)*width+j-width/L2];
		  sum3=sum3+image[(i+iheight/L2)*width+j];
	  }
	  mean1=sum1/((height/L2)*(width/L2));
	  mean2=sum2/((height/L2)*(width/L2));
	  mean3=sum3/((height/L2)*(width/L2));
	  
	  
	  for(i=0;i<(height/L2);i++)
	    for(j=(width/L2);j<2*(width/L2);j++)
	  {
		var4=var4+(image[i*width+j]-mean1)*(image[i*width+j]-mean1);
		var5=var5+(image[(i+height/L2)*width+j-width/L2]-mean2)*(image[(i+height/L2)*width+j-width/L2]-mean2);
		var6=var6+(image[(i+iheight/L2)*width+j]-mean3)*(image[(i+height/L2)*width+j]-mean3);
	  }
     
	  var4=var4/((height/L2)*(width/L2));
      var5=var5/((height/L2)*(width/L2));
      var6=var6/((height/L2)*(width/L2));

    //计算分辨率为1的各子带的均值和方差
  
      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="Denoisedimage";

   	//输出图象
	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::OnVideodenoising() 
{
	// TODO: Add your command handler code here


	
}