rxxpowellmethod.cpp

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

	lfSz = m_lfOptimizedScaleZ;
	lfSHx = m_lfShearX;
	lfSHy = m_lfShearY;
	lfSHz = m_lfShearZ;

	switch(nParameterNum)
	{
	case 0:
		lfTx = 0.0;
		lfTy = 0.0;
		lfTz = - (g_iFloatBinaryTransZ - g_iRefBinaryTransZ);
		lfRadX = 0.0;
		lfRadY = 0.0;
		lfRadZ = 0.0;
		lfSx = m_lfOptimizedScaleX;
		lfSy = m_lfOptimizedScaleY;
		lfSz = m_lfOptimizedScaleZ;
		lfSHx = 0.0;
		lfSHy = 0.0;
		lfSHz = 0.0;
	case 1:
		lfTx = lfValue;
		break;
	case 2:
		lfTy = lfValue;
		break;
	case 3:
		lfTz = lfValue; 
		break;
	case 4:
		lfRadX = lfValue; 
		break;
	case 5:
		lfRadY = lfValue; 
		break;
	case 6:
		lfRadZ = lfValue; 
		break;
	case 7:
		lfSHx = lfValue;
		break;
	case 8:
		lfSHy = lfValue;
		break;
	case 9:
		lfSHz = lfValue;
		break;
	case 10:
		lfTx = m_lfBoundingBoxX;
		lfTy = m_lfBoundingBoxY;
		lfTz = 0.0;
		lfRadX = 0.0;
		lfRadY = 0.0;
		lfRadZ = 0.0;
		lfSx = m_lfOptimizedScaleX;
		lfSy = m_lfOptimizedScaleY;
		lfSz = m_lfOptimizedScaleZ;
		lfSHx = 0.0;
		lfSHy = 0.0;
		lfSHz = 0.0;
		break;
	}

	double lfCorrelation = 0.0;
	_FPoint ptRot;
	ptRot.x = lfRadX;		ptRot.y = lfRadY;		ptRot.z = lfRadZ;
	RxMatrix4D TrRot = MakeRotationMatrix(ptRot);
	RxMatrix4D SH(1.0, lfSHx, lfSHx, 0.0, 
		lfSHy, 1.0, lfSHy, 0.0,
		lfSHz, lfSHz, 1.0, 0.0,
		0.0, 0.0, 0.0, 1.0);

	double lfFltCx = g_FloatCenX;
	double lfFltCy = g_FloatCenY;
	double lfFltCz = g_FloatCenZ;

	double lfRefCx = g_RefCenX;
	double lfRefCy = g_RefCenY;
	double lfRefCz = g_RefCenZ;

	int nSamplingRate = 1;
	int nTz, nTy, nTx;
	int nCnt=0;
	int i, j, k, walk;
	double lfSumOfDistance = 0.0;
	double lfNumOfDistance = 0.0;

	for (walk = 0; walk < m_nNumEdge; walk++)
	{
		i = m_pEdgeListX[walk];
		j = m_pEdgeListY[walk];
		k = m_pEdgeListZ[walk];
		RxVect4D v(((double)i - lfFltCx)*m_lfOptimizedScaleX, 
			((double)j - lfFltCy)*m_lfOptimizedScaleY, 
			((double)k - lfFltCz)*m_lfOptimizedScaleZ, 1.);

		v = TrRot*v;
		v = SH*v;
		
		v.m[0] += (lfRefCx + lfTx);
		v.m[1] += (lfRefCy + lfTy);
		v.m[2] += (lfRefCz + lfTz);
				
		if(v.m[0]>=0 && v.m[0]<m_pRefDistanceMap->m_VolX-1 &&
			v.m[1]>=0 && v.m[1]<m_pRefDistanceMap->m_VolY-1 &&
			v.m[2]>=0 && v.m[2]<m_pRefDistanceMap->m_VolZ-1){
			switch (nDistanceMode)
			{
			case 0:
				lfSumOfDistance += m_pRefDistanceMap->GetDistance(v.m[0], v.m[1], v.m[2]);
				break;
			case 1:
				lfSumOfDistance += m_pRefDistanceMap->GetPropagationDistance(v.m[0], v.m[1], v.m[2]);
//				lfSumOfDistance += m_pRefDistanceMap->GetPropagationDistanceZRatio(v.m[0], v.m[1], v.m[2]);
				break;
			case 2:
				lfSumOfDistance += m_pRefDistanceMap->GetEuclideanDistance(v.m[0], v.m[1], v.m[2]);
				break;
			}

			lfNumOfDistance++;
//				lfCorrelation += pRefMap[nTz*m_pRefer->nVolX*m_pRefer->nVolY+nTy*m_pRefer->nVolX+nTx];
//				lfCorrelation += PartialVolumeInterpolation(v.m[0], v.m[1], v.m[2], pRefMap);
//				lfCorrelation += TLInterPolation(v.m[0], v.m[1], v.m[2], m_pRefer->nVolX, m_pRefer->nVolY, m_pRefer->nVolZ, pRefMap);
//				lfCorrelation += NearestInterpolation(v.m[0], v.m[1], v.m[2], pRefMap);
		}
	}

	if (lfNumOfDistance == 0)
		return 1000.0f;

	return lfSumOfDistance/lfNumOfDistance;
}

