fmwnd3dmpr.cpp

3D reconstruction, medical image processing from colons, using intel image processing for based clas

}

unsigned short *RxFMWnd3DMPR::GetRawMPRImage(RxMatrix4D mxInvViewing, double uMin, double uMax, double vMin, double vMax, int nZMin, int nZMax, 
											 double fRatioZ, int nThickness, int nResetMode, CSize &szImage, int nModeMPR, int iLevel, 
											 int &xStride, int &yStride, BOOL &bRotated, BOOL &bXMirror, BOOL &bYMirror, BOOL b2DMPRMode)
{
	int			nSlope;
	int			nSizeX = szImage.cx;
	int			nSizeY = szImage.cy;
	double		vecX, vecY, vecZ, dDistance;

	_3DMPR_INFO	*pStMPR = theFMDocVR.MPR_Get_Struct();

	RxVect4D	StartPt, EndPt;	
	CPoint		CenterPt = pStMPR->ptCenter;	
	

	if (pStMPR->iSelected == HT_3DMPR_HORIZONTAL)
	{
		StartPt	= mxInvViewing*RxVect4D(uMin, CenterPt.y, nZMin);
		EndPt	= mxInvViewing*RxVect4D(uMax, CenterPt.y, nZMin);

		StartPt.m[2]	*= fRatioZ;
		EndPt.m[2]		*= fRatioZ;
		nSlope			= HORIZONTAL;
	}
	else if(pStMPR->iSelected == HT_3DMPR_VERTICAL)
	{
		StartPt	= mxInvViewing*RxVect4D(CenterPt.x, vMin, nZMin);
		EndPt	= mxInvViewing*RxVect4D(CenterPt.x, vMax, nZMin);

		StartPt.m[2]	*= fRatioZ;
		EndPt.m[2]		*= fRatioZ;
		nSlope			= VERTICAL;
	}
	else if(pStMPR->iSelected == HT_3DMPR_SLANT || pStMPR->iSelected == HT_3DMPR_CIRCLE)
	{
		RxFMWndVR* pFMWndVR = RxGetFrameMain()->m_pFMWndVR[m_nSeries];
		if(pFMWndVR == NULL)				return NULL;
		
		SetStartEndPt(pStMPR->ptCenterSlant, pStMPR->iCircleAngle, pFMWndVR->m_szWindow);

		int nDX = m_StartPt.x-m_EndPt.x;
		int nDY = m_StartPt.y-m_EndPt.y;
		if(nDX == 0)			nSlope = VERTICAL;
		else if(nDY == 0)		nSlope = HORIZONTAL;
		else
		{
			double dSlope = ((double)nDY/(double)nDX);
			if(fabs(dSlope)>1.0)		nSlope = VERTICAL;
			else						nSlope = HORIZONTAL;
		}

		double x1 = (double)m_StartPt.x, y1 = (double)m_StartPt.y;
		double x2 = (double)m_EndPt.x, y2 = (double)m_EndPt.y;
		double slope;

		if (x1==x2) {
			StartPt = mxInvViewing*RxVect4D(x1, vMin, nZMin);
			EndPt	= mxInvViewing*RxVect4D(x1, vMax, nZMin);
		}
		else if (nSlope==HORIZONTAL) {
			slope	= (y2-y1)/(x2-x1);
			StartPt = mxInvViewing*RxVect4D(uMin, slope*(uMin-x1)+y1, nZMin);
			EndPt	= mxInvViewing*RxVect4D(uMax, slope*(uMax-x1)+y1, nZMin);
		}
		else {
			slope	= (y2-y1)/(x2-x1);
			StartPt = mxInvViewing*RxVect4D((vMin-y1)/slope+x1, vMin, nZMin);
			EndPt	= mxInvViewing*RxVect4D((vMax-y1)/slope+x1, vMax, nZMin);
		}
		// Determine the slope of curve
		StartPt.m[2]	*= fRatioZ;
		EndPt.m[2]		*= fRatioZ;
	}

	//Compute Eigenvector
	dDistance = sqrt((EndPt.m[0]-StartPt.m[0])*(EndPt.m[0]-StartPt.m[0])+
					 (EndPt.m[1]-StartPt.m[1])*(EndPt.m[1]-StartPt.m[1])+
					 (EndPt.m[2]-StartPt.m[2])*(EndPt.m[2]-StartPt.m[2]));
 	vecX = (EndPt.m[0]-StartPt.m[0])/dDistance;
	vecY = (EndPt.m[1]-StartPt.m[1])/dDistance;
	vecZ = (EndPt.m[2]-StartPt.m[2])/dDistance;	

	RemoveSamplePt();
	for(int i=0; i<dDistance; i++) {
		RxPoint3D <double> *pAddPtD = new RxPoint3D <double> 
			(StartPt.m[0]+i*vecX, StartPt.m[1]+i*vecY, StartPt.m[2]+i*vecZ);
		m_apSamplePt.Add((RxPoint3D <double> *)pAddPtD);
	}	

	return GetMPRData(mxInvViewing, &m_apSamplePt, nZMin, nZMax, fRatioZ, nThickness, 
		szImage, nModeMPR, iLevel, NULL, xStride, yStride, bRotated, bXMirror, bYMirror);
}

