sample.c

可调节自适应滤波器-代码-ERDAS版本

	}
	


	//**************************************************************************************
	
	//向输出图层写入数据
	eimg_LayerWrite(md->outlayer,0,0,md->outlayer->width,md->outlayer->height,md->pixelblock,&err);
	if (err) 
	{
		return 8;
    }



	//关闭文件,做清理工作.
	eimg_PixelRectDelete(md->pixelblock,&err);
	if (err) 
	{
		return 9;
    }
	eimg_LayerClose(md->outlayer,&err);
	if (err) 
	{
		return 10;
    }
	estr_StringListDelete(md->layerNames, &err);
	if (err) 
	{
		return 11;
    }

	//tempstring = estr_Sprintf(NULL,"Done",&err);
	emet_MeterInfoChangeTaskName(meterinfo,"Done",&err);

	emet_MeterInfoDelete(meterinfo,&err);

	return 0;
}




static int
FilterImage(md)
MyData *md;
{
	Eerr_ErrorReport *err=NULL;
	int i,j,k;//temp variable for loop
	long addwindow;
	double temp;
	Eimg_PixelRect *pixeltmp;
	char *tempstring;


	addwindow = md->windowsize - 1 ;


	

//**************************************************************************************************************************************
	//将图像复制到pixeltmp,补齐边界
	pixeltmp = eimg_PixelRectCreate(md->npixx + addwindow, md->npixy + addwindow, EGDA_TYPE_F64, &err);
	if(!pixeltmp || err)return 1;

	for(i = addwindow / 2 ; i < md->npixy + addwindow / 2; i++)
	{
		for(j = addwindow / 2 ; j < md->npixx + addwindow / 2; j++)
		{
			temp = *((Egda_F64 *)(EIMG_DATAPTR(md->pixelblock, j - addwindow / 2, i - addwindow / 2 )));
			*((Egda_F64 *)(EIMG_DATAPTR(pixeltmp, j, i))) = temp;
		}
	}

	//补边界
	for(i=0;i<addwindow/2;i++)//上
	{
		for(j = addwindow / 2 ; j < md->npixx+addwindow / 2; j++)
		{
			temp = *((Egda_F64 *)(EIMG_DATAPTR(md->pixelblock, j - addwindow / 2, 0)));
			*((Egda_F64 *)(EIMG_DATAPTR(pixeltmp, j, i))) = temp;
		}
	}

	for(i=md->npixy+addwindow/2;i<md->npixy+addwindow;i++)//下
	{
		for(j = addwindow / 2 ; j < md->npixx+addwindow / 2; j++)
		{
			temp = *((Egda_F64 *)(EIMG_DATAPTR(md->pixelblock, j - addwindow / 2, md->npixy-1)));
			*((Egda_F64 *)(EIMG_DATAPTR(pixeltmp, j, i))) = temp;
		}
	}

	for(i = addwindow / 2;i < md->npixy+addwindow / 2 ; i++)//左
	{
		for(j = 0 ; j < addwindow / 2; j++)
		{
			temp = *((Egda_F64 *)(EIMG_DATAPTR(md->pixelblock, 0, i-addwindow / 2)));
			*((Egda_F64 *)(EIMG_DATAPTR(pixeltmp, j, i))) = temp;
		}
	}

	for(i = addwindow / 2;i < md->npixy+addwindow / 2 ; i++)//右
	{
		for(j = md->npixx + addwindow / 2; j < md->npixx + addwindow ; j++)
		{
			temp = *((Egda_F64 *)(EIMG_DATAPTR(md->pixelblock, md->npixx-1, i-addwindow / 2)));
			*((Egda_F64 *)(EIMG_DATAPTR(pixeltmp, j, i))) = temp;
		}
	}

	//4个角落
	for(i=0;i<addwindow / 2;i++)//左上
	{
		for(j=0;j<addwindow/2;j++)
		{
			temp = *((Egda_F64 *)(EIMG_DATAPTR(md->pixelblock, 0, 0)));
			*((Egda_F64 *)(EIMG_DATAPTR(pixeltmp, j, i))) = temp;
		}
	}

	for(i=md->npixy+addwindow/2;i<md->npixy+addwindow;i++)//左下
	{
		for(j=0;j<addwindow/2;j++)
		{
			temp = *((Egda_F64 *)(EIMG_DATAPTR(md->pixelblock, 0, md->npixy-1)));
			*((Egda_F64 *)(EIMG_DATAPTR(pixeltmp, j, i))) = temp;
		}
	}

	for(i=0;i<addwindow/2;i++)//右上
	{
		for(j=md->npixx+addwindow/2;j<md->npixx+addwindow;j++)
		{
			temp = *((Egda_F64 *)(EIMG_DATAPTR(md->pixelblock, md->npixx-1, 0)));
			*((Egda_F64 *)(EIMG_DATAPTR(pixeltmp, j, i))) = temp;
		}
	}

	for(i=md->npixy+addwindow/2;i<md->npixy+addwindow;i++)//右下
	{
		for(j=md->npixx+addwindow/2;j<md->npixx+addwindow;j++)
		{
			temp = *((Egda_F64 *)(EIMG_DATAPTR(md->pixelblock, md->npixx-1, md->npixy-1)));
			*((Egda_F64 *)(EIMG_DATAPTR(pixeltmp, j, i))) = temp;
		}
	}


	

//**************************************************************************************************************************************
//**************************************************************************************************************************************

	//自己的滤波器

	//在pixeltemp中进行滤波处理,将滤波后的像素点写入md->pixelblock 中.
	AdaptiveFilter(md,pixeltmp,md->pixelblock);


//**************************************************************************************************************************************
//**************************************************************************************************************************************

	
	//清理现场
	emet_MeterInfoPrint(meterinfo,50,&err);
	

	eimg_LayerChangeOptions(md->outlayer,&err,EIMG_LAYER_OPTION_PROGRESS_METER_ON_CLOSE,meterinfo,EIMG_LAYER_OPTION_END);

	tempstring = estr_Sprintf(NULL,"Stage 3 of 3 : Computing statistics ... ...",&err);
	emet_MeterInfoChangeTaskName(meterinfo,tempstring,&err);

	emet_MeterInfoSet(meterinfo,0,100,&err);
	emsc_Free(tempstring);

	eimg_PixelRectDelete(pixeltmp,&err);
	if (err) 
	{
		return 100;
    }

	return 0;
}


