rigidreg3dagent.cpp

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

// RigidReg3DAgent.cpp: implementation of the RxRigidReg3DAgent class.
//
//////////////////////////////////////////////////////////////////////
//
//	Title: 3D Rigid-body Registration based on Mutual Information
//
//		1. Rigid-body transformation
//		   - 6-DOF (3-translation, 3-rotation)
//		2. Automatic Sampling
//		   - Random sampling
//		   - Uniform sampling
//		   - 3D Levoy-Sobel operator
//		3. Interpolation
//		   - Nearest neighbor
//		   - Trilinear interpolation
//		   - Original Trilinear Interpolation
//		4. Cost function
//		   - Derivative of the mutual information
//		5. Stochastic optimization
//		   - steepest descent method
//		6. Automatic convergence
//
//////////////////////////////////////////////////////////////////////
//
//	Author: Helen Hong, 3DMed co. LTD
//	138-dong 417-ho, Seoul National Univ.
//	Shinlim-dong, Kwanak-gu, Seoul, Korea
//	(Email. [e-mail])
//
//	Date	: 2002. 9. 5.
//	Update	: 2003. 1. 14.
//
//////////////////////////////////////////////////////////////////////

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

#include "stdafx.h"
#include "RigidReg3DAgent.h"

#include "RigidReg3DInfo.h"
#include "Transform3DInfo.h"
#include "Matrix3D.h"
#include "AutoSample3D.h"
#include "AutoSample3DInfo.h"
#include "DerivateMI3D.h"
#include "RxxMutualInfoAgent.h"

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

#define	fabs(a)						((a < 0) ? -a : a)
#define	DESTROY_CLASS_INSTANCE(ci)	{if(ci) {delete ci; ci = NULL;}}

#define	NEAREST_NEIGHBOR

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

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

RxRigidReg3DAgent::RxRigidReg3DAgent()
{
	m_pnReferVolume		= NULL;
	m_pnFloatVolume		= NULL;

	m_pAutoSample3D		= new RxAutoSample3D;
	m_pDerivateMI3D		= new RxDerivateMI3D;
	m_pMIAnalysis		= new RxxMutualInfoAgent;
}

RxRigidReg3DAgent::~RxRigidReg3DAgent()
{
	DESTROY_CLASS_INSTANCE(m_pAutoSample3D);
	DESTROY_CLASS_INSTANCE(m_pDerivateMI3D);
	DESTROY_CLASS_INSTANCE(m_pMIAnalysis);
}

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

//********************************************************************
//	InitRegistration(): initialize registration
//********************************************************************
BOOL RxRigidReg3DAgent::InitRegistration(RxRigidReg3DInfo *RegistInfos,RxTransform3DInfo *TransformInfos)
{
	// Set reference and float volume data
	m_pnReferVolume		= RegistInfos->m_pnReferVolume;
	m_pnFloatVolume		= RegistInfos->m_pnFloatVolume;

	// Set reference and float volume resolution and voxel size
	m_pRigidReg3DInfo	= RegistInfos;

	// Set transformation 
	m_pTransform		= TransformInfos;

	return TRUE;
}

//********************************************************************
//	SetInitialParameters(): set initial parameters for MI analysis
//********************************************************************
BOOL RxRigidReg3DAgent::SetInitialParameters(_MODAL_INFO *Refer, _MODAL_INFO *Float)
{
	// Set reference volume information
	Refer->lfVoxelX	= m_pRigidReg3DInfo->m_ReferVolInfo.Su;
	Refer->lfVoxelY	= m_pRigidReg3DInfo->m_ReferVolInfo.Sv;
	Refer->lfVoxelZ	= m_pRigidReg3DInfo->m_ReferVolInfo.Sw;

	Refer->nVolX		= m_pRigidReg3DInfo->m_ReferVolInfo.Ru;
	Refer->nVolY		= m_pRigidReg3DInfo->m_ReferVolInfo.Rv;
	Refer->nVolZ		= m_pRigidReg3DInfo->m_ReferVolInfo.Rw;

	Refer->nCenterX	= Refer->nVolX / 2.0;
	Refer->nCenterY	= Refer->nVolY / 2.0;
	Refer->nCenterZ	= Refer->nVolZ / 2.0;

	Refer->pVol		= m_pnReferVolume;
	Refer->pWindowingTable = RxGetFrameTuning()->GetWindowingTable(0);

	// Set float volume information
	Float->lfVoxelX	= m_pRigidReg3DInfo->m_FloatVolInfo.Su;
	Float->lfVoxelY	= m_pRigidReg3DInfo->m_FloatVolInfo.Sv;
	Float->lfVoxelZ	= m_pRigidReg3DInfo->m_FloatVolInfo.Sw;
	
	Float->nVolX		= m_pRigidReg3DInfo->m_FloatVolInfo.Ru;
	Float->nVolY		= m_pRigidReg3DInfo->m_FloatVolInfo.Rv;
	Float->nVolZ		= m_pRigidReg3DInfo->m_FloatVolInfo.Rw;

	Float->nCenterX	= Float->nVolX / 2.0;
	Float->nCenterY	= Float->nVolY / 2.0;
	Float->nCenterZ	= Float->nVolZ / 2.0;

	Float->pVol		= m_pnFloatVolume;
	Float->pWindowingTable = RxGetFrameTuning()->GetWindowingTable(1);

	return TRUE;
}

