rxxmutualinfoagent.cpp

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

// RxxMutualInfoAgent.cpp: implementation of the RxxMutualInfoAgent class.
//
//////////////////////////////////////////////////////////////////////
//	Title: Mutual Information based registration
//
//////////////////////////////////////////////////////////////////////
//
//	Author: Ho Lee
//	302-dong 314-2-ho Seoul National Univ.
//	Shinlim-dong, Kwanak-gu, Seoul, 151-742, Korea
//	(Email. [e-mail])
//	Update		: 
//	Last update	: 2004/06/07
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "RxxMutualInfoAgent.h"
#include "ProgressWnd.h"
#include "Matrix4D.h"
#include "Filters.h"

#include<math.h>

#define SAFE_DELETE_ARRAY(p) { if(p) { delete[] (p);   (p)=NULL; } }
#define REF_DEPTH			256
#define FLT_DEPTH			256
#define DEGREETORADIAN		3.141592653589793238462643383279502884197169399375105820975/180

#define PV	0
#define TRI	1
#define NN	2

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
RxxMutualInfoAgent::RxxMutualInfoAgent()
{	
	m_plfJointHisto		= NULL;
	m_plfProbF			= NULL;
	m_plfProbR			= NULL;
	m_lfMaxMI			= 0.0;
}

RxxMutualInfoAgent::RxxMutualInfoAgent(_MODAL_INFO* pFloat, _MODAL_INFO* pRefer, UINT nInterpID, double lfGlobalSx, double lfGlobalSy, double lfGlobalSz)
{
	m_plfJointHisto		= NULL;
	m_lfMaxMI			= 0.0;
	m_lfScaleX			= lfGlobalSx;
	m_lfScaleY			= lfGlobalSy;
	m_lfScaleZ			= lfGlobalSz;

	m_lfTransX			= 0.0;
	m_lfTransY			= 0.0;
	m_lfTransZ			= 0.0;
	m_lfRadX			= 0.0;		
	m_lfRadY			= 0.0;		
	m_lfRadZ			= 0.0;
	m_lfRotX			= 0.0;
	m_lfRotY			= 0.0;
	m_lfRotZ			= 0.0;
	
	m_pFloat			= pFloat;
	m_pRefer			= pRefer;

	m_lfRatioFltZ			= pFloat->lfVoxelZ / pFloat->lfVoxelX;
	m_lfRatioRefZ		= pRefer->lfVoxelZ / pRefer->lfVoxelX;		

	m_lfRegionInPointCnt = 0.0;

	m_lfGradientInfo	= 0.0;
	m_nLabelIdx			= 0;

	m_nInterpID			= nInterpID;

	m_plfProbF			= NULL;
	m_plfProbR			= NULL;	
}

RxxMutualInfoAgent::~RxxMutualInfoAgent()
{
	Empty();
}

void RxxMutualInfoAgent::Empty()
{	
	SAFE_DELETE_ARRAY(m_plfJointHisto);

	if(m_plfProbF){ delete[] m_plfProbF;	m_plfProbF = NULL; }
	if(m_plfProbR){ delete[] m_plfProbR;	m_plfProbR = NULL; }
}

double RxxMutualInfoAgent::ComputeMI(double ScaleX, double ScaleY, double ScaleZ, double fTransX, double fTransY, double fTransZ, double thetaX, double thetaY, double thetaZ)
{	
	InitialJointHistogram();
	ComputeTransformedVector(ScaleX, ScaleY, ScaleZ, fTransX, fTransY, fTransZ, thetaX, thetaY, thetaZ);
	
	double MI = ComputeMISub();	

//	double NewMI = m_lfGradientInfo*MI;
	double NewMI = MI;
//	double NewMI = m_lfGradientInfo;

	#ifdef _DEBUG
		TRACE(_T("Scale=(%lf,%lf,%lf), Translate=(%lf,%lf,%lf), Theta=(%lf,%lf,%lf), MI=%lf\n"),ScaleX, ScaleY, ScaleZ, fTransX, fTransY, fTransZ, thetaX, thetaY, thetaZ, NewMI);
	#endif // _DEBUG

	//Reset
	m_lfRegionInPointCnt = 0.0;
	m_lfGradientInfo = 0.0;

	return NewMI;
}

void RxxMutualInfoAgent::InitialJointHistogram()
{		
	m_plfJointHisto = new double[REF_DEPTH*FLT_DEPTH];
	m_plfProbF	= new double[FLT_DEPTH];
	m_plfProbR	= new double[REF_DEPTH];
}

