ip.cpp

转化为DIB位图再显示出来的dicom文件C++代码

	/*调整频域点，写入存储器*/
	for(i=0;i<count;i++)
	{
		X[i].re=F[i]*cos(i*PI/(count*2));
		X[i].im=F[i]*sin(i*PI/(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 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 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;
	}
}

//////////////////////////////////////////////////////////////////////////////////
// internal functions

#define Point(x,y) lpPoints[(x)+(y)*nWidth]

void GetPoints(int nWidth,int nHeight,BYTE *lpBits,BYTE *lpPoints)
{
	int x,y,p;
	int nByteWidth = WIDTHBYTES(nWidth*24); //nWidth*3;
	//if (nByteWidth%4) nByteWidth+=4-(nByteWidth%4);

	for(y=0;y<nHeight;y++)
	{
		for(x=0;x<nWidth;x++)
		{
			p=x*3+y*nByteWidth;
			lpPoints[x+y*nWidth]=(BYTE)(0.299*(float)lpBits[p+2]+0.587*(float)lpBits[p+1]+0.114*(float)lpBits[p]+0.1);
		}
	}
}

void PutPoints(int nWidth,int nHeight,BYTE *lpBits,BYTE *lpPoints)
{
	int x,y,p,p1;
	int nByteWidth = WIDTHBYTES(nWidth*24); //nWidth*3;
	//if (nByteWidth%4) nByteWidth+=4-(nByteWidth%4);

	for(y=0;y<nHeight;y++)
	{
		for(x=0;x<nWidth;x++)
		{
			p=x*3+y*nByteWidth;
			p1=x+y*nWidth;
			lpBits[p]=lpPoints[p1];
			lpBits[p+1]=lpPoints[p1];
			lpBits[p+2]=lpPoints[p1];
		}
	}
}
//////////////////////////////////////////////////////////////////////////////

/****************************************************
	FFTDIB()

	参数：

		hDIB为输入的DIB句柄

	返回值：

		成功为TRUE；失败为FALSE

	说明：

		本函数实现DIB位图的快速傅立叶变换
****************************************************/
BOOL FFTDIB(HDIB hDIB)
{
	if (hDIB == NULL)
		return FALSE;

	// start wait cursor
	WaitCursorBegin();

	HDIB hDib = NULL;
	HDIB hNewDib = NULL;
	// we only convolute 24bpp DIB, so first convert DIB to 24bpp
	WORD wBitCount = DIBBitCount(hDIB);
	if (wBitCount != 24)
	{
		hNewDib = ConvertDIBFormat(hDIB, 24, NULL);
		hDib = CopyHandle(hNewDib);
	}
	else
	{
		hNewDib = CopyHandle(hDIB);
		hDib = CopyHandle(hDIB);
	}

	if (hNewDib == NULL && hDib == NULL)
	{
		WaitCursorEnd();
		return FALSE;
	}

	// process!
	LPBYTE lpSrcDIB = (LPBYTE)GlobalLock(hDib);
	LPBYTE lpDIB = (LPBYTE)GlobalLock(hNewDib);
	LPBYTE lpInput = FindDIBBits(lpSrcDIB);
	LPBYTE lpOutput = FindDIBBits(lpDIB);
	int nWidth = DIBWidth(lpSrcDIB);
	int nHeight = DIBHeight(lpSrcDIB);
	
	int w=1,h=1,wp=0,hp=0;
	while(w*2<=nWidth)
	{
		w*=2;
		wp++;
	}
	while(h*2<=nHeight)
	{
		h*=2;
		hp++;
	}
	int x,y;
	BYTE *lpPoints=new BYTE[nWidth*nHeight];
	GetPoints(nWidth,nHeight,lpInput,lpPoints);

	COMPLEX *TD=new COMPLEX[w*h];
	COMPLEX *FD=new COMPLEX[w*h];

	for(y=0;y<h;y++)
	{
		for(x=0;x<w;x++)
		{
			TD[x+w*y].re=Point(x,y);
			TD[x+w*y].im=0;
		}
	}

	for(y=0;y<h;y++)
	{
		FFT(&TD[w*y],&FD[w*y],wp);
	}
	for(y=0;y<h;y++)
	{
		for(x=0;x<w;x++)
		{
			TD[y+h*x]=FD[x+w*y];
//			TD[x+w*y]=FD[x*h+y];
		}
	}
	for(x=0;x<w;x++)
	{
		FFT(&TD[x*h],&FD[x*h],hp);
	}

	memset(lpPoints,0,nWidth*nHeight);
	double m;
	for(y=0;y<h;y++)
	{
		for(x=0;x<w;x++)
		{
			m=sqrt(FD[x*h+y].re*FD[x*h+y].re+FD[x*h+y].im*FD[x*h+y].im)/100;
			if (m>255) m=255;
			Point((x<w/2?x+w/2:x-w/2),nHeight-1-(y<h/2?y+h/2:y-h/2))=(BYTE)(m);
		}
	}
	delete TD;
	delete FD;
	PutPoints(nWidth,nHeight,lpOutput,lpPoints);
	delete lpPoints;

	// recover
	DWORD dwSize = GlobalSize(hDib);
	memcpy(lpSrcDIB, lpDIB, dwSize);
	GlobalUnlock(hDib);
	GlobalUnlock(hNewDib);
	if (wBitCount != 24)
	{
		hNewDib = ConvertDIBFormat(hDib, wBitCount, NULL);

		lpSrcDIB = (LPBYTE)GlobalLock(hDIB);
		lpDIB = (LPBYTE)GlobalLock(hNewDib);
		dwSize = GlobalSize(hNewDib);
		memcpy(lpSrcDIB, lpDIB, dwSize);
		GlobalUnlock(hDIB);
		GlobalUnlock(hNewDib);
	}
	else
	{
		lpSrcDIB = (LPBYTE)GlobalLock(hDIB);
		lpDIB = (LPBYTE)GlobalLock(hDib);
		dwSize = GlobalSize(hDib);
		memcpy(lpSrcDIB, lpDIB, dwSize);
		GlobalUnlock(hDIB);
		GlobalUnlock(hDib);
	}

	// cleanup
	GlobalFree(hDib);
	GlobalFree(hNewDib);

	// return
	WaitCursorEnd();
	return TRUE;
}

/****************************************************
	DCTDIB()

	参数：

		hDIB为输入的DIB句柄

	返回值：

		成功为TRUE；失败为FALSE

	说明：

		本函数实现DIB位图的快速余弦变换
****************************************************/
BOOL DCTDIB(HDIB hDIB)
{
	if (hDIB == NULL)
		return FALSE;

	// start wait cursor
	WaitCursorBegin();

	HDIB hDib = NULL;
	HDIB hNewDib = NULL;
	// we only convolute 24bpp DIB, so first convert DIB to 24bpp
	WORD wBitCount = DIBBitCount(hDIB);
	if (wBitCount != 24)
	{
		hNewDib = ConvertDIBFormat(hDIB, 24, NULL);
		hDib = CopyHandle(hNewDib);
	}
	else
	{
		hNewDib = CopyHandle(hDIB);
		hDib = CopyHandle(hDIB);
	}

	if (hNewDib == NULL && hDib == NULL)
	{
		WaitCursorEnd();
		return FALSE;
	}

	// process!
	LPBYTE lpSrcDIB = (LPBYTE)GlobalLock(hDib);
	LPBYTE lpDIB = (LPBYTE)GlobalLock(hNewDib);
	LPBYTE lpInput = FindDIBBits(lpSrcDIB);
	LPBYTE lpOutput = FindDIBBits(lpDIB);
	int nWidth = DIBWidth(lpSrcDIB);
	int nHeight = DIBHeight(lpSrcDIB);

	int w=1,h=1,wp=0,hp=0;
	while(w*2<=nWidth)
	{
		w*=2;
		wp++;
	}
	while(h*2<=nHeight)
	{
		h*=2;
		hp++;
	}
	int x,y;
	BYTE *lpPoints=new BYTE[nWidth*nHeight];
	GetPoints(nWidth,nHeight,lpInput,lpPoints);

	double *f=new double[w*h];
	double *W=new double[w*h];

	for(y=0;y<h;y++)
	{
		for(x=0;x<w;x++)
		{
			f[x+y*w]=Point(x,y);
		}
	}

	for(y=0;y<h;y++)
	{
		DCT(&f[w*y],&W[w*y],wp);
	}
	for(y=0;y<h;y++)
	{
		for(x=0;x<w;x++)
		{
			f[x*h+y]=W[x+w*y];
		}
	}
	for(x=0;x<w;x++)
	{
		DCT(&f[x*h],&W[x*h],hp);
	}
	double a;
	memset(lpPoints,0,nWidth*nHeight);

	for(y=0;y<h;y++)
	{
		for(x=0;x<w;x++)
		{
			a=fabs(W[x*h+y]);
			if (a>255) a=255;
			Point(x,nHeight-y-1)=(BYTE)(a);
		}
	}
	delete f;
	delete W;
	PutPoints(nWidth,nHeight,lpOutput,lpPoints);
	delete lpPoints;

	// recover
	DWORD dwSize = GlobalSize(hDib);
	memcpy(lpSrcDIB, lpDIB, dwSize);
	GlobalUnlock(hDib);
	GlobalUnlock(hNewDib);
	if (wBitCount != 24)
	{
		hNewDib = ConvertDIBFormat(hDib, wBitCount, NULL);

		lpSrcDIB = (LPBYTE)GlobalLock(hDIB);
		lpDIB = (LPBYTE)GlobalLock(hNewDib);
		dwSize = GlobalSize(hNewDib);
		memcpy(lpSrcDIB, lpDIB, dwSize);
		GlobalUnlock(hDIB);
		GlobalUnlock(hNewDib);
	}
	else
	{
		lpSrcDIB = (LPBYTE)GlobalLock(hDIB);
		lpDIB = (LPBYTE)GlobalLock(hDib);
		dwSize = GlobalSize(hDib);
		memcpy(lpSrcDIB, lpDIB, dwSize);
		GlobalUnlock(hDIB);
		GlobalUnlock(hDib);
	}

	// cleanup
	GlobalFree(hDib);
	GlobalFree(hNewDib);

	// return
	WaitCursorEnd();
	return TRUE;
}

/****************************************************
	WALhDIB()

	参数：

		hDIB为输入的DIB句柄

	返回值：

		成功为TRUE；失败为FALSE

	说明：

		本函数实现DIB位图的快速沃尔什-哈达玛变换
****************************************************/
BOOL WALhDIB(HDIB hDIB)
{
	if (hDIB == NULL)
		return FALSE;

	// start wait cursor
	WaitCursorBegin();

	HDIB hDib = NULL;
	HDIB hNewDib = NULL;
	// we only convolute 24bpp DIB, so first convert DIB to 24bpp
	WORD wBitCount = DIBBitCount(hDIB);
	if (wBitCount != 24)
	{
		hNewDib = ConvertDIBFormat(hDIB, 24, NULL);
		hDib = CopyHandle(hNewDib);
	}
	else
	{
		hNewDib = CopyHandle(hDIB);
		hDib = CopyHandle(hDIB);
	}

	if (hNewDib == NULL && hDib == NULL)
	{
		WaitCursorEnd();
		return FALSE;
	}

	// process!
	LPBYTE lpSrcDIB = (LPBYTE)GlobalLock(hDib);
	LPBYTE lpDIB = (LPBYTE)GlobalLock(hNewDib);
	LPBYTE lpInput = FindDIBBits(lpSrcDIB);
	LPBYTE lpOutput = FindDIBBits(lpDIB);
	int nWidth = DIBWidth(lpSrcDIB);
	int nHeight = DIBHeight(lpSrcDIB);

	int w=1,h=1,wp=0,hp=0;
	while(w*2<=nWidth)
	{
		w*=2;
		wp++;
	}
	while(h*2<=nHeight)
	{
		h*=2;
		hp++;
	}
	int x,y;
	BYTE *lpPoints=new BYTE[nWidth*nHeight];
	GetPoints(nWidth,nHeight,lpInput,lpPoints);

	double *f=new double[w*h];
	double *W=new double[w*h];

	for(y=0;y<h;y++)
	{
		for(x=0;x<w;x++)
		{
			f[x+y*w]=Point(x,y);
		}
	}

	for(y=0;y<h;y++)
	{
		WALh(f+w*y,W+w*y,wp);
	}
	for(y=0;y<h;y++)
	{
		for(x=0;x<w;x++)
		{
			f[x*h+y]=W[x+w*y];
		}
	}
	for(x=0;x<w;x++)
	{
		WALh(f+x*h,W+x*h,hp);
	}
	double a;
	memset(lpPoints,0,nWidth*nHeight);

	for(y=0;y<h;y++)
	{
		for(x=0;x<w;x++)
		{
			a=fabs(W[x*h+y]*1000);
			if (a>255) a=255;
			Point(x,nHeight-y-1)=(BYTE)a;
		}
	}
	delete f;
	delete W;
	PutPoints(nWidth,nHeight,lpOutput,lpPoints);
	delete lpPoints;

	// recover
	DWORD dwSize = GlobalSize(hDib);
	memcpy(lpSrcDIB, lpDIB, dwSize);
	GlobalUnlock(hDib);
	GlobalUnlock(hNewDib);
	if (wBitCount != 24)
	{
		hNewDib = ConvertDIBFormat(hDib, wBitCount, NULL);

		lpSrcDIB = (LPBYTE)GlobalLock(hDIB);
		lpDIB = (LPBYTE)GlobalLock(hNewDib);
		dwSize = GlobalSize(hNewDib);
		memcpy(lpSrcDIB, lpDIB, dwSize);
		GlobalUnlock(hDIB);
		GlobalUnlock(hNewDib);
	}
	else
	{
		lpSrcDIB = (LPBYTE)GlobalLock(hDIB);
		lpDIB = (LPBYTE)GlobalLock(hDib);
		dwSize = GlobalSize(hDib);
		memcpy(lpSrcDIB, lpDIB, dwSize);
		GlobalUnlock(hDIB);
		GlobalUnlock(hDib);
	}

	// cleanup
	GlobalFree(hDib);
	GlobalFree(hNewDib);

	// return
	WaitCursorEnd();
	return TRUE;
}