rxxmutualinfoagent.cpp

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

		pFloatGaussian[nFloatNextZ+nXofFlt+1]- pFloatGaussian[nFloatNextY+1])/4.0;

	RefP = (pRefGaussian[nRefNextX]- pRefGaussian[nRefNowIdx] + 
		pRefGaussian[nRefNextY+1]- pRefGaussian[nRefNextY] +
		pRefGaussian[nRefNextZ+1]- pRefGaussian[nRefNextZ] +
		pRefGaussian[nRefNextZ+nXofRef+1]- pRefGaussian[nRefNextZ+nXofRef])/4.0;

	RefQ = (pRefGaussian[nRefNextY]- pRefGaussian[nRefNowIdx] + 
		pRefGaussian[nRefNextY+1]- pRefGaussian[nRefNextX] +
		pRefGaussian[nRefNextZ+nXofRef]- pRefGaussian[nRefNextZ] +
		pRefGaussian[nRefNextZ+nXofRef+1]- pRefGaussian[nRefNextZ+1])/4.0;

	RefR = (pRefGaussian[nRefNextZ]- pRefGaussian[nRefNowIdx] + 
		pRefGaussian[nRefNextZ+1]- pRefGaussian[nRefNextX] +
		pRefGaussian[nRefNextZ+nXofRef]- pRefGaussian[nRefNextY] +
		pRefGaussian[nRefNextZ+nXofRef+1]- pRefGaussian[nRefNextY+1])/4.0;

	Float_Mag = sqrt(FltP*FltP + FltQ*FltQ + FltR*FltR);
	Ref_Mag = sqrt(RefP*RefP + RefQ*RefQ + RefR*RefR);

	if(Float_Mag !=0 && Ref_Mag !=0){
		double degree = (FltP*RefP+FltQ*RefQ+FltR*RefR)/(Float_Mag*Ref_Mag);

		double Radangle = (degree==-1) ? 3.141592:acos(degree);
		double w = (cos(2*Radangle)+1)/2.0;

		UpdateGradientInfo(w, Float_Mag, Ref_Mag);
	}	
}

void RxxMutualInfoAgent::NearestInterpolation(double lfTransformedX, double lfTransformedY, double lfTransformedZ, int nFltValue)
{
	int nTransX = (int)(lfTransformedX+0.5);
	int nTransY = (int)(lfTransformedY+0.5);
	int nTransZ = (int)(lfTransformedZ+0.5);			
	
	unsigned short nRefValue = (unsigned short)((double)GetRefValue(nTransX,nTransY,nTransZ)/16.0f);
	m_lfRegionInPointCnt++;	
	UpdateJointHistogram(nRefValue, nFltValue, 1);
}

void RxxMutualInfoAgent::TrilinearVolumeInterpolation(double lfTransformedX, double lfTransformedY, double lfTransformedZ, int nFltValue)
{
	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 lfSum=0.0;
	for(int nIndex =0 ; nIndex<8 ; nIndex++){
		int nXRef = nTransformedXInt + nOffsetX[nIndex];
		int nYRef = nTransformedYInt + nOffsetY[nIndex];
		int nZRef = nTransformedZInt + nOffsetZ[nIndex];
		unsigned short nRefValue = (unsigned short)((double)GetRefValue(nXRef,nYRef,nZRef)/16.0f);
		lfSum=lfSum+nRefValue*w[nIndex];		
	}
	unsigned short nInterpolatedRefValue = (unsigned short)(lfSum+0.5);
	m_lfRegionInPointCnt++;
	UpdateJointHistogram(nInterpolatedRefValue, nFltValue, 1);
}

void RxxMutualInfoAgent::PartialVolumeInterpolation(double lfTransformedX, double lfTransformedY, double lfTransformedZ, int nFltValue)
{
	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);
	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);
	for(int nIndex =0 ; nIndex<8 ; nIndex++){
		int nXRef = nTransformedXInt + nOffsetX[nIndex];
		int nYRef = nTransformedYInt + nOffsetY[nIndex];
		int nZRef = nTransformedZInt + nOffsetZ[nIndex];
		unsigned short nRefValue = (unsigned short)((double)(GetRefValue(nXRef,nYRef,nZRef))/16.0f);

		m_lfRegionInPointCnt+=w[nIndex];
		UpdateJointHistogram(nRefValue, nFltValue, w[nIndex]);
	}
}

