mathex.cpp

这是一款蛮COOL的图像处理系统

// MathEx.cpp: implementation of the CMathEx class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "QuickImage.h"
#include "MathEx.h"
#include <math.h>

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
/*
COMPLEX::COMPLEX()
{
	re = 0.0;
	im = 0.0;
}

COMPLEX::COMPLEX(double real, double image)
{
	re = real;
	im = image;
}

COMPLEX::COMPLEX(const COMPLEX &c)
{
	re = c.re;
	im = c.im;
}

COMPLEX::~COMPLEX()
{
}

COMPLEX COMPLEX::operator +(const COMPLEX &c)
{
	COMPLEX ret;
	ret.re = re + c.re;
	ret.im = im + c.im;
	return ret;
}


COMPLEX COMPLEX::operator -(const COMPLEX &c)
{
	COMPLEX ret;
	ret.re = re - c.re;
	ret.im = im - c.im;
	return ret;
}

COMPLEX COMPLEX::operator *(const COMPLEX &c)
{
	COMPLEX ret;
	ret.re = re * c.re - im * c.im;
	ret.im = re * c.im + im * c.re;
	return ret;
}

COMPLEX& COMPLEX::operator =(const COMPLEX &c)
{
	re = c.re;
	im = c.im;
	return *this;
}

BOOL COMPLEX::operator ==(const COMPLEX &c)
{
	return ((re == c.re ) && (im == c.im));
}

COMPLEX COMPLEX::operator *(double d)
{
	COMPLEX ret;
	ret.re = re * d;
	ret.im = im * d;
	return ret;
}

COMPLEX operator *(double d, const COMPLEX &c)
{
	COMPLEX ret;
	ret.re = c.re * d;
	ret.im = c.im * d;
	return ret;
}
*/
/*complex add*/
COMPLEX Add(COMPLEX c1, COMPLEX c2)
{
	COMPLEX c;
	c.re=c1.re+c2.re;
	c.im=c1.im+c2.im;
	return c;
}

/*complex substract*/
COMPLEX Sub(COMPLEX c1, COMPLEX c2)
{
	COMPLEX c;
	c.re=c1.re-c2.re;
	c.im=c1.im-c2.im;
	return c;
}

/*complex multiple*/
COMPLEX Mul(COMPLEX c1, COMPLEX c2)
{
	COMPLEX c;
	c.re=c1.re*c2.re-c1.im*c2.im;
	c.im=c1.re*c2.im+c2.re*c1.im;
	return c;
}

CMathEx::CMathEx()
{

}

CMathEx::~CMathEx()
{

}

void CMathEx::MatrixRotate(double *after, const double *before, int row, int col)
{
	ASSERT( NULL != after);
	ASSERT( NULL != before);

	for(int i=0; i<col; i++)
		for(int j=0; j<row; j++)
			after[i * row + j] = before[j * col + i];

}

void CMathEx::MatrixMutiply
	(double *result, const double *left, const double *right, int row, int coll, int colr)
{
	ASSERT(NULL != result);
	ASSERT(NULL != left);
	ASSERT(NULL != right);
	int i,j,k;
	double *pRes = result;
	
	for(i = 0; i < row * colr; i++, pRes++)
		*pRes = 0.0;
	
	for(i = 0; i < row; i++)
		for(j = 0; j< colr;j++)
			for(k = 0;k< coll ; k++)
				result[colr*i+j] += (left[coll*i+k] * right[colr*k+j]);			
}