static void
AdaptiveFilter(md,src,tar)
MyData *md;
Eimg_PixelRect *src;
Eimg_PixelRect *tar;
{
	Eerr_ErrorReport *err=NULL;

	long i,j;
	long addwindow;
	double Xmean,Dmean;
	double dT,dTn;
	double minT,maxT;
	double p;
	double k;
	double X0;
	Eimg_PixelRect *pixTemp;
	char *mystring;


	addwindow = md->windowsize - 1;
	minT = 999999;
	maxT = -999999;
	pixTemp = eimg_PixelRectCreate(md->npixx, md->npixy, EGDA_TYPE_F64, &err );

	emet_MeterInfoPrint(meterinfo,0.0,&err);

	for(i=addwindow/2;i<md->npixy+addwindow/2;i++)
	{
		mystring = estr_Sprintf(NULL,"Stage 1 of 3 : Computing pixels : line %d of %d",&err,i+1,md->npixy);
		emet_MeterInfoChangeTaskName(meterinfo,mystring,&err);
		emet_MeterInfoSet(meterinfo,(100.0*i)/md->npixy,(100.0*(i+0.5))/md->npixy,&err);
		emet_MeterInfoPrint(meterinfo,(100.0*i)/md->npixy,&err);

		for(j=addwindow/2;j<md->npixx+addwindow/2;j++)
		{
			//得到以j,i点为中心点的窗口内的均值和方差
			GetXmeanDmean(src,j,i,md->windowsize,&Xmean,&Dmean);

			//将该点的均值Xmean保存在pixTemp中
			*((Egda_F64 *)(EIMG_DATAPTR(pixTemp, j-addwindow/2, i-addwindow/2))) = Xmean;


			//得到以j,i点为中心点的窗口内的T统计量
			dT = GetT((*((Egda_F64 *)(EIMG_DATAPTR(tar, j-addwindow/2, i-addwindow/2)))),Xmean,Dmean);

			//得到最大和最小的T统计量
			if(dT > maxT)
			{
				maxT = dT;
			}
			if(dT < minT)
			{
				minT = dT;
			}


			//将T统计量写入目标区tar.
			*((Egda_F64 *)(EIMG_DATAPTR(tar, j-addwindow/2, i-addwindow/2))) = dT;
		}
	}

	emet_MeterInfoPrint(meterinfo,100.0,&err);

	emet_MeterInfoPrint(meterinfo,0.0,&err);

	//归一化目标区tar,并将最终结果x=k*x0 + (1-k)*Xmean写入目标区tar
	for(i=0;i<md->npixy;i++)
	{
		mystring = estr_Sprintf(NULL,"Stage 2 of 3 : Writing pixels to line : %d of %d",&err,i+1,md->npixy);
		emet_MeterInfoChangeTaskName(meterinfo,mystring,&err);
		emet_MeterInfoSet(meterinfo,(100.0*i)/md->npixy,(100.0*(i+0.5))/md->npixy,&err);
		emet_MeterInfoPrint(meterinfo,(100.0*i)/md->npixy,&err);

		for(j=0;j<md->npixx;j++)
		{
			dT = *((Egda_F64 *)(EIMG_DATAPTR(tar, j, i)));
			dTn = GetTn(dT,minT,maxT);//得到j,i点的归一化T值
			p = 1 - dTn;

			//得到k的值
			k = getK(md->a,md->b,p);

			//利用公式计算新的x值
			X0 = *((Egda_F64 *)(EIMG_DATAPTR(src, j+addwindow/2, i+addwindow/2)));
			
			dT = k * X0 + (1-k) * (*((Egda_F64 *)(EIMG_DATAPTR(pixTemp, j, i))));

			//写入结果
			*((Egda_F64 *)(EIMG_DATAPTR(tar, j, i))) = dT;

		}
	}

	emet_MeterInfoPrint(meterinfo,100.0,&err);

	emsc_Free(mystring);
	eimg_PixelRectDelete(pixTemp,&err);
	return;
}