void RxxMutualInfoAgent::ComputeTransformedVector(double lfScaleX, double lfScaleY, double lfScaleZ, double lfTransX, double lfTransY, double lfTransZ, double thetaX, double thetaY, double thetaZ)
{	
	unsigned short* pFltVol	= m_pFloat->pVol;
//	unsigned short* pFltVol	= m_pFloat->pGaussianVol;
	BYTE* pWT = m_pFloat->pWindowingTable;

	double lfCosThetaX = cos(thetaX);
	double lfSinThetaX = sin(thetaX);

	double lfCosThetaY = cos(thetaY);
	double lfSinThetaY = sin(thetaY);

	double lfCosThetaZ = cos(thetaZ);
	double lfSinThetaZ = sin(thetaZ);	

	double lfFltCenterX = m_pFloat->nCenterX;
	double lfFltCenterY = m_pFloat->nCenterY;
	double lfFltCenterZ = m_pFloat->nCenterZ;

	double lfRefCenterX = m_pRefer->nCenterX;
	double lfRefCenterY = m_pRefer->nCenterY;
	double lfRefCenterZ = m_pRefer->nCenterZ;
	
	//Modify!!
	_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;		
	
	int nIdx = 0;	
	int nSamplingNum = 10;

	//Search for the region in the rectangle
	for(int k = szLT.z ; k<szRB.z ; k+=nSamplingNum){
		for(int j=szLT.y ; j<szRB.y ; j+=nSamplingNum){
			for(int i=szLT.x ; i<szRB.x ; i+=nSamplingNum){
//				unsigned short nFltVol = pFltVol[nIdx++];
				
				unsigned short nFltVol = pWT[pFltVol[nIdx++]];
				
				//Calculate Transformed location
				double lfShiftOriginX = (i - lfFltCenterX) * lfScaleX;
				double lfShiftOriginY = (j - lfFltCenterY) * lfScaleY;
				double lfShiftOriginZ = (k - lfFltCenterZ) * lfScaleZ;

				double lfTransformedX = lfCosThetaY*lfCosThetaZ*lfShiftOriginX+(lfSinThetaX*lfSinThetaY*lfCosThetaZ-lfCosThetaX*lfSinThetaZ)*lfShiftOriginY+(lfCosThetaX*lfSinThetaY*lfCosThetaZ+lfSinThetaX*lfSinThetaZ)*lfShiftOriginZ+lfTransX+lfRefCenterX;
				double lfTransformedY = lfCosThetaY*lfSinThetaZ*lfShiftOriginX+(lfSinThetaX*lfSinThetaY*lfSinThetaZ+lfCosThetaX*lfCosThetaZ)*lfShiftOriginY+(lfCosThetaX*lfSinThetaY*lfSinThetaZ+lfSinThetaX*lfCosThetaZ)*lfShiftOriginZ+lfTransY+lfRefCenterY;
				double lfTransformedZ = -lfSinThetaY*lfShiftOriginX+lfSinThetaX*lfCosThetaY*lfShiftOriginY+lfCosThetaX*lfCosThetaY*lfShiftOriginZ+lfTransZ+lfRefCenterZ;


				if(IsValid((double)i,(double)j,(double)k, lfTransformedX, lfTransformedY, lfTransformedZ)){
					/////////////////////////
					//Classify interpolation
					switch(m_nInterpID)
					{
						case PV :
							PartialVolumeInterpolation(lfTransformedX, lfTransformedY, lfTransformedZ, nFltVol);							
							break;
						case TRI :							
							TrilinearVolumeInterpolation(lfTransformedX, lfTransformedY, lfTransformedZ, nFltVol);							
							break;
						case NN :							
							NearestInterpolation(lfTransformedX, lfTransformedY, lfTransformedZ, nFltVol);
							break;		
					}					
					/////////////////////////

					//Implement Gradient Information				
					BOOL bGradientInfo = FALSE;
					if(bGradientInfo){
						if(IsValid((double)i+1,(double)j+1,(double)k+1, lfTransformedX+1, lfTransformedY+1, lfTransformedZ+1))
							ComputeGradientInfo(i,j,k, (int)(lfTransformedX+0.5), (int)(lfTransformedY+0.5), (int)(lfTransformedZ+0.5));
					}
				}
			}
		}
	}
}