BOOL CMathEx::MatrixConvert(double *result, const double *before, int size)
{
	ASSERT(NULL != result);
	ASSERT(NULL != before);

	int i = size * size * sizeof(double);
	HANDLE hHeap = HeapCreate(HEAP_NO_SERIALIZE, i, 0);
	if(hHeap == NULL)
	{
		return FALSE;
	}
	HeapLock(hHeap);
	double *pSubthis = NULL;
	pSubthis=(double*)HeapAlloc(hHeap, HEAP_NO_SERIALIZE, i);
	if(pSubthis == NULL)
	{
		HeapUnlock(hHeap);
		HeapDestroy(hHeap);
		return FALSE;
	}
	memcpy((void*)pSubthis, (void*)before, i);

	double *pRes = result;
	for(i=0; i<size * size; i++, pRes++)
		*pRes = 1.0;

	int k, k2;
	int imax;
	double dmax;
	for(i=0;i<size;i++)
	{
//////////////////////////////////////////////////////////////////
		imax=i;
		dmax=pSubthis[i*size+i];
		for(k=i+1;k<size;k++)
			if(fabs(pSubthis[k*size+i])>fabs(dmax))
			{
				dmax=pSubthis[k*size+i];
				imax=k;
			}
		if(fabs(dmax)<1e-30)
		{
			AfxMessageBox("Can not convert a ill matrix!",MB_ICONERROR);
			HeapFree(hHeap,HEAP_NO_SERIALIZE | HEAP_GENERATE_EXCEPTIONS,pSubthis);
			pSubthis = NULL;
			HeapUnlock(hHeap);
			HeapDestroy(hHeap);
			return FALSE;
		}
		if(imax!=i)
			for(k=0;k<size;k++)
			{
				dmax=pSubthis[i*size+k];
				pSubthis[i*size+k]=pSubthis[imax*size+k];
				pSubthis[imax*size+k]=dmax;

				dmax = result[i*size+k];
				result[i*size+k] = result[imax*size+k];
				result[imax*size+k] = dmax;
			}
	//此间为列选主元
////////////////////////////////////////////////////////////////////
		dmax= pSubthis[i*size+i];
		for(k=i;k<size;k++)
		{
			imax = i*size+k;
			pSubthis[imax] /= dmax;
			result[imax] /= dmax;
		}
		for(k=0; k<i; k++)
		{
			dmax=pSubthis[k*size+i];
			for(k2 = i; k2 < size; k2++)
			{
				pSubthis[k*size+k2] -= (pSubthis[i*size+k2]*dmax);
				result[k*size+k2] -= (result[i*size+k2]*dmax);
			}
		}
		for(k=i+1; k<size; k++)
		{
			dmax=pSubthis[k*size+i];
			for(k2 = i; k2 < size; k2++)
			{
				pSubthis[k*size+k2] -= (pSubthis[i*size+k2]*dmax);
				result[k*size+k2] -= (result[i*size+k2]*dmax);
			}
		}
	}

	try
	{
		HeapFree(hHeap,HEAP_NO_SERIALIZE | HEAP_GENERATE_EXCEPTIONS,pSubthis);
		pSubthis = NULL;
		HeapUnlock(hHeap);
		HeapDestroy(hHeap);
	}
	catch(CMemoryException *e)
	{
		char msg[256];
		e->GetErrorMessage(msg,255);
		e->ReportError();
		e->Delete();
		return FALSE;
	}

	return TRUE;
}

void CMathEx::JieFC(double *xsz, int row, int col, double *result)
{
	double b,ab1;
	int i,j,k,j1,n1,n2,l;
	for( i=0; i<row; i++)
	{
		result[i]=0.0;
		b=0.0;
		for(j=i; j<row; j++)
		{
			if(fabs(b)<=fabs(xsz[j * col + i]))
			{
				b=xsz[j * col + i];
				j1=j;
			}
		}
		for(k=i; k<col; k++)
		{
			ab1 = xsz[j1 * col + k]/b;
			xsz[j1 * col + k] = xsz[i * col + k];
			xsz[i * col + k] = ab1;
		}
		n1=i+1;
		if(n1<row)
		{
			for(j=n1; j<row;j++)
			{
				for(k=n1;k<col;k++)
				{
					xsz[j * col + k] -= (xsz[j * col + i]*xsz[i * col + k]);
				}
			}
		}
	}
	result[row-1] = xsz[row * col -1];//xsz[(row-1) * col + col-1];
	n2=row-1;
	for(i=0;i<n2;i++)
	{
		l=n2- 1 - i;
		result[l]=xsz[l * col -1];
		for(j=l + 1; j<row; j++)
		{
			result[l] -= result[j]*xsz[l * col + j];
		}
	}
}