unsigned short* RxFMWnd3DMPR::GetMPRData(RxMatrix4D mxViewing, CPtrArray *pSamplePt, int nZMin, int nZMax, double fRatioZ,
							 int nMPRThickness, CSize &szImage, int nModeMPR, int level, BYTE *pbMPRDirection, 
							 int &xStride, int &yStride, BOOL &bRotated, BOOL &bXMirror, BOOL &bYMirror)
{
	int nSeries = m_nSeries;
	if(m_nSeries == RXSERIES_COM){
		nSeries = RXSERIES_REF;
	}

	// Fake thickness for odd value effect
	int nOriginalThick = nMPRThickness;
	if(nOriginalThick%2==1) nMPRThickness = nOriginalThick+1;

	int nSampleNum = pSamplePt->GetSize();
	// If sample points too small, return NULL
	const int nLimit = 2;
	if(nSampleNum<nLimit) return NULL;	

	RxVect4D vecTemp1 = mxViewing*RxVect4D(0, 0, nZMax);
	vecTemp1.m[2] *= fRatioZ;
	RxVect4D vecTemp2 = mxViewing*RxVect4D(0, 0, nZMin);
	vecTemp2.m[2] *= fRatioZ;

	int iSample = (int)ceil(sqrt((vecTemp2.m[0]-vecTemp1.m[0])*(vecTemp2.m[0]-vecTemp1.m[0]) + 
		(vecTemp2.m[1]-vecTemp1.m[1])*(vecTemp2.m[1]-vecTemp1.m[1]) +
		(vecTemp2.m[2]-vecTemp1.m[2])*(vecTemp2.m[2]-vecTemp1.m[2])));

	szImage.cx = iSample;
	szImage.cy = nSampleNum;

	// 捞固瘤甫 掘阑 规氢 vector 拌魂
	RxPoint3D <double> tmpV1 = *((RxPoint3D <double> *)pSamplePt->GetAt(0));
	RxPoint3D <double> tmpV2 = *((RxPoint3D <double> *)pSamplePt->GetAt(pSamplePt->GetSize()-1));
	RxPoint3D <double> tmpV3 = tmpV2-tmpV1;
	RxVect4D vecYTmp(tmpV3.x, tmpV3.y, tmpV3.z);

	double dTempX, dTempY, dTempZ;
	mxViewing.TransformVect3D(0, 0, 1, &dTempX, &dTempY, &dTempZ);
	RxVect4D Dir(dTempX, dTempY, dTempZ*fRatioZ, 1);
	Dir.Normalize();
	

	int iXSign = Dir.GetPrincipalAxisWithSign();
	int iYSign = vecYTmp.GetPrincipalAxisWithSign();

	if (iXSign==iYSign || iXSign==-iYSign)
		iYSign = vecYTmp.GetSecondaryAxisWithSign(iXSign);

//	CSize szOriginal;
	if (level==1) {
		m_szOriginal = szImage;
		szImage.cx/=2;
		szImage.cy/=2;
	}

	int iSkip;
	if (level==1) iSkip = 2;
	else iSkip = 1;

	xStride = 1;
	yStride = szImage.cx;
	if (!GetMPRDirection(iXSign, iYSign, xStride, yStride, bRotated, bXMirror, bYMirror, szImage))
		return NULL;

	// 2D curve indicator 曼绊..
	if(pbMPRDirection!=NULL) {
		if (bRotated) {
			// x 绵俊 弊妨具摆瘤夸
			if (!bYMirror)
				*pbMPRDirection = 0;
			else
				*pbMPRDirection = 1;
		}
		else {
			// y 绵俊 弊府绰 疤聪促.
			if (bYMirror)
				*pbMPRDirection = 2;
			else
				*pbMPRDirection = 3;
		}
	}

	// Make image be returned
	unsigned short *pImage = new unsigned short[szImage.cx*szImage.cy];
	memset(pImage, 0x00, szImage.cx*szImage.cy*sizeof(unsigned short));

	// Make Image
	if(nOriginalThick==1) {	// MPR thickness啊 绝绰 版快 
		RxVect4D pixLoc;
		for(int y=0; y<nSampleNum/iSkip; y++) {
			// y绵捞 弊篮 急			
			if (level==0 && ::WaitForSingleObject(g_MPRThreadKillEvent.m_hObject, 0) == WAIT_OBJECT_0) {
				delete[] pImage;
				return NULL;
			}
			
			RxPoint3D <double> TempPtY = *((RxPoint3D <double> *)pSamplePt->GetAt(y*iSkip));

			double dLocX = TempPtY.x;
			double dLocY = TempPtY.y;
			double dLocZ = TempPtY.z;

			for(int x=0; x<iSample/iSkip; x++) {
				// x 绵篮 表槛洒 ..
				dLocX += Dir[0]*iSkip;
				dLocY += Dir[1]*iSkip;
				dLocZ += Dir[2]*iSkip;

				if (!bRotated) {
					if (bXMirror && !bYMirror)
						pImage[(szImage.cx-x-1)*xStride+y*yStride] = trilinear(dLocX, dLocY, dLocZ/fRatioZ, nSeries);
					else if (bXMirror && bYMirror)
						pImage[(szImage.cx-x-1)*xStride+(szImage.cy-y-1)*yStride] = trilinear(dLocX, dLocY, dLocZ/fRatioZ, nSeries);
					else if (!bXMirror && !bYMirror)
						//OK
						pImage[x*xStride+y*yStride] = trilinear(dLocX, dLocY, dLocZ/fRatioZ, nSeries);
					else
						//OK
						pImage[x*xStride+(szImage.cy-y-1)*yStride] = trilinear(dLocX, dLocY, dLocZ/fRatioZ, nSeries);
				}
				else {
					if (bXMirror && bYMirror)
						//OK
						pImage[(szImage.cy-x-1)*xStride+(szImage.cx-y-1)*yStride] = trilinear(dLocX, dLocY, dLocZ/fRatioZ, nSeries);
					else if (bXMirror && !bYMirror)
						//OK
						pImage[(szImage.cy-x-1)*xStride+y*yStride] = trilinear(dLocX, dLocY, dLocZ/fRatioZ, nSeries);
					else if (!bXMirror && bYMirror)
						//OK
						pImage[x*xStride+(szImage.cx-y-1)*yStride] = trilinear(dLocX, dLocY, dLocZ/fRatioZ, nSeries);
					else
						//OK
						pImage[x*xStride+y*yStride] = trilinear(dLocX, dLocY, dLocZ/fRatioZ, nSeries);
				}
			}
		}
	}
	else {	// MPR thickness啊 乐绰 版快 
		RxVect4D *pixLoc = new RxVect4D[nMPRThickness];
		int nHalfThick = (int)nMPRThickness/2;

		for(int y=0; y<nSampleNum/iSkip; y++) {
			// MPR thread kill			
			if (level==0 && ::WaitForSingleObject(g_MPRThreadKillEvent.m_hObject, 0) == WAIT_OBJECT_0) {
				delete[] pixLoc;
				delete[] pImage;
				return NULL;
			}
			

			RxPoint3D <double> TempPtY = *((RxPoint3D <double> *)pSamplePt->GetAt(y*iSkip));

			// Thickness父怒 气阑 滴绊 sampling且 捞固瘤 谅钎甫 固府 拌魂
			for(int i=-nHalfThick; i<nHalfThick; i++) {
				RxVect4D thickDir;
				if (y < nSampleNum/iSkip-1) {
					RxPoint3D <double> nextPt = *((RxPoint3D <double> *)pSamplePt->GetAt(y*iSkip+1));
					RxPoint3D <double> slope = nextPt-TempPtY;
					RxVect4D tmpV = RxVect4D(slope.x, slope.y, slope.z);
					thickDir = tmpV.CrossProduct(Dir);
					thickDir.Normalize();
				}
				else {
					RxPoint3D <double> nextPt = *((RxPoint3D <double> *)pSamplePt->GetAt(y*iSkip-1));
					RxPoint3D <double> slope = nextPt-TempPtY;
					RxVect4D tmpV = RxVect4D(slope.x, slope.y, slope.z);
					thickDir = tmpV.CrossProduct(Dir);
					thickDir.Normalize();
				}
				pixLoc[i+nHalfThick] = RxVect4D(TempPtY.x+i*thickDir.m[0], TempPtY.y+i*thickDir.m[1], TempPtY.z+i*thickDir.m[2]);
			}

			for(int x=0;x<iSample/iSkip;x++) {
				long nGrey = -1;
				
				// For MIP and MINIP modes
				int nTrash = -1;
				short nVal = -1;
				short nMin = 10000, nMax = -1;
				for(int i=0; i<nMPRThickness; i++) {	// Thickness父怒 sampling秦辑 averaging
					nVal = trilinear((pixLoc[i])[0]+x*Dir[0]*iSkip, (pixLoc[i])[1]+x*Dir[1]*iSkip, ((pixLoc[i])[2]+x*Dir[2]*iSkip)/fRatioZ, nSeries);

					if(nVal==-1) { nTrash++; continue; }
					else if (level==1) nVal++;
					
					nGrey += nVal;

					if(nMin>nVal) nMin = nVal;
					if(nMax<nVal) nMax = nVal;
				}

				if (nMin==10000) nMin = -1;

				if(nOriginalThick%2==1)	{ if(nVal!=-1) nGrey -= nVal; }

				if(nOriginalThick-nTrash>0) nGrey /= (nOriginalThick-nTrash);
				else nGrey = -1;

				short nResult;
				switch(nModeMPR) {
				case MPR_MODE_RAYSUM:
					nResult = (short)nGrey;
					break;
				case MPR_MODE_MIP: 
					nResult = nMax;
					break;
				case MPR_MODE_MINIP:
					nResult = nMin;
					break;
				}

				if (!bRotated) {
					if (bXMirror && !bYMirror)
						pImage[(szImage.cx-x-1)*xStride+y*yStride] = nResult;
					else if (bXMirror && bYMirror)
						pImage[(szImage.cx-x-1)*xStride+(szImage.cy-y-1)*yStride] = nResult;
					else if (!bXMirror && !bYMirror)
						//OK
						pImage[x*xStride+y*yStride] = nResult;
					else
						//OK
						pImage[x*xStride+(szImage.cy-y-1)*yStride] = nResult;
				}
				else {
					if (bXMirror && bYMirror)
						//OK
						pImage[(szImage.cy-x-1)*xStride+(szImage.cx-y-1)*yStride] = nResult;
					else if (bXMirror && !bYMirror)
						//OK
						pImage[(szImage.cy-x-1)*xStride+y*yStride] = nResult;
					else if (!bXMirror && bYMirror)
						//OK
						pImage[x*xStride+(szImage.cx-y-1)*yStride] = nResult;
					else
						//OK
						pImage[x*xStride+y*yStride] = nResult;
				}
			}
		}
		delete []pixLoc;
	}

	if (level == 1)
	{
		short *pRetImage = new short[m_szOriginal.cx*m_szOriginal.cy*sizeof(unsigned short)];
		memset(pRetImage, 0x00, m_szOriginal.cx*m_szOriginal.cy*sizeof(unsigned short));
		short *pTempImage = (short*)pImage;

		int x, y;
		for (y=0;y<szImage.cy;y++)
			for (x=0;x<szImage.cx;x++)
				pRetImage[2*x+2*y*szImage.cx*2] = pTempImage[x+y*szImage.cx];

		for (y=0;y<szImage.cy*2;y+=2)
			for (x=1;x<szImage.cx*2-1;x+=2)
				pRetImage[x+y*szImage.cx*2] = (pRetImage[(x-1)+y*szImage.cx*2]+pRetImage[(x+1)+y*szImage.cx*2])/2;

		for (y=1;y<szImage.cy*2-1;y+=2)
			for (x=0;x<szImage.cx*2;x++)
				pRetImage[x+y*szImage.cx*2] = (pRetImage[x+(y-1)*szImage.cx*2]+pRetImage[x+(y+1)*szImage.cx*2])/2;

		CSize szTemp = szImage;
		szImage.cx = szTemp.cx*2;
		szImage.cy = szTemp.cy*2;

		delete []pImage;

		return (unsigned short*)pRetImage;	
	}
	else
		return (unsigned short*)pImage;
}

