fmwnd3dmpr.cpp

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

			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::GetFittedFltData(RxRayCastingCommonInfo *pFltInfo, RxRayCastingCommonInfo *pRefInfo, CPtrArray *pSamplePt, int nZMin, int nZMax, double fRatioZ, double fRefRatioZ,
							 int nMPRThickness, CSize szImage, int nModeMPR, int level)
{	
	// 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;
	
	CPtrArray apSamplePt;
	for(int i=0 ; i<nSampleNum ; i++){
		RxPoint3D <double> *pVecPtD = (RxPoint3D <double> *)m_apSamplePt.GetAt(i);
		RxVect4D ViewingVec = pRefInfo->m_mxViewing*RxVect4D(pVecPtD->x, pVecPtD->y, pVecPtD->z);
		RxVect4D Fltvec = pFltInfo->m_mxInvViewing*ViewingVec;

		RxPoint3D <double> *pAddPtD = new RxPoint3D <double> (Fltvec.m[0], Fltvec.m[1], Fltvec.m[2]);
		apSamplePt.Add((RxPoint3D <double> *)pAddPtD);
	}

	
	RxVect4D vecTemp1 = pRefInfo->m_mxInvViewing*RxVect4D(0, 0, nZMax);
	vecTemp1.m[2] *= fRatioZ;
	RxVect4D vecTemp2 = pRefInfo->m_mxInvViewing*RxVect4D(0, 0, nZMin);
	vecTemp2.m[2] *= fRatioZ;
	
	int nDis = (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])));

	double vecX = (vecTemp2.m[0]-vecTemp1.m[0])/nDis;
	double vecY = (vecTemp2.m[1]-vecTemp1.m[1])/nDis;
	double vecZ = (vecTemp2.m[2]-vecTemp1.m[2])/nDis;
	RxVect4D Viewing1 = pRefInfo->m_mxViewing*RxVect4D(vecTemp2.m[0]+vecX, vecTemp2.m[1]+vecY, vecTemp2.m[2]+vecZ);
	RxVect4D Viewing2 = pRefInfo->m_mxViewing*RxVect4D(vecTemp2.m[0]+vecX*2, vecTemp2.m[1]+vecY*2, vecTemp2.m[2]+vecZ*2);
	
	RxVect4D Dir1 = pFltInfo->m_mxInvViewing*Viewing1;
	RxVect4D Dir2 = pFltInfo->m_mxInvViewing*Viewing2;		
	RxVect4D Dir = Dir1-Dir2;		

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

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

	
	double dTempX, dTempY, dTempZ;
	pFltInfo->m_mxInvViewing.TransformVect3D(0, 0, 1, &dTempX, &dTempY, &dTempZ);
	RxVect4D RDir(dTempX, dTempY, dTempZ*fRatioZ, 1);
	RDir.Normalize();
	

	