BOOL CMathEx::GSXQ(double *x, double *a, double *b, int size)
{
	int i, j, k;
	double dMax;
	int iMax;
	BOOL bIll = TRUE;
	for(i = 0; i< size; i++)
	{
		iMax = i;
		dMax = x[i * size + i];
		for(j = i +1; j< size; j++)
		{
			if(fabs(dMax) < fabs(a[j * size + i]))
			{
				dMax = a[j * size + i];
				iMax = j;
			}
		}
		if(dMax > -0.00000000001 && dMax < 0.00000000001)
			bIll = FALSE;
		if(iMax != i)
		{
			for(j = i; j < size; j++)
			{
				dMax = a[i * size + j];
				a[i * size + j] = a[iMax * size + j];
				a[iMax * size + j] = dMax;
			}
			dMax = b[i];
			b[i] = b[iMax];
			b[iMax] = dMax;
		}
		for(j = i +1; j < size; j++)
		{
			dMax = a[j * size + i] / a[i * size + i];
			a[j * size + i] = 0.0;
			for(k = i + 1; k< size; k++)
			{
				a[j * size + k] -= dMax * a[i * size + k];
			}
			b[j] -= dMax * b[i];
		}

	}
	x[size -1] = b[size -1] / a[size * size -1];
	for(k = size - 2; k > -1; k--)
	{
		dMax = 0.0;
		for(j = k + 1; j< size; j++)
		{
			dMax += a[k * size + j] * x[j];
		}
		x[k] = (b[k] - dMax) / a[k * size + k];
	}
	return bIll;
}

BOOL CMathEx::FFT(const COMPLEX *TD, COMPLEX *FD, int power)
{
	ASSERT((NULL != TD) && (NULL != FD));
	int i, j, k, bfsize, p, iIndex;;
	double angle;
	int count = 1 << power;//变换点数
	double PAI2 = PAI * 2.0;
	COMPLEX *W, *X1, *X2, *X = NULL;
	W = new COMPLEX[count / 2];
	if(NULL == W)
	{
		printf("Sorry, insufficent memory available!");
		return FALSE;
	}
	X1 = new COMPLEX[count];
	if(NULL == X1)
	{
		delete W;
		printf("Sorry, insufficent memory available!");
		return FALSE;
	}
	X2 = new COMPLEX[count];
	if(NULL == X2)
	{
		printf("Sorry, insufficent memory available!");
		delete X1;
		delete W;
		return FALSE;
	}
//计算加权系数
	for(i = 0; i < count / 2 ; i++)
	{
		angle = -i * PAI2 / count;
		W[i].re = cos(angle);
		W[i].im = sin(angle);
	}
	memcpy(X1, TD, sizeof(COMPLEX) * count);
//蝶形运算
	for(k = 0; k < power; k++)
	{
		for(j = 0; j < 1<<k; j++)
		{
			bfsize = (1<<(power - k)) / 2;
			for(i = 0; i < bfsize; i++)
			{
				p = j * bfsize * 2;
				iIndex = i + p;
				X2[iIndex] = Add(X1[iIndex], X1[iIndex + bfsize]);
				X2[iIndex + bfsize] = Mul(Sub(X1[iIndex], X1[iIndex + bfsize]),
					(W[i * (1 << k)]));
			}
		}
		X = X1;
		X1 = X2;
		X2 = X;
	}
//重新排续
	for(j =0 ; j < count; j++)
	{
		p = 0;
		for(i = 0; i < power; i++)
		{
			if(j & (1<<i))
			{
				p += 1<<(power - i -1);
			}
			FD[j] = X1[p];
		}
	}
	delete W;
	delete X1;
	delete X2;
	return TRUE;
}