BOOL RxxMutualInfoAgent::FindOptimalLoc(unsigned short *pFltVol, unsigned short *pRefVol)
{
	DWORD CurrentTime = GetTickCount();
	int nIteration = 10;
	RxProgressWnd*	pWndProgress;
	pWndProgress = new RxProgressWnd;
	int nSetRange = ((nIteration)*6);
	pWndProgress->GoModal();
	pWndProgress->SetWindowSize(0);
	pWndProgress->SetRange(0, nSetRange);
	pWndProgress->SetText(_T("Registering"));
	
	double lfStep = 1.0f;
	double lfTransRange = 30.0f;
	double lfRotRange = 30.0f; //Degree

	double lfAttenuation = 0.9;
	double lfMI;

	InitialJointHistogram();
	double lfLastTransX, lfLastTransY, lfLastTransZ, lfLastRotX, lfLastRotY, lfLastRotZ;
	CString strOut;

	for(int i=0 ; i<nIteration ; i++){
		lfLastTransX = m_lfTransX;
		lfLastTransY = m_lfTransY; 
		lfLastTransZ = m_lfTransZ;
		lfLastRotX = m_lfRotX;
		lfLastRotY = m_lfRotY;
		lfLastRotZ = m_lfRotZ;

		//x trans
		for(double Tx=-lfTransRange/2.0f ; Tx<=lfTransRange/2.0f ; Tx+=lfStep){
			if(Tx == 0 && i!=0)	continue;
			lfMI = VoxelBasedMI(pFltVol, pRefVol, lfLastTransX+Tx, m_lfTransY, m_lfTransZ, m_lfScaleX, m_lfScaleY, m_lfScaleZ, m_lfRotX, m_lfRotY, m_lfRotZ);
			strOut.Format(_T("Translate=(%lf,%lf,%lf), Theta=(%lf,%lf,%lf), MI=%lf\r\n"),lfLastTransX+Tx,m_lfTransY,m_lfTransZ, m_lfRotX, m_lfRotY, m_lfRotZ, lfMI);
			SaveResultsToFile(strOut);
//			TRACE(_T("Translate=(%lf,%lf,%lf), Theta=(%lf,%lf,%lf), MI=%lf\n"),lfLastTransX+Tx,m_lfTransY,m_lfTransZ, m_lfRotX, m_lfRotY, m_lfRotZ, lfMI);
			if(lfMI> m_lfMaxMI){
				KeepParameter(m_lfScaleX, m_lfScaleY, m_lfScaleZ, lfLastTransX+Tx, m_lfTransY, m_lfTransZ, m_lfRotX, m_lfRotY, m_lfRotZ);
				m_lfMaxMI = lfMI;
				TRACE(_T("Change\n"));				
				SaveResultsToFile(_T("Change\r\n"));
			}
		}
		if(!ControlProgress(pWndProgress))	break;
		//y trans
		for(double Ty=-lfTransRange/2.0f ; Ty<=lfTransRange/2.0f ; Ty+=lfStep){
			if(Ty == 0 && i!=0)	continue;
			lfMI = VoxelBasedMI(pFltVol, pRefVol, m_lfTransX, lfLastTransY+Ty, m_lfTransZ, m_lfScaleX, m_lfScaleY, m_lfScaleZ, m_lfRotX, m_lfRotY, m_lfRotZ);	
			strOut.Format(_T("Translate=(%lf,%lf,%lf), Theta=(%lf,%lf,%lf), MI=%lf\r\n"),m_lfTransX,lfLastTransY+Ty,m_lfTransZ, m_lfRotX, m_lfRotY, m_lfRotZ, lfMI);
			SaveResultsToFile(strOut);
//			TRACE(_T("Translate=(%lf,%lf,%lf), Theta=(%lf,%lf,%lf), MI=%lf\n"),m_lfTransX,lfLastTransY+Ty,m_lfTransZ, m_lfRotX, m_lfRotY, m_lfRotZ, lfMI);
			if(lfMI> m_lfMaxMI){
				KeepParameter(m_lfScaleX, m_lfScaleY, m_lfScaleZ, m_lfTransX, lfLastTransY+Ty, m_lfTransZ, m_lfRotX, m_lfRotY, m_lfRotZ);
				m_lfMaxMI = lfMI;
				TRACE(_T("Change\n"));
				SaveResultsToFile(_T("Change\r\n"));
			}
		}
		if(!ControlProgress(pWndProgress))	break;
		//z Rot		
		for(double Rz=-lfRotRange/2.0f ; Rz<=lfRotRange/2.0f ; Rz+=lfStep){
			if(Rz == 0 && i!=0)	continue;
			lfMI = VoxelBasedMI(pFltVol, pRefVol, m_lfTransX, m_lfTransY, m_lfTransZ, m_lfScaleX, m_lfScaleY, m_lfScaleZ, m_lfRotX, m_lfRotY, lfLastRotZ+Rz);
			strOut.Format(_T("Translate=(%lf,%lf,%lf), Theta=(%lf,%lf,%lf), MI=%lf\r\n"),m_lfTransX,m_lfTransY,m_lfTransZ, m_lfRotX, m_lfRotY, lfLastRotZ+Rz, lfMI);
			SaveResultsToFile(strOut);

//			TRACE(_T("Translate=(%lf,%lf,%lf), Theta=(%lf,%lf,%lf), MI=%lf\n"),m_lfTransX,m_lfTransY,m_lfTransZ, m_lfRotX, m_lfRotY, lfLastRotZ+Rz, lfMI);
			if(lfMI> m_lfMaxMI){
				KeepParameter(m_lfScaleX, m_lfScaleY, m_lfScaleZ, m_lfTransX, m_lfTransY, m_lfTransZ, m_lfRotX, m_lfRotY, lfLastRotZ+Rz);
				m_lfMaxMI = lfMI;
				TRACE(_T("Change\n"));
				SaveResultsToFile(_T("Change\r\n"));
			}		
		}
		if(!ControlProgress(pWndProgress))	break;
		//x Rot
		for(double Rx=-lfRotRange/2.0f ; Rx<=lfRotRange/2.0f ; Rx+=lfStep){
			if(Rx == 0 && i!=0)	continue;
			lfMI = VoxelBasedMI(pFltVol, pRefVol, m_lfTransX, m_lfTransY, m_lfTransZ, m_lfScaleX, m_lfScaleY, m_lfScaleZ, lfLastRotX+Rx, m_lfRotY, m_lfRotZ);
			strOut.Format(_T("Translate=(%lf,%lf,%lf), Theta=(%lf,%lf,%lf), MI=%lf\r\n"),m_lfTransX,m_lfTransY,m_lfTransZ, Rx+lfLastRotX, m_lfRotY, m_lfRotZ, lfMI);
			SaveResultsToFile(strOut);
//			TRACE(_T("Translate=(%lf,%lf,%lf), Theta=(%lf,%lf,%lf), MI=%lf\n"),m_lfTransX,m_lfTransY,m_lfTransZ, Rx+lfLastRotX, m_lfRotY, m_lfRotZ, lfMI);
			if(lfMI> m_lfMaxMI){
				KeepParameter(m_lfScaleX, m_lfScaleY, m_lfScaleZ, m_lfTransX, m_lfTransY, m_lfTransZ, Rx+lfLastRotX, m_lfRotY, m_lfRotZ);
				m_lfMaxMI = lfMI;
				TRACE(_T("Change\n"));
				SaveResultsToFile(_T("Change\r\n"));
			}
		}
		if(!ControlProgress(pWndProgress))	break;
		//y Rot
		for(double Ry=-lfRotRange/2.0f ; Ry<=lfRotRange/2.0f ; Ry+=lfStep){
			if(Ry == 0 && i!=0)	continue;
			lfMI = VoxelBasedMI(pFltVol, pRefVol, m_lfTransX, m_lfTransY, m_lfTransZ, m_lfScaleX, m_lfScaleY, m_lfScaleZ, m_lfRotX, lfLastRotY+Ry, m_lfRotZ);
			strOut.Format(_T("Translate=(%lf,%lf,%lf), Theta=(%lf,%lf,%lf), MI=%lf\r\n"),m_lfTransX,m_lfTransY,m_lfTransZ, m_lfRotX, lfLastRotY+Ry, m_lfRotZ, lfMI);
			SaveResultsToFile(strOut);
//			TRACE(_T("Translate=(%lf,%lf,%lf), Theta=(%lf,%lf,%lf), MI=%lf\n"),m_lfTransX,m_lfTransY,m_lfTransZ, m_lfRotX, lfLastRotY+Ry, m_lfRotZ, lfMI);
			if(lfMI> m_lfMaxMI){
				KeepParameter(m_lfScaleX, m_lfScaleY, m_lfScaleZ, m_lfTransX, m_lfTransY, m_lfTransZ, m_lfRotX, lfLastRotY+Ry, m_lfRotZ);
				m_lfMaxMI = lfMI;
				TRACE(_T("Change\n"));
				SaveResultsToFile(_T("Change\r\n"));
			}
		}
		if(!ControlProgress(pWndProgress))	break;
		//z trans
		for(double Tz=-lfTransRange/2.0f ; Tz<=lfTransRange/2.0f ; Tz+=lfStep){
			if(Tz == 0 && i!=0)	continue;
			lfMI = VoxelBasedMI(pFltVol, pRefVol, m_lfTransX, m_lfTransY, lfLastTransZ+Tz, m_lfScaleX, m_lfScaleY, m_lfScaleZ, m_lfRotX, m_lfRotY, m_lfRotZ);
			strOut.Format(_T("Translate=(%lf,%lf,%lf), Theta=(%lf,%lf,%lf), MI=%lf\r\n"), m_lfTransX,m_lfTransY, lfLastTransZ+Tz, m_lfRotX, m_lfRotY, m_lfRotZ, lfMI);
			SaveResultsToFile(strOut);
//			TRACE(_T("Translate=(%lf,%lf,%lf), Theta=(%lf,%lf,%lf), MI=%lf\n"), m_lfTransX,m_lfTransY, lfLastTransZ+Tz, m_lfRotX, m_lfRotY, m_lfRotZ, lfMI);
			if(lfMI> m_lfMaxMI){
				KeepParameter(m_lfScaleX, m_lfScaleY, m_lfScaleZ, m_lfTransX, m_lfTransY, lfLastTransZ+Tz, m_lfRotX, m_lfRotY, m_lfRotZ);
				m_lfMaxMI = lfMI;
				TRACE(_T("Change\n"));
				SaveResultsToFile(_T("Change\r\n"));
			}
		}
		if(!ControlProgress(pWndProgress))	break;
		lfTransRange *= lfAttenuation;
		lfRotRange *= lfAttenuation;
		lfStep *= lfAttenuation;
	}
	
	if (pWndProgress) {	delete pWndProgress;	pWndProgress = NULL;	}

	
	DWORD ElapsedTime = GetTickCount() - CurrentTime;	
	CString str;
	str.Format(_T("S = (%0.3lf, %0.3lf, %0.3lf)\r\nT = (%0.3lf, %0.3lf, %0.3lf)\r\nR(Rad) = (%0.3lf, %0.3lf, %0.3lf)\r\nR(Degree) = (%0.3lf, %0.3lf, %0.3lf)\r\n\r\nMaxMI = %lf\r\n\r\nTotal Processing Time = %dms"), 
				m_lfScaleX, m_lfScaleY, m_lfScaleZ, 
				m_lfTransX, m_lfTransY, m_lfTransZ,				
				m_lfRadX, m_lfRadY, m_lfRadZ, 
				m_lfRotX, m_lfRotY, m_lfRotZ, m_lfMaxMI, ElapsedTime);
	
	AfxMessageBox(str);
	
	return TRUE;
}

