rxxpowellmethod.cpp

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

	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	n8	n7
	// w3 w4	n5	n6

	//Bottom
	// w6 w5	n4	n3
	// w7 w8	n1	n2

	// Direction : Bottom to Top

	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 lfValue = 0;
	for(int nIndex =0 ; nIndex<8 ; nIndex++){
		int nXRef = nTransformedXInt + nOffsetX[nIndex];
		int nYRef = nTransformedYInt + nOffsetY[nIndex];
		int nZRef = nTransformedZInt + nOffsetZ[nIndex];
		int nRefValue = pRefMap[nZRef*m_pRefer->nVolX*m_pRefer->nVolY+nYRef*m_pRefer->nVolX+nXRef];		
		lfValue += nRefValue*w[nIndex];
	}
	return lfValue;
}

float RxxPowellMethod::TLInterPolation(float x, float y, float z, int iX, int iY, int iZ, BYTE* pRefImg)
{
	float tx, ty, tz, txc, tyc;
	int ix, iy, iz;	

	ix = x;
	iy = y;
	iz = z;	

	unsigned short a,b,c,d,e,f,g,h;
	a = pRefImg[iz*iX*iY + iy*iX + ix];
	b = pRefImg[(iz+1)*iX*iY + iy*iX + ix];
	c = pRefImg[iz*iX*iY + (iy+1)*iX + ix];
	d = pRefImg[iz*iX*iY + iy*iX + ix+1];
	e = pRefImg[(iz+1)*iX*iY + (iy+1)*iX + ix];
	f = pRefImg[(iz+1)*iX*iY + iy*iX + ix];
	g = pRefImg[iz*iX*iY + (iy+1)*iX + ix+1];
	h = pRefImg[(iz+1)*iX*iY + (iy+1)*iX + ix+1];	

	tx = x-ix;
	ty = y-iy;
	tz = z-iz;

	txc = 1.f-tx;
	tyc = 1.f-ty;


	return
		(
			(
				((float)(a) * txc + (float)(b) * tx) * tyc +
				((float)(c) * txc + (float)(e) * tx) * ty
			) * (1 - tz) + 
			(
				((float)(d) * txc + (float)(f) * tx) * tyc +
				((float)(g) * txc + (float)(h) * tx) * ty
			) * tz
		);
}

void RxxPowellMethod::ComputeSquareWeight(double *w, double xFract, double yFract, double zFract)
{
	//compute weight
	//Top
	// w1 w0
	// w2 w3

	//Bottom
	// w5 w4
	// w6 w7

	// w2 w1	n8	n7
	// w3 w4	n5	n6

	//Bottom
	// w6 w5	n4	n3
	// w7 w8	n1	n2

	double xFractComplement = (1-xFract);
	double yFractComplement = (1-yFract);
	double zFractComplement = (1-zFract);

	w[0] = xFractComplement * yFract * zFractComplement;
	w[1] = xFract * yFract * zFractComplement;
	w[2] = xFract * yFractComplement * zFractComplement;
	w[3] = xFractComplement * yFractComplement* zFractComplement;
	w[4] = xFractComplement * yFract*zFract;
	w[5] = xFract * yFract * zFract;
	w[6] = xFract * yFractComplement * zFract;
	w[7] = xFractComplement * yFractComplement * zFract;

	for(int i=0 ; i<8 ; i++){
		if(w[i] < 0)	w[i] = 0.0f;		
	}
}