BOOL CMathEx::IFFT(const COMPLEX *FD, COMPLEX *TD, int power)
{
	ASSERT((NULL != TD) && (NULL != FD));
	int count = 1 << power;
	COMPLEX *x = new COMPLEX[count];
	if(NULL == x)
	{
		printf("Sorry, insufficent memory available!");
		return FALSE;
	}
	memcpy(x, FD, sizeof(COMPLEX) * count);
	for(int i =0; i < count; i ++)
	{
		x[i].im = - x[i].im;
	}
	FFT(x, TD, power);

	for(i =0; i < count; i++)
	{
		TD[i].re /= count;
		TD[i].im = -TD[i].im / count;
	}
	delete x;
	return TRUE;
}


/*******************************************************
	DCT()

	参数：

		f为时域值
		F为频域值
		power为2的幂数

	返回值：

		无

	说明：

		本函数利用快速傅立叶变换实现快速离散余弦变换
********************************************************/
void CMathEx::DCT(double *f, double *F, int power)
{
	int i,count;
	COMPLEX *X;
	double s;

	/*计算离散余弦变换点数*/
	count=1<<power;
	
	/*分配运算所需存储器*/
	X=(COMPLEX *)malloc(sizeof(COMPLEX)*count*2);
	
	/*延拓*/
	memset(X,0,sizeof(COMPLEX)*count*2);
	
	/*将时域点写入存储器*/
	for(i=0;i<count;i++)
	{
		X[i].re=f[i];
	}
	
	/*调用快速傅立叶变换*/
	FFT(X,X,power+1);
	
	/*调整系数*/
	s=1/sqrt(count);
	F[0]=X[0].re*s;
	s*=sqrt(2);
	for(i=1;i<count;i++)
	{
		F[i]=(X[i].re*cos(i*PAI/(count*2))+X[i].im*sin(i*PAI/(count*2)))*s;
	}
	
	/*释放存储器*/
	free(X);
}

/************************************************************
	IDCT()

	参数：

		F为频域值
		f为时域值
		power为2的幂数

	返回值：

		无

	说明：

		本函数利用快速傅立叶反变换实现快速离散反余弦变换
*************************************************************/
void CMathEx::IDCT(double *F, double *f, int power)
{
	int i,count;
	COMPLEX *X;
	double s;

	/*计算离散反余弦变换点数*/
	count=1<<power;
	
	/*分配运算所需存储器*/
	X=(COMPLEX *)malloc(sizeof(COMPLEX)*count*2);
	
	/*延拓*/
	memset(X,0,sizeof(COMPLEX)*count*2);
	
	/*调整频域点，写入存储器*/
	for(i=0;i<count;i++)
	{
		X[i].re=F[i]*cos(i*PAI/(count*2));
		X[i].im=F[i]*sin(i*PAI/(count*2));
	}
	
	/*调用快速傅立叶反变换*/
	IFFT(X,X,power+1);
	
	/*调整系数*/
	s=1/sqrt(count);
	for(i=1;i<count;i++)
	{
		f[i]=(1-sqrt(2))*s*F[0]+sqrt(2)*s*X[i].re*count*2;
	}
	
	/*释放存储器*/
	free(X);
}

/**********************************************************
	WALh()

	参数：

		f为时域值
		W为频域值
		power为2的幂数

	返回值：

		无

	说明：

		本函数利用快速傅立叶变换实现快速沃尔什-哈达玛变换
*************************************************************/
void CMathEx::WALh(double *f, double *W, int power)
{
	int count;
	int i,j,k,bfsize,p;
	double *X1,*X2,*X;

	/*计算快速沃尔什变换点数*/

	count=1<<power;
	
	/*分配运算所需存储器*/
	X1=(double *)malloc(sizeof(double)*count);
	X2=(double *)malloc(sizeof(double)*count);
	
	/*将时域点写入存储器*/
	memcpy(X1,f,sizeof(double)*count);
	
	/*蝶形运算*/
	for(k=0;k<power;k++)
	{
		for(j=0;j<1<<k;j++)
		{
			bfsize=1<<(power-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];
			}
		}
		X=X1;
		X1=X2;
		X2=X;
	}
	/*调整系数*/
