liftingschemeview.cpp

正交变换.rar文件压缩的基础知识无毒good

					var=var+g[k]*tmp[(abs(i+k-(g0-1))-(ln_g-1)%2)*ww+j];
				else
					var=var+g[k]*tmp[(i+k-(g0-1))*ww+j];
			}
			tmp[i*ww+j+ww/2]=var;
		}
		for(i=g0-1;i<hh-(ln_g-g0);i++)
		{
			var=0;
			for(k=0;k<ln_g;k++)
				var=var+g[k]*tmp[(i+k-(g0-1))*ww+j];
			tmp[i*ww+j+ww/2]=var;
		}
		for(i=hh-(ln_g-g0);i<hh;i++)
		{
			var=0;
			for(k=0;k<ln_g;k++)
			//	var=var+g[k]*tmp[(2*hh-(i+k-(g0-1))-2+(ln_g-1)%2)*ww+j];
			{
				if((i+k-g0+1)<hh)
				var=var+g[k]*tmp[(i+k-g0+1)*ww+j];
				else
					var=var+g[k]*tmp[(2*hh-(i+k-g0+1)-2+(ln_g-1)%2)*ww+j];}
			tmp[i*ww+j+ww/2]=var;
		}
	}//return tmp;

	for(i=0;i<hh/2;i++)
		for(j=0;j<ww/2;j++)
		{
			if(Low)
				resultImage[i*(ww/2)+j]=tmp[2*i*ww+j+ww/2];
			else
				resultImage[(i+hh/2)*(ww/2)+j]=tmp[2*i*ww+j+ww/2];
		}
//计算LL HH
	for(j=0;j<ww/2;j++)
	{
		for(i=0;i<h0-1;i++)
		{
			var=0;
			for(k=0;k<ln_h;k++)
			//	var=var+h[k]*tmp[(abs(i+k-(h0-1))-(ln_h-1)%2)*ww+j];
			{
					if((i+k-(h0-1)<0))
					var=var+h[k]*tmp[(abs(i+k-(h0-1))-(ln_h-1)%2)*ww+j];
				else
					var=var+h[k]*tmp[(i+k-(h0-1))*ww+j];
			}
			tmp[i*ww+j+ww/2]=var;
		}
		for(i=h0-1;i<hh-(ln_h-h0);i++)
		{
			var=0;
			for(k=0;k<ln_h;k++)
				var=var+h[k]*tmp[(i+k-(h0-1))*ww+j];
			tmp[i*ww+j+ww/2]=var;
		}
		for(i=hh-(ln_h-h0);i<hh;i++)
		{
			var=0;
			for(k=0;k<ln_h;k++)
			{
				if((i+k-h0+1)<hh)
				var=var+h[k]*tmp[(i+k-h0+1)*ww+j];//(2*hh-(i+k-(h0-1))-2+(ln_h-1)%2)
				else
				var=var+h[k]*tmp[(2*hh-(i+k-h0+1)-2+(ln_h-1)%2)*ww+j];
			}
			tmp[i*ww+j+ww/2]=var;
		}
	}//return tmp;
   
	for(i=hh/2;i<hh;i++)
		for(j=0;j<ww/2;j++)
		{
			if(Low)
				resultImage[i*(ww/2)+j]=tmp[2*(i-hh/2)*ww+j+ww/2];
			else
				resultImage[(i-hh/2)*(ww/2)+j]=tmp[2*(i-hh/2)*ww+j+ww/2];
		}
//	DisplayImage(m_orgWidth/2,m_orgHeight,0,resultImage);
	return resultImage;
}

