bip.cpp

变化检测源程序

	    originimage[i][j*bonum+bons-1]=0;
	
	delete [] originimage;
	return TRUE;

}
//////求图像平均灰度值：图像数据，波段数，宽度，高度，处理波段
double CBip::MeanData(BYTE *lpbuf, int bonum, int w, int h, int bons)
{
	BYTE **originimage;
	int i,j;
	double mean;
	long sum;
	sum=0;
	originimage=new BYTE*[h];
	for(i=0;i<h;i++)
			originimage[i]=lpbuf+i*bonum*w;
	for(i=0;i<h;i++)
	for(j=0;j<w;j++)
			sum+=originimage[i][j*bonum+bons-1];
	mean=sum/w/h;
	delete [] originimage;
	return mean;

}





void CBip::DCT(double *f, double *F, int r)
{
	// 离散余弦变换点数
/*	LONG	count;
	
	// 循环变量
	int		i;
	
	// 中间变量
	double	dTemp;
	
	complex<double> *X;
	
	// 计算离散余弦变换点数
	count = 1<<r;
	
	// 分配内存
	X = new complex<double>[count*2];
	
	// 赋初值为0
	memset(X, 0, sizeof(complex<double>) * count * 2);
	
	// 将时域点写入数组X
	for(i=0;i<count;i++)
	{
		X[i] = complex<double> (f[i], 0);
	}
	
	// 调用快速付立叶变换
	FFT(X,X,r+1);
	
	// 调整系数
	dTemp = 1/sqrt(count);
	
	// 求F[0]
	F[0] = X[0].real() * dTemp;
	
	dTemp *= sqrt(2);
	
	// 求F[u]	
	for(i = 1; i < count; i++)
	{
		F[i]=(X[i].real() * cos(i*PI/(count*2)) + X[i].imag() * sin(i*PI/(count*2))) * dTemp;
	}
	
	// 释放内存
	delete X;*/

}

void CBip::IDCT(double *F, double *f, int r)
{
	// 离散余弦反变换点数
/*	LONG	count;
	
	// 循环变量
	int		i;
	
	// 中间变量
	double	dTemp, d0;
	
	complex<double> *X;
	
	// 计算离散余弦变换点数
	count = 1<<r;
	
	// 分配内存
	X = new complex<double>[count*2];
	
	// 赋初值为0
	memset(X, 0, sizeof(complex<double>) * count * 2);
	
	// 将频域变换后点写入数组X
	for(i=0;i<count;i++)
	{
		X[i] = complex<double> (F[i] * cos(i*PI/(count*2)), F[i] * sin(i*PI/(count*2)));
	}
	
	// 调用快速付立叶反变换
	IFFT(X,X,r+1);
	
	// 调整系数
	dTemp = sqrt(2.0/count);
	d0 = (sqrt(1.0/count) - dTemp) * F[0];
	
	// 计算f(x)
	for(i = 0; i < count; i++)
	{
		f[i] = d0 + X[i].real()* dTemp * 2 * count;
	}
	
	// 释放内存
	delete X;*/
}

void CBip::WALSH(double *f, double *F, int r)
{
	// 沃尔什-哈达玛变换点数
	LONG	count;
	
	// 循环变量
	int		i,j,k;
	
	// 中间变量
	int		bfsize,p;
	
	double *X1,*X2,*X;
	
	// 计算快速沃尔什变换点数
	count = 1 << r;
	
	// 分配运算所需的数组
	X1 = new double[count];
	X2 = new double[count];
	
	// 将时域点写入数组X1
	memcpy(X1, f, sizeof(double) * count);
	
	// 蝶形运算
	for(k = 0; k < r; k++)
	{
		for(j = 0; j < 1<<k; j++)
		{
			bfsize = 1 << (r-k);
			for(i = 0; i < bfsize / 2; i++)
			{
				p = j * bfsize;
				X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2];
				X2[i + p + bfsize / 2] = X1[i + p] - X1[i + p + bfsize / 2];
			}
		}
		
		// 互换X1和X2  
		X = X1;
		X1 = X2;
		X2 = X;
	}
	
	// 调整系数
	for(j = 0; j < count; j++)
	{
		p = 0;
		for(i = 0; i < r; i++)
		{
			if (j & (1<<i))
			{
				p += 1 << (r-i-1);
			}
		}

		F[j] = X1[p] / count;
	}
	
	// 释放内存
	delete X1;
	delete X2;

}

void CBip::IWALSH(double *F, double *f, int r)
{
	// 变换点数
	LONG	count;
	
	// 循环变量
	int		i;
	
	// 计算变换点数
	count = 1 << r;
	
	// 调用快速沃尔什－哈达玛变换进行反变换
	WALSH(F, f, r);
	
	// 调整系数
	for(i = 0; i < count; i++)
	{
		f[i] *= count;
	}

}

