rxxpowellmethod.cpp

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

	for(double lfThetaY=m_lfThetaY-lfIncR ; ; lfThetaY -=lfIncR){
		MI = RxxMutualInfoAgent::ComputeMI(m_lfGlobalScaleX, m_lfGlobalScaleY, m_lfGlobalScaleZ, m_lfTranslateXofFlt, m_lfTranslateYofFlt, m_lfTranslateZofFlt, m_lfThetaX, lfThetaY, m_lfThetaZ);
//		TRACE(_T("Scale=(%lf,%lf,%lf), Translate=(%lf,%lf,%lf), Theta=(%lf,%lf,%lf), MI=%d\n"),m_lfGlobalScaleX, m_lfGlobalScaleY, m_lfGlobalScaleZ,m_lfTranslateXofFlt, m_lfTranslateYofFlt, m_lfTranslateZofFlt, m_lfThetaX, lfThetaY, m_lfThetaZ, MI);
		if(MI>lfPresentMI){
			RxxMutualInfoAgent::KeepParameter(m_lfGlobalScaleX, m_lfGlobalScaleY, m_lfGlobalScaleZ, m_lfTranslateXofFlt, m_lfTranslateYofFlt, m_lfTranslateZofFlt, m_lfThetaX, lfThetaY, m_lfThetaZ);
			lfPresentMI = MI;
			nIsThetaCOUNTERState = EMPTYCOUNTER;
		}
		else if(nIsThetaCOUNTERState!=FULLCOUNTER)	nIsThetaCOUNTERState++;		
		else	break;		
	}
	return lfPresentMI;
}

double RxxPowellMethod::CCWRotZvector(double lfPresentMI, double lfIncR)
{
	double MI =0.0;
	int nIsThetaCOUNTERState = EMPTYCOUNTER;
	for(double lfThetaZ=m_lfThetaZ-lfIncR ; ; lfThetaZ -=lfIncR){
		MI = RxxMutualInfoAgent::ComputeMI(m_lfGlobalScaleX, m_lfGlobalScaleY, m_lfGlobalScaleZ, m_lfTranslateXofFlt, m_lfTranslateYofFlt, m_lfTranslateZofFlt, m_lfThetaX, m_lfThetaY, lfThetaZ);
//		TRACE(_T("Scale=(%lf,%lf,%lf), Translate=(%lf,%lf,%lf), Theta=(%lf,%lf,%lf), MI=%d\n"),m_lfGlobalScaleX, m_lfGlobalScaleY, m_lfGlobalScaleZ,m_lfTranslateXofFlt, m_lfTranslateYofFlt, m_lfTranslateZofFlt, m_lfThetaX, m_lfThetaY, lfThetaZ, MI);
		if(MI>lfPresentMI){
			RxxMutualInfoAgent::KeepParameter(m_lfGlobalScaleX, m_lfGlobalScaleY, m_lfGlobalScaleZ, m_lfTranslateXofFlt, m_lfTranslateYofFlt, m_lfTranslateZofFlt, m_lfThetaX, m_lfThetaY, lfThetaZ);
			lfPresentMI = MI;
			nIsThetaCOUNTERState = EMPTYCOUNTER;
		}
		else if(nIsThetaCOUNTERState!=FULLCOUNTER)	nIsThetaCOUNTERState++;		
		else	break;		
	}
	return lfPresentMI;	
}