BOOL RxxMutualInfoAgent::IsValid(double lfOriginX, double lfOriginY, double lfOriginZ, double lfTranformedX, double lfTransformedY, double lfTransformedZ)
{
	if(lfOriginX<0 || lfOriginX>=m_pFloat->nVolX-1 ||
		lfOriginY<0 || lfOriginY>=m_pFloat->nVolY-1 ||
		lfOriginZ<0 || lfOriginZ>=m_pFloat->nVolZ-1)
		return FALSE;

	if(lfTranformedX<0 || lfTranformedX>=m_pRefer->nVolX-1 ||
		lfTransformedY<0 || lfTransformedY>=m_pRefer->nVolY-1 || 
		lfTransformedZ<0 || lfTransformedZ>=m_pRefer->nVolZ-1)
		return FALSE;
	else	return TRUE;
}

unsigned short RxxMutualInfoAgent::GetRefValue(int nXRefIdx, int nYRefIdx, int nZRefIdx)
{
	int nRefIdx = nZRefIdx*m_pRefer->nVolX*m_pRefer->nVolY+nYRefIdx*m_pRefer->nVolX+nXRefIdx;
//	unsigned short byRefValue = m_pRefer->pGaussianVol[nRefIdx];
//	unsigned short byRefValue = m_pRefer->pVol[nRefIdx];
	
	unsigned short nRefValue = m_pRefer->pVol[nRefIdx];
	return nRefValue;
}

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

	//Bottom
	// w5 w4
	// w6 w7

	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;
}

void RxxMutualInfoAgent::UpdateJointHistogram(unsigned short nRefValue, unsigned short nFloatValue, double lfWeight)
{	
//	m_plfJointHisto[nFloatValue][nRefValue] += lfWeight;
	int nIndex = nRefValue*FLT_DEPTH+nFloatValue;
	m_plfJointHisto[nIndex] +=lfWeight;
}

double RxxMutualInfoAgent::ComputeMISub()
{	
	double lfMapSize = m_lfRegionInPointCnt;
	double lfInfo = ComputeInfo(lfMapSize);
	return lfInfo;
}

/*
double RxxMutualInfoAgent::ComputeInfo(double lfMapSize)
{
	double lfInfo = 0.0;
	double T = 0.0;
	double pfProbF=0.0;
	double pfProbR=0.0;	
	double pfProbFR = 0.0;

	double Hf = 0.0;
	double Hr = 0.0;
	double Hfr = 0.0;
	

	for(int f=0 ; f<PETDEPTH ; f++){
		pfProbF = Prob_F(f,lfMapSize);
		if(pfProbF!=0)
			Hf -= (pfProbF*log(pfProbF));
	}
	
	for(int r=0 ; r<MRDEPTH ; r++){
		pfProbR = Prob_R(r,lfMapSize);
		if(pfProbR!=0)
			Hr -= (pfProbR*log(pfProbR));
	}
		
	for(f=0 ; f<CTDEPTH ; f++){
		for(int r=0 ; r<COLORDEPTH ; r++){
			pfProbFR = Prob_FR(f, r, lfMapSize);
			if(pfProbFR!=0)
				Hfr -= (pfProbFR*log(pfProbFR));
		}
	}
	

//	lfInfo = (Hf + Hr)/Hfr;
	lfInfo = Hf + Hr - Hfr;
	return lfInfo;
}
*/

double RxxMutualInfoAgent::ComputeInfo(double lfMapSize)
{
	double lfInfo = 0.0;
	double T = 0.0;	
	
	memset(m_plfProbF, 0x00, sizeof(double)*FLT_DEPTH);
	memset(m_plfProbR, 0x00, sizeof(double)*REF_DEPTH);

	for(int f=0 ; f<FLT_DEPTH ; f++){
		m_plfProbF[f] = Prob_F(f,lfMapSize);
//		if(pfProbF[f]!=0)
//			TRACE(_T("F=%lf\n"), pfProbF[f]);
	}
	
	for(int r=0 ; r<REF_DEPTH ; r++){
		m_plfProbR[r] = Prob_R(r,lfMapSize);
//		if(pfProbR[r]!=0)
//			TRACE(_T("R=%lf\n"), pfProbR[r]);
	}
		
	for(f=0 ; f<FLT_DEPTH ; f++){
		for(int r=0 ; r<REF_DEPTH ; r++){			
			if(m_plfProbF[f]>0 && m_plfProbR[r]>0){
				double lfProbFR = Prob_FR(f,r,lfMapSize);
				T = (double)log((lfProbFR/(m_plfProbF[f]*m_plfProbR[r])));
				T /= log(2);
				if(T>0)					
					lfInfo += lfProbFR*T;
			}			
		}
	}
	return lfInfo;
}