void CBip::DctBip(BYTE *lpbuf, int bonum, int width, int height, int bons)
{

}

void CBip::WalshBip(BYTE *lpbuf, int bonum, int width, int height, int bons)
{
	BYTE **originimage;
	int i,j;
	double dTemp;
	originimage=new BYTE*[height];
	for(i=0;i<height;i++)
			originimage[i]=lpbuf+i*bonum*width;
	// 进行付立叶变换的宽度和高度（2的整数次方）
	long	w;
	long	h;
	int		wp;
	int		hp;
	// 赋初值
	w = 1;
	h = 1;
	wp = 0;
	hp = 0;
	
	// 计算进行付立叶变换的宽度和高度（2的整数次方）
	while(w * 2 <= width)
	{
		w *= 2;
		wp++;
	}
	
	while(h * 2 <= height)
	{
		h *= 2;
		hp++;
	}
	double *f = new double[w * h];
	double *F = new double[w * h];
	for(i = 0; i < h; i++)
	for(j = 0; j < w; j++)
	
		f[j + w * i]= originimage[i][j*bonum+bons-1];

	for(i = 0; i < h; i++)// 对y方向进行快速付立叶变换
		WALSH(f + w * i, F + w * i, wp);

		// 保存变换结果
	for(i = 0; i < h; i++)
	for(j = 0; j < w; j++)
		f[j * h + i] = F[j + w * i];
	for(i = 0; i < w; i++)// 对x方向进行快速付立叶变换
		WALSH(f + j * h, F + j * h, hp);

	for(i = 0; i < h; i++)
	for(j = 0; j < w; j++)
	{
		dTemp = fabs(F[j + i*w] * 1000);
		if(dTemp>255)
			dTemp=255;
		originimage[i][j*bonum+bons-1]=(int)(dTemp);
	}
	//释放内存
	delete f;
	delete F;

}
//////求图像一个波段的方差：图像数据，波段数，宽度，高度，处理波段

double CBip::SS(BYTE *lpbuf, int bonum, int w, int h, int bons)
{
	BYTE **originimage;
	int i,j;
	double mean,ss;
	long sum;
	sum=0;
	ss=0;
	originimage=new BYTE*[h];
	for(i=0;i<h;i++)
			originimage[i]=lpbuf+i*bonum*w;
	for(i=0;i<h;i++)
	for(j=0;j<w;j++)
			sum+=originimage[i][j*bonum+bons-1];
	mean=sum/w/h;
	for(i=0;i<h;i++)
	for(j=0;j<w;j++)
			ss+=(originimage[i][j*bonum+bons-1]-mean)*(originimage[i][j*bonum+bons-1]-mean)/w/h;
	delete [] originimage;
	return ss;

}
//////求图像两个波段的协方差：图像数据，波段数，宽度，高度，处理波段

double CBip::SSmn(BYTE *lpbuf, int bonum, int w, int h, int bons1,int bons2)
{
	BYTE **originimage;
	int i,j;
	double mean1,mean2,ss;
	long sum1,sum2;
	sum1=0;
	sum2=0;
	ss=0;
	originimage=new BYTE*[h];
	for(i=0;i<h;i++)
			originimage[i]=lpbuf+i*bonum*w;
	for(i=0;i<h;i++)
	for(j=0;j<w;j++)
			sum1+=originimage[i][j*bonum+bons1-1];
	mean1=sum1/w/h;
	for(i=0;i<h;i++)
	for(j=0;j<w;j++)
			sum2+=originimage[i][j*bonum+bons2-1];
	mean2=sum2/w/h;

	for(i=0;i<h;i++)
	for(j=0;j<w;j++)
			ss+=(originimage[i][j*bonum+bons1-1]-mean1)*(originimage[i][j*bonum+bons2-1]-mean2)/w/h;
	delete [] originimage;
	return ss;

}
//////求图像的模板匹配：图像数据，模板数据，波段数，宽度，高度，处理波段，模板宽度，模板高度

POINT CBip::TemplateMatch(BYTE *lpbuf, BYTE *lptempate, int bonum, int w, int h, int bons, int tw, int th)
{
	BYTE **originimage;
	originimage=new BYTE*[h];
	double dSigmaST;
	double dSigmaS;
	double dSigmaT;
	POINT maxp;
	//相似性测度
	double R;
	//最大相似性测度
	double MaxR;
    int j,i,m,n;
	dSigmaT = 0;
	for(i=0;i<h;i++)
			originimage[i]=lpbuf+i*bonum*w;

	for (n = 0;n < tw ;n++)
	for(m = 0;m < th ;m++)
		dSigmaT+=(*(lptempate+m*tw+n))*(*(lptempate+m*tw+n));

	MaxR = 0.0;
	for(i = 0;i < h-th+ 1;i++)
	for (j = 0;j < w-tw+1 ;j++)
		{
			dSigmaST = 0;
			dSigmaS = 0;
	
			for (n = 0;n < th ;n++)
			for(m = 0;m < tw ;m++)
				{

					dSigmaS += (double)originimage[i][j*bonum+bons-1]*originimage[i][j*bonum+bons-1];
					dSigmaST += (double)originimage[i][j*bonum+bons-1]*(*(lptempate+n*tw+m));
				}
			
			//计算相似性
			R = dSigmaST / ( sqrt(dSigmaS)*sqrt(dSigmaT));
			//与最大相似性比较
			if (R > MaxR)
			{
				MaxR = R;
				maxp.x = i;
				maxp.y= j;
			}
		
		}
	return maxp;
}
//////求图像的协方差矩阵：图像数据，波段数，宽度，高度