//得到一个窗口内的均值和标准差
static void 
GetXmeanDmean(src, x,y,windowsize,Xmean,Dmean)
Eimg_PixelRect *src;//图像的范围
long x;
long y;//窗口中心点的坐标
int windowsize;//窗口的大小
double *Xmean;//是返回值:均值
double *Dmean;//是返回值:方差
{
	int i,j;
	double Xall,Dall;

	Xall = 0;
	Dall = 0;

	for(i=x-(windowsize-1)/2;i<x+(windowsize-1)/2;i++)
	{
		for(j=y-(windowsize-1)/2;j<y+(windowsize-1)/2;j++)
		{
			Xall = Xall + (*((Egda_F64 *)(EIMG_DATAPTR(src, i, j))));
		}
	}

	*Xmean = Xall / (windowsize * windowsize);

	for(i=x-(windowsize-1)/2;i<x+(windowsize-1)/2;i++)
	{
		for(j=y-(windowsize-1)/2;j<y+(windowsize-1)/2;j++)
		{
			Dall = Dall + ((*((Egda_F64 *)(EIMG_DATAPTR(src, i, j))) - *Xmean)*(*((Egda_F64 *)(EIMG_DATAPTR(src, i, j))) - *Xmean));
		}
	}

	*Dmean = sqrt(Dall / (windowsize * windowsize));


	return;
}



//得到窗口内的T统计量
static double
GetT(X0,Xmean,Dmean)
double X0;
double Xmean;
double Dmean;
{
	double ret;

	ret = 0;
	
	ret = (X0 - Xmean)/Dmean;

	return ret;
}


//得到归一化的Tn
static double 
GetTn(dT,minT,maxT)
double dT;
double minT;
double maxT;
{
	double ret;


	ret = 0;

	dT = fabs(dT);
	//minT = fabs(minT);
	//maxT = fabs(maxT);

	ret = (dT - minT) / (maxT - minT);


	return ret;
}


static double 
getK(a,b,p)
double a;
double b;
double p;
{
	double ret;


	ret = 0;
	
	if(p>=0 && p<=a)
	{
		ret = 0;
	}
	else if(p>a && p<b)
	{
		ret = (p-a)/(b-a);
	}
	else if(p>=b && p<1)
	{
		ret = 1;
	}
	else
	{
		ret = -999;
	}


	return ret;
}