void RxxPowellMethod::FindOptimalTransformationForNoduleRegistration(int nDistanceMode)
{
	DWORD CurrentTime = GetTickCount();

	m_nCount = 0;

	// for paper resolution
	int nIteration = 4;
	double lfTransRange = 20.0f;
	double lfRotRange = 5.0f; //Degree
	double lfAttenuation = 0.5;
	double lfStep = 0.5f;

	g_outFile << "* Nodule registration optimization parameters\n";
	g_outFile << "(nIteration, lfTransRange, lfRotRange, lfAttenuation, lfStep) = "
		<< nIteration << ", " << lfTransRange << ", " << lfRotRange << ", " << lfAttenuation << ", " <<  lfStep << "\n";

	RxProgressWnd*	pWndProgress;
	pWndProgress = new RxProgressWnd;
	int nSetRange = ((nIteration)*6);
	pWndProgress->GoModal();
	pWndProgress->SetWindowSize(0);
	pWndProgress->SetRange(0, nSetRange);
	pWndProgress->SetText(_T("Nodule Registering"));

//	CalculateNoduleCorrespondence(m_lfTransX + m_lfBoundingBoxX, m_lfTransY + m_lfBoundingBoxY, m_lfTransZ,
//		m_lfOptimizedScaleX, m_lfOptimizedScaleY, m_lfOptimizedScaleZ,
//		m_lfRotX, m_lfRotY, m_lfRotZ);
	
	m_lfMaxNoduleDistance = AverageNoduleDistance(m_lfTransX + m_lfBoundingBoxX, m_lfTransY + m_lfBoundingBoxY, m_lfTransZ,
		m_lfOptimizedScaleX, m_lfOptimizedScaleY, m_lfOptimizedScaleZ,
		m_lfRotX, m_lfRotY, m_lfRotZ);

	m_lfTempNoduleDistance = m_lfMaxNoduleDistance;

	double lfCC;

	for(int i=0 ; i<nIteration ; i++){
		//x trans
		for(double Tx=-lfTransRange/2.0f ; Tx<lfTransRange/2.0f ; Tx+=lfStep){
			m_nCount++;
			lfCC = AverageNoduleDistance(m_lfBoundingBoxX + m_lfTransX+Tx, m_lfBoundingBoxY + m_lfTransY, m_lfTransZ, m_lfOptimizedScaleX, m_lfOptimizedScaleY, m_lfOptimizedScaleZ, m_lfRotX, m_lfRotY, m_lfRotZ);
			if(lfCC < m_lfMaxNoduleDistance){
				m_lfMaxNoduleDistance = lfCC;
				KeepParameterForNR(m_lfOptimizedScaleX, m_lfOptimizedScaleY, m_lfOptimizedScaleZ, m_lfTransX+Tx, m_lfTransY, m_lfTransZ, m_lfRotX, m_lfRotY, m_lfRotZ);
			}
		}
 		ControlProgress(pWndProgress);
		//y trans
		for(double Ty=-lfTransRange/2.0f ; Ty<lfTransRange/2.0f ; Ty+=lfStep){
			m_nCount++;
			lfCC = AverageNoduleDistance(m_lfBoundingBoxX + m_lfTransX, m_lfBoundingBoxY + m_lfTransY+Ty, m_lfTransZ, m_lfOptimizedScaleX, m_lfOptimizedScaleY, m_lfOptimizedScaleZ, m_lfRotX, m_lfRotY, m_lfRotZ);	
			if(lfCC < m_lfMaxNoduleDistance){
				m_lfMaxNoduleDistance = lfCC;
				KeepParameterForNR(m_lfOptimizedScaleX, m_lfOptimizedScaleY, m_lfOptimizedScaleZ, m_lfTransX, m_lfTransY+Ty, m_lfTransZ, m_lfRotX, m_lfRotY, m_lfRotZ);
			}
		}
		ControlProgress(pWndProgress);
		//z Rot		
		for(double Rz=-lfRotRange/2.0f ; Rz<lfRotRange/2.0f ; Rz+=lfStep){
			m_nCount++;
			lfCC = AverageNoduleDistance(m_lfBoundingBoxX + m_lfTransX, m_lfBoundingBoxY + m_lfTransY, m_lfTransZ, m_lfOptimizedScaleX, m_lfOptimizedScaleY, m_lfOptimizedScaleZ, m_lfRotX, m_lfRotY, m_lfRotZ + Rz);
			if(lfCC < m_lfMaxNoduleDistance){
				m_lfMaxNoduleDistance = lfCC;
				KeepParameterForNR(m_lfOptimizedScaleX, m_lfOptimizedScaleY, m_lfOptimizedScaleZ, m_lfTransX, m_lfTransY, m_lfTransZ, m_lfRotX, m_lfRotY, m_lfRotZ+Rz);
			}		
		}
		ControlProgress(pWndProgress);
		//x Rot
		for(double Rx=-lfRotRange/2.0f ; Rx<lfRotRange/2.0f ; Rx+=lfStep){
			m_nCount++;
			lfCC = AverageNoduleDistance(m_lfBoundingBoxX + m_lfTransX, m_lfBoundingBoxY + m_lfTransY, m_lfTransZ, m_lfOptimizedScaleX, m_lfOptimizedScaleY, m_lfOptimizedScaleZ, m_lfRotX + Rx, m_lfRotY, m_lfRotZ);
			if(lfCC < m_lfMaxNoduleDistance){
				m_lfMaxNoduleDistance = lfCC;	
				KeepParameterForNR(m_lfOptimizedScaleX, m_lfOptimizedScaleY, m_lfOptimizedScaleZ, m_lfTransX, m_lfTransY, m_lfTransZ, Rx+m_lfRotX, m_lfRotY, m_lfRotZ);
			}
		}
		ControlProgress(pWndProgress);
		//y Rot
		for(double Ry=-lfRotRange/2.0f ; Ry<lfRotRange/2.0f ; Ry+=lfStep){
			m_nCount++;
			lfCC = AverageNoduleDistance(m_lfBoundingBoxX + m_lfTransX, m_lfBoundingBoxY + m_lfTransY, m_lfTransZ, m_lfOptimizedScaleX, m_lfOptimizedScaleY, m_lfOptimizedScaleZ, m_lfRotX, m_lfRotY + Ry, m_lfRotZ);
			if(lfCC < m_lfMaxNoduleDistance){
				m_lfMaxNoduleDistance = lfCC;
				KeepParameterForNR(m_lfOptimizedScaleX, m_lfOptimizedScaleY, m_lfOptimizedScaleZ, m_lfTransX, m_lfTransY, m_lfTransZ, m_lfRotX, m_lfRotY+Ry, m_lfRotZ);
			}
		}
		ControlProgress(pWndProgress);
		//z trans
		for(double Tz=-lfTransRange/2.0f ; Tz<lfTransRange/2.0f ; Tz+=lfStep){
			m_nCount++;
			lfCC = AverageNoduleDistance(m_lfBoundingBoxX + m_lfTransX, m_lfBoundingBoxY + m_lfTransY, m_lfTransZ+Tz, m_lfOptimizedScaleX, m_lfOptimizedScaleY, m_lfOptimizedScaleZ, m_lfRotX, m_lfRotY, m_lfRotZ);
			if(lfCC < m_lfMaxNoduleDistance){
				m_lfMaxNoduleDistance = lfCC;
				KeepParameterForNR(m_lfOptimizedScaleX, m_lfOptimizedScaleY, m_lfOptimizedScaleZ, m_lfTransX, m_lfTransY, m_lfTransZ+Tz, m_lfRotX, m_lfRotY, m_lfRotZ);
			}
		}
		ControlProgress(pWndProgress);
		lfTransRange *= lfAttenuation;
		lfRotRange *= lfAttenuation;
		lfStep *= lfAttenuation;
	}

	if (pWndProgress) {	delete pWndProgress;	pWndProgress = NULL;	}

	m_lfTempCC = AverageDistanceDifference(nDistanceMode, 0, 0);
	m_lfSurfaceAccuracyIR = AverageDistanceDifference(nDistanceMode, 10, 0);
	m_lfNoduleAccuracyIR = AverageNoduleDistance(m_lfBoundingBoxX, m_lfBoundingBoxY, 0.0, m_lfOptimizedScaleX, m_lfOptimizedScaleY, m_lfOptimizedScaleZ, 0.0, 0.0, 0.0);

	DWORD ElapsedTime = GetTickCount() - CurrentTime;	

	g_outFile << "\n* Nodule registration result\n";
	g_outFile << "Scale( " << m_lfOptimizedScaleX << ", " << m_lfOptimizedScaleY << ", " << m_lfOptimizedScaleZ << " )\n";
	g_outFile << "Trans( " << m_lfBoundingBoxX + m_lfTransX << ", " << m_lfBoundingBoxY + m_lfTransY << ", " << m_lfTransZ << " )\n";
	g_outFile << "RotationDeg( " << m_lfRotX << ", " << m_lfRotY << ", " << m_lfRotZ << " )\n";
	g_outFile << "MinAverageDistance = " << m_lfMaxCC << '\n';
	g_outFile << "MinAverageNoduleDistance = " << m_lfMaxNoduleDistance << '\n';
	g_outFile << "Registration time = " << ElapsedTime / 1000.0 << '\n';
	m_lfTempTime = ElapsedTime / 1000.0;
}

