mutualinfosm.cpp

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

	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];
		BYTE nRefValue = m_pbyRefVol[nZRef*m_pRefInfo->nVolX*m_pRefInfo->nVolY+nYRef*m_pRefInfo->nVolX+nXRef];
		lfSum=lfSum+nRefValue*w[nIndex];		
	}
	BYTE nInterpolatedRefValue = (BYTE)(lfSum+0.5);
	m_lfJointHistoSize++;
	float fVote = 1.0f;
	UpdateJointHistogram(nFltValue, nInterpolatedRefValue, fVote);
}
*/

void LHMutualInfoSM::PartialVolumeInterpolation(double fTransformedX, double fTransformedY, double fTransformedZ, BYTE nFltValue)
{	
	int ix, iy, iz;	
//	ix = FloatToInt(fTransformedX);
//	iy = FloatToInt(fTransformedY);
//	iz = FloatToInt(fTransformedZ);
	ix = (int)(fTransformedX);
	iy = (int)(fTransformedY);
	iz = (int)(fTransformedZ);
	double xFract = fTransformedX - ix;
	double yFract = fTransformedY - iy;
	double zFract = fTransformedZ - iz;
	double xFractComplement = (1.0f-xFract);
	double yFractComplement = (1.0f-yFract);
	double zFractComplement = (1.0f-zFract);	

	// Direction : Bottom to Top
	//compute weight
	//Top
	// w2 w1	n8	n7
	// w3 w4	n5	n6

	//Bottom
	// w6 w5	n4	n3
	// w7 w8	n1	n2		
	m_fWeight[0] = float(xFractComplement * yFractComplement*zFractComplement);
	m_fWeight[1] = float(xFract * yFractComplement*zFractComplement);
	m_fWeight[2] = float(xFract * yFract*zFractComplement);
	m_fWeight[3] = float(xFractComplement * yFract* zFractComplement);
	m_fWeight[4] = float(xFractComplement * yFractComplement*zFract);
	m_fWeight[5] = float(xFract * yFractComplement*zFract);
	m_fWeight[6] = float(xFract * yFract*zFract);
	m_fWeight[7] = float(xFractComplement * yFract* zFract);

	int iRslice = m_pRefInfo->nVolX*m_pRefInfo->nVolY;
	for(int nIndex = 0 ; nIndex<8 ; nIndex++){
		int nXRef = ix + OFFSET_X[nIndex];
		int nYRef = iy + OFFSET_Y[nIndex];
		int nZRef = iz + OFFSET_Z[nIndex];
		if(m_fWeight[nIndex]==0/* || !IsValid(nXRef, nYRef, nZRef)*/)	continue;
		BYTE nRefValue = m_pbyRefVol[nZRef*iRslice+nYRef*m_pRefInfo->nVolX+nXRef];
//		m_lfJointHistoSize+=m_fWeight[nIndex];
		UpdateJointHistogram(nFltValue, nRefValue, m_fWeight[nIndex]);
	}
//	m_lfJointHistoSize++;
}

/*
void LHMutualInfoSM::PartialVolumeInterpolation(float fTransformedX, float fTransformedY, float fTransformedZ, BYTE nFltValue)
{	
	int ix, iy, iz;	

//	ix = FloatToInt(fTransformedX);
//	iy = FloatToInt(fTransformedY);
//	iz = FloatToInt(fTransformedZ);
	ix = (int)(fTransformedX);
	iy = (int)(fTransformedY);
	iz = (int)(fTransformedZ);
	float xFract = fTransformedX - ix;
	float yFract = fTransformedY - iy;
	float zFract = fTransformedZ - iz;
	float xFractComplement = (1.0f-xFract);
	float yFractComplement = (1.0f-yFract);
	float zFractComplement = (1.0f-zFract);	

	// Direction : Bottom to Top
	//compute weight
	//Top
	// w2 w1	n8	n7
	// w3 w4	n5	n6

	//Bottom
	// w6 w5	n4	n3
	// w7 w8	n1	n2		
	m_fWeight[0] = xFractComplement * yFractComplement*zFractComplement;
	m_fWeight[1] = xFract * yFractComplement*zFractComplement;
	m_fWeight[2] = xFract * yFract*zFractComplement;
	m_fWeight[3] = xFractComplement * yFract* zFractComplement;
	m_fWeight[4] = xFractComplement * yFractComplement*zFract;
	m_fWeight[5] = xFract * yFractComplement*zFract;
	m_fWeight[6] = xFract * yFract*zFract;
	m_fWeight[7] = xFractComplement * yFract* zFract;

	int iRslice = m_pRefInfo->nVolX*m_pRefInfo->nVolY;
	for(int nIndex = 0 ; nIndex<8 ; nIndex++){
		int nXRef = ix + OFFSET_X[nIndex];
		int nYRef = iy + OFFSET_Y[nIndex];
		int nZRef = iz + OFFSET_Z[nIndex];
		if(m_fWeight[nIndex]==0 || IsValid(nXRef, nYRef, nZRef)==FALSE)	continue;
		BYTE nRefValue = m_pbyRefVol[nZRef*iRslice+nYRef*m_pRefInfo->nVolX+nXRef];
//		if(nRefValue<0 || nRefValue>=m_iJointHistoY || nFltValue<0 || nFltValue>=m_iJointHistoX) continue;
		m_lfJointHistoSize+=m_fWeight[nIndex];
		UpdateJointHistogram(nFltValue, nRefValue, m_fWeight[nIndex]);
	}
}
*/