//********************************************************************
//	Registration(): 3D rigid-body registration
//********************************************************************
BOOL RxRigidReg3DAgent::Registration(RxTransform3DInfo **FloatTrans)
{
	DWORD CurrentTime = GetTickCount();

	// Calculate real volume size on voxel size
	SetViewingMatrix();
	if(!SameTransformFactor()) {
		m_pMatrix3D->make_viewmat3D(m_obj2eye,m_pTransform);
		m_pMatrix3D->invert_mat3D(m_obj2eye,m_inverse_obj2eye);

		// Generate sample points in reference volume
		m_pAutoSample3D->SetVolumeData(m_pnReferVolume);
		m_pAutoSample3D->SetVolumeResolution(m_pRigidReg3DInfo->m_ReferVolInfo);
		m_pAutoSample3D->SetSampleInfo(m_pRigidReg3DInfo->m_SampleNum);
		switch (m_pRigidReg3DInfo->m_SampleMode) {
		case 0:										// uniform sampling
			m_pAutoSample3D->UniformSampling(&m_sampleA);
			break;
		case 1:										// random sampling
			m_pAutoSample3D->RandomSampling(&m_sampleA);
			break;
		case 2:										// sobel operator
			m_pAutoSample3D->SobelSampling(&m_sampleA);
			break;
		case 3:										// sobel-levoy sampling
			m_pAutoSample3D->SobelLevoySampling(&m_sampleA);
			break;
		}

		// Automatic convergence
		m_Repeat = 0;
		AutomaticConverge();

		delete [] m_sampleA;
	}
	
	DWORD ElapsedTime = GetTickCount() - CurrentTime;
	CString	str;
	str.Format(_T("\n Total Processing Time : %d ms"), ElapsedTime);
	AfxMessageBox(str);

	// calculate mutual information for analysis
	if(IsValid_MIAnalysis(m_pRigidReg3DInfo->m_MIAnalysis)) {
		_MODAL_INFO* pRefer = new _MODAL_INFO;
		_MODAL_INFO* pFloat = new _MODAL_INFO;

		SetInitialParameters(pRefer, pFloat);
		int nInterpID = RxGetFrameMain()->m_pFMWndVR[2]->m_MIInterpolation;

		RxxMutualInfoAgent *pPM = new RxxMutualInfoAgent(pRefer, pFloat, nInterpID, m_pTransform->xscale, m_pTransform->yscale, m_pTransform->zscale);
		double MI = pPM->ComputeMI(m_pTransform->xscale,m_pTransform->yscale,m_pTransform->zscale,
												  m_pTransform->xmove,m_pTransform->ymove,m_pTransform->zmove,
												  m_pTransform->xrot,m_pTransform->yrot,m_pTransform->zrot);
		TRACE(_T("\n MI analysis: %f"), MI);

		delete pPM;
		delete pRefer;
		delete pFloat;
	}

	*FloatTrans = m_pTransform;

	return TRUE;
}

//********************************************************************
//	SetViewingMatrix(): set viewing matrix 
//********************************************************************
BOOL RxRigidReg3DAgent::SetViewingMatrix()
{
	// calculate real volume size depending on voxel size
	m_pTransform->lengthX1 = (double)m_pRigidReg3DInfo->m_ReferVolInfo.Ru * m_pRigidReg3DInfo->m_ReferVolInfo.Su;
	m_pTransform->lengthX2 = (double)m_pRigidReg3DInfo->m_FloatVolInfo.Ru * m_pRigidReg3DInfo->m_FloatVolInfo.Su;
	m_pTransform->lengthY1 = (double)m_pRigidReg3DInfo->m_ReferVolInfo.Rv * m_pRigidReg3DInfo->m_ReferVolInfo.Sv;
	m_pTransform->lengthY2 = (double)m_pRigidReg3DInfo->m_FloatVolInfo.Rv * m_pRigidReg3DInfo->m_FloatVolInfo.Sv;
	m_pTransform->lengthZ1 = (double)m_pRigidReg3DInfo->m_ReferVolInfo.Rw * m_pRigidReg3DInfo->m_ReferVolInfo.Sw;
	m_pTransform->lengthZ2 = (double)m_pRigidReg3DInfo->m_FloatVolInfo.Rw * m_pRigidReg3DInfo->m_FloatVolInfo.Sw;

	// movement in float volume (shift factor)
	m_pTransform->xmove		= ((m_pTransform->lengthX2-m_pTransform->lengthX1)/2.0)/m_pTransform->lengthX2;
	m_pTransform->ymove		= ((m_pTransform->lengthY2-m_pTransform->lengthY1)/2.0)/m_pTransform->lengthY2;
	m_pTransform->zmove		= ((m_pTransform->lengthZ2-m_pTransform->lengthZ1)/2.0)/m_pTransform->lengthZ2;

	// scaling in float volume (scaling factor)
	m_pTransform->xscale	= m_pRigidReg3DInfo->m_FloatVolInfo.Su / m_pRigidReg3DInfo->m_ReferVolInfo.Su;
	m_pTransform->yscale	= m_pRigidReg3DInfo->m_FloatVolInfo.Sv / m_pRigidReg3DInfo->m_ReferVolInfo.Sv;
	if(m_pRigidReg3DInfo->m_ReferVolInfo.Ru == 512)
		m_pRigidReg3DInfo->m_ReferVolInfo.Sw = 0.6;
	if(m_pRigidReg3DInfo->m_FloatVolInfo.Ru == 512)
		m_pRigidReg3DInfo->m_FloatVolInfo.Sw = 0.6;
	m_pTransform->zscale	= m_pRigidReg3DInfo->m_FloatVolInfo.Sw / m_pRigidReg3DInfo->m_ReferVolInfo.Sw;

	// rotating in float volume
	m_pTransform->xrot		= 0.0;
	m_pTransform->yrot		= 0.0;
	m_pTransform->zrot		= 0.0;

	return TRUE;
}