mianalysis.cpp

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

// MIAnalysis.cpp: implementation of the RxMIAnalysis class.
//
//////////////////////////////////////////////////////////////////////
//
//	Title: Calculate mutual information for analysis
//
//////////////////////////////////////////////////////////////////////
//
//	Author: Helen Hong, 3DMed co. LTD
//	138-dong 417-ho, Seoul National Univ.
//	Shinlim-dong, Kwanak-gu, Seoul, Korea
//	(Email. [e-mail])
//
//	Date	: 2002. 10. 10.
//	Update	: 2002. 11. 5.
//
//////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////
//	include
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "fusion.h"
#include "MIAnalysis.h"

#include "Transform3DInfo.h"
#include "Matrix3D.h"

#include <math.h>

//////////////////////////////////////////////////////////////////////
//	define
//////////////////////////////////////////////////////////////////////

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

#define	COLORDEPTH		255
#define	CTDEPTH			4096
#define	MRDEPTH			4096
#define	PETDEPTH		4096

//////////////////////////////////////////////////////////////////////
//	declaration
//////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

RxMIAnalysis::RxMIAnalysis()
{
	m_pnReferVolume		= NULL;
	m_pnFloatVolume		= NULL;
	m_JointHistogram	= NULL;
}

RxMIAnalysis::~RxMIAnalysis()
{
}

//////////////////////////////////////////////////////////////////////

//*********************************************************************
//	SetReferVolume(): Set reference volume and resolution
//*********************************************************************
BOOL RxMIAnalysis::SetReferVolume(unsigned short *Buf,short iVolX,short iVolY,short iVolZ)
{
	m_pnReferVolume	= Buf;

	m_ReferVolX		= iVolX;
	m_ReferVolY		= iVolY;
	m_ReferVolZ		= iVolZ;

	return TRUE;
}

//*********************************************************************
//	SetFloatVolume(): Set float volume and resolution
//*********************************************************************
BOOL RxMIAnalysis::SetFloatVolume(unsigned short *Buf,short iVolX,short iVolY,short iVolZ)
{
	m_pnFloatVolume	= Buf;

	m_FloatVolX		= iVolX;
	m_FloatVolY		= iVolY;
	m_FloatVolZ		= iVolZ;

	return TRUE;
}

//*********************************************************************
//	SetTransformation(): Set geometrical transformation
//*********************************************************************
BOOL RxMIAnalysis::SetTransformation(RxTransform3DInfo *Trans)
{
	m_pTrans = Trans;

	return TRUE;
}

//*********************************************************************
//	ComputeMI(): Calculate mutual information
//*********************************************************************
void RxMIAnalysis::ComputeMI()
{
	/*
	InitialJointHistogram();
	ComputeTransformedVector();
	
	double	MI = ComputeMISub();

	return MI;
	*/
}

