//////////////////////////////////////////////////////////////////////////
//得到Moment特征
//lpDIBBits是經過二值化處理以后的圖像指針
//lWidth不必是4的倍數
//pFeature是長度為8的實數數組
BOOL WINAPI GetMomentFeatures (LPSTR lpDIBBits, LONG lWidth, LONG lHeight, double *pFeature)
{
	//**********************************************************************
	//下面這一段是標準代碼，可以用在每個函數的開頭
	LONG lLineBytes;	//每行的字節數
	lLineBytes = WIDTHBYTES(lWidth * 8);
	unsigned char **cTemp;
	cTemp = cmatrix(0, lHeight-1, 0, lLineBytes-1);//提取信息可以用這個數組
	unsigned char* pmid = (unsigned char *) lpDIBBits;
	long i,j;
	
	for(i = 0; i < lHeight; i++)
	{
		for(j = 0; j < lLineBytes; j++)
		{
			cTemp[i][j] = *pmid;
			pmid++;
		}
	}
	//標準代碼結束
	//***********************************************************************

	double pi = 3.1415926535;
	double fg1,fg2,fg3,fg4,fg5,fg6,fg7,fg8;
	int x,y;
	double theta; //angle
	double lm1,lm2;
	double u00=0,u02=0,u20=0,u11=0;
	double u03a=0,u03s=0,u03;//u03a>0,u03s<0
	double u30a=0,u30s=0,u30;//u30a>0,u30s<0
	double u12a=0,u12s=0,u12;//u12a>0,u12s<0
	double u21a=0,u21s=0,u21;//u21a>0,u21s<0
	double size;
	double xmax=0,xmin=10000000.0,ymax=0,ymin=100000000.0;
	double xav=0,yav=0;

	int ndots = 0;
	for( y = 0; y < lHeight; y++ )
	{
		for( x = 0; x < lWidth; x++ )
		{
			if(cTemp[y][x] == 1)
			{
				ndots++;
				xav+= x;
				yav+= y;
			}
		}
	}

	xav/=(double)ndots;
	yav/=(double)ndots;

	for( y = 0; y < lHeight; y++ )
	{
		for( x = 0; x < lWidth; x++ )
		{
			if(cTemp[y][x] == 1)
			{
				if(x>xmax)
					xmax=x;
				if(x<xmin)
					xmin=x;
				if(y>ymax)
					ymax=y;
				if(y<ymin)
					ymin=y;
			
				//add u02,u20......
				u02+=(y-yav)*(y-yav);
				u20+=(x-xav)*(x-xav);
				u11+=(x-xav)*(y-yav);
		
				u30=(x-xav)*(x-xav)*(x-xav);
				if(u30>0)
					u30a+=u30;
				else
					u30s-=u30;
			
				u03=(y-yav)*(y-yav)*(y-yav);
				if(u03>0)
					u03a+=u03;
				else
					u03s-=u03;
			
				u12=(x-xav)*(y-yav)*(y-yav);
				if(u12>0)
					u12a+=u12;
				else
					u12s-=u12;
			
				u21=(x-xav)*(x-xav)*(y-yav);
				if(u21>0)
					u21a+=u21;
				else
					u21s-=u21;
			}
		}
	}
	
	u30 = u30a - u30s;
	u03 = u03a - u03s;
	u21 = u21a - u21s;
	u12 = u12a - u12s;
	
	double yt02 = u02 / (double) ndots; 
	double yt20 = u20 / (double) ndots;
	double yt03 = u03 / (double) ndots;
	double yt30 = u30 / (double) ndots;
	double yt12 = u12 / (double) ndots;
	double yt21 = u21 / (double) ndots;
	double yt11 = u11 / (double) nodts;
	
	fg1 = yt20 + yt02;
	
	double temp, temp1;
	temp = yt20 - yt02;
	fg2 = temp * temp + 4 * yt11 * yt11;
	
	temp = yt30 - 3 * yt12;
	temp1 = 3 * yt21 + yt03;
	fg3 = temp * temp + temp1 * temp1;
	
	temp = yt30 + yt12;
	temp1 = yt21 + yt03;
	fg4 = temp * temp + temp1 * temp1;
	
	temp = yt30 + yt12;
	temp1 = yt21 + yt03;
	fg5 = (yt30 - 3 * yt12)(yt30 + yt12)(temp * temp - 3 * temp1 * temp1)
		+ (3 * yt21 - yt03)(yt21 + yt03)(3 * temp * temp - temp1 * temp1);
		
	fg6 = (yt20 - yt02)(temp * temp - temp1 * temp1)
		+ 4 * yt11 * temp * temp1;
		
	fg7 = (3 * yt12 - yt30) + temp * (temp * temp - 3 * temp1 * temp1)
		+ (3 * yt11 - yt03) * (yt21 + yt03) * (3 * temp * temp - temp1 * temp1)
	
	pFeature[0] = fg1;
	pFeature[1] = fg2;
	pFeature[2] = fg3;
	pFeature[3] = fg4;
	pFeature[4] = fg5;
	pFeature[5] = fg6;
	pFeature[6] = fg7;

	return true;

	free_cmatrix(cTemp, 0, lHeight-1, 0, lLineBytes-1);
}

double MomentFeagureDistance(double* feagure1, double* feagure2)
{
	double dist=0.0;
	
	dist+=fabs(feagure1[0] - feagure2[0]);
	dist+=fabs(feagure1[1] - feagure2[1]);
	dist+=fabs(feagure1[2] - feagure2[2]);
	dist+=fabs(feagure1[3] - feagure2[3]);
	dist+=fabs(feagure1[4] - feagure2[4]);
	dist+=fabs(feagure1[5] - feagure2[5]);
	dist+=fabs(feagure1[6] - feagure2[6]);
	
	return dist;
}

/*****************************************************************************
 * 函數名稱：cmatrix
 * 功能：申請一個無符號字符型二維數組（矩陣）
 * 參數：nrl：數組的行最低下標
 *       nrh：數組的行最高下標
 *       ncl：數組的列最低下標
 *       nch：數組的列最高下標
 * 返回：數組地址
 */
unsigned char **cmatrix(long nrl,long nrh,long ncl,long nch)
{
	int i;
	unsigned char **m;

	m=(unsigned char **) calloc((unsigned) (nrh-nrl+1),sizeof(unsigned char*));
	if (!m) nrerror("allocation failure 1 in dmatrix()");
	m -= nrl;

	for(i=nrl;i<=nrh;i++) {
		m[i]=(unsigned char *) calloc((unsigned) (nch-ncl+1),sizeof(unsigned char));
		if (!m[i]) nrerror("allocation failure 2 in dmatrix()");
		m[i] -= ncl;
	}
	return m;
}


/*****************************************************************************
 * 函數名稱：free_cmatrix
 * 功能：釋放一個無符號字符二維數組（矩陣）
 * 參數：m：原數組地址
 *       nrl：數組的行最低下標
 *       nrh：數組的行最高下標
 *       ncl：數組的列最低下標
 *       nch：數組的列最高下標
 * 無返回值
 */
void free_cmatrix(unsigned char **m,long nrl,long nrh,long ncl,long nch)
{
	int i;

	for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
	free((char*) (m+nrl));
}