/*
void LHMutualInfoSM::PartialVolumeInterpolation(double lfTransformedX, double lfTransformedY, double lfTransformedZ, BYTE 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	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];
		BYTE nRefValue = m_pbyRefVol[nZRef*m_pRefInfo->nVolX*m_pRefInfo->nVolY+nYRef*m_pRefInfo->nVolX+nXRef];

		m_lfJointHistoSize+=w[nIndex];
		UpdateJointHistogram(nFltValue, nRefValue, w[nIndex]);
	}
}
*/

#include<direct.h>
BOOL LHMutualInfoSM::SaveResultsToFile(CString strResult)
{
	CFile file;
	CString strPath = _T("D:\\Storage\\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 LHMutualInfoSM::SaveResultsToFile(CString strFile, CString strResult)
{	
	CFile file;
	CString strPath = _T("D:\\Storage\\MIresults\\");
	CString strFilename = strPath+strFile+_T(".xls");
	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;
}

RxMatrix4D LHMutualInfoSM::MakeTransformMatrix(_TransformParameters _OptTransformPar)
{
	double lfFltCx = m_pFltInfo->nCenterX;
	double lfFltCy = m_pFltInfo->nCenterY;
	double lfFltCz = m_pFltInfo->nCenterZ;

	double lfRefCx = m_pRefInfo->nCenterX;
	double lfRefCy = m_pRefInfo->nCenterY;
	double lfRefCz = m_pRefInfo->nCenterZ;
	
	RxMatrix4D matTR;
	matTR.LoadIdentity();
	matTR.Translate(-lfFltCx, -lfFltCy, -lfFltCz);
	matTR.Scale(_OptTransformPar.lfSx, _OptTransformPar.lfSy, m_lfRatioRefZ*_OptTransformPar.lfSz);
	matTR.Rotate(0, _OptTransformPar.lfRotX);
	matTR.Rotate(1, _OptTransformPar.lfRotY);
	matTR.Rotate(2, _OptTransformPar.lfRotZ);		
	matTR.Translate(_OptTransformPar.lfTx, _OptTransformPar.lfTy, _OptTransformPar.lfTz);	
	matTR.Scale(1.0f,1.0f, 1.0f/m_lfRatioRefZ);
	matTR.Translate(lfRefCx, lfRefCy, lfRefCz);	

	return matTR;
}

BOOL LHMutualInfoSM::ComputeRMSE(CString strFilename, int nCnt)
{
	double Tx = 10.5;
	double Ty = 9.4;
	double Tz = -10.2;
	double Rx = -10.8;
	double Ry = -9.2;
	double Rz = 11.3;

	double lfMSE = 				
		(
		(m_OptTransformPar.lfTx-Tx)*(m_OptTransformPar.lfTx-Tx)
		+(m_OptTransformPar.lfTy-Ty)*(m_OptTransformPar.lfTy-Ty)
		+(m_OptTransformPar.lfTz-Tz)*(m_OptTransformPar.lfTz-Tz)
		+(m_OptTransformPar.lfRotX-Rx)*(m_OptTransformPar.lfRotX-Rx)
		+(m_OptTransformPar.lfRotY-Ry)*(m_OptTransformPar.lfRotY-Ry)
		+(m_OptTransformPar.lfRotZ-Rz)*(m_OptTransformPar.lfRotZ-Rz)
		)/6;

	double lfRMSE = sqrt(lfMSE);
	CString strResult;
	strResult.Format(_T("%d\t%lf\n"), nCnt, lfRMSE);
	SaveResultsToFile(strFilename, strResult);

	return TRUE;
}

BYTE* LHMutualInfoSM::ConversionData(unsigned short *pusVol, int iX, int iY, int iZ, int& nMaxSize)
{	
	int iSize = iX*iY*iZ;
	BYTE* pbyVol = new BYTE[iSize];
	int itemp;

	int iMax =0;
	for(int i=0 ; i<iSize ; i++){
		if(iMax<pusVol[i])	iMax = pusVol[i];
	}

//	LHRegData RegData;
//	int iMax = RegData.GetMaxValueUsingMMX(pusVol, iSize);

	int iBin = 255; //bin size 256*256
//	int iBin = 127; //bin size 128*128

	if(iMax<=iBin){
		for(int i=0 ; i<iSize ; i++){
			pbyVol[i] = pusVol[i];
		}
		nMaxSize = iMax+1;
	}
	else{
		double lfRatio = (double)iMax/(double)iBin;
		for(int i=0 ; i<iSize ; i++){
			itemp = (int)((double)(pusVol[i])/lfRatio+0.5);
			if(itemp<=iBin)	pbyVol[i] = itemp;
			else pbyVol[i] = iBin;			
		}
		nMaxSize = (int)((double)(iMax)/lfRatio+0.5)+1;
	}
	return pbyVol;	
}

void LHMutualInfoSM::GetOptTransformPar(double &lfSx, double &lfSy, double &lfSz, 										 
										double &lfRx, double &lfRy, double &lfRz,
										double &lfTx, double &lfTy, double &lfTz)
{
	lfSx	= m_OptTransformPar.lfSx;
	lfSy	= m_OptTransformPar.lfSy;	
	lfSz	= m_OptTransformPar.lfSz;	
	lfRx	= m_OptTransformPar.lfRotX;
	lfRy	= m_OptTransformPar.lfRotY;
	lfRz	= m_OptTransformPar.lfRotZ;
	lfTx	= m_OptTransformPar.lfTx;	
	lfTy	= m_OptTransformPar.lfTy;	
	lfTz	= m_OptTransformPar.lfTz;	
}

BOOL LHMutualInfoSM::IsValid(int nX, int nY, int nZ)
{
	if(nX>=0 && nX<m_pRefInfo->nVolX &&
		nY>=0 && nY<m_pRefInfo->nVolY && 
		nZ>=0 && nZ<m_pRefInfo->nVolZ)
		return TRUE;
	
	return FALSE;
}