void RxxPowellMethod::FindOptimalTransformation(short nDistanceMode)
{
	DWORD CurrentTime = GetTickCount();

	m_nCount = 0;

	int nIteration = 3;
	double lfTransRange = 40.0f;
	double lfRotRange = 6.0f; //Degree
	double lfAttenuation = 0.5;
	double lfStep = 1.0f;
	double lfSHMax = 0.01f, lfSHStep = 0.0005f;

	g_outFile << "\n* Optimization Parameters\n";
	g_outFile << "(nIteration, lfTransRange, lfRotRange, lfAttenuation, lfStep, lfSHMax, lfSHStep) = "
		<< nIteration << ", " << lfTransRange << ", " << lfRotRange << ", " << lfAttenuation << ", " <<  lfStep << ", " << lfSHMax << ", " << lfSHStep << "\n";

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

	double lfCC;
	double lfLocalOptimum, lfLocalOptimumDeg;

	for(int i=0 ; i<nIteration ; i++){
		//x trans
		for(double Tx=-lfTransRange/2.0f ; Tx<lfTransRange/2.0f ; Tx+=lfStep){
			m_nCount++;
			lfLocalOptimum = m_lfTransX;
			lfCC = AverageDistanceDifference(nDistanceMode, 1, m_lfBoundingBoxX + lfLocalOptimum + Tx);
			if(lfCC < m_lfMaxCC){
				m_lfMaxCC = lfCC;
				KeepParameter(1, lfLocalOptimum + Tx);
			}
		}
 		ControlProgress(pWndProgress);
		//y trans
		for(double Ty=-lfTransRange/2.0f ; Ty<lfTransRange/2.0f ; Ty+=lfStep){
			m_nCount++;
			lfLocalOptimum = m_lfTransY;
			lfCC = AverageDistanceDifference(nDistanceMode, 2, m_lfBoundingBoxY + lfLocalOptimum + Ty);
			if(lfCC< m_lfMaxCC){
				m_lfMaxCC = lfCC;
				KeepParameter(2, lfLocalOptimum + Ty);
			}
		}
		ControlProgress(pWndProgress);
		//z Rot		
		for(double Rz=-lfRotRange/2.0f ; Rz<lfRotRange/2.0f ; Rz+=lfStep){
			m_nCount++;
			lfLocalOptimum = m_lfRadZ;
			lfLocalOptimumDeg = m_lfRotZ;
			lfCC = AverageDistanceDifference(nDistanceMode, 6, lfLocalOptimum + Rz*DEGREETORADIAN);
			if(lfCC< m_lfMaxCC){
				m_lfMaxCC = lfCC;
				KeepParameter(6, lfLocalOptimumDeg + Rz);
			}		
		}
		ControlProgress(pWndProgress);
		//x Rot
		for(double Rx=-lfRotRange/2.0f ; Rx<lfRotRange/2.0f ; Rx+=lfStep){
			m_nCount++;
			lfLocalOptimum = m_lfRadX;
			lfLocalOptimumDeg = m_lfRotX;
			lfCC = AverageDistanceDifference(nDistanceMode, 4, lfLocalOptimum + Rx*DEGREETORADIAN);
			if(lfCC< m_lfMaxCC){
				m_lfMaxCC = lfCC;
				KeepParameter(4, lfLocalOptimumDeg + Rx);
			}		
		}
		ControlProgress(pWndProgress);
		//y Rot
		for(double Ry=-lfRotRange/2.0f ; Ry<lfRotRange/2.0f ; Ry+=lfStep){
			m_nCount++;
			lfLocalOptimum = m_lfRadY;
			lfLocalOptimumDeg = m_lfRotY;
			lfCC = AverageDistanceDifference(nDistanceMode, 5, lfLocalOptimum + Ry*DEGREETORADIAN);
			if(lfCC< m_lfMaxCC){
				m_lfMaxCC = lfCC;
				KeepParameter(5, lfLocalOptimumDeg + Ry);
			}		
		}
		ControlProgress(pWndProgress);
		//z trans
		for(double Tz=-lfTransRange/2.0f ; Tz<lfTransRange/2.0f ; Tz+=lfStep){
			m_nCount++;
			lfLocalOptimum = m_lfTransZ;
			lfCC = AverageDistanceDifference(nDistanceMode, 3, m_lfBoundingBoxZ + lfLocalOptimum + Tz);
			if(lfCC< m_lfMaxCC){
				m_lfMaxCC = lfCC;
				KeepParameter(3, lfLocalOptimum + Tz);
			}
		}
		ControlProgress(pWndProgress);
		lfTransRange *= lfAttenuation;
		lfRotRange *= lfAttenuation;
		lfStep *= lfAttenuation;
	}

	DWORD ElapsedTimeRigid = GetTickCount() - CurrentTime;	

	if (g_nTransMode == 1)
	{
		for (double SHx = -lfSHMax ; SHx < lfSHMax ; SHx += lfSHStep)
		{
			m_nCount++;
			lfLocalOptimum = m_lfShearX;
			lfCC = AverageDistanceDifference(nDistanceMode, 7, lfLocalOptimum + SHx);
			if(lfCC< m_lfMaxCC){
				m_lfMaxCC = lfCC;
				KeepParameter(7, lfLocalOptimum + SHx);
			}
		}

		for (double SHy = -lfSHMax ; SHy < lfSHMax ; SHy += lfSHStep)
		{
			m_nCount++;
			lfLocalOptimum = m_lfShearY;
			lfCC = AverageDistanceDifference(nDistanceMode, 8, lfLocalOptimum + SHy);
			if(lfCC< m_lfMaxCC){
				m_lfMaxCC = lfCC;
				KeepParameter(8, lfLocalOptimum + SHy);
			}
		}

		for (double SHz = -lfSHMax ; SHz < lfSHMax ; SHz += lfSHStep)
		{
			m_nCount++;
			lfLocalOptimum = m_lfShearZ;
			lfCC = AverageDistanceDifference(nDistanceMode, 9, lfLocalOptimum + SHz);
			if(lfCC< m_lfMaxCC){
				m_lfMaxCC = lfCC;
				KeepParameter(9, lfLocalOptimum + SHz);
			}
		}
	}

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

	DWORD ElapsedTime = GetTickCount() - CurrentTime;	

	g_outFile << "\n* 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_lfBoundingBoxZ + m_lfTransZ << " )\n";
	g_outFile << "RotationDeg( " << m_lfRotX << ", " << m_lfRotY << ", " << m_lfRotZ << " )\n";
	g_outFile << "Shear( " << m_lfShearX << ", " << m_lfShearY << ", " << m_lfShearZ << " )\n";
	g_outFile << "MinAverageDistance = " << m_lfMaxCC << '\n';
	g_outFile << "Registration time = " << ElapsedTimeRigid / 1000.0 << '\t' << ElapsedTime / 1000.0 << '\n';

	m_lfRegistrationTime = ElapsedTime / 1000.0;
}