void CLiftingschemeView::DisplayInfoBox(char tp[], int level, char method[], long timespan,BOOL ReCon)
{
		CInfo InfoDlg;
		char *buffer=new char[10];
		InfoDlg.m_Info_wavelet=tp;
		_itoa(level,buffer,10);
		InfoDlg.m_Info_level=buffer;
		InfoDlg.m_Info_method="Mallat";
		_itoa(timespan,buffer,10);
		InfoDlg.m_Info_time=buffer;
		InfoDlg.m_Info_Image=m_szFileName;
		_itoa(m_orgWidth,buffer,10);
		InfoDlg.m_Info_width=buffer;
		_itoa(m_orgHeight,buffer,10);
		InfoDlg.m_Info_height=buffer;
		if(ReCon)
		{
			double Pnsr1;
			Pnsr1=PNSR(m_pImageData,ReconImage,m_orgWidth,m_orgHeight);
			char TempString[10];
			sprintf(TempString,"%6.2f",Pnsr1);
				
			InfoDlg.m_Info_pnsr=TempString;
			
		}
		
		InfoDlg.DoModal();
}

void CLiftingschemeView::Reconstruct2(double Hi_R[], int Hi_R0, int Hi_R_l, double Lo_R[], int Lo_R0, int Lo_R_l, double *Image, long www, long hhh, int level, char tp[], BOOL disp)
{
	long ww,hh;int i,j,loop;
	int shrink=0;
	ReconImage=new double[(www-2*shrink)*(hhh-2*shrink)];
	double *tmp;//=new double[www*hhh];
	double *tmp2=new double[www*hhh];
	long timespan=0;
	loop=level;
	for(loop;loop>0;loop--)
	{
		ww=www>>(loop-1);
		hh=hhh>>(loop-1);
		long t1=GetTickCount();
		tmp=new double[ww*hh];
		tmp=RealReconstruct2(Hi_R,Hi_R0,Hi_R_l,Lo_R,Lo_R0,Lo_R_l,Image,www,hhh,loop);
		t1=GetTickCount()-t1;
		timespan=timespan+t1;
		for(i=hhh-hh;i<hhh;i++)
			for(j=0;j<ww;j++)
				Image[i*www+j]=tmp[(i+hh-hhh)*ww+j];//[i*www+j]=tmp[i*ww+j];ReconImagetmp2
		delete[] tmp;
		}
//	delete[]tmp;

//////动态范围调整//////////
	double mx=-3000;double mn=3000;
	for(i=0+shrink;i<www-shrink;i++)
		for(j=0+shrink;j<hhh-shrink;j++)
		{
			if(Image[i*www+j]>mx)
				mx=Image[i*www+j];
			if(Image[i*www+j]<mn)
				mn=Image[i*www+j];
		}
	for(i=0+shrink;i<www-shrink;i++)
		for(j=0+shrink;j<hhh-shrink;j++)
		ReconImage[(i-shrink)*(www-2*shrink)+j-shrink]=(Image[i*www+j]-mn)*255/(mx-mn);
		//delete[]tmp2;
///////////////动态范围调整结束////////////////

	DisplayImage(www-2*shrink,hhh-2*shrink,0,FALSE,ReconImage);//ReconImage
	if(disp)
		DisplayInfoBox(tp,level,"Mallat",timespan,TRUE);
}

