wavelettransform.cpp

基于小波的SAR斑点处理

			}

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

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

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

// 小波逆变换,要求遥感图像的宽度和高度均能被2^times整除.
void InvWTBIOTS2D2(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;
		}
		
		// 先进行列变换
		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;
			}

			// 先求出下标为奇数的信号
			for(j=0; j<height-2; j=j+2)
			{
				pY[j+1] = -data[imgWidth*(height+j)/2+i] + (data[imgWidth*j/2+i] + data[imgWidth*(j/2+1)+i])/2;
			}
			pY[height-1] = -data[(height-1)*imgWidth+i] + data[(height/2-1)*imgWidth+i];
			
			// 再求出下标为偶数的信号
			pY[0] = 2*data[i] - pY[1];
			for(j=2; j<height; j=j+2)
			{
				pY[j] = (4*data[(j/2)*imgWidth+i]-pY[j-1]-pY[j+1])/2;
			}

			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;
			}

			// 先求出下标为奇数的信号
			for(j=0; j<width-2; j=j+2)
			{
				pX[j+1] = -data[i*imgWidth+(width+j)/2] + (data[i*imgWidth+j/2] + data[i*imgWidth+j/2+1])/2;
			}
			pX[width-1] = -data[i*imgWidth+width-1] + data[i*imgWidth+width/2-1];

			// 再求出下标为偶数的信号
			pX[0] = 2*data[i*imgWidth]-pX[1];
			for(j=2; j<width; j=j+2)
			{
				pX[j] = (4*data[i*imgWidth+j/2]-pX[j-1]-pX[j+1])/2;
			}
						
			for(j=0; j<width; j++)
				data[i*imgWidth+j] = pX[j];
		}

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

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

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

// 小波变换,要求遥感图像的宽度和高度均能被2^times整除.
// 对边界作了对称延拓处理,比如 3 2 1 0 1 2 3.......
void WTBIOTS3D1(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] = ( 3*data[i*imgWidth] + 4*data[i*imgWidth+1] + data[i*imgWidth+2] )/8;
			pX[width/2] = ( 3*data[i*imgWidth] - 2*data[i*imgWidth+1] - data[i*imgWidth+2] )/4;

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

			// 更新原始数据
			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] = ( 3*data[i] + 4*data[imgWidth+i] + data[2*imgWidth+i] )/8;
			pY[height/2] = ( 3*data[i] - 2*data[imgWidth+i] - data[2*imgWidth+i] )/4;
			
			// 中间部分处理
			for(j=2; j<height-2; j=j+2)
			{
				pY[j/2] = ( data[(j-1)*imgWidth+i] + 3*data[j*imgWidth+i]
							+ 3*data[(j+1)*imgWidth+i] + data[(j+2)*imgWidth+i] )/8;
				pY[(height+j)/2] = ( data[(j-1)*imgWidth+i] + 3*data[j*imgWidth+i]
							- 3*data[(j+1)*imgWidth+i] - data[(j+2)*imgWidth+i] )/4;
			}

			// 下边界处理
			pY[height/2-1] = ( data[(height-3)*imgWidth+i] + 4*data[(height-2)*imgWidth+i]
							+ 3*data[(height-1)*imgWidth+i] )/8;
			pY[height-1] = ( data[(height-3)*imgWidth+i] + 2*data[(height-2)*imgWidth+i]
							- 3*data[(height-1)*imgWidth+i] )/4;
			
			// 更新原始数据
			for(j=0; j<height;j++)
				data[j*imgWidth+i] = pY[j];
		}

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

// 小波逆变换,要求遥感图像的宽度和高度均能被2^times整除.
void InvWTBIOTS3D1(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 * pXP  = (float *)new float[width/2];
		float * pXM  = (float *)new float[width/2];
		float * pYP  = (float *)new float[height/2];
		float * pYM  = (float *)new float[height/2];
		if( pXP==NULL || pXM==NULL ||
			pYP==NULL || pYM==NULL )
		{
			AfxMessageBox("申请内存失败!");
			if(pXP != NULL) delete pXP;
			if(pXM != NULL) delete pXM;
			if(pYP != NULL) delete pYP;
			if(pYM != NULL) delete pYM;
			return;
		}

		// 先进行列变换
		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;
			}

			// 获取中间值,取出全部有用信息,以便后面可以抹去工作数据原有信息
			for(j=0; j<height; j=j+2)
			{
				pYP[j/2] = 2*data[j/2*imgWidth+i] + data[(height+j)/2*imgWidth+i];
				pYM[j/2] = 2*data[j/2*imgWidth+i] - data[(height+j)/2*imgWidth+i];
			}

			// 先求出下标序号为1、2、3、4、5至2N-2的信号
			for(j=0; j<height-2; j=j+2)
			{
				data[(j+1)*imgWidth+i] = ( 3*pYM[j/2] - pYP[j/2+1] )/4;
				data[(j+2)*imgWidth+i] = ( 3*pYP[j/2+1] - pYM[j/2] )/4;
			}

			// 再求出下标序号为0、2N-1的信号
			data[i] = ( 2*pYP[0] - data[imgWidth+i] )/3;
			data[(height-1)*imgWidth+i] = ( 2*pYM[height/2-1] - data[(height-2)*imgWidth+i] )/3;
		}
			
		// 再进行行变换
		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)
			{
				pXP[j/2] = 2*data[i*imgWidth+j/2] + data[i*imgWidth+(width+j)/2];
				pXM[j/2] = 2*data[i*imgWidth+j/2] - data[i*imgWidth+(width+j)/2];
			}

			// 先求出下标序号为1、2、3、4、5、2N-2的信号
			for(j=0; j<width-2; j=j+2)
			{
				data[i*imgWidth+j+1] = ( 3*pXM[j/2] - pXP[j/2+1] )/4;
				data[i*imgWidth+j+2] = ( 3*pXP[j/2+1] - pXM[j/2] )/4;
			}

			// 再求出下标序号为0、2N-1的信号
			data[i*imgWidth] = ( 2*pXP[0] - data[i*imgWidth+1] )/3;
			data[i*imgWidth+width-1] = ( 2*pXM[width/2-1] - data[i*imgWidth+width-2] )/3;
		}

		// 释放临时工作内存
		delete pXP;
		delete pXM;
		delete pYP;
		delete pYM;
	}
}