comput.cpp

基于小波的SAR斑点处理

#include "stdafx.h"
#include <math.h>
#include "Comput.h"

#define EPS 0.000001
#define PI 3.1415926

// 分布表
static double F_DISTR[31][9] ={
	{161.4,199.5,215.7,224.6,230.2,234.0,236.8,238.9,243.9},
	{18.5, 19.0, 19.2, 19.2, 19.3, 19.3, 19.4, 19.4, 19.4},
	{10.1, 9.55, 9.28, 9.12, 9.01, 8.94, 8.89, 8.85, 8.74},
	{7.71, 6.94, 6.59, 6.39, 6.26, 6.16, 6.09, 6.04, 5.91},
	{6.61, 5.79, 5.41, 5.19, 5.05, 4.95, 4.88, 4.82, 4.68},
	{5.99, 5.14, 4.76, 4.53, 4.39, 4.28, 4.21, 4.15, 4.00},
	{5.59, 4.74, 4.35, 4.12, 3.97, 3.87, 3.79, 3.73, 3.57},
	{5.32, 4.46, 4.07, 3.84, 3.69, 3.58, 3.50, 3.44, 3.28},
	{5.12, 4.26, 3.86, 3.63, 3.48, 3.37, 3.29, 3.23, 3.04},
	{4.96, 4.10, 3.71, 3.48, 3.33, 3.22, 3.14, 3.07, 2.91},
	{4.84, 3.98, 3.59, 3.36, 3.20, 3.09, 3.01, 2.95, 2.79},
	{4.75, 3.89, 3.49, 3.26, 3.11, 3.00, 2.91, 2.85, 2.69},
	{4.67, 3.81, 3.41, 3.18, 3.03, 2.92, 2.83, 2.77, 2.60},
	{4.60, 3.74, 3.34, 3.11, 2.96, 2.85, 2.76, 2.70, 2.53},
	{4.54, 3.68, 3.29, 3.06, 2.90, 2.79, 2.71, 2.64, 2.48},
	{4.49, 3.63, 3.24, 3.01, 2.85, 2.74, 2.66, 2.59, 2.42},
	{4.45, 3.59, 3.20, 2.96, 2.81, 2.70, 2.61, 2.55, 2.38},
	{4.41, 3.55, 3.16, 2.93, 2.77, 2.66, 2.58, 2.51, 2.34},
	{4.38, 4.52, 3.13, 2.90, 2.74, 2.63, 2.54, 2.48, 2.31},
	{4.35, 3.49, 3.10, 2.87, 2.71, 2.60, 2.51, 2.45, 2.28},
	{4.32, 3.47, 3.07, 2.84, 2.68, 2.57, 2.49, 2.42, 2.25},
	{4.30, 3.44, 3.05, 2.82, 2.66, 2.55, 2.46, 2.40, 2.34},
	{4.28, 3.42, 3.03, 2.80, 2.64, 2.53, 2.44, 2.37, 2.20},
	{4.26, 3.40, 3.01, 2.78, 2.62, 2.51, 2.42, 2.36, 2.29},
	{4.24, 3.39, 2.99, 2.76, 2.60, 2.49, 2.40, 2.34, 2.16},
	{4.23, 3.37, 2.98, 2.74, 2.59, 2.47, 2.39, 2.32, 2.15},
	{4.21, 3.35, 2.96, 2.73, 2.57, 2.46, 2.37, 2.31, 2.13},
	{4.20, 3.34, 2.95, 2.71, 2.56, 2.43, 2.36, 2.29, 2.12},
	{4.18, 3.33, 2.93, 2.70, 2.55, 2.42, 2.35, 2.28, 2.10},
	{4.17, 3.32, 2.92, 2.69, 2.53, 2.34, 2.33, 2.27, 2.09},
	{4.08, 3.23, 2.84, 2.61, 2.45, 2.25, 2.25, 2.18, 2.00}
};

// 查找分布表
double look_up(int n,int m)
{
 if(n>8)
	n = 9;
 if(m>30)
	m = 31;
 return (double)F_DISTR[m-1][n-1];
}

// 查找最大的双精度数，并返回其在数组中的位置
int max_i(double *x,int n)
{
	int i;
	int l=0;
	for(i=1; i<n; i++)
	{
		if(*(x+l)<*(x+i))
		l=i;
	}
	return l;
}

// 查找最大的双精度数，并返回其值
double  max_f(double *x,int n)
{
	return *(x+max_i(x,n));
}