#include<direct.h>
BOOL RxxMutualInfoAgent::SaveResultsToFile(CString strResult)
{
	CFile file;
	CString strPath = _T("d:\\MIresults\\");	
	CString strFilename = strPath+_T("MIresult.txt");
	int nLength = strResult.GetLength();

	char* chResult= new char[nLength];
	for(int i=0 ; i<nLength ; i++){
		chResult[i] = strResult.GetAt(i);
	}

	CFileException fe;
	if(file.Open(strFilename, CFile::modeCreate|CFile::modeNoTruncate|CFile::modeReadWrite, &fe)){	
		file.SeekToEnd();
		file.Write(chResult, nLength);
	}
	else return FALSE;
	file.Close();
	delete[] chResult;
	return TRUE;
}

BOOL RxxMutualInfoAgent::ControlProgress(RxProgressWnd *pProWnd)
{
	if(pProWnd && pProWnd->Cancelled()){		
		return FALSE;
	}			
	else if(pProWnd){
		pProWnd->StepIt();
		pProWnd->PeekAndPump();		
	}
	return TRUE;
}

RxMatrix4D RxxMutualInfoAgent::MakeTransformMatrix(double lfTx, double lfTy, double lfTz,
										double lfSx, double lfSy, double lfSz,
										 double lfRotX, double lfRotY, double lfRotZ)
{
	double lfFltCx = m_pFloat->nCenterX;
	double lfFltCy = m_pFloat->nCenterY;
	double lfFltCz = m_pFloat->nCenterZ;

	double lfRefCx = m_pRefer->nCenterX;
	double lfRefCy = m_pRefer->nCenterY;
	double lfRefCz = m_pRefer->nCenterZ;

	RxMatrix4D Tr;
	Tr.LoadIdentity();
	Tr.Translate(-lfFltCx, -lfFltCy, -lfFltCz);
	Tr.Scale(lfSx, lfSy, m_lfRatioFltZ*lfSz);
	Tr.Rotate(0, lfRotX);
	Tr.Rotate(1, lfRotY);
	Tr.Rotate(2, lfRotZ);		
	Tr.Translate(lfTx, lfTy, lfTz);	
	Tr.Scale(1.0f,1.0f, 1/m_lfRatioRefZ);
	Tr.Translate(lfRefCx, lfRefCy, lfRefCz);
	
	return Tr;
}