void RxxPowellMethod::FindOptimalTransformationWithSDM(void)
{
	DWORD CurrentTime = GetTickCount();

	m_nCount = 0;

	int nIteration = 3;
	double lfTransRange = 40.0f;
	double lfRotRange = 6.0f; //Degree
	double lfAttenuation = 0.5;
	double lfStep = 1.0f;
	double lfSHMax = 0.01f, lfSHStep = 0.0005f;

//  faster debug
/*	int nIteration = 1;
	double lfTransRange = 1.0f;
	double lfRotRange = 1.0f; //Degree
	double lfAttenuation = 0.5;
	double lfStep = 1.0f;
	double lfSHMax = 0.01f, lfSHStep = 0.01f;*/

	g_outFile << "* Optimization parameters\n";
	g_outFile << "\n* Optimization Parameters\n";
	g_outFile << "(nIteration, lfTransRange, lfRotRange, lfAttenuation, lfStep, lfSHMax, lfSHStep) = "
		<< nIteration << ", " << lfTransRange << ", " << lfRotRange << ", " << lfAttenuation << ", " <<  lfStep << ", " << lfSHMax << ", " << lfSHStep << "\n";

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

	double lfCC, lfLocalOptimum, lfLocalOptimumDeg;

	for(int i=0 ; i<nIteration ; i++){
		//x trans
		for(double Tx=-lfTransRange/2.0f ; Tx<lfTransRange/2.0f ; Tx+=lfStep){
			lfLocalOptimum = m_lfTransX;
			if (fabs(m_lfTransX+Tx) <= 10.0)	{
				m_nCount++;
				lfCC = AverageDistanceDifference(1, 1, m_lfBoundingBoxX + lfLocalOptimum + Tx);
				if(lfCC < m_lfMaxCC){
					m_lfMaxCC = lfCC;
					KeepParameter(1, lfLocalOptimum + Tx);
				}
			}
		}
 		ControlProgress(pWndProgress);
		//y trans
		for(double Ty=-lfTransRange/2.0f ; Ty<lfTransRange/2.0f ; Ty+=lfStep){
			lfLocalOptimum = m_lfTransY;
			if (fabs(m_lfTransY+Ty) <= 10.0) {
				m_nCount++;
				lfCC = AverageDistanceDifference(1, 2, m_lfBoundingBoxY + lfLocalOptimum + Ty);
				if(lfCC< m_lfMaxCC){
					m_lfMaxCC = lfCC;
					KeepParameter(2, lfLocalOptimum + Ty);
				}
			}
		}
		ControlProgress(pWndProgress);
		//z Rot		
		for(double Rz=-lfRotRange/2.0f ; Rz<lfRotRange/2.0f ; Rz+=lfStep){
			m_nCount++;
			lfLocalOptimum = m_lfRadZ;
			lfLocalOptimumDeg = m_lfRotZ;
			lfCC = AverageDistanceDifference(1, 6, lfLocalOptimum + Rz*DEGREETORADIAN);
			if(lfCC< m_lfMaxCC){
				m_lfMaxCC = lfCC;
				KeepParameter(6, lfLocalOptimumDeg + Rz);
			}		
		}
		ControlProgress(pWndProgress);
		//x Rot
		for(double Rx=-lfRotRange/2.0f ; Rx<lfRotRange/2.0f ; Rx+=lfStep){
			m_nCount++;
			lfLocalOptimum = m_lfRadX;
			lfLocalOptimumDeg = m_lfRotX;
			lfCC = AverageDistanceDifference(1, 4, lfLocalOptimum + Rx*DEGREETORADIAN);
			if(lfCC< m_lfMaxCC){
				m_lfMaxCC = lfCC;
				KeepParameter(4, lfLocalOptimumDeg + Rx);
			}		
		}
		ControlProgress(pWndProgress);
		//y Rot
		for(double Ry=-lfRotRange/2.0f ; Ry<lfRotRange/2.0f ; Ry+=lfStep){
			m_nCount++;
			lfLocalOptimum = m_lfRadY;
			lfLocalOptimumDeg = m_lfRotY;
			lfCC = AverageDistanceDifference(1, 5, lfLocalOptimum + Ry*DEGREETORADIAN);
			if(lfCC< m_lfMaxCC){
				m_lfMaxCC = lfCC;
				KeepParameter(5, lfLocalOptimumDeg + Ry);
			}		
		}
		ControlProgress(pWndProgress);
		//z trans
		for(double Tz=-lfTransRange/2.0f ; Tz<lfTransRange/2.0f ; Tz+=lfStep){
			lfLocalOptimum = m_lfTransZ;
			if (fabs(m_lfTransZ+Tz) <= 10.0) {
				m_nCount++;
				lfCC = AverageDistanceDifference(1, 3, m_lfBoundingBoxZ + lfLocalOptimum + Tz);
				if(lfCC< m_lfMaxCC){
					m_lfMaxCC = lfCC;
					KeepParameter(3, lfLocalOptimum + Tz);
				}
			}
		}
		ControlProgress(pWndProgress);
		lfTransRange *= lfAttenuation;
		lfRotRange *= lfAttenuation;
		lfStep *= lfAttenuation;
	}

	if (g_nTransMode == 1)
	{
		for (double SHx = -lfSHMax ; SHx < lfSHMax ; SHx += lfSHStep)
		{
			m_nCount++;
			lfLocalOptimum = m_lfShearX;
			lfCC = AverageDistanceDifference(1, 7, lfLocalOptimum + SHx);
			if(lfCC< m_lfMaxCC){
				m_lfMaxCC = lfCC;
				KeepParameter(7, lfLocalOptimum + SHx);
			}
		}

		for (double SHy = -lfSHMax ; SHy < lfSHMax ; SHy += lfSHStep)
		{
			m_nCount++;
			lfLocalOptimum = m_lfShearY;
			lfCC = AverageDistanceDifference(1, 8, lfLocalOptimum + SHy);
			if(lfCC< m_lfMaxCC){
				m_lfMaxCC = lfCC;
				KeepParameter(8, lfLocalOptimum + SHy);
			}
		}

		for (double SHz = -lfSHMax ; SHz < lfSHMax ; SHz += lfSHStep)
		{
			m_nCount++;
			lfLocalOptimum = m_lfShearZ;
			lfCC = AverageDistanceDifference(1, 9, lfLocalOptimum + SHz);
			if(lfCC< m_lfMaxCC){
				m_lfMaxCC = lfCC;
				KeepParameter(9, lfLocalOptimum + SHz);
			}
		}
	}

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

	DWORD ElapsedTime = GetTickCount() - CurrentTime;	

	m_lfTempCC = AverageDistanceDifference(1, 0, 0);
	m_lfSurfaceAccuracyIR = AverageDistanceDifference(1, 10, 0);

	g_outFile << "\n* 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_lfBoundingBoxZ + m_lfTransZ << " )\n";
	g_outFile << "RotationDeg( " << m_lfRotX << ", " << m_lfRotY << ", " << m_lfRotZ << " )\n";
	g_outFile << "Shear( " << m_lfShearX << ", " << m_lfShearY << ", " << m_lfShearZ << " )\n";
	g_outFile << "MinAverageDistance = " << m_lfMaxCC << '\n';
	g_outFile << "Registration time = " << ElapsedTime / 1000.0 << '\n';
	m_lfRegistrationTime = ElapsedTime / 1000.0;
}

double RxxPowellMethod::AverageDistanceDifference(short nDistanceMode, short nParameterNum, double lfValue)
{
	double lfTx, lfTy, lfTz, lfRadX, lfRadY, lfRadZ, lfSHx, lfSHy, lfSHz, lfSx, lfSy, lfSz;

	lfTx = m_lfBoundingBoxX + m_lfTransX;
	lfTy = m_lfBoundingBoxY + m_lfTransY;
	lfTz = m_lfBoundingBoxZ + m_lfTransZ;
	lfRadX = m_lfRadX;
	lfRadY = m_lfRadY;
	lfRadZ = m_lfRadZ;
	lfSx = m_lfOptimizedScaleX;
	lfSy = m_lfOptimizedScaleY;