double* CLiftingschemeView::RealReconstruct2(double Hi_R[], int Hi_R0, int Hi_R_l, double Lo_R[], int Lo_R0, int Lo_R_l, double *Image, long www, long hhh, int loop)
{
		long ww,hh;
	ww=www>>loop;
	hh=hhh>>loop;
	//BYTE *tmp=new BYTE[ww*hh*4];
	double *tmp1=new double[ww*hh*4];//=;NULL
	double *tmp2=new double[ww*hh*4];
	double *tmp3=new double[ww*hh*4];
	double *tmp4=new double[ww*hh*4];
	double *subImage=new double[ww*hh*4];
	int i,j;
	tmp1=SubReconstruct2(Lo_R,Lo_R0,Lo_R_l,Lo_R,Lo_R0,Lo_R_l,Image,ww,hh,LL);

	tmp2=SubReconstruct2(Lo_R,Lo_R0,Lo_R_l,Hi_R,Hi_R0,Hi_R_l,Image,ww,hh,LH);//
	tmp3=SubReconstruct2(Hi_R,Hi_R0,Hi_R_l,Lo_R,Lo_R0,Lo_R_l,Image,ww,hh,HL);//
	tmp4=SubReconstruct2(Hi_R,Hi_R0,Hi_R_l,Hi_R,Hi_R0,Hi_R_l,Image,ww,hh,HH);
/*double maxx2=-2000,minn2=2000;
	for(i=0;i<2*hh;i++)
		for(j=0;j<2*ww;j++)
		{
			if(tmp2[i*(2*ww)+j]<minn2)
				minn2=tmp2[i*(2*ww)+j];
			if(tmp2[i*(2*ww)+j]>maxx2)
				maxx2=tmp2[i*(2*ww)+j];
		}
		double maxx1=-2000, minn1=2000;
	for(i=0;i<2*hh;i++)
		for(j=0;j<2*ww;j++)
		{
			if(tmp1[i*(2*ww)+j]<minn1)
				minn1=tmp1[i*(2*ww)+j];
			if(tmp1[i*(2*ww)+j]>maxx1)
				maxx1=tmp1[i*(2*ww)+j];
		}*/

	for(i=0;i<2*hh;i++)
		for(j=0;j<2*ww;j++)
            subImage[i*(2*ww)+j]=tmp1[i*(2*ww)+j]+tmp4[i*(2*ww)+j]+tmp2[i*(2*ww)+j]+tmp3[i*(2*ww)+j];//tmp1[i*(2*ww)+j]+;//+tmp2[i*(2*ww)+j];//+tmp3[i*(2*ww)+j]+tmp4[i*(2*ww)+j];
	/*double maxx,minn;
		maxx=-2000;minn=2000;
	for(i=0;i<2*hh;i++)
		for(j=0;j<2*ww;j++)
		{
			if(subImage[i*(2*ww)+j]>maxx)
				maxx=subImage[i*(2*ww)+j];
			if(subImage[i*(2*ww)+j]<minn)
				minn=subImage[i*(2*ww)+j];
		}
	for(i=0;i<2*hh;i++)
		for(j=0;j<2*ww;j++)
            tmp1[i*(2*ww)+j]=(subImage[i*(2*ww)+j]-minn)*255.0/(maxx-minn);
		delete[]tmp2;
		delete[]tmp3;
		delete[]tmp4;
*/
//	DisplayImage(ww*2,hh*2,256,FALSE,subImage);//tmp1
	//MessageBox("subRec");

	return subImage;//tmp1

}