double RxxMutualInfoAgent::VoxelBasedMI(unsigned short *pFltVol, unsigned short *pRefVol, double lfTx, double lfTy, double lfTz,
										 double lfSx, double lfSy, double lfSz, double lfRotX, double lfRotY, double lfRotZ)
{
	int nSampIntval = 1;
	int nTz, nTy, nTx;
	int nCnt=0;

	RxMatrix4D matTR = MakeTransformMatrix(lfTx, lfTy, lfTz, lfSx, lfSy, lfSz, lfRotX, lfRotY, lfRotZ);	
	
	_sz3DVol szLT;		szLT.x =0;	szLT.y =0;	szLT.z =0;
	_sz3DVol szRB;		szRB.x =m_pFloat->nVolX;	szRB.y =m_pFloat->nVolY;	szRB.z =m_pFloat->nVolZ;

	memset(m_plfJointHisto, 0x00, sizeof(double)*FLT_DEPTH*REF_DEPTH);

	int nFltSliceSize = m_pFloat->nVolX*m_pFloat->nVolY;

	RxxFilters filters;

	BOOL* pTemplate = new BOOL[m_pFloat->nVolX*m_pFloat->nVolY];
	ZeroMemory(pTemplate, sizeof(BOOL)*m_pFloat->nVolX*m_pFloat->nVolY);
	BOOL bIsTemplate = FALSE;

	for(int k = szLT.z ; k<szRB.z ; k+=nSampIntval){
		for(int j=szLT.y ; j<szRB.y ; j+=nSampIntval){
			for(int i=szLT.x ; i<szRB.x ; i+=nSampIntval){

//				if(!filters.IsInCircle2D(&i, &j, 60, &m_pFloat->nVolX, &m_pFloat->nVolY, pTemplate, bIsTemplate))	continue;				

				unsigned short nFltVol = (unsigned short)((double)(pFltVol[nFltSliceSize*k+j*m_pFloat->nVolX+i])/16.0f);

				RxVect4D v( (double)i, j, k, 1.0f);
				v = matTR*v;			
				
				if(v.m[0]>=0 && v.m[0]<m_pRefer->nVolX-1 &&
					v.m[1]>=0 && v.m[1]<m_pRefer->nVolY-1 &&
						v.m[2]>=0 && v.m[2]<m_pRefer->nVolZ-1){
					switch(m_nInterpID)
					{
						case PV :
							PartialVolumeInterpolation(v.m[0], v.m[1], v.m[2], nFltVol);							
							break;
						case TRI :							
							TrilinearVolumeInterpolation(v.m[0], v.m[1], v.m[2], nFltVol);							
							break;
						case NN :							
							NearestInterpolation(v.m[0], v.m[1], v.m[2], nFltVol);
							break;
					}		
				}
			}
		}
		bIsTemplate = TRUE;
	}

	delete[] pTemplate;
	
	double MI = ComputeMISub();	
//	double NewMI = m_lfGradientInfo*MI;
	double NewMI = MI;
//	double NewMI = m_lfGradientInfo;	

	//Reset
	m_lfRegionInPointCnt = 0.0;
	m_lfGradientInfo = 0.0;

	return NewMI;
}

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

	Cx = cos(Rad.x);
	Sx = sin(Rad.x);
	Cy = cos(Rad.y);
	Sy = sin(Rad.y);
	Cz = cos(Rad.z);
	Sz = sin(Rad.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;
}