mutualinfosm.cpp

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

// MutualInfoSM.cpp: implementation of the LHMutualInfoSM 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 "fusion.h"
#include "MutualInfoSM.h"

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

#include "FusionGlobal.h"
#include "FastInlineFuncs.h"
#include "RegData.h"
#include<math.h>


#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
LHMutualInfoSM::LHMutualInfoSM()
{	
	m_plfJointHisto		= NULL;	
	m_iJointHistoX = 256;
	m_iJointHistoY = 256;
	m_iSubSampLevel = 1;
}

LHMutualInfoSM::LHMutualInfoSM(_MODAL_INFO* pFltInfo, _MODAL_INFO* pRefInfo)
{
	m_pbyFltVol		= ConversionData(pFltInfo->pVol, pFltInfo->nVolX, pFltInfo->nVolY, pFltInfo->nVolZ, m_iJointHistoX);
	m_pbyRefVol		= ConversionData(pRefInfo->pVol, pRefInfo->nVolX, pRefInfo->nVolY, pRefInfo->nVolZ, m_iJointHistoY);
	m_pFltInfo		= pFltInfo;
	m_pRefInfo		= pRefInfo;

	m_lfRatioFltZ	= m_pFltInfo->lfVoxelZ / m_pFltInfo->lfVoxelX;
	m_lfRatioRefZ	= m_pRefInfo->lfVoxelZ / m_pRefInfo->lfVoxelX;	

	int	nSeries		= RxGetCurSeries();
	m_nInterpID		= RxGetFrameMain()->m_pFMWndVR[nSeries]->m_MIInterpolation;	

	m_iDirection	= XTRANS;
//	m_mxTR.LoadIdentity();

	m_lfMaxMI = 0.0f;
	m_iSubSampLevel = 1;

	m_plfJointHisto = new float[m_iJointHistoX*m_iJointHistoY];	
	InitialTransformParameters();	
}

double LHMutualInfoSM::VoxelBasedMI( double lfT)
{
	m_lfJointHistoSize = 0.0f;
	memset(m_plfJointHisto, 0x00, sizeof(float)*m_iJointHistoX*m_iJointHistoY);
	
	switch(m_iDirection){
		case XTRANS:{	m_OptTransformPar.lfTx = lfT;	break; }
		case YTRANS:{	m_OptTransformPar.lfTy = lfT;	break; }
		case ZTRANS:{	m_OptTransformPar.lfTz = lfT;	break; }
		case XROT:	{	m_OptTransformPar.lfRotX = lfT;	break; }
		case YROT:	{	m_OptTransformPar.lfRotY = lfT;	break; }
		case ZROT:	{	m_OptTransformPar.lfRotZ = lfT;	break; }
		defaut: break;
	}
	RxMatrix4D matTR = MakeTransformMatrix(m_OptTransformPar);
	
	int nSampIntval = 1;
	int iFslice = m_pFltInfo->nVolX*m_pFltInfo->nVolY;	
//	int nOptFltTh = 120; //SPECT phantom1 tolerance
//	int nOptFltTh = 200; //SPECT phantom2 tolerance
	
	for(int z = 0 ; z<m_pFltInfo->nVolZ ; z+=nSampIntval){
		if(z%m_iSubSampLevel!=0)	continue;
		for(int y=0 ; y<m_pFltInfo->nVolY ; y+=nSampIntval){
			if(y%m_iSubSampLevel!=0)	continue;
			for(int x=0 ; x<m_pFltInfo->nVolX ; x+=nSampIntval){
				if(x%m_iSubSampLevel!=0)	continue;

//				if(m_pFltInfo->pVol[z*iFslice+y*m_pFltInfo->nVolX+x]<nOptFltTh)	continue;

				RxVect4D v(x, y, z, 1.0f);
				v = matTR*v;

				if(v.m[0]>=0 && v.m[0]<=m_pRefInfo->nVolX-1 &&
					v.m[1]>=0 && v.m[1]<=m_pRefInfo->nVolY-1 &&
						v.m[2]>=0 && v.m[2]<=m_pRefInfo->nVolZ-1){
					switch(m_nInterpID)
					{
						case PV :
							PartialVolumeInterpolation(v.m[0], v.m[1], v.m[2], m_pbyFltVol[iFslice*z+y*m_pFltInfo->nVolX+x]);
							break;
						case TRI :							
							TrilinearVolumeInterpolation(v.m[0], v.m[1], v.m[2], m_pbyFltVol[iFslice*z+y*m_pFltInfo->nVolX+x]);							
							break;
						case NN :							
							NearestInterpolation(v.m[0], v.m[1], v.m[2], m_pbyFltVol[iFslice*z+y*m_pFltInfo->nVolX+x]);
							break;
					}
				}
			}
		}
	}
	double MI = ComputeMI();	
	return 1.0f/MI;
}

CString LHMutualInfoSM::OutputToString(_TransformParameters _OptTransformPar, double lfHf, double lfHr, double lfHfr, double lfMI, double lfMin)
{
	BOOL bMax=FALSE;
	if(m_lfMaxMI<lfMI){
		m_lfMaxMI = lfMI;
		bMax = TRUE;
	}
	CString strResults= _T("");	
	strResults.Format(_T("(Tx,Ty,Tz)=( %lf, %lf, %lf)\t(Rx,Ry,Rz)=( %lf, %lf, %lf )\tHf=%lf\tHr=%lf\tHfr=%lf\tMI=%lf\tMin=%lf\t%d\r\n"),
		_OptTransformPar.lfTx, _OptTransformPar.lfTy, _OptTransformPar.lfTz,
		_OptTransformPar.lfRotX, _OptTransformPar.lfRotY, _OptTransformPar.lfRotZ,
		lfHf, lfHr, lfHfr, lfMI, lfMin, bMax);
		
	return strResults;
}