double RxxMutualInfoAgent::Prob_FR(int f, int r, double lfJointSum)
{
	int nIndex = r*FLT_DEPTH+f;
	return (m_plfJointHisto[nIndex]/lfJointSum);
//	return (m_plfJointHisto[f][r]/lfJointSum);
}

double RxxMutualInfoAgent::Prob_F(int f, double lfJointSum)
{	
	double lfSum = 0.0;
	for(int r=0 ; r<REF_DEPTH ; r++){
		int nIndex = r*FLT_DEPTH+f;
		lfSum += m_plfJointHisto[nIndex];
	}
	return (lfSum/lfJointSum);
}

double RxxMutualInfoAgent::Prob_R(int r, double lfJointSum)
{	
	double lfSum = 0.0;
	int nRIndex = r*FLT_DEPTH;
	for(int f=0 ; f<FLT_DEPTH ; f++){
		lfSum += m_plfJointHisto[nRIndex+f];
	}
	return (lfSum/lfJointSum);
}

void RxxMutualInfoAgent::KeepParameter(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_lfScaleX = lfScaleVectorX;
	m_lfScaleY = lfScaleVectorY;
	m_lfScaleZ = 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;
}

void RxxMutualInfoAgent::UpdateGradientInfo(double w, double Flt_Mag, double Ref_Mag)
{
	double info = w * (Flt_Mag > Ref_Mag) ? Ref_Mag : Flt_Mag;
	m_lfGradientInfo += info;
//	TRACE(_T("m_lfGradientInfo=%lf\n"), m_lfGradientInfo);
}

void RxxMutualInfoAgent::ComputeGradientInfo(int nX, int nY, int nZ, int nTransX, int nTransY, int nTransZ)
{	
	double FltP =0.0, FltQ = 0.0, FltR = 0.0;
	double RefP = 0.0, RefQ = 0.0, RefR = 0.0;
	double Float_Mag = 0.0, Ref_Mag=0.0;	
	
	int nXofFlt = m_pFloat->nVolX;	int nYofFlt = m_pFloat->nVolY;	int nZofFlt = m_pFloat->nVolZ;
	int nXofRef = m_pRefer->nVolX;	int nYofRef = m_pRefer->nVolY;	int nZofRef = m_pRefer->nVolZ;	

	unsigned short* pFloatGaussian = m_pFloat->pGaussianVol;
	unsigned short* pRefGaussian = m_pRefer->pGaussianVol;

	int nFloatNowIdx = nZ*nXofFlt*nYofFlt + nY*nXofFlt + nX;
	int nFloatNextX = nZ*nXofFlt*nYofFlt + nY*nXofFlt + nX + 1;
	int nFloatNextY = nZ*nXofFlt*nYofFlt + (nY+1)*nXofFlt + nX;
	int nFloatNextZ = (nZ+1)*nXofFlt*nYofFlt + nY*nXofFlt + nX;

	int nRefNowIdx = nTransZ*nXofRef*nYofRef + nTransY*nXofRef + nTransX;
	int nRefNextX = nTransZ*nXofRef*nYofRef + nTransY*nXofRef + nTransX + 1;
	int nRefNextY = nTransZ*nXofRef*nYofRef + (nTransY+1)*nXofRef + nTransX;
	int nRefNextZ = (nTransZ+1)*nXofRef*nYofRef + nTransY*nXofRef + nTransX;
	
	FltP = (pFloatGaussian[nFloatNextX]- pFloatGaussian[nFloatNowIdx] + 
		pFloatGaussian[nFloatNextY+1]- pFloatGaussian[nFloatNextY] +
		pFloatGaussian[nFloatNextZ+1]- pFloatGaussian[nFloatNextZ] +
		pFloatGaussian[nFloatNextZ+nXofFlt+1]- pFloatGaussian[nFloatNextZ+nXofFlt])/4.0;

	FltQ = (pFloatGaussian[nFloatNextY]- pFloatGaussian[nFloatNowIdx] + 
		pFloatGaussian[nFloatNextY+1]- pFloatGaussian[nFloatNextX] +
		pFloatGaussian[nFloatNextZ+nXofFlt]- pFloatGaussian[nFloatNextZ] +
		pFloatGaussian[nFloatNextZ+nXofFlt+1]- pFloatGaussian[nFloatNextZ+1])/4.0;

	FltR = (pFloatGaussian[nFloatNextZ]- pFloatGaussian[nFloatNowIdx] + 
		pFloatGaussian[nFloatNextZ+1]- pFloatGaussian[nFloatNextX] +
		pFloatGaussian[nFloatNextZ+nXofFlt]- pFloatGaussian[nFloatNextY] +