//	double dTempX, dTempY, dTempZ;
//	pFltInfo->m_mxInvViewing.TransformVect3D(0, 0, 1, &dTempX, &dTempY, &dTempZ);
//	RxVect4D Dir(dTempX, dTempY, dTempZ*fRatioZ, 1);
//	Dir.Normalize();
	

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

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

	CSize szOriginal;
	if (level==1) {
		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;
	
	// 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> *)apSamplePt.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, RXSERIES_FLT);
					else if (bXMirror && bYMirror)
						pImage[(szImage.cx-x-1)*xStride+(szImage.cy-y-1)*yStride] = trilinear(dLocX, dLocY, dLocZ/fRatioZ, RXSERIES_FLT);
					else if (!bXMirror && !bYMirror)
						//OK
						pImage[x*xStride+y*yStride] = trilinear(dLocX, dLocY, dLocZ/fRatioZ, RXSERIES_FLT);
					else
						//OK
						pImage[x*xStride+(szImage.cy-y-1)*yStride] = trilinear(dLocX, dLocY, dLocZ/fRatioZ, RXSERIES_FLT);
				}
				else {
					if (bXMirror && bYMirror)
						//OK
						pImage[(szImage.cy-x-1)*xStride+(szImage.cx-y-1)*yStride] = trilinear(dLocX, dLocY, dLocZ/fRatioZ, RXSERIES_FLT);
					else if (bXMirror && !bYMirror)
						//OK
						pImage[(szImage.cy-x-1)*xStride+y*yStride] = trilinear(dLocX, dLocY, dLocZ/fRatioZ, RXSERIES_FLT);
					else if (!bXMirror && bYMirror)
						//OK
						pImage[x*xStride+(szImage.cx-y-1)*yStride] = trilinear(dLocX, dLocY, dLocZ/fRatioZ, RXSERIES_FLT);
					else
						//OK
						pImage[x*xStride+y*yStride] = trilinear(dLocX, dLocY, dLocZ/fRatioZ, RXSERIES_FLT);
				}
			}
		}
	}		

	//Remove apSamplePt
	int nNewSampleNum = apSamplePt.GetSize();
	for(i=0; i<nNewSampleNum; i++)
	{
		RxPoint3D <double> *pPt = (RxPoint3D <double> *)apSamplePt.GetAt(i);
		delete pPt;
	}
	apSamplePt.RemoveAll();


	if (level == 1)
	{
		short *pRetImage = new short[szOriginal.cx*szOriginal.cy*sizeof(unsigned short)];
		memset(pRetImage, 0x00, szOriginal.cx*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::GetFittedFltData(RxRayCastingCommonInfo *pFltInfo, RxRayCastingCommonInfo *pRefInfo, CPtrArray *pSamplePt, int nZMin, int nZMax, double fRatioZ, double fRefRatioZ,
							 int nMPRThickness, CSize szImage, CSize szFittedImage, int nModeMPR, int level,
							 int &xStride, int &yStride, BOOL &bRotated, BOOL &bXMirror, BOOL &bYMirror)
{	
	// 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;
	
	CPtrArray apSamplePt;
	for(int i=0 ; i<nSampleNum ; i++){
		RxPoint3D <double> *pVecPtD = (RxPoint3D <double> *)m_apSamplePt.GetAt(i);
		RxVect4D ViewingVec = pRefInfo->m_mxViewing*RxVect4D(pVecPtD->x, pVecPtD->y, pVecPtD->z);
		RxVect4D Fltvec = pFltInfo->m_mxInvViewing*ViewingVec;

		RxPoint3D <double> *pAddPtD = new RxPoint3D <double> (Fltvec.m[0], Fltvec.m[1], Fltvec.m[2]);
		apSamplePt.Add((RxPoint3D <double> *)pAddPtD);
	}
	
	RxVect4D vecTemp1 = pRefInfo->m_mxInvViewing*RxVect4D(0, 0, nZMax);
	vecTemp1.m[2] *= fRatioZ;
	RxVect4D vecTemp2 = pRefInfo->m_mxInvViewing*RxVect4D(0, 0, nZMin);
	vecTemp2.m[2] *= fRatioZ;
	
	int nDis = (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])));

	double vecX = (vecTemp1.m[0]-vecTemp2.m[0])/nDis;
	double vecY = (vecTemp1.m[1]-vecTemp2.m[1])/nDis;
	double vecZ = (vecTemp1.m[2]-vecTemp2.m[2])/nDis;
	RxVect4D Viewing1 = pRefInfo->m_mxViewing*RxVect4D(vecTemp2.m[0]+vecX, vecTemp2.m[1]+vecY, vecTemp2.m[2]+vecZ);
	RxVect4D Viewing2 = pRefInfo->m_mxViewing*RxVect4D(vecTemp2.m[0]+vecX*2, vecTemp2.m[1]+vecY*2, vecTemp2.m[2]+vecZ*2);
	
	RxVect4D Dir1 = pFltInfo->m_mxInvViewing*Viewing1;
	RxVect4D Dir2 = pFltInfo->m_mxInvViewing*Viewing2;		
	RxVect4D Dir = Dir2-Dir1;

	int iSample = szImage.cx;
	// 捞固瘤甫 掘阑 规氢 vector 拌魂

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

	int iSkip;
	if (level==1) iSkip = 2;
	else iSkip = 1;
	
	// Make image be returned
	unsigned short *pImage = new unsigned short[szFittedImage.cx*szFittedImage.cy];
	memset(pImage, 0x00, szFittedImage.cx*szFittedImage.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> *)apSamplePt.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[(szFittedImage.cx-x-1)*xStride+y*yStride] = trilinear(dLocX, dLocY, dLocZ/fRatioZ, RXSERIES_FLT);
					else if (bXMirror && bYMirror)
						pImage[(szFittedImage.cx-x-1)*xStride+(szFittedImage.cy-y-1)*yStride] = trilinear(dLocX, dLocY, dLocZ/fRatioZ, RXSERIES_FLT);
					else if (!bXMirror && !bYMirror)
						//OK
						pImage[x*xStride+y*yStride] = trilinear(dLocX, dLocY, dLocZ/fRatioZ, RXSERIES_FLT);
					else
						//OK
						pImage[x*xStride+(szFittedImage.cy-y-1)*yStride] = trilinear(dLocX, dLocY, dLocZ/fRatioZ, RXSERIES_FLT);
				}
				else {
					if (bXMirror && bYMirror)
						//OK
						pImage[(szFittedImage.cy-x-1)*xStride+(szFittedImage.cx-y-1)*yStride] = trilinear(dLocX, dLocY, dLocZ/fRatioZ, RXSERIES_FLT);
					else if (bXMirror && !bYMirror)
						//OK
						pImage[(szFittedImage.cy-x-1)*xStride+y*yStride] = trilinear(dLocX, dLocY, dLocZ/fRatioZ, RXSERIES_FLT);
					else if (!bXMirror && bYMirror)
						//OK
						pImage[x*xStride+(szFittedImage.cx-y-1)*yStride] = trilinear(dLocX, dLocY, dLocZ/fRatioZ, RXSERIES_FLT);
					else
						//OK
						pImage[x*xStride+y*yStride] = trilinear(dLocX, dLocY, dLocZ/fRatioZ, RXSERIES_FLT);