// 查找最小的双精度数，并返回其在数组中的位置
int  min_i(double *x,int n)
{
	int i;
	int l=0;
	for(i=1;i<n;i++)
	{
		if(*(x+l)>*(x+i))
		l=i;
	}
	return l;
}

// 查找最小的双精度数，并返回其值
double  min_f(double *x,int n)
{
	return *(x+min_i(x,n));
}

// 求某双精度数的平方
double  sqr(double x)
{
	return x*x;
}

// 求1-n的自然数列的和
int  add(int n)
{
	return n*(n+1)/2;
}

// 求均值 
BYTE  mean(BYTE *p,int m)
{
	long int s=0;
	int i;
	for(i=0;i<m;i++)
		s+=*(p+i);
	return (BYTE)(s/m);
}

// 求中值
BYTE  media(BYTE *data,int num)
{
	BYTE temp;
	int i,j;
	for(i=0;i<num;i++)
	for(j=i;j<num;j++)
		if(*(data+j)>*(data+i))
		{
			temp=*(data+j);
			*(data+j)=*(data+i);
			*(data+i)=temp;
		}
	if(num%2 == 1)
		return data[num/2];
	else
		return (BYTE)(((double)data[num/2]+(double)data[num/2-1])/2);
}

// 求中值
float  mediaF(float *data,int num)
{
	float temp;
	int i,j;
	for(i=0;i<num;i++)
	for(j=i;j<num;j++)
		if(*(data+j)>*(data+i))
		{
			temp=*(data+j);
			*(data+j)=*(data+i);
			*(data+i)=temp;
		}
	if(num%2 == 1)
		return data[num/2];
	else
		return (float)(((double)data[num/2]+(double)data[num/2-1])/2);
}


// 求众数
BYTE  many(BYTE *data,int num)
{
	int t,w;
	BYTE md,oldmode;
	int count,oldcount;
	oldmode=0;
	oldcount=0;
	for(t=0; t<num; t++)
	{
		md=data[t];
		count=1;
		for(w=t+1; w<num; w++)
			if(md == data[w])
				count++;
		if(count>oldcount)
		{
			oldmode=md;
			oldcount=count;
		}
	}
	return oldmode;
}

// 数据块转赋值
void  assign(double *x1,double *x2,int n,int m)
{
	int i,j;
	for(i=0; i<n; i++)
	{
		for(j=0; j<m; j++)
			x2[i*m+j] = x1[i*m+j];
	}
}

// 返回数的符号：1—正数，0—零，-1—负数
double sign(double x)
{
	if(x>0)
		return 1.0;
	else if(x<0)
		return -1.0;
	else
		return 0.0;
}

// 求双精度数组的和
double sum(float *x,int n)
{
	double s=0;
	int i;
	for(i=0;i<n;i++)
		s+=(*(x+i));
	return s;
}

// 求协方差值*n
double l(float *x1,float *x2,int n)
{
	double medx1,medx2;
	int i;
	double s=0;
	medx1=sum(x1,n)/n;
	medx2=sum(x2,n)/n;
	for(i=0;i<n;i++)
		s+=(*(x1+i)-medx1)*(*(x2+i)-medx2);
	return s;
}

// 求相关系数
double relate(float *x1,float * x2,int n)
{
	return l(x1,x2,n)/sqrt(l(x1,x1,n)*l(x2,x2,n));
}

// 求双精度数组的和
double sum(double *x,int n)
{
	double s=0;
	int i;
	for(i=0;i<n;i++)
		s+=(*(x+i));
	return s;
}

// 求协方差值*n
double l(double *x1,double *x2,int n)
{
	double medx1,medx2;
	int i;
	double s=0;
	medx1=sum(x1,n)/n;
	medx2=sum(x2,n)/n;
	for(i=0;i<n;i++)
		s+=(*(x1+i)-medx1)*(*(x2+i)-medx2);
	return s;
}

// 求相关系数
double relate(double *x1,double * x2,int n)
{
	return l(x1,x2,n)/sqrt(l(x1,x1,n)*l(x2,x2,n));
}

// 求双精度数组的和
double sum(BYTE *x,int n)
{
	double s=0;
	int i;
	for(i=0;i<n;i++)
		s+=(*(x+i));
	return s;
}

// 求协方差值*n
double l(BYTE *x1,BYTE *x2,int n)
{
	double medx1,medx2;
	int i;
	double s=0;
	medx1=sum(x1,n)/n;
	medx2=sum(x2,n)/n;
	for(i=0;i<n;i++)
		s+=(*(x1+i)-medx1)*(*(x2+i)-medx2);
	return s;
}

// 求相关系数
double relate(BYTE *x1,BYTE * x2,int n)
{
	return l(x1,x2,n)/sqrt(l(x1,x1,n)*l(x2,x2,n));
}

