comput.cpp

基于小波的SAR斑点处理

void Qsm32(double *ResX,
			  double *ResY,
			  double *Src,
			  int VarNumber,
			  double *ResB)
{
 int t,time1,time2;
 int a_i,b_i;
 int i,j,k,l;
 double *a,*b,*r;
 double *r2;
 //float temp;
 //float s1,s2,s3;
 //float sum_z;
 t = 2;

 time1=add(t+1);
 time2=add(2*t+1);
 a = new double[time2];
 b = new double[time1];
 r =  new double[time1*time1];
 r2 = new double[time1*time1];
 a_i=0;
 for(i=0;i<=2*t;i++)
  {
	for(j=0;j<=i;j++)
	 {
	 *(a+a_i)=0;
	 for(k=0;k<VarNumber;k++)
	  *(a+a_i)+=new_pow(*(ResX+k),*(ResY+k),j,i-j);
	 a_i++;
	 }
  }

 b_i=0;
 for(i=0;i<=t;i++) /*Z*/
  {
	for(j=0;j<=i;j++)
	 {
	 *(b+b_i)=0;
	 for(k=0;k<VarNumber;k++)
	  *(b+b_i)+=new_pow(*(ResX+k),*(ResY+k),j,i-j)*(*(Src+k));
	 b_i++;
	 }
  }
 a_i=0;
 for(i=0;i<=t;i++)
  for(j=0;j<=i;j++)
	{
	 b_i=0;
	 for(k=i;k<=i+t;k++)
	  for(l=0;l<=k-i;l++)
	  {
		*(ResB+b_i)=*(a+k*(k+1)/2+j+l);
		b_i++;
		}
	 assign(ResB,r+a_i*time1,1,time1);
	 a_i++;
	 }
  inv(r,r2,time1);
  mul(r2,b,ResB,time1,time1,1);
  delete a;
  delete b;
  delete r;
  delete r2;
}

int TrendSurface2(double *x,
						 double *y,
						 double *z,
						 int VarNum,
						 int Times,
						 double *bb)
{
	int i;
	int j;
	int k;
	int l;

	int Time1 = add(Times+1);
	int Time2 = add(2*Times+1);
	double *a = new double[Time2];
	double *b = new double[Time1];
	double *r = new double[Time1*Time1];
	double *r2 = new double[Time1*Time1];
	int a_i;
	int b_i;

	a_i=0;
	for(i=0;i<=2*Times;i++)
	{
		for(j=0;j<=i;j++)
		{
			*(a+a_i)=0;
			for(k=0;k<VarNum;k++)
			*(a+a_i)+=new_pow(*(x+k),*(y+k),j,i-j);
			a_i++;
		}
	}

	b_i=0;
	for(i=0;i<=Times;i++) /*Z*/
	{
		for(j=0;j<=i;j++)
		{
			*(b+b_i)=0;
			for(k=0;k<VarNum;k++)
			*(b+b_i)+=new_pow(*(x+k),*(y+k),j,i-j)*(*(z+k));
			b_i++;
		}
	}

	a_i=0;
	for(i=0;i<=Times;i++)
		for(j=0;j<=i;j++)
		{
			b_i=0;
			for(k=i;k<=i+Times;k++)
				for(l=0;l<=k-i;l++)
				{
					*(bb+b_i)=*(a+k*(k+1)/2+j+l);
					b_i++;
				}
			assign(bb,r+a_i*Time1,1,Time1);
			a_i++;
		}
	if(inv(r,r2,Time1)==-1)
		return -1;
	mul(r2,b,bb,Time1,Time1,1);
	delete []a;
	delete []b;
	delete []r;
	delete []r2;
	return 0;
}

double GetMean(float * fp,long Size)
{
	double Sum=0;
	for(long order=0;order<Size;order++)
		Sum+=(double)fp[order];
	return (double)(Sum/Size);
}

double GetStdE(float *fp,double Mean,long Size)
{
	double Sum=0;
	for(long order=0;order<Size;order++)
		Sum+=(fp[order]-Mean)*(fp[order]-Mean);
	return (double)sqrt(Sum/(Size-1));
}