//	for(i=0;i<count;i++)
//	{
//		W[i]=X1[i]/count;
//	}
	for(j=0;j<count;j++)
	{
		p=0;
		for(i=0;i<power;i++)
		{
			if (j&(1<<i)) p+=1<<(power-i-1);
		}
		W[j]=X1[p]/count;
	}
	
	/*释放存储器*/
	free(X1);
	free(X2);
}

/*********************************************************************
	IWALh()

	参数：

		W为频域值
		f为时域值
		power为2的幂数

	返回值：

		无

	说明：

		本函数利用快速沃尔什-哈达玛变换实现快速沃尔什-哈达玛反变换
**********************************************************************/
void CMathEx::IWALh(double *W, double *f, int power)
{
	int i, count;

	/*计算快速沃尔什反变换点数*/
	count=1<<power;

	/*调用快速沃尔什-哈达玛变换*/
	WALh(W, f, power);

	/*调整系数*/
	for(i=0;i<count;i++)
	{
		f[i] *= count;
	}
}

double CMathEx::GetArea(POINT vert, POINT from, POINT to)
{
	int dx, dy;
	dx = from.x - vert.x;
	dy = from.y - vert.y;
	double sf = sqrt(double(dx * dx + dy * dy));
	double sinf = (0 == dx && 0 == dy) ? 0 : dy / sf;
	double cosf = (0 == dx && 0 == dy) ? 0 : dx / sf;
	dx = to.x - vert.x;
	dy = to.y - vert.y;
	double st = sqrt(double(dx * dx + dy * dy));
	double sint = (0 == dx && 0 == dy) ? 0 : dy / st;
	double cost = (0 == dx && 0 == dy) ? 0 : dx / st;

	return sf * (sinf * cost - cosf * sint) * st * 0.5;
}

int CMathEx::compare(const void *arg1, const void *arg2)
{
	if(*((float*)arg1) < *((float*)arg2))
		return -1;
	if(*((float*)arg1) > *((float*)arg2))
		return 1;
	return 0;
}

BOOL CMathEx::MedianFilter(float *pData, int iWidth, int iHeight)
{
	float *pSub = new float[iWidth * 2];
	if(NULL == pSub)
	{
		return FALSE;
	}
	float *pWin = new float[9];
	if(NULL == pWin)
	{
		delete pWin;
		return FALSE;
	}
	int x, y;
	int iRowSize = iWidth * sizeof(float);
	int row;

	memcpy(&pSub[iWidth], pData, iRowSize);
	for(y = 1; y < iHeight -1; y++)
	{
		memmove(pSub, &pSub[iWidth], iRowSize);
		row = y * iWidth;
		memcpy(&pSub[iWidth], &pData[row], iRowSize);

		for(x = 1; x <  iWidth-1; x++)
		{
			pWin[0] = pSub[x -1];
			pWin[1] = pSub[x];
			pWin[2] = pSub[x +1];
			pWin[3] = pSub[iWidth + x -1];
			pWin[4] = pSub[iWidth + x];
			pWin[5] = pSub[iWidth + x +1];
			pWin[6] = pData[row + iWidth + x -1];
			pWin[7] = pData[row + iWidth + x];
			pWin[8] = pData[row + iWidth + x +1];
			qsort(pWin, 9, sizeof(float), compare);

			pData[row + x] = pWin[4];
		}
	}

	delete pWin;
	delete pSub;
	return TRUE;
}