LHMutualInfoSM::~LHMutualInfoSM()
{
	SAFE_DELETE_ARRAY(m_plfJointHisto);
	SAFE_DELETE_ARRAY(m_pbyFltVol);
	SAFE_DELETE_ARRAY(m_pbyRefVol);
}

void LHMutualInfoSM::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 LHMutualInfoSM::UpdateJointHistogram(BYTE nFltValue, BYTE nRefValue, float& fWeight)
{	
//	m_plfJointHisto[nFloatValue][nRefValue] += fWeight;
	int nIndex = nRefValue*m_iJointHistoX+nFltValue;
	m_plfJointHisto[nIndex] +=fWeight;
}

double LHMutualInfoSM::ComputeMI()
{
	double lfInfo = 0.0;
	double Pf = 0.0;
	double Pr = 0.0;	
	double Pfr = 0.0;

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

	for(int iIdx=0 ; iIdx<m_iJointHistoX*m_iJointHistoY ; iIdx++)
		m_lfJointHistoSize+=m_plfJointHisto[iIdx];

	for(int f=0 ; f<m_iJointHistoX ; f++){
		Pf = Prob_F(f,m_lfJointHistoSize);
		if(Pf!=0)	Hf -= (Pf*(log(Pf)/log(2)));
	}
	
	for(int r=0 ; r<m_iJointHistoY ; r++){
		Pr = Prob_R(r,m_lfJointHistoSize);
		if(Pr!=0)	Hr -= (Pr*(log(Pr)/log(2)));
	}
		
	for(r=0 ; r<m_iJointHistoY ; r++){
		for(int f=0 ; f<m_iJointHistoX ; f++){
			Pfr = Prob_FR(f, r, m_lfJointHistoSize);
			if(Pfr!=0)
				Hfr -= (Pfr*(log(Pfr)/log(2)));
		}
	}	
	//NMI
	lfInfo = (Hf + Hr)/Hfr;
	//MI
//	lfInfo = Hf + Hr - Hfr;

//	SaveResultsToFile(_T("MI"), OutputToString(m_OptTransformPar, Hf, Hr, Hfr, lfInfo, 1.0f/lfInfo));
	return lfInfo;
}

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

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

void LHMutualInfoSM::InitialTransformParameters()
{
	m_OptTransformPar.lfSx			= m_pFltInfo->lfVoxelX/m_pRefInfo->lfVoxelX;
	m_OptTransformPar.lfSy			= m_pFltInfo->lfVoxelY/m_pRefInfo->lfVoxelY;
	m_OptTransformPar.lfSz			= m_pFltInfo->lfVoxelZ/m_pRefInfo->lfVoxelZ;
	m_OptTransformPar.lfRotX		= 0.0;
	m_OptTransformPar.lfRotY		= 0.0;
	m_OptTransformPar.lfRotZ		= 0.0;
	m_OptTransformPar.lfTx			= 0.0;
	m_OptTransformPar.lfTy			= 0.0;
	m_OptTransformPar.lfTz			= 0.0;
}

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

void LHMutualInfoSM::NearestInterpolation(double lfTransformedX, double lfTransformedY, double lfTransformedZ, BYTE nFltValue)
{
	int nXRef = (int)(lfTransformedX+0.5);
	int nYRef = (int)(lfTransformedY+0.5);
	int nZRef = (int)(lfTransformedZ+0.5);
//	int nXRef = FloatToInt(lfTransformedX+0.5f);
//	int nYRef = FloatToInt(lfTransformedY+0.5f);
//	int nZRef = FloatToInt(lfTransformedZ+0.5f);
	
	BYTE nRefValue = m_pbyRefVol[nZRef*m_pRefInfo->nVolX*m_pRefInfo->nVolY+nYRef*m_pRefInfo->nVolX+nXRef];
//	m_lfJointHistoSize++;	
	float fVote = 1.0f;
	UpdateJointHistogram(nRefValue, nFltValue, fVote);
}

void LHMutualInfoSM::TrilinearVolumeInterpolation(double lfTransformedX, double lfTransformedY, double lfTransformedZ, BYTE nFltValue)
{
	int ix, iy, iz;

//	ix = FloatToInt(fTransformedX);
//	iy = FloatToInt(fTransformedY);
//	iz = FloatToInt(fTransformedZ);
	ix = (int)(lfTransformedX);
	iy = (int)(lfTransformedY);
	iz = (int)(lfTransformedZ);
	double xFract = lfTransformedX - ix;
	double yFract = lfTransformedY - iy;
	double zFract = lfTransformedZ - 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;
	double lfSum=0.0;
	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];		
		lfSum=lfSum+nRefValue*m_fWeight[nIndex];
	}
	BYTE nInterpolatedRefValue = (BYTE)(lfSum+0.5);
//	m_lfJointHistoSize++;
	float fVote = 1.0f;
	UpdateJointHistogram(nFltValue, nInterpolatedRefValue, fVote);	
}
/*
void LHMutualInfoSM::TrilinearVolumeInterpolation(double lfTransformedX, double lfTransformedY, double lfTransformedZ, BYTE nFltValue)
{
	int nTransformedXInt = (int)(lfTransformedX);
	int nTransformedYInt = (int)(lfTransformedY);