double GetMean(BYTE * fp,long Size)
{
	double Sum=0;
	for(long order=0;order<Size;order++)
		Sum+=(double)fp[order];
	return (double)(Sum/Size);
}

double GetStdE(BYTE *fp,double Mean,long Size)
{
	double Sum=0;
	for(long order=0;order<Size;order++)
		Sum+=(fp[order]-Mean)*(fp[order]-Mean);
	return (double)sqrt(Sum/(Size-1));
}

//FFT
// compute x to the power of y
int powerof2(int x, int y)
{
	int z,i;
	z=1;
	i=0;
	while( z<y )
	{
		z=z*x;
		i++;
	}
	return i;
}

// One-dimensional fast Fourier transform
void FFT1D(float xr[],float xi[],int LENGTH,int inve,int POWER)
{
	int n,nv2,nm1,k,i,ii,j,ie,iei,ip;
	float treal,timag,ur,ui,wr,wi,yr,yi,pi=3.141593F;

	n=(int)pow((double)2,(double)POWER);
	nv2=n/2;
	nm1=n-1;
	j=0;
	for (i=0;i<nm1;++i) 
	{
		if (i<=j) 
		{
			treal=xr[j]; timag=xi[j];
			xr[j]=xr[i]; xi[j]=xi[i];
			xr[i]=treal; xi[i]=timag;
		}
		k=nv2;
		while (k<=j) 
		{
			j -= k;
			k /=2;
		}
		j += k;
	}

	for (ii=1;ii<=POWER;ii++) 
	{
		ie=(int)pow((double)2,(double)ii);
		iei=ie/2;
		ur=1.0F; ui=0.0F;
		wr=(float)cos(pi/iei);
		if (inve==1)
			wi=(float)sin(pi/iei);
		else
			wi=-(float)sin(pi/iei);
		for (j=0;j<iei;++j) 
		{
			for (i=j;i<LENGTH;i+=ie)
			{
				ip=i+iei;
				treal=xr[ip]*ur-xi[ip]*ui;
				timag=xi[ip]*ur+xr[ip]*ui;
				xr[ip]=xr[i]-treal;
				xi[ip]=xi[i]-timag;
				xr[i] +=treal;
				xi[i] +=timag;
			}
			yr=ur*wr-ui*wi;
			yi=ui*wr+ur*wi;
			ur=yr;
			ui=yi;
		}
	}
	return;
} // FFT1D
                                                      
// two-dimensional Fourier transform and highpass filtering
void FFT2D( float *lProfilereal,float *lProfileimag,int LENGTH,int inve,int POWER)
{
	float *lpreal;
	float *lpimag;

	float *xreal = new float[LENGTH];
	ASSERT(xreal != NULL);
	float *ximag = new float[LENGTH];		// used for the FFT
	ASSERT(ximag != NULL);
	
	int i,j;
	float fLENGTH;
	DWORD lLENGTH;

	lLENGTH = LENGTH;
	fLENGTH = (float)LENGTH*LENGTH;
	//POWER = powerof2( 2, LENGTH );
	//MAXF = (int)ceil((double)( 1.414F * LENGTH * 0.5F ));
	for (j=0;j<LENGTH;j++)
	{
		for (i=0;i<LENGTH;i++)
		{
			xreal[i]=*( lProfilereal+(DWORD)LENGTH*i+j );
			ximag[i]=*( lProfileimag+(DWORD)LENGTH*i+j );
		}
		FFT1D(xreal, ximag, LENGTH, inve,POWER);
		for (i=0;i<LENGTH;i++)
		{
			*( lProfilereal+(DWORD)LENGTH*i+j )=xreal[i];
			*( lProfileimag+(DWORD)LENGTH*i+j )=ximag[i];
		}
	}

	for (i=0;i<LENGTH;i++)
	{
		lpreal = lProfilereal + i*lLENGTH;
		lpimag = lProfileimag + i*lLENGTH;
		for (j=0;j<LENGTH;j++, lpreal++, lpimag++)
		{
			xreal[j]=*lpreal;
			ximag[j]=*lpimag;
		}
		FFT1D( xreal, ximag, LENGTH, inve,POWER);
		lpreal = lProfilereal + i*lLENGTH;
		lpimag = lProfileimag + i*lLENGTH;
		for (j=0;j<LENGTH;j++, lpreal++, lpimag++)
		{
			*lpreal =xreal[j];
			*lpimag =ximag[j];
		}
	}

	if ( inve == 0 )
	{
		lpreal = lProfilereal;
		lpimag = lProfileimag;
		for (i=0;i<LENGTH;i++)
		{
			for (j=0;j<LENGTH;j++, lpreal++, lpimag++)
			{
				*lpreal /= fLENGTH;
				*lpimag /= fLENGTH;
			}
		}
	}
	delete []xreal;
	delete []ximag;

	return;
} // FFT2D

