wavelettransform.cpp

基于小波的SAR斑点处理

//WaveletTransform.cpp------WaveletTransform.h的执行文件
//对遥感图像进行对称紧支集双正交小波变换和逆变换
//开发者:		余家忠
//开发单位:		北京大学遥感与地理信息系统研究所GeoSIS实验室
//开发时间:		2000年3月22日

#include "stdafx.h"
#include "WaveletTransform.h"
#include "Comput.h"
#include "Progress.h"

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

// 说明:以下每组小波变换的参数均应乘以2^0.5,为加快运算速度,将该因子去掉.

// ****************************S=1,D=3的情形****************************

// 变换系数,其中H[]、G[]用于正变换,^H[]、^G[]用于逆变换
// 其中G[]和^H[]中的变换值均为应有值的2倍,在做有关统计时应注意.
//           -2      -1      0     1    2        3
// H[] = {    0       0     0.5,  0.5   0        0  };
// G[] = { -0.125, -0.125,   1,   -1,  0.125,  0.125};
//^H[] = { -0.125,  0.125,   1     1   0.125  -0.125};
//^G[] = {    0       0      0.5 -0.5   0        0  }

// 小波变换,要求遥感图像的宽度和高度必须能被2^times整除.
void WTBIOTS1D3(float * data,int imgHeight,int imgWidth,int times)
{
	int i,j,scale,time;
	int width,height;

	//进行m次小波分解
	for(time=0; time<times; time++)
	{
		// 为临时工作区申请内存
		scale = 1;
		for(i=0; i<time;i++)
			scale <<= 1;
		width  = imgWidth/scale;
		height = imgHeight/scale;
		
		float * pX = (float *)new float[width];
		float * pY = (float *)new float[height];
		if( pX==NULL || pY==NULL )
		{
			AfxMessageBox("申请内存失败!");
			if( pX != NULL )delete pX;
			if( pY != NULL )delete pY;
			return;
		}

		// 先进行行变换
		WORD pro,process = 0;
		for(i=0; i<height; i++)
		{
			// 进度指示
			pro=(int)(100.0 * (i+1) / height);
			if(pro>process)
			{
				for(j=0; j<pro-process; j++)
					UpdateStatusBar();
				process=pro;
			}

			// 求低频部分
			for(j=0; j<width; j=j+2)
			{
				pX[j/2] = float( (data[i*imgWidth+j]+data[i*imgWidth+j+1])/2 );
			}
			
			// 求高频部分
			pX[width/2] = data[i*imgWidth] - data[i*imgWidth+1];
			pX[width/2] = float( pX[width/2]-(pX[0]-pX[1])/2 );
			
			for(j=2; j<width-2; j=j+2)
			{
				pX[(j+width)/2] = data[i*imgWidth+j]-data[i*imgWidth+j+1];
				pX[(j+width)/2] = float( pX[(j+width)/2]-(pX[j/2-1]-pX[j/2+1])/4 );
			}

			pX[width-1] = data[i*imgWidth+width-2] - data[i*imgWidth+width-1];
			pX[width-1] = float( pX[width-1] - (pX[width/2-2]-pX[width/2-1])/2 );

			// 更新原始数据
			for(j=0; j<width; j++)
				data[i*imgWidth+j] = pX[j];
		}
			
		// 再进行列变换
		for(i=0; i<width; i++)
		{
			// 进度指示
			pro=(int)(100.0 * (i+1) / width);
			if(pro>process)
			{
				for(j=0; j<pro-process; j++)
					UpdateStatusBar();
				process=pro;
			}

			// 求低频部分
			for(j=0; j<height; j=j+2)
			{
				pY[j/2] = float( (data[j*imgWidth+i]+data[(j+1)*imgWidth+i])/2 );
			}

			// 求高频部分
			pY[height/2] = data[i]-data[imgWidth+i];
			pY[height/2] = float( pY[height/2] - (pY[0]-pY[1])/2 );

			for(j=2; j<height-2; j=j+2)
			{
				pY[(height+j)/2] = data[j*imgWidth+i]-data[(j+1)*imgWidth+i];
				pY[(height+j)/2] = float( pY[(height+j)/2]	- (pY[j/2-1]-pY[j/2+1])/4 );
			}
			
			pY[height-1] = data[(height-2)*imgWidth+i]-data[(height-1)*imgWidth+i];
			pY[height-1] = float( pY[height-1] - (pY[height/2-2]-pY[height/2-1])/2 );

			// 更新原始数据
			for(j=0; j<height;j++)
				data[j*imgWidth+i] = pY[j];
		}

		// 释放临时工作内存
		delete pX;
		delete pY;
	}
}