// 转赋数据块的值
void trn(double *x1,double *x2,int n1,int m1)
{
	int t1,t2;
	for(t1=0;t1<m1;t1++)
	{
		for(t2=0;t2<n1;t2++)
		{
			x2[t1*n1+t2] = x1[t2*m1+t1];
		}
	}
}

void mul(double *x1,double *y1,double *result,int n1,int i,int m1)
{
	int t1,t2,t3;
	double s;
	for(t1=0;t1<n1;t1++)
	for(t2=0;t2<m1;t2++)
	{
		s=0;
		for(t3=0;t3<i;t3++)
		{
			s=s+x1[t1*i+t3]*y1[t3*m1+t2];
		}
		result[t1*m1+t2]=s;
	}
}

void swap(double * x,int m,int i,int j) /*swamp two lines in the matrix*/
{
	int t;
	double temp;
	for(t=0;t<m;t++)
	{
		temp=x[i*m+t];
		x[i*m+t]=x[j*m+t];
		x[j*m+t]=temp;
	}
}

void div1(double *x,int m,int i,int j)
{
	int t;
	double a;
	a=x[i*m+j];
	for(t=0;t<m;t++)
		x[i*m+t]=x[i*m+t]/a;
}

void div(double *x,int m,int m1,int i,int j)
{
	 double a;
	 int t;
	 a=-*(x+j*m+m1)/(*(x+i*m+m1));
	 for(t=0;t<m;t++)
		x[j*m+t]=x[j*m+t]+x[i*m+t]*a;
}

// 求逆矩阵
BOOL inv(double *x,double *y,int n)
{
	int i,j;
	double *temp;
	int line;
	temp=new double[n*n*2];
	
	// 将原始矩阵和单位矩阵放在一起
	for(i=0;i<n;i++) /*y=E*/
	{
		for(j=0;j<n;j++)
			temp[2*i*n+j]=x[i*n+j];
		for(j=n;j<2*n;j++)
		{
			if(j-n==i)
				temp[2*n*i+j]=1;
			else
				temp[2*i*n+j]=0;
		}
	}

	// 使下三角矩阵为0
	for(i=0;i<n;i++)
	{
		line=i+1;
		while(fabs(*(temp+i*n*2+i)-0)<EPS)
		{
			if(line==n)
				return FALSE;
			swap(temp,2*n,i,line);
			line++;
		}
		for(j=line;j<n;j++)
			div(temp,2*n,i,i,j);
	}

	// 判断是否可逆
	if(fabs(*(temp+(n-1)*2*n+n-1)-0)<EPS)
		return FALSE;//Fail

	// 使上三角矩阵为0
	for(i=n-1;i>0;i--)
	{
		for(j=i-1;j>=0;j--)
			div(temp,2*n,i,i,j);
	}

	// 使对角线元素为1
	for(i=0;i<n;i++)
		div1(temp,2*n,i,i);
	
	// 转存结果
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			*(y+i*n+j)=*(temp+i*2*n+j+n);
	
	delete temp;
	return TRUE;
 }

double Row_Col_V(double *x,int n)
{
 double *r1=new double[n*n];
 int line;
 double s=1.0;
 assign(x,r1,n,n);
 for(int i=0;i<n;i++)
  {
	line=i+1;
	while(fabs(*(r1+i*n+i)-0)<EPS)
	 {
	 swap(r1,n,i,line);
	 line++;
	 if(line==n)
	  return 0.0;
	 }
	for(int j=line;j<n;j++)
	 div(r1,n,i,i,j);
  }
  for(i=0;i<n;i++)
	s*=r1[i*n+i];
  return s;
}

void yakbi1(double *x,double *x1,int i,int j,double thita,int n)
{
 int m;
 assign(x,x1,n,n);

 *(x1+i*n+i)=*(x+i*n+i)*sqr(cos(thita))+
		*(x+j*n+j)*sqr(sin(thita))+2*(*(x+i*n+j))*sin(thita)*cos(thita);
 *(x1+j*n+j)=*(x+i*n+i)*sqr(sin(thita))+
		*(x+j*n+j)*sqr(cos(thita))-2*(*(x+i*n+j))*sin(thita)*cos(thita);

 *(x1+i*n+j)=0;
 *(x1+j*n+i)=0;
 for(m=0;m<n;m++)
	  if(m!=i&&m!=j)
	 {
	  *(x1+i*n+m)=*(x+i*n+m)*cos(thita)+*(x+j*n+m)*sin(thita);
	  *(x1+i+m*n)=*(x1+i*n+m);
	  }
 for(m=0;m<n;m++)
	  if(m!=j&&m!=i)
	 {
	  *(x1+j*n+m)=-*(x+i*n+m)*sin(thita)+*(x+j*n+m)*cos(thita);
	  *(x1+j+m*n)=*(x1+j*n+m);
	 }
}