double* CLiftingschemeView::SubReconstruct2(double Lo_R[], int Lo_R0, int Lo_R_l, double Hi_R[], int Hi_R0, int Hi_R_l, double *Image, long ww, long hh, int p)
{

	int x,y,i,j,k;double var;
	if(p==LL){x=1;y=0;}
	if(p==LH){x=0;y=0;}
	if(p==HL){x=1;y=1;}
	if(p==HH){x=0;y=1;}
	double *tmp1=new double[ww*hh*2];
	double *tmp11=new double[ww*hh*2];
	double *tmp2=new double[ww*hh*4];
	double *result=new double[ww*hh*4];
//	BYTE *tp=new BYTE[ww*hh];
/////////////////	
	for(i=0;i<hh;i++)
		for(j=0;j<ww;j++)
		{///要重构的区域
		//	tp[i*ww+j]=Image[(i+m_orgHeight-2*hh+x*hh)*m_orgWidth+j+y*ww];
			tmp1[2*i*ww+j]=Image[(i+m_orgHeight-2*hh+x*hh)*m_orgWidth+j+y*ww];
			tmp1[(2*i+1)*ww+j]=0;
		}
		//DisplayImage(ww,hh,0,FALSE,tp);
	/*	MessageBox("special");
		for(i=0;i<hh;i++)
		for(j=0;j<ww;j++)
		{
			Image[(i+m_orgHeight-2*hh+x*hh)*m_orgWidth+j+y*ww]=255;
		}
		DisplayImage(m_orgWidth,m_orgWidth,0,FALSE,Image);
	MessageBox("upsample");
		DisplayImage(ww,hh*2,0,FALSE,tmp1);*/
		
//////////////第一次上采样结束，列////
//////列卷积

	for(j=0;j<ww;j++)
	{
	
		for(i=0;i<Hi_R0-1;i++)
		{
			var=0;
			for(k=0;k<Hi_R_l;k++)
				//var=var+Hi_R[k]*tmp1[(abs(i+k-(Hi_R0-1))-(Hi_R_l-1)%2)*ww+j];
			{	if( (i+k-(Hi_R0-1))<0 )
					var=var+Hi_R[k]*tmp1[(abs(i+k-(Hi_R0-1))-(Hi_R_l-1)%2)*ww+j];
				else
					var=var+Hi_R[k]*tmp1[(i+k-(Hi_R0-1))*ww+j];
			}
					
			tmp11[i*ww+j]=var;
		}
	
		for(i=Hi_R0-1;i<2*hh-(Hi_R_l-Hi_R0);i++)
		{
			var=0;
			for(k=0;k<Hi_R_l;k++)
				var=var+Hi_R[k]*tmp1[(i+k-(Hi_R0-1))*ww+j];
		
			tmp11[i*ww+j]=var;
		}
		for(i=2*hh-(Hi_R_l-Hi_R0);i<2*hh;i++)
		{
			var=0;
			for(k=0;k<Hi_R_l;k++)
				//var=var+Hi_R[k]*tmp1[(4*hh-(i+k-(Hi_R0-1))-2+(Hi_R_l-1)%2)*ww+j];
			{
				if((i+k-Hi_R0+1)<2*hh)
					var=var+Hi_R[k]*tmp1[(i+k-Hi_R0+1)*ww+j];//(2*hh-(i+k-(h0-1))-2+(ln_h-1)%2)
				else
					var=var+Hi_R[k]*tmp1[(4*hh-(i+k-Hi_R0+1)-2+(Hi_R_l-1)%2)*ww+j];
			}
			
			
			tmp11[i*ww+j]=var;
		}
	}
	////第二次上采样，行////////
		for(j=0;j<ww;j++)
			for(i=0;i<(2*hh);i++)
			{
				tmp2[i*2*ww+2*j]=tmp11[i*ww+j];
				tmp2[i*2*ww+2*j+1]=0;
			}
		//	return tmp2;
///////////////////////////////////////
////////行卷积/////
		for(i=0;i<2*hh;i++)
		{
			for(j=0;j<Lo_R0-1;j++)
			{
				var=0;
				for(k=0;k<Lo_R_l;k++)
				//	var=var+Lo_R[k]*tmp2[i*2*ww+abs(j+k-(Lo_R0-1))-(Lo_R_l-1)%2];
				{
					if( (j+k-(Lo_R0-1))<0 )
					var=var+Lo_R[k]*tmp2[i*2*ww+abs(j+k-(Lo_R0-1))-(Lo_R_l-1)%2];
				else
					var=var+Lo_R[k]*tmp2[i*ww+j+k-(Lo_R0-1)];
			
				}
				result[i*ww*2+j]=var;
			}
		
			for(j=Lo_R0-1;j<2*ww-(Lo_R_l-Lo_R0);j++)
			{
				var=0;
				for(k=0;k<Lo_R_l;k++)
					var=var+Lo_R[k]*tmp2[i*2*ww+j+k-(Lo_R0-1)];
				result[i*ww*2+j]=var;
			}
			for(j=2*ww-(Lo_R_l-Lo_R0);j<2*ww;j++)
			{
				var=0;
				for(k=0;k<Lo_R_l;k++)
				//	var=var+Lo_R[k]*tmp2[i*2*ww+4*ww-(j+k-(Lo_R0-1))-2+(Lo_R_l-1)%2];
				{if((j+k-Lo_R0+1)<2*ww)
				var=var+Lo_R[k]*tmp2[(j+k-Lo_R0+1)*2*ww+j];//(2*hh-(i+k-(h0-1))-2+(ln_h-1)%2)
				else
					var=var+Lo_R[k]*tmp2[(4*ww-(j+k-Lo_R0+1)-2+(Lo_R_l-1)%2)*2*ww+j];}
			
				result[i*ww*2+j]=var;
			}
		}	

	delete[]tmp2;
	delete[]tmp1;
	delete[]tmp11;
	return result;
}