double RxxPowellMethod::StdDistanceDifference(short nDistanceMode, double lfTx, double lfTy, double lfTz,
										 double lfSx, double lfSy, double lfSz, double lfRadX, double lfRadY, double lfRadZ)
{
	double lfCorrelation = 0.0;
	_FPoint ptRot;
	ptRot.x = lfRadX;		ptRot.y = lfRadY;		ptRot.z = lfRadZ;
	RxMatrix4D TrRot = MakeRotationMatrix(ptRot);

	double lfFltCx = g_FloatCenX;
	double lfFltCy = g_FloatCenY;
	double lfFltCz = g_FloatCenZ;

	double lfRefCx = g_RefCenX;
	double lfRefCy = g_RefCenY;
	double lfRefCz = g_RefCenZ;

	int nSamplingRate = 1;
	int nTz, nTy, nTx;
	int nCnt=0;
	int i, j, k, walk;
	double lfSumOfDistance = 0.0;
	double lfNumOfDistance = 0.0;

	for (walk = 0; walk < m_nNumEdge; walk++)
	{
			i = m_pEdgeListX[walk];
			j = m_pEdgeListY[walk];
			k = m_pEdgeListZ[walk];
			RxVect4D v(((double)i - lfFltCx)*m_lfOptimizedScaleX, 
				((double)j - lfFltCy)*m_lfOptimizedScaleY, 
				((double)k - lfFltCz)*m_lfOptimizedScaleZ, 1.);

			v = TrRot*v;
			
			v.m[0] += (lfRefCx + lfTx);
			v.m[1] += (lfRefCy + lfTy);
			v.m[2] += (lfRefCz + lfTz);
			
			if(v.m[0]>=0 && v.m[0]<m_pRefDistanceMap->m_VolX-1 &&
				v.m[1]>=0 && v.m[1]<m_pRefDistanceMap->m_VolY-1 &&
				v.m[2]>=0 && v.m[2]<m_pRefDistanceMap->m_VolZ-1){
				switch (nDistanceMode)
				{
				case 0:
					lfSumOfDistance += m_pRefDistanceMap->GetDistance(v.m[0], v.m[1], v.m[2]);
					break;
				case 1:
					lfSumOfDistance += m_pRefDistanceMap->GetPropagationDistance(v.m[0], v.m[1], v.m[2]);
//					lfSumOfDistance += m_pRefDistanceMap->GetPropagationDistanceZRatio(v.m[0], v.m[1], v.m[2]);
					break;
				case 2:
					lfSumOfDistance += m_pRefDistanceMap->GetEuclideanDistance(v.m[0], v.m[1], v.m[2]);
					break;
				}
				lfNumOfDistance++;
			}
	}

	double average = lfSumOfDistance/lfNumOfDistance;
	double std;
	double num = 0;
	
	for (walk = 0; walk < m_nNumEdge; walk++)
	{
			i = m_pEdgeListX[walk];
			j = m_pEdgeListY[walk];
			k = m_pEdgeListZ[walk];
			RxVect4D v(((double)i - lfFltCx)*m_lfOptimizedScaleX, 
				((double)j - lfFltCy)*m_lfOptimizedScaleY, 
				((double)k - lfFltCz)*m_lfOptimizedScaleZ, 1.);

			v = TrRot*v;
			
			v.m[0] += (lfRefCx + lfTx);
			v.m[1] += (lfRefCy + lfTy);
			v.m[2] += (lfRefCz + lfTz);
			
			if(v.m[0]>=0 && v.m[0]<m_pRefDistanceMap->m_VolX-1 &&
				v.m[1]>=0 && v.m[1]<m_pRefDistanceMap->m_VolY-1 &&
				v.m[2]>=0 && v.m[2]<m_pRefDistanceMap->m_VolZ-1){
				switch (nDistanceMode)
				{
				case 0:
					lfSumOfDistance += (m_pRefDistanceMap->GetDistance(v.m[0], v.m[1], v.m[2]) - average)*(m_pRefDistanceMap->GetDistance(v.m[0], v.m[1], v.m[2]) - average);
					break;
				case 1:
					lfSumOfDistance += (m_pRefDistanceMap->GetPropagationDistance(v.m[0], v.m[1], v.m[2]) - average)*(m_pRefDistanceMap->GetPropagationDistance(v.m[0], v.m[1], v.m[2]) - average);
//					lfSumOfDistance += (m_pRefDistanceMap->GetPropagationDistanceZRatio(v.m[0], v.m[1], v.m[2]) - average)*(m_pRefDistanceMap->GetPropagationDistanceZRatio(v.m[0], v.m[1], v.m[2]) - average);
					break;
				case 2:
					lfSumOfDistance += (m_pRefDistanceMap->GetEuclideanDistance(v.m[0], v.m[1], v.m[2]) - average)*(m_pRefDistanceMap->GetEuclideanDistance(v.m[0], v.m[1], v.m[2]) - average);
					break;
				}
			}
	}

	return sqrt(lfSumOfDistance/lfNumOfDistance);
}