void DrawRect(CDC * pDC, POINT * pnt1, POINT * pnt2)
{
	int minX,maxX,minY,maxY;
	if(pnt1->x < pnt2->x)
	{
		minX = pnt1->x;
		maxX = pnt2->x;
	}
	else
	{
		minX = pnt2->x;
		maxX = pnt1->x;
	}
	if(pnt1->y < pnt2->y)
	{
		minY = pnt1->y;
		maxY = pnt2->y;
	}
	else
	{
		minY = pnt2->y;
		maxY = pnt1->y;
	}
	pDC->MoveTo(minX,minY);
	pDC->LineTo(minX,maxY);
	
	//pDC->MoveTo(point2.x,point1.y);
	pDC->LineTo(maxX,maxY);

	//pDC->MoveTo(point2);
	pDC->LineTo(maxX,minY);

	//pDC->MoveTo(point1.x,point2.y);
	pDC->LineTo(minX,minY);
}

int intAbs(int nTemp)
{
	if(nTemp > 0)
		return nTemp;
	else
		return -nTemp;
}

float fAbs(float fTemp)
{
	if(fTemp>0)
		return fTemp;
	else
		return -fTemp;
}

double doubleAbs(double dTemp)
{
	if(dTemp >0)
		return dTemp;
	else
		return -dTemp;
}

int SymFold(int x,int xSize)
{
	if(x < 0)
		x = -x;
	if(x >= xSize)
		x = 2*xSize-x-2;

	return x;
}

void Sort(BYTE * data, int size)
{	// 采用简单的冒泡排序
	BYTE temp;
	for(int flag=0; flag<size-1; flag++)
		for(int i=flag; i<size-1; i++)
			if(data[i]>data[i+1])
			{
				temp = data[i];
				data[i] = data[i+1];
				data[i+1] = temp;
			}
}

float Mean(BYTE * a,int m)
{
	double sum = 0;
	for(int i=0; i<m; i++)
		sum = sum + a[i];
	return float(sum/m);
}

float STD(BYTE * a,int m,float fMean)
{
	if(fMean != 65535)
		fMean = Mean(a,m);
	double v = 0;
	for(int i=0; i<m; i++)
		v = v + (a[i]-fMean)*(a[i]-fMean);
	v = sqrt(v/m);

	return float(v);
}

// 融合规则一
// 采用能量作为"显著性因子",对某尺度某方向的小波系数整体作选择估计
// 结果保存在 "pATemp" 中;
void WholeSelect(float * pETemp,float * pATemp,int width,int height)
{
	int i,j;
	double eSum,aSum;
	float eTemp,aTemp;

	// 计算测度
	eSum = 0;
	aSum = 0;
	for(i=0; i<height; i++)
		for(j=0; j<width; j++)
		{
			eTemp = pETemp[i*width+j];
			aTemp = pATemp[i*width+j];
			
			eSum += eTemp * eTemp;
			aSum += aTemp * aTemp;
		}

	// 根据测度,选择系数
	if( eSum>aSum )
	{
		for(i=0; i<height; i++)
			for(j=0; j<width; j++)
				pATemp[i*width+j] = pETemp[i*width+j];
	}
}