double  CLiftingschemeView::PNSR(BYTE *orgImage, double *reImage, long w, long h)
{
	int i,j;
	double sum=0;
	double pnsr=0;
	double maxx=0;
	for(i=0;i<h;i++)
		for(j=0;j<w;j++)
		{
			if(orgImage[i*w+j]>maxx)
				maxx=orgImage[i*w+j];
			sum=sum+(orgImage[i*w+j]-reImage[i*w+j])*(orgImage[i*w+j]-reImage[i*w+j]);
		}
		if(sum==0)
			return 0;
		else
		{
			pnsr=10*log10(w*h*maxx*maxx/sum);
			return pnsr;
		}

}

void CLiftingschemeView::OnTransFFT() 
{
	CChoose dlg;
	BOOL reconstruct;
	if(dlg.DoModal()==IDOK)
		reconstruct=TRUE;
	else
		reconstruct=FALSE;

    FFTtransform(m_pImageData,m_orgWidth,m_orgHeight,reconstruct);
	
}

void CLiftingschemeView::OnTransDCT() 
{
	CChoose dlg;
	BOOL reconstruct;
	if(dlg.DoModal()==IDOK)
		reconstruct=TRUE;
	else
		reconstruct=FALSE;
    DCTtransform(m_pImageData,m_orgWidth,m_orgHeight,reconstruct,0);

}

void CLiftingschemeView::FFTtransform(BYTE *Image, long ww, long hh, BOOL reconstruct)
{
	finalImage=new double[ww*hh];
	ImageToRe=new double[ww*hh];
	memset(finalImage, 0, sizeof(double) * ww* hh);
	memset(ImageToRe, 0, sizeof(double) * ww* hh);
	unsigned char *lpSrc;
	// 中间变量
	double	dTemp;
	// 循环变量
	LONG	i;
	LONG	j;
	LONG	w=1;
	LONG	h=1;
	int		wp=0;
	int		hp=0;
	LONG	lLineBytes;
	lLineBytes = WIDTHBYTES(ww * 8);
	
	// 计算进行付立叶变换的宽度和高度（2的整数次方）
	while(w * 2 <= ww)
	{
		w *= 2;
		wp++;
	}
	while(h * 2 <= hh)
	{
		h *= 2;
		hp++;
	}
	// 分配内存
	complex<double> *TD = new complex<double>[w * h];
	complex<double> *FD = new complex<double>[w * h];
	// 行
	for(i = 0; i < h; i++)
	{
		// 列
		for(j = 0; j < w; j++)
		{
			lpSrc = (unsigned char*)Image + lLineBytes * (hh - 1 - i) + j;
			TD[j + w * i] = complex<double>(*(lpSrc), 0);
		}
	}
	for(i = 0; i < h; i++)
	{
	    // 对y方向进行快速付立叶变换
		FFT(&TD[w * i], &FD[w * i], wp);
	}
	// 保存变换结果
	for(i = 0; i < h; i++)
	{
		for(j = 0; j < w; j++)
		{
			TD[i + h * j] = FD[j + w * i];
		}
	}
	for(i = 0; i < w; i++)
	{
		// 对x方向进行快速付立叶变换
		FFT(&TD[i * h], &FD[i * h], hp);
	}
	// 行
	for(i = 0; i < h; i++)
	{
		// 列
		for(j = 0; j < w; j++)
		{
			// 计算频谱
			dTemp = sqrt(FD[j * h + i].real() * FD[j * h + i].real() + 
				         FD[j * h + i].imag() * FD[j * h + i].imag()) / 100;
			
			dTemp=(dTemp > 255)?255:dTemp;
			ImageToRe[lLineBytes * (hh - 1 - (i<h/2 ? i+h/2 : i-h/2))+ (j<w/2 ? j+w/2 : j-w/2)]=dTemp;
		}
	}

    //反变换
	for(i = 0; i < w; i++)
	{
		// 对x方向进行快速付立叶反变换
		IFFT(&FD[i * h], &TD[i * h], hp);
	}
	for(i = 0; i < h; i++)
	{
		for(j = 0; j < w; j++)
		{
			FD[i + h * j]= TD[j + w * i];
		}
	}