double RxxPowellMethod::AverageNoduleDistance(double lfTx, double lfTy, double lfTz, double lfSx, double lfSy, double lfSz, double lfDegX, double lfDegY, double lfDegZ)
{
	RxVolumeData* pRefVol = RxGetVolumeData(0);	
	RxVolumeData* pFloatVol = RxGetVolumeData(1);	

	double lfRealX, lfRealY, lfRealZ;
	pRefVol->GetVoxelSize(&lfRealX, &lfRealY, &lfRealZ);	

	double DistAfter = 0.0;

	RxMatrix4D m_mxRegistration;

	double lfSHx = m_lfShearX;
	double lfSHy = m_lfShearY;
	double lfSHz = m_lfShearZ;

	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);
	RxMatrix4D Temp;

	m_mxRegistration.LoadIdentity();
	m_mxRegistration.Translate(-g_FloatCenX, -g_FloatCenY, -g_FloatCenZ + g_iFloatBinaryTransZ);
	m_mxRegistration.Scale(lfSx, lfSy, lfSz);
	m_mxRegistration.Rotate(0, lfDegX);
	m_mxRegistration.Rotate(1, lfDegY);
	m_mxRegistration.Rotate(2, lfDegZ);

	Temp = m_mxRegistration;
	m_mxRegistration = SH * Temp;

	m_mxRegistration.Translate(g_RefCenX + lfTx, g_RefCenY + lfTy, g_RefCenZ + lfTz - g_iRefBinaryTransZ);

	for (int num = 0; num < pRefVol->m_nNumNodule; num++)
	{
		RxVect4D registeredNodule = m_mxRegistration * RxVect4D(pFloatVol->m_nNoduleCenterX[FloatNoduleNumber(num)], pFloatVol->m_nNoduleCenterY[FloatNoduleNumber(num)], pFloatVol->m_nNoduleCenterZ[FloatNoduleNumber(num)], 1);		

		DistAfter += sqrt(pow(pRefVol->m_nNoduleCenterX[num]*lfRealX - registeredNodule.m[0]*lfRealX, 2.0) + 
			pow(pRefVol->m_nNoduleCenterY[num]*lfRealY - registeredNodule.m[1]*lfRealY, 2.0) + 
			pow(pRefVol->m_nNoduleCenterZ[num]*lfRealZ - registeredNodule.m[2]*lfRealZ, 2.0));
	}

	DistAfter /= (double)pRefVol->m_nNumNodule;
	
	return DistAfter;
}