/*
unsigned short* RxFMWnd3DMPR::GetMPRData(RxMatrix4D mxViewing, CPtrArray *pSamplePt, int nZMin, int nZMax, double fRatioZ,
							 int nMPRThickness, CSize &szImage, int nModeMPR, int level, BYTE *pbMPRDirection)
{
	// Get volume global functions
//	RxVolumeData *pVolumeData = RxGetVolumeData(m_nSelectedSeries);
//	unsigned short * Volume = pVolumeData->GetBigVolume();

//	int xsize, ysize, zsize;
//	pVolumeData->GetBigVolumeDmsn(&xsize, &ysize, &zsize);

	int nSeries = m_nSeries;
	if(m_nSeries == RXSERIES_COM){
		nSeries = RXSERIES_REF;
	}

	// Fake thickness for odd value effect
	int nOriginalThick = nMPRThickness;
	if(nOriginalThick%2==1) nMPRThickness = nOriginalThick+1;

	int nSampleNum = pSamplePt->GetSize();
	// If sample points too small, return NULL
	const int nLimit = 2;
	if(nSampleNum<nLimit) return NULL;	

	RxVect4D vecTemp1 = mxViewing*RxVect4D(0, 0, nZMax);
	vecTemp1.m[2] *= fRatioZ;
	RxVect4D vecTemp2 = mxViewing*RxVect4D(0, 0, nZMin);
	vecTemp2.m[2] *= fRatioZ;

	int iSample = (int)ceil(sqrt((vecTemp2.m[0]-vecTemp1.m[0])*(vecTemp2.m[0]-vecTemp1.m[0]) + 
		(vecTemp2.m[1]-vecTemp1.m[1])*(vecTemp2.m[1]-vecTemp1.m[1]) +
		(vecTemp2.m[2]-vecTemp1.m[2])*(vecTemp2.m[2]-vecTemp1.m[2])));

	szImage.cx = iSample;
	szImage.cy = nSampleNum;

	// 捞固瘤甫 掘阑 规氢 vector 拌魂
	RxPoint3D <double> tmpV1 = *((RxPoint3D <double> *)pSamplePt->GetAt(0));
	RxPoint3D <double> tmpV2 = *((RxPoint3D <double> *)pSamplePt->GetAt(pSamplePt->GetSize()-1));
	RxPoint3D <double> tmpV3 = tmpV2-tmpV1;
	RxVect4D vecYTmp(tmpV3.x, tmpV3.y, tmpV3.z);

	double dTempX, dTempY, dTempZ;
	mxViewing.TransformVect3D(0, 0, 1, &dTempX, &dTempY, &dTempZ);
	RxVect4D Dir(dTempX, dTempY, dTempZ*fRatioZ, 1);
	Dir.Normalize();

	int xStride, yStride;
	BOOL bRotated = FALSE;
	BOOL bXMirror = FALSE, bYMirror = FALSE;

	int iXSign = Dir.GetPrincipalAxisWithSign();
	int iYSign = vecYTmp.GetPrincipalAxisWithSign();

	if (iXSign==iYSign || iXSign==-iYSign)
		iYSign = vecYTmp.GetSecondaryAxisWithSign(iXSign);

//	CSize szOriginal;
	if (level==1) {
		m_szOriginal = szImage;
		szImage.cx/=2;
		szImage.cy/=2;
	}

	int iSkip;
	if (level==1) iSkip = 2;
	else iSkip = 1;

	xStride = 1;
	yStride = szImage.cx;
	if (!GetMPRDirection(iXSign, iYSign, xStride, yStride, bRotated, bXMirror, bYMirror, szImage))
		return NULL;

	// 2D curve indicator 曼绊..
	if(pbMPRDirection!=NULL) {
		if (bRotated) {
			// x 绵俊 弊妨具摆瘤夸
			if (!bYMirror)
				*pbMPRDirection = 0;
			else
				*pbMPRDirection = 1;
		}
		else {
			// y 绵俊 弊府绰 疤聪促.
			if (bYMirror)
				*pbMPRDirection = 2;
			else
				*pbMPRDirection = 3;
		}
	}

	// Make image be returned
	unsigned short *pImage = new unsigned short[szImage.cx*szImage.cy];
	memset(pImage, 0x00, szImage.cx*szImage.cy*sizeof(unsigned short));

	// Make Image
	if(nOriginalThick==1) {	// MPR thickness啊 绝绰 版快 
		RxVect4D pixLoc;
		for(int y=0; y<nSampleNum/iSkip; y++) {
			// y绵捞 弊篮 急			
			if (level==0 && ::WaitForSingleObject(g_MPRThreadKillEvent.m_hObject, 0) == WAIT_OBJECT_0) {
				delete[] pImage;
				return NULL;
			}