/*
//*********************************************************************
//	InitialJointHistogram(): Initialize joint histogram
//*********************************************************************
BOOL RxMIAnalysis::InitialJointHistogram()
{
	if(m_JointHistogram == NULL) {
		m_JointHistogram = new double*[PETDEPTH];
		for (int j = 0; j < PETDEPTH; j++)
			m_JointHistogram[j] = new double[MRDEPTH];
	}
	for(int j = 0; j < PETDEPTH; j++) {
		memset(m_JointHistogram[j], 0, sizeof(double)*MRDEPTH);
	}

	return TRUE;
}

//*********************************************************************
//	ComputeTransformedVector(): Compute transformed vector
//*********************************************************************
BOOL RxMIAnalysis::ComputeTransformedVector()
{
	// Transform object coordinate to eye coordinate
	m_pMatrix->make_viewmat3D(m_obj2eye, m_pTrans);

	// Transform object coordinate to eye coordinate
	m_pMatrix->make_viewmat3D(m_obj2eye, m_pTrans);

	// Generate transformed float volume data
	double	vector[4], transform[4];
	for (int z = 0; z < m_ReferVolZ; z++) {
		for (int y = 0; y < m_ReferVolY; y++) {
			for (int x = 0; x < m_ReferVolX; x++) {
				vector[0] = (double)x;
				vector[1] = (double)y;
				vector[2] = (double)z;
				vector[3] = 1.0;

				m_pMatrix->matrix3D_vector4(m_obj2eye,vector,transform);

				int xCoord = (int)(transform[0]+0.5);
				int yCoord = (int)(transform[1]+0.5);
				int	zCoord = (int)(transform[2]+0.5);

				int	index = ((z*m_ReferVolY)+y)*m_ReferVolX+x;
				m_pnTransVolume[index] = GetIntensity3D(m_pnFloatVolume,m_FloatVolX,m_FloatVolY,m_FloatVolZ,xCoord,yCoord,zCoord);
				TRACE(_T("\n m_pnTransVolume[%d]: %d"), index, m_pnTransVolume[index]);
			} // end of x for
		} // end of y for
	} // end of z for

	return TRUE;


	int	nIdx = 0;
	for(int k = 0; k < m_FloatVolZ; k++) {
		for(int j = 0; j < m_FloatVolY; j++) {
			for(int i = 0; i < m_FloatVolX; i++) {
				unsigned short nFltVol = m_pnFloatVolume[nIdx++];

				// Calculate transformed location
				double	lfShiftOriginX = (i - m_FloatCenterX) * Trans->xscale;
				double	lfShiftOriginY = (j - m_FloatCenterY) * Trans->yscale;
				double	lfShiftOriginZ = (k - m_FloatCenterZ) * Trans->zscale;
	
				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
					NearestInterplation(lfTransformedX, lfTransformedY, lfTransformedZ, nFltVol);

//					TrilinearVolumeInterplation(lfTransformedX, lfTransformedY, lfTransformedZ, nFltVol);

//					PartialVolumeInterplation(lfTransformedX, lfTransformedY, lfTransformedZ, nFltVol);					
					////////////
					//Implement Gradient Information				
					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));
				}
			}
		}
	}




//*********************************************************************
//	ReferMarginalEntropy(): Calculate reference marginal entropy
//*********************************************************************
double RxMIAnalysis::ReferMarginalEntropy()
{
	double	result = 0.0;
	double	rp = 0.0;
	for(int row = 0; row < CTDEPTH; row++) {
		rp = ReferMarginalProbability(row);
		if(rp != 0.0) {
			result += (rp * log(rp));
		}
	}

	return -result;
}

//*********************************************************************
//	ReferMarginalProbaility(): Calculate reference marginal probability
//*********************************************************************
double RxMIAnalysis::ReferMarginalProbability(int index)
{
	double	result = 0.0;
	for(int col = 0; col < CTDEPTH; col++) {
		if(m_BiPDF[index][col] != 0.0) {
			result += m_BiPDF[index][col];
		}
	}

	return result;
}

//*********************************************************************
//	FloatMarginalEntropy(): Calculate float marginal entropy
//*********************************************************************
double RxMIAnalysis::FloatMarginalEntropy()
{
	double	result = 0.0;
	double	fp = 0.0;
	for (int col = 0; col < CTDEPTH; col++) {
		fp = FloatMarginalProbability(col);
		if(fp != 0.0) {
			result += (fp * log(fp));
		}
	}

	return -result;
}

//*********************************************************************
//	FloatMarginalProbability(): Calculate float marginal probability
//*********************************************************************
double RxMIAnalysis::FloatMarginalProbability(int index)
{
	double	result = 0.0;
	for(int row = 0; row < CTDEPTH; row++) {
		if(m_BiPDF[row][index] != 0.0) {
			result += m_BiPDF[row][index];
		}
	}

	return result;
}

//*********************************************************************
//	MutualInformation(): Calculate mutual information
//*********************************************************************
BOOL RxMIAnalysis::MutualInformation()
{
	GetTransVolumeData();

	m_JointEntropy			= JointEntropy();
	m_ReferMarginalEntropy	= ReferMarginalEntropy();
	m_FloatMarginalEntropy	= FloatMarginalEntropy();
	m_MutualInfo			= m_ReferMarginalEntropy + m_FloatMarginalEntropy - m_JointEntropy;
	
	TRACE(_T("\n %f+%f-%f=%f"), m_ReferMarginalEntropy,m_FloatMarginalEntropy,m_JointEntropy,m_MutualInfo);

	ReleaseVolumeData(m_pnTransVolume,m_ReferVolX,m_ReferVolY,m_ReferVolZ);

	return TRUE;
}

//*********************************************************************
//	GetTransVolumeData(): Generate transformed test volume data
//*********************************************************************
BOOL RxMIAnalysis::GetTransVolumeData()
{
	AllocateVolumeData(&m_pnTransVolume,m_ReferVolX,m_ReferVolY,m_ReferVolZ);

	// Transform object coordinate to eye coordinate
	m_pMatrix->make_viewmat3D(m_obj2eye, m_pTrans);

	// Generate transformed float volume data
	double	vector[4], transform[4];
	for (int z = 0; z < m_ReferVolZ; z++) {
		for (int y = 0; y < m_ReferVolY; y++) {
			for (int x = 0; x < m_ReferVolX; x++) {
				vector[0] = (double)x;
				vector[1] = (double)y;
				vector[2] = (double)z;
				vector[3] = 1.0;

				m_pMatrix->matrix3D_vector4(m_obj2eye,vector,transform);

				int xCoord = (int)(transform[0]+0.5);
				int yCoord = (int)(transform[1]+0.5);
				int	zCoord = (int)(transform[2]+0.5);

				int	index = ((z*m_ReferVolY)+y)*m_ReferVolX+x;
				m_pnTransVolume[index] = GetIntensity3D(m_pnFloatVolume,m_FloatVolX,m_FloatVolY,m_FloatVolZ,xCoord,yCoord,zCoord);
				TRACE(_T("\n m_pnTransVolume[%d]: %d"), index, m_pnTransVolume[index]);
			} // end of x for
		} // end of y for
	} // end of z for

	return TRUE;
}

//*********************************************************************
//	GetIntensity3D(): Get intensity value from float volume
//*********************************************************************
BOOL RxMIAnalysis::GetIntensity3D(unsigned short *Buf,short iVolX,short iVolY,short iVolZ,short x,short y,short z)
{
	if(x < 0 || y < 0 || z < 0 || x > iVolX-1 || y > iVolY-1 || z > iVolZ-1)
		return 0;
	return Buf[((z*iVolY)+y)*iVolX+x];
}
*/