void RxxPowellMethod::CalculateNoduleCorrespondence(double lfTx, double lfTy, double lfTz, double lfSx, double lfSy, double lfSz, double lfDegX, double lfDegY, double lfDegZ)
{
	RxVolumeData* pRefVol = RxGetVolumeData(0);	
	RxVolumeData* pFloatVol = RxGetVolumeData(1);	
	
	RxMatrix4D m_mxRegistration;

	double lfSHx = m_lfShearX;
	double lfSHy = m_lfShearY;
	double lfSHz = m_lfShearZ;

	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);
	RxMatrix4D Temp;

	m_mxRegistration.LoadIdentity();
	m_mxRegistration.Translate(-g_FloatCenX, -g_FloatCenY, -g_FloatCenZ + g_iFloatBinaryTransZ);
	m_mxRegistration.Scale(lfSx, lfSy, lfSz);
	m_mxRegistration.Rotate(0, lfDegX);
	m_mxRegistration.Rotate(1, lfDegY);
	m_mxRegistration.Rotate(2, lfDegZ);

	Temp = m_mxRegistration;
	m_mxRegistration = SH * Temp;

	m_mxRegistration.Translate(g_RefCenX + lfTx, g_RefCenY + lfTy, g_RefCenZ + lfTz - g_iRefBinaryTransZ);

	double MinDist, Dist;

	for (int numRef = 0; numRef < pRefVol->m_nNumNodule; numRef++)
	{
		RxVect4D registeredNodule = m_mxRegistration * RxVect4D(pFloatVol->m_nNoduleCenterX[FloatNoduleNumber(numRef)], 
			pFloatVol->m_nNoduleCenterY[FloatNoduleNumber(numRef)], 
			pFloatVol->m_nNoduleCenterZ[FloatNoduleNumber(numRef)], 1);		

		MinDist = sqrt(pow(pRefVol->m_nNoduleCenterX[numRef] - registeredNodule.m[0], 2.0) + 
			pow(pRefVol->m_nNoduleCenterY[numRef] - registeredNodule.m[1], 2.0) + 
			pow(pRefVol->m_nNoduleCenterZ[numRef] - registeredNodule.m[2], 2.0));

		for (int num = 0; num < pRefVol->m_nNumNodule; num++)
		{
			RxVect4D registeredNodule = m_mxRegistration * RxVect4D(pFloatVol->m_nNoduleCenterX[num], pFloatVol->m_nNoduleCenterY[num], pFloatVol->m_nNoduleCenterZ[num], 1);		

			Dist = sqrt(pow(pRefVol->m_nNoduleCenterX[numRef] - registeredNodule.m[0], 2.0) + 
				pow(pRefVol->m_nNoduleCenterY[numRef] - registeredNodule.m[1], 2.0) + 
				pow(pRefVol->m_nNoduleCenterZ[numRef] - registeredNodule.m[2], 2.0));

			if (num != numRef && Dist < MinDist)
			{
				g_outFile << '\n' << numRef << "*****ERROR*****" << num;
				m_arrNoduleMapping[numRef] = num;
				MinDist = Dist;
			}
		}
	}
}

int RxxPowellMethod::FloatNoduleNumber(int RefNoduleNumber)
{
	return m_arrNoduleMapping[RefNoduleNumber];
}

void RxxPowellMethod::KeepParameterForNR(double lfScaleVectorX, double lfScaleVectorY, double lfScaleVectorZ, 
											 double lfTransX, double lfTransY, double lfTransZ, 
											 double RotX, double RotY, double RotZ)
{	
	m_lfTransX = lfTransX;
	m_lfTransY = lfTransY;
	m_lfTransZ = lfTransZ;

	m_lfOptimizedScaleX = lfScaleVectorX;
	m_lfOptimizedScaleY = lfScaleVectorY;
	m_lfOptimizedScaleZ = lfScaleVectorZ;	

	//degree
	m_lfRotX = RotX;
	m_lfRotY = RotY;
	m_lfRotZ = RotZ;

	//radian
	m_lfRadX = m_lfRotX*DEGREETORADIAN;
	m_lfRadY = m_lfRotY*DEGREETORADIAN;
	m_lfRadZ = m_lfRotZ*DEGREETORADIAN;

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