// 融合规则二
// 高频融合:系数选择,结果保存在toP中
// 采用能量测度作为系数的估计依据,同时考虑了整体性
void CoeffSelect(float * fromP,float * toP,int width,int height,int size)
{
	//指示,系数取fromP置为T,取toP置为F;
	BOOL * pTemp = (BOOL *)new BOOL[width*height];
	BOOL * pT = (BOOL *)new BOOL[width*height];
	ASSERT(pTemp != NULL && pT != NULL);

	int i,j,m,n,num,length;
	length = size*size;
	BOOL * L = (BOOL *)new BOOL[length];
	ASSERT( L != NULL);
	float temp;
	double fromEnergy,toEnergy;

	for(i=0; i<height; i++)
		for(j=0; j<width; j++)
		{
			// 计算能量测度
			fromEnergy = 0;
			toEnergy = 0;
			for(m=-size/2; m<(size+1)/2; m++)
				for(n=-size/2; n<(size+1)/2; n++)
				{
					temp = fromP[SymFold(i+m,height)*width+SymFold(j+n,width)];
					fromEnergy += temp*temp;
					temp = toP[SymFold(i+m,height)*width+SymFold(j+n,width)];
					toEnergy += temp*temp;
				}
			
			// 判断中心点的选取
			if( fromEnergy > toEnergy )
				pTemp[i*width+j] = TRUE;
			else
				pTemp[i*width+j] = FALSE;
		}
	
	// "众数一致性"调整
	for(i=0; i<height; i++)
		for(j=0; j<width; j++)
		{
			// 取指示
			for(m=-size/2; m<(size+1)/2; m++)
				for(n=-size/2; n<(size+1)/2; n++)
					L[(m+size/2)*size+n+size/2] = pTemp[SymFold(i+m,height)*width+SymFold(j+n,width)];

			// 计算中心周围取法与自己不同的数目
			num = 0;
			for(m=0; m<size; m++)
				for(n=0; n<size; n++)
				{
					if(L[m*size+n] != L[size/2*size+size/2])
						num++;
				}

			// 如该数目>=6,变更中心的取法
			if(num>=6)
				pT[i*width+j] = !pTemp[i*width+j];
			else
				pT[i*width+j] = pTemp[i*width+j];
		}

	// 根据指示确定系数的选择
	for(i=0; i<height; i++)
		for(j=0; j<width; j++)
		{
			if(pT[i*width+j] == TRUE)
				toP[i*width+j] = fromP[i*width+j];
		}
	
	delete pTemp;
	delete pT;
	delete L;
}

// 融合规则三
// 高频融合,结果保存在pA中
// 采用 "相似因子" 和 "显著因子" 综合判断
void SelectCoeff(float * pE,float * pA,int width,int height,float threshold,int size)
{
	// 申请内存,存放融合结果
	float * pTemp = (float *)new float[width*height];
	ASSERT(pTemp != NULL);

	int i,j,m,n;
	float eTemp,aTemp;
	double eEnergy,aEnergy,cEnergy,match;
	float eWeight,aWeight;

	for(i=0; i<height; i++)
		for(j=0; j<width; j++)
		{
			// 初始化
			eEnergy = 0;
			aEnergy = 0;
			cEnergy = 0;
			match = 0;

			// 计算各个模板内的能量
			for(m=-size/2; m<(size+1)/2; m++)
				for(n=-size/2; n<(size+1)/2; n++)
				{
					eTemp = pE[SymFold(i+m,height)*width+SymFold(j+n,width)];
					aTemp = pA[SymFold(i+m,height)*width+SymFold(j+n,width)];

					eEnergy += eTemp * eTemp;								
					aEnergy += aTemp * cEnergy;
					cEnergy += eTemp * aTemp;
				}

			// 计算匹配度
			match = 2*cEnergy/(eEnergy+aEnergy);

			// 判断中心点的选取
			if(match<threshold)
			{
				if(eEnergy>aEnergy)
					pTemp[i*width+j] = pE[i*width+j];
				else
					pTemp[i*width+j] = pA[i*width+j];
			}
			else
			{
				if(eEnergy>aEnergy)
					eWeight = float( 1+(1-match)/match )/2;
				else
					eWeight = float( 1-(1-match)/match )/2;
				aWeight = 1-eWeight;

				pTemp[i*width+j] = eWeight*pE[i*width+j] + aWeight*pA[i*width+j];
			}
		}
	
	// 转存结果
	for(i=0; i<height; i++)
		for(j=0; j<width; j++)
			pA[i*width+j] = pTemp[i*width+j];

	delete pTemp;
}