// 小波逆变换,要求遥感图像的宽度和高度均能被2^times整除.
void InvWTBIOTS1D3(float * data,int imgHeight,int imgWidth,int times)
{
	int i,j,scale,time;
	int width,height;
	
	//进行m次小波分解
	for(time=times-1; time>=0; time--)
	{
		// 为临时工作区申请内存
		scale = 1;
		for(i=0; i<time;i++)
			scale <<= 1;
		width  = imgWidth/scale;
		height = imgHeight/scale;
		float * pX  = (float *)new float[width];
		float * pY  = (float *)new float[height];
		if( pX==NULL || pY==NULL )
		{
			AfxMessageBox("申请内存失败!");
			if(pX != NULL) delete pX;
			if(pY != NULL) delete pY;
			return;
		}
		
		float temp;

		// 先进行列变换
		WORD pro,process = 0;
		for(i=0; i<width; i++)
		{
			// 进度指示
			pro=(int)(100.0 * (i+1) / width);
			if(pro>process)
			{
				for(j=0; j<pro-process; j++)
					UpdateStatusBar();
				process=pro;
			}

			// 上边界处理
			temp  = float( data[imgWidth*height/2+i]+(data[i]-data[imgWidth+i])/2 );
			temp  = temp/2;
			pY[0] = data[i] + temp;
			pY[1] = data[i] - temp;
						
			// 中间部分处理
			for(j=2; j<height-2; j=j+2)
			{
				temp    = data[imgWidth*(height+j)/2+i]+(data[imgWidth*(j/2-1)+i]-data[imgWidth*(j/2+1)+i])/4;
				temp    = temp/2;
				pY[j]   = data[imgWidth*j/2+i] + temp;
				pY[j+1] = data[imgWidth*j/2+i] - temp;
			}

			// 下边界处理
			temp = data[imgWidth*(height-1)+i]+(data[imgWidth*(height/2-2)+i]-data[imgWidth*(height/2-1)+i])/2;
			temp = temp/2;
			pY[height-2] = data[imgWidth*(height/2-1)+i] + temp;
			pY[height-1] = data[imgWidth*(height/2-1)+i] - temp;
						
			for(j=0; j<height;j++)
				data[j*imgWidth+i] = pY[j];
		}
			
		// 再进行行变换
		for(i=0; i<height; i++)
		{
			// 进度指示
			pro=(int)(100.0 * (i+1) / height);
			if(pro>process)
			{
				for(j=0; j<pro-process; j++)
					UpdateStatusBar();
				process=pro;
			}

			// 左边界处理
			temp  = data[i*imgWidth+width/2] + (data[i*imgWidth]-data[i*imgWidth+1]) / 2;
			temp  = temp/2;
			pX[0] = data[i*imgWidth] + temp;
			pX[1] = data[i*imgWidth] - temp;
			
			// 中间部分处理
			for(j=2; j<width-2; j=j+2)
			{
				temp    = data[i*imgWidth+(width+j)/2] + (data[i*imgWidth+j/2-1]-data[i*imgWidth+j/2+1]) / 4;
				temp    = temp/2;
				pX[j]   = data[i*imgWidth+j/2] + temp;
				pX[j+1] = data[i*imgWidth+j/2] - temp;
			}

			// 右边界处理
			temp = data[i*imgWidth+width-1] + (data[i*imgWidth+width/2-2]-data[i*imgWidth+width/2-1])/2;
			temp = temp/2;
			pX[width-2] = data[i*imgWidth+width/2-1] + temp;
			pX[width-1] = data[i*imgWidth+width/2-1] - temp;
			
			for(j=0; j<width; j++)
				data[i*imgWidth+j] = pX[j];
		}

		// 释放临时工作内存
		delete pX;
		delete pY;
	}
}

// ****************************S=2,D=2的情形****************************

//           -2     -1     0     1      2  	   3
// H[] = {    0    0.25   0.5   0.25    0      0    };
// G[] = {	  0    0.125  0.25 -0.75   0.25   0.125 };
//^H[] = { -0.125  0.25   0.75  0.25  -0.125   0    };
//^G[] = {    0      0    0.25  -0.5   0.25    0    };

// 小波变换,要求遥感图像的宽度和高度均能被2^times整除.
// 对边界作了对称延拓处理,比如 3 2 1 0 1 2 3.......
void WTBIOTS2D2(float * data,int imgHeight,int imgWidth,int times)
{
	int i,j,scale,time;
	int width,height;

	//进行m次小波分解
	for(time=0; time<times; time++)
	{
		// 为临时工作区申请内存
		scale = 1;
		for(i=0; i<time;i++)
			scale <<= 1;
		width  = imgWidth/scale;
		height = imgHeight/scale;
		
		float * pX = (float *)new float[width];
		float * pY = (float *)new float[height];
		if( pX==NULL || pY==NULL )
		{
			AfxMessageBox("申请内存失败!");
			if( pX != NULL )delete pX;
			if( pY != NULL )delete pY;
			return;
		}

		// 先进行行变换
		WORD pro,process = 0;
		for(i=0; i<height; i++)
		{
			// 进度指示
			pro=(int)(100.0 * (i+1) / height);
			if(pro>process)
			{
				for(j=0; j<pro-process; j++)
					UpdateStatusBar();
				process=pro;
			}

			// 求低频部分
			pX[0] = ( data[i*imgWidth] + data[i*imgWidth+1] )/2;
			for(j=2; j<width; j=j+2)
			{
				pX[j/2] = (data[i*imgWidth+j-1] + 2*data[i*imgWidth+j] + data[i*imgWidth+j+1])/4;
			}
			
			// 求高频部分
			for(j=0; j<width-2; j=j+2)
			{
				pX[(j+width)/2] = -data[i*imgWidth+j+1]+(pX[j/2]+pX[j/2+1])/2;
			}
			pX[width-1] = -data[i*imgWidth+width-1] + pX[width/2-1];

			// 更新原始数据
			for(j=0; j<width; j++)
				data[i*imgWidth+j] = pX[j];
		}
			
		// 再进行列变换
		for(i=0; i<width; i++)
		{
			// 进度指示
			pro=(int)(100.0 * (i+1) / width);
			if(pro>process)
			{
				for(j=0; j<pro-process; j++)
					UpdateStatusBar();
				process=pro;
			}

			// 求低频部分
			pY[0] = ( data[i] + data[imgWidth+i] )/2;
			for(j=2; j<height; j=j+2)
			{
				pY[j/2] = (data[(j-1)*imgWidth+i] + 2*data[j*imgWidth+i] + data[(j+1)*imgWidth+i])/4;