derivatemi3d.cpp

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

// DerivateMI3D.cpp: implementation of the RxDerivateMI3D class.
//
//////////////////////////////////////////////////////////////////////
//
//	Title: Calculate of the derivaitve of mutual information
//
//	1. Factor1: intensity difference (1x1 vector)
//	2. Factor2: weighting factor on image intensity (constant)
//	3. Factor3: weighting factor on square distance (constant)
//	4. Factor4: difference of weighting factors (constant)
//	5. Factor5: gradient difference (3x1 vector)
//			  : image gradient and derivative of transformed coordinate
//
//////////////////////////////////////////////////////////////////////
//
//	Author: Helen Hong, 3DMed co. LTD
//	138-dong 417-ho Seoul National Univ.
//	San 56-1 Shinlim-dong Kwanak-gu, Seoul, Korea
//	Email. [e-mail]
//
//	Date	: 2002. 9. 10.
//	Update	: 2003. 1. 8.
//	
//////////////////////////////////////////////////////////////////////

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

#include "stdafx.h"
#include "DerivateMI3D.h"

#include "VolumeInfo.h"
#include "Matrix3D.h"
#include "AutoSample3DInfo.h"

#include <cmath>

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

#define	PI							3.1415926535897931160
#define	PI_2						6.2831853071795862320
#define	pi_180						PI / 180.0
#define	pi_360						PI / 360.0

#define	DESTROY_CLASS_INSTANCE(ci)	{if(ci) {delete ci; ci = NULL;}}

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

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

RxDerivateMI3D::RxDerivateMI3D()
{
	m_pnVolume	= NULL;

	m_pMatrix3D	= new RxMatrix3D;
	m_pResol	= new RxVolumeInfo;
}

RxDerivateMI3D::~RxDerivateMI3D()
{
	DESTROY_CLASS_INSTANCE(m_pMatrix3D);
	DESTROY_CLASS_INSTANCE(m_pResol);
}

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

//*********************************************************************
//	SetVolumeData(): set reference volume data
//*********************************************************************
BOOL RxDerivateMI3D::SetVolumeData(unsigned short *Buf)
{
	m_pnVolume		= Buf;

	return TRUE;
}

//*********************************************************************
//	SetReferVolInfo(): set reference volume information
//*********************************************************************
BOOL RxDerivateMI3D::SetReferVolInfo(RxVolumeInfo ReferInfo)
{
	m_pResol->Ru	= ReferInfo.Ru;
	m_pResol->Rv	= ReferInfo.Rv;
	m_pResol->Rw	= ReferInfo.Rw;

	return TRUE;
}

//*********************************************************************
//	SetSampleNumber(): set sample number
//*********************************************************************
BOOL RxDerivateMI3D::SetSampleNumber(int number)
{
	m_SampleNum		= number;

	return TRUE;
}

//*********************************************************************
//	SetSampleA(): set sampleA
//*********************************************************************s
BOOL RxDerivateMI3D::SetSampleA(RxAutoSample3DInfo *sample)
{
	sampleA			= sample;

	return TRUE;
}

//*********************************************************************
//	SetTransform(): set geometrical transformation
//*********************************************************************
BOOL RxDerivateMI3D::SetTransform(RxTransform3DInfo *Trans)
{
	m_Trans			= *Trans;

	return TRUE;
}

//*********************************************************************
//	SetMIPParameter(): set parameters for mutual information
//*********************************************************************
BOOL RxDerivateMI3D::SetMIParameters(double varFloat, double covarRefer, double covarFloat)
{
	m_varFloat		= varFloat;		// variance of float volume
	m_covarRefer	= covarRefer;	// covariance of reference volume
	m_covarFloat	= covarFloat;	// covariance of float volume

	return TRUE;
}

//*********************************************************************
//	Calculate(): calculate the derivative of mutual information
//*********************************************************************
BOOL RxDerivateMI3D::Calculate(RxAutoSample3DInfo *sampleB,double (*totalfactor)[6])
{
	double	factor1, factor2, factor3, factor4, factor5[6],
			gauss_sum, SumofFactor[6];

	SumofFactor[0] = SumofFactor[1] = SumofFactor[2] = SumofFactor[3] = SumofFactor[4] = SumofFactor[5] = 0.0;
			
	for (int i = 0; i < m_SampleNum; i++) {
		for (int j = 0; j < m_SampleNum; j++) {
			if(i != j) {
				// Initialize factor values
				factor1 = factor2 = factor3 = factor4 = 0.0;
				factor5[0] = factor5[1] = factor5[2] = factor5[3] = factor5[4] = factor5[5] = 0.0;
				
				// FACTOR1: calculate intensity difference between vi and vj
				factor1 = (double)(sampleB[i].val-sampleB[j].val);

				if(factor1 != 0.0) {
					// FACTOR2: calculate weighting factor for vi & vj relation
					factor2 = Gauss_point(sampleB[i].val,sampleB[j].val,m_varFloat);
					gauss_sum = 0.0;
					if(factor2 != 0.0) {
						for (int k = 0; k < m_SampleNum; k++) {
							if(k != i) {
								gauss_sum += Gauss_point(sampleB[i].val,sampleB[k].val,m_varFloat);
							}
						}
						if(gauss_sum != 0.0) {
							factor2 /= gauss_sum;
							factor2 /= m_varFloat;
						}
						else {
							factor2 = 0.0;
						}
					}
					else {
						factor2 = 0.0;
					}

					// FACTOR3: calculate weighting factor wi & wj relation
					factor3 = Gauss_pair(sampleA[i].val,sampleB[i].val,sampleA[j].val,sampleB[j].val,m_covarRefer,m_covarFloat);
					gauss_sum = 0.0;
					if(factor3 != 0.0) {
						for (int k = 0; k < m_SampleNum; k++) {
							if(k != i) {
								gauss_sum += Gauss_pair(sampleA[i].val,sampleB[i].val,sampleA[k].val,sampleB[k].val,m_covarRefer,m_covarFloat);
							}
						}

						if(gauss_sum != 0.0) {
							factor3 /= gauss_sum;
							factor3 /= m_covarFloat;
						}
						else {
							factor3 = 0.0;
						}
					}
					else {
						factor3 = 0.0;
					}

					// FACTOR4: calculate difference of weighting function
					factor4 = factor2 - factor3;

					// FACTOR5: involve image gradient and the derivative of transformed coordinate
					double	result[6];
					result[0]=result[1]=result[2]=result[3]=result[4]=result[5]=0.0;
					if(factor4 != 0.0) {
						DerivateTrans3D(sampleB[i],sampleB[j],&factor5);