double * CBip::SS(BYTE *lpbuf, int bonum, int w, int h)
{
	BYTE **originimage;
	int i,j,k,l;
	double *mean=new double[bonum];
	double *ss;
	ss=new double[bonum*bonum];
	long* sum=new long[bonum];
	for(i=0;i<bonum;i++)
	{
	sum[i]=0;
	for(j=0;j<bonum;j++)
	ss[i*bonum+j]=0;
	}
	originimage=new BYTE*[h];
	for(i=0;i<h;i++)
			originimage[i]=lpbuf+i*bonum*w;
	for(k=0;k<bonum;k++)
	for(i=0;i<h;i++)
	for(j=0;j<w;j++)
			sum[k]+=originimage[i][j*bonum+k];

	for(k=0;k<bonum;k++)
	mean[k]=sum[k]/w/h;

	for(k=0;k<bonum;k++)
	for(i=0;i<h;i++)
	for(j=0;j<w;j++)
	{
		ss[k*bonum+k]+=(originimage[i][j*bonum+k]-mean[k])*(originimage[i][j*bonum+k]-mean[k])/w/h;
		for(l=k+1;l<bonum;l++)
		ss[k*bonum+l]+=(originimage[i][j*bonum+k]-mean[k])*(originimage[i][j*bonum+l]-mean[l])/w/h;
	}
	delete [] originimage;
	return ss;

}
////矩阵相乘：A*B=C
BOOL  CBip::Mulmatri(int *a, int *b, int *c,int m, int n, int k)
{
	int i,j,l;
	for(i=0;i<m;i++)
	for(l=0;l<k;l++)
	*(c+i*m+l)=0;
	for(i=0;i<m;i++)
	for(j=0;j<n;j++)
	for(l=0;l<k;l++)
		*(c+m*i+l)+=(*(a+i*m+j))*(*(b+j*n+k));
	return TRUE;
}
////矩阵相乘：A*B=C

int * CBip::Mulmatri(int *a, int *b, int m, int n, int k)
{
	int i,j,l;
	int *c=new int[m*k];
	for(i=0;i<m;i++)
	for(l=0;l<k;l++)
	*(c+i*m+l)=0;
	for(i=0;i<m;i++)
	for(j=0;j<n;j++)
	for(l=0;l<k;l++)
		*(c+m*i+l)+=(*(a+i*m+j))*(*(b+j*n+k));
	return c;

}



long CBip::pow(int a, int r)
{
	long c;
	c=0;
	for(int i=0;i<r;i++)
		c+=a*c;
	return c;
}

long* CBip::Texture10(BYTE *lpbuf, int bonum, int w, int h, int bons)
{
	long *c=new long[256*256];
	long *p=c;
	BYTE **originimage;
	int i,j,m,n;
	originimage=new BYTE*[h];
	for(i=0;i<h;i++)
			originimage[i]=lpbuf+i*bonum*w;
	for(i=0;i<256;i++)
	for(j=0;j<256;j++)
		*(c+i*256+j)=0;

	for(i=0;i<256;i++)
	for(j=0;j<256;j++)
	for(n=0;n<h;n++)
	for(m=0;m<w-1;m++)
	{
		if(originimage[n][m*bonum+bons-1]==i&&originimage[n][(m+1)*bonum+bons-1]==j)
		*(c+i*256+j)=*(c+i*256+j)+1;
	}
	for(i=0;i<256;i++)
	for(j=i;j<256;j++)
		*(c+i*256+j)+=*(c+j*256+i);
	for(i=1;i<256;i++)
	for(j=0;j<i;j++)
		*(c+j*256+i)=*(c+i*256+j);

	delete [] originimage;
	return p;


}



void CBip::kernelsusan_corners_quick(BYTE *lpbuf, int bonum, int w, int h, int *r, int max_no, long &CornerNum)
{
	BYTE **img;
	int i,j;
	
	img=new BYTE*[h];
	for(i=0;i<h;i++)
			img[i]=lpbuf+i*bonum*w;



}