void Yakbi(double *CoMatrix,double *FMatrix,double *Lamda,int n)
{
 double max;
 double thita;
 double *r1,*r2,*t1,*t2,*t3;
 int time=1;
 int row_num;
 int col_num;
 r1=new double[n*n];
 r2=new double[n*n];
 t1=new double[n*n];
 t2=new double[n*n];
 t3=new double[n*n];
 assign(CoMatrix,r1,n,n);
 while(1)
 {
  max=*(r1+1);
  row_num=0;
  col_num=1;
  /*find the max number of the matrix*/
  for(int i=0;i<n;i++)
	for(int j=0;j<n;j++)
	if(fabs(*(r1+i*n+j))>max&&i!=j)
	  {
		max=fabs(*(r1+i*n+j));
		row_num=i;
		col_num=j;
	  }
  if(max<EPS)
	 break;
  if(fabs(*(r1+row_num*(n+1))-*(r1+col_num*(n+1)))<EPS)
	 thita=PI/4;
  else
	 thita=atan(2*(*(r1+row_num*n+col_num))/
		 (*(r1+row_num*(n+1))-*(r1+col_num*(n+1))))/2;
  for(i=0;i<n;i++)
	for(int j=0;j<n;j++)
  {
	if(i==j)
	 *(t2+i*n+j)=1;
	else
	 *(t2+i*n+j)=0;
	}
  *(t2+row_num*(n+1))=cos(thita);
  *(t2+col_num*(n+1))=cos(thita);
  *(t2+row_num*n+col_num)=-sin(thita);
  *(t2+col_num*n+row_num)=sin(thita);

  if(time==1)
	assign(t2,t1,n,n);
  else
	{
	mul(t1,t2,t3,n,n,n);
	assign(t3,t1,n,n);
  }
  yakbi1(r1,r2,row_num,col_num,thita,n);
  assign(r2,r1,n,n);
  time++;
 }
 assign(t1,FMatrix,n,n);

 for(int i=0;i<n;i++)
  *(Lamda+i)=*(r1+i*(n+1));
 delete r1;
 delete r2;
 delete t1;
 delete t2;
 delete t3;
}

// 降序排序
void sortd(double *x,int n,int *m)
{
	int i,j,tmp;
	double temp;
	for(i=0;i<n;i++)
		*(m+i)=i;
	for(i=0;i<n;i++)
		for(j=i;j<n;j++)
			if(*(x+j)>*(x+i))
			{
				 temp=*(x+j);
				 *(x+j)=*(x+i);
				 *(x+i)=temp;
				 tmp=*(m+j);
				 *(m+j)=*(m+i);
				 *(m+i)=tmp;
			}
}

// 升序排序(冒泡)
void sortA(double *x,int n,int *m)
{
	int i,j,tmp;
	double temp;
	for(i=0;i<n;i++)
		*(m+i)=i;
	for(i=0;i<n;i++)
		for(j=i;j<n;j++)
			if(*(x+j)>*(x+i))
			{
				 temp=*(x+j);
				 *(x+j)=*(x+i);
				 *(x+i)=temp;
				 tmp=*(m+j);
				 *(m+j)=*(m+i);
				 *(m+i)=tmp;
			}
}

double new_pow(double x,double y,int i,int j)
{
 double s=1;
 int n;
 for(n=0;n<i;n++)
  s*=x;
 for(n=0;n<j;n++)
  s*=y;
 return s;
}
void qjqn(double *x1,double *x2,int n,int k)
{
 int t1,t2;
 for(t1=0;t1<n;t1++)
  for(t2=0;t2<n;t2++)
  {
	if(t1!=k&&t2!=k)
	 *(x2+t1*n+t2)=*(x1+t1*n+t2)-(*(x1+t1*n+k))*(*(x1+k*n+t2))/(*(x1+k*n+k));
	else if(t1!=k&&t2==k)
	 *(x2+t1*n+t2)=-(*(x1+t1*n+t2))/(*(x1+k*n+k));
	else if(t2!=k&&t1==k)
		 *(x2+t1*n+t2)=(*(x1+t1*n+t2))/(*(x1+k*n+k));
			else
			 *(x2+t1*n+t2)=1/(*(x1+t1*n+t2));
	}
}