void RxxPowellMethod::ControlProgress(RxProgressWnd *pProWnd)
{
	if(pProWnd && pProWnd->Cancelled()){
		delete pProWnd;
		pProWnd = NULL;
		return;
	}			
	else if(pProWnd){
		pProWnd->StepIt();
		pProWnd->PeekAndPump();
	}
}

RxMatrix4D RxxPowellMethod::MakeRotationMatrix(_FPoint Rot)
{
	double	Cx, Sx, 
			Cy, Sy,
			Cz, Sz;

	Cx = cos(Rot.x);
	Sx = sin(Rot.x);
	Cy = cos(Rot.y);
	Sy = sin(Rot.y);
	Cz = cos(Rot.z);
	Sz = sin(Rot.z);


	RxMatrix4D Rx(1,	0,		0,		0,
				  0,	Cx,		-Sx,	0,
				  0,	Sx,		Cx,		0,
				  0,	0,		0,		1);
	
	RxMatrix4D Ry(Cy,	0,		Sy,		0,
				  0,	1,		0,		0,
				  -Sy,	0,		Cy,		0,
				  0,	0,		0,		1);

	RxMatrix4D Rz(Cz,	-Sz,	0,		0,
				  Sz,	Cz,		0,		0,
				  0,	0,		1,		0,
				  0,	0,		0,		1);
		
	RxMatrix4D Tr = Rz*Ry*Rx;
	
	return Tr;
}

void RxxPowellMethod::KeepParameter(short nParameterNum, double lfValue)
{
	switch (nParameterNum)
	{
	case 1:
		m_lfTransX = lfValue;
		break;
	case 2:
		m_lfTransY = lfValue;
		break;
	case 3:
		m_lfTransZ = lfValue;
		break;
	case 4:
		m_lfRotX = lfValue;
		m_lfRadX = m_lfRotX*DEGREETORADIAN;
		break;
	case 5:
		m_lfRotY = lfValue;
		m_lfRadY = m_lfRotY*DEGREETORADIAN;
		break;
	case 6:
		m_lfRotZ = lfValue;
		m_lfRadZ = m_lfRotZ*DEGREETORADIAN;
		break;
	case 7:
		m_lfShearX = lfValue;
		break;
	case 8:
		m_lfShearY = lfValue;
		break;
	case 9:
		m_lfShearZ = lfValue;
		break;
	}

	g_outFile << m_nCount << '\t' << m_lfMaxCC << '\n';
}

double RxxPowellMethod::NearestInterpolation(double lfTransformedX, double lfTransformedY, double lfTransformedZ, unsigned char* pRefMap)
{
	int nTransformedXInt = (int)lfTransformedX;
	int nTransformedYInt = (int)lfTransformedY;
	int nTransformedZInt = (int)lfTransformedZ;
	double XFract = lfTransformedX - (double)nTransformedXInt;
	double YFract = lfTransformedY - (double)nTransformedYInt;
	double ZFract = lfTransformedZ - (double)nTransformedZInt;
	double w[8];
	//compute weight
	//Top
	// w2 w1
	// w3 w4

	//Bottom
	// w6 w5
	// w7 w8

	int nOffsetX[8] = {0,1,1,0,0,1,1,0};
	int nOffsetY[8] = {1,1,0,0,1,1,0,0};
	int nOffsetZ[8] = {0,0,0,0,1,1,1,1};		

	ComputeSquareWeight(w, XFract, YFract, ZFract);		
	double lfMaxWeight = 0;
	int nMinDistIdx = 0;
	for(int nIndex =0 ; nIndex<8 ; nIndex++){		
		if(lfMaxWeight < w[nIndex]){
			lfMaxWeight = w[nIndex];
			nMinDistIdx = nIndex;
		}
	}
		
	int nXRef = nTransformedXInt + nOffsetX[nMinDistIdx];
	int nYRef = nTransformedYInt + nOffsetY[nMinDistIdx];
	int nZRef = nTransformedZInt + nOffsetZ[nMinDistIdx];

	double lfRefValue = pRefMap[nZRef*m_pRefer->nVolX*m_pRefer->nVolY+nYRef*m_pRefer->nVolX+nXRef];		
	return lfRefValue;
}

double RxxPowellMethod::PartialVolumeInterpolation(double lfTransformedX, double lfTransformedY, double lfTransformedZ, unsigned char* pRefMap)
{
	double lfT = 1000000.0;
	UINT lfFixedX = (UINT)((lfTransformedX+0.00000000001)*lfT);
	UINT lfFixedY = (UINT)((lfTransformedY+0.00000000001)*lfT);
	UINT lfFixedZ = (UINT)((lfTransformedZ+0.00000000001)*lfT);
	lfTransformedX = (double)lfFixedX/lfT;
	lfTransformedY = (double)lfFixedY/lfT;
	lfTransformedZ = (double)lfFixedZ/lfT;

	int nTransformedXInt = (int)(lfTransformedX);