crosscorrelationagent.cpp

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

// CrossCorrelationAgent.cpp: implementation of the RxxCrossCorrelationAgent class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "fusion.h"
#include "CrossCorrelationAgent.h"

#include <math.h>

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

#define BOUNDARY 0
#define EDGE 1

#define	TRANSLATE_VECTOR	7.0
#define CIRCLEDEGREE		3.141592/180
#define CIRCLERADIAN		180/3.141592
#define THETA				2.0*CIRCLEDEGREE
#define ANGLE				2.0
#define	SCALE_VECTOR		0.0

#define FULLCOUNTER		3
#define EMPTYCOUNTER		0

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

RxxCrossCorrelationAgent::RxxCrossCorrelationAgent()
{

}

RxxCrossCorrelationAgent::~RxxCrossCorrelationAgent()
{
	if(m_ptFloatTrace){	delete[] m_ptFloatTrace;	m_ptFloatTrace=NULL; }
	if(m_ptReferTrace){	delete[] m_ptReferTrace;	m_ptReferTrace=NULL; }
	
	
	for(int z=0 ; z<m_pFloat->nVolZ ; z++){		
		delete[]	m_pFltTraceInfo[z].ptTraceContour;		
		m_pFltTraceInfo[z].nTraceCnt = 0;
	}	
	if(m_pFltTraceInfo){ delete[] m_pFltTraceInfo;	m_pFltTraceInfo = NULL;}
	
	for(z=0 ; z<m_pRefer->nVolZ ; z++){		
		delete[]	m_pRefTraceInfo[z].ptTraceContour;		
		m_pRefTraceInfo[z].nTraceCnt = 0;
	}	
	if(m_pRefTraceInfo){ delete[] m_pRefTraceInfo;	m_pRefTraceInfo = NULL;}
}

RxxCrossCorrelationAgent::RxxCrossCorrelationAgent(_MODAL_INFO* pFloat, _MODAL_INFO* pRefer, double lfGlobalScaleX,double lfGlobalScaleY,double lfGlobalScaleZ)
{
	m_lfOptimizedScaleX = 0.0;
	m_lfOptimizedScaleY = 0.0;
	m_lfOptimizedScaleZ = 0.0;

	m_lfTranslateXofFlt = 0.0;
	m_lfTranslateYofFlt = 0.0;
	m_lfTranslateZofFlt = 0.0;

	m_lfThetaX = 0.0;		
	m_lfThetaY = 0.0;		
	m_lfThetaZ = 0.0;

	m_lfRotX = 0.0;
	m_lfRotY = 0.0;
	m_lfRotZ = 0.0;
	
	m_pFloat = pFloat;
	m_pRefer = pRefer;

	m_pFltTraceInfo = NULL;
	m_pRefTraceInfo = NULL;

	m_ptFloatTrace = NULL;
	m_ptReferTrace = NULL;
	m_nFloatTraceCnt = 0;
	m_nReferTraceCnt = 0;

	m_nLabelIdx = 0;
	
	m_lfGlobalScaleX = lfGlobalScaleX;
	m_lfGlobalScaleY = lfGlobalScaleY;
	m_lfGlobalScaleZ = lfGlobalScaleZ;

	m_nFltContourAve = 0;
	m_nRefContourAve = 0;
}

// Implement rigid model
void RxxCrossCorrelationAgent::ComputeMaxCorrelation(int nRelation)
{ 		
	DWORD CurrentTime = GetTickCount();

//	int nSearchSize = 1.5;
	int nSearchSize = 4;
 	double lfCorrelation = 0.0; 
//	int nMax = 0;
	double lfMax = 0.0;
//	int nMin = 1000000;	
 
 	double fScaleVectorX = m_lfGlobalScaleX;
 	double fScaleVectorY = m_lfGlobalScaleY; 
	double fScaleVectorZ = m_lfGlobalScaleZ;	

	double lfThetaVector = 1.0;	
 
 	int nTransStartX = - nSearchSize;	int nTransEndX = nSearchSize;
 	int nTransStartY = - nSearchSize;	int nTransEndY = nSearchSize; 
	int nTransStartZ = - nSearchSize;	int nTransEndZ= nSearchSize;
 	
 	double lfThetaIncrementValue = 0.2;
//	double lfThetaIncrementValue = 0.7;
	double lfTransIncrementValue = 0.5;
//	double lfTransIncrementValue = 0.5;

	/*
	int nReferVolSize = m_pRefer->nVolX*m_pRefer->nVolY*m_pRefer->nVolZ;
	BYTE *pMap = new BYTE[nReferVolSize];
	memset(pMap, 0x00, nReferVolSize);

	CPoint ptReferTC;

//	int nNeighborX[26] = {-1,0,1,-1,0,1,-1,0,1,-1,0,1,-1,1,-1,0,1,-1,0,1,-1,0,1,-1,0,1};
//	int nNeighborY[26] = {-1,-1,-1,0,0,0,1,1,1,-1,-1,-1,0,0,1,1,1,-1,-1,-1,0,0,0,1,1,1};
//	int nNeighborZ[26] = {-1,-1,-1,-1,-1,-1,-1,-1,-1,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1};

	int nNeighborX[8] = {-1,0,1,-1,1,-1,0,1};
	int nNeighborY[8] = {-1,-1,-1,0,0,1,1,1};

	for(int z=0 ; z<m_pRefer->nVolZ ; z++){
		for(int nIdx=0 ; nIdx< m_pRefTraceInfo[z].nTraceCnt ; nIdx++){
			ptReferTC = m_pRefTraceInfo[z].ptTraceContour[nIdx];			
			pMap[z*m_pRefer->nVolX*m_pRefer->nVolY+ptReferTC.y*m_pRefer->nVolX+ptReferTC.x] = 2;

			for(int i=0 ; i<8; i++){
//				pMap[(z+nNeighborZ[i])*m_pRefer->nVolX*m_pRefer->nVolY+(ptReferTC.y+nNeighborY[i])*m_pRefer->nVolX+ptReferTC.x+nNeighborX[i]] = 1;
				pMap[z*m_pRefer->nVolX*m_pRefer->nVolY+(ptReferTC.y+nNeighborY[i])*m_pRefer->nVolX+ptReferTC.x+nNeighborX[i]] = 1;
			}
		}
	}
	*/
	int nReferVolSize = m_pRefer->nVolX*m_pRefer->nVolY*m_pRefer->nVolZ;
	BYTE *pMap = new BYTE[nReferVolSize];
	memset(pMap, 0x00, nReferVolSize);

	CPoint ptReferTC;	
	for(int z=0 ; z<m_pRefer->nVolZ ; z++){
		for(int nIdx=0 ; nIdx< m_pRefTraceInfo[z].nTraceCnt ; nIdx++){
			ptReferTC = m_pRefTraceInfo[z].ptTraceContour[nIdx];			
			pMap[z*m_pRefer->nVolX*m_pRefer->nVolY+ptReferTC.y*m_pRefer->nVolX+ptReferTC.x] = 5;

			int nIdx;
			for(int y=-1 ; y<1 ; y++){
				for(int x=-1 ; x<1 ; x++){
					nIdx = z*m_pRefer->nVolX*m_pRefer->nVolY+(ptReferTC.y+y)*m_pRefer->nVolX+ptReferTC.x+x;
					if(x==0 && y==0 || pMap[nIdx]>=4)	continue;
					pMap[nIdx] = 4;

					for(int yy=-1 ; yy<1 ; yy++){
						for(int xx=-1 ; xx<1 ; xx++){
							nIdx = z*m_pRefer->nVolX*m_pRefer->nVolY+(ptReferTC.y+y+yy)*m_pRefer->nVolX+ptReferTC.x+x+xx;
							if(x==0 && y==0 || pMap[nIdx]>=3)	continue;
							pMap[nIdx] = 3;

							for(int yyy=-1 ; yyy<1 ; yyy++){
								for(int xxx=-1 ; xxx<1 ; xxx++){
									nIdx = z*m_pRefer->nVolX*m_pRefer->nVolY+(ptReferTC.y+y+yy+yyy)*m_pRefer->nVolX+ptReferTC.x+x+xx+xxx;
									if(x==0 && y==0 || pMap[nIdx]>=2)	continue;
									pMap[nIdx] = 2;
								}
							}
						}
					}
				}
			}
		}
	}
 	 
	//Full search!!
	for(double thetaZ = -lfThetaVector ; thetaZ<=lfThetaVector ; thetaZ += lfThetaIncrementValue){
		for(double thetaY = -lfThetaVector ; thetaY<=lfThetaVector ; thetaY += lfThetaIncrementValue){
			for(double thetaX = -lfThetaVector ; thetaX<=lfThetaVector ; thetaX += lfThetaIncrementValue){
				for(double k= (double)nTransStartZ ; k <= nTransEndZ; k = k+lfTransIncrementValue){
 					for(double j= (double)nTransStartY ; j <= nTransEndY; j = j+lfTransIncrementValue){
 						for(double i= (double)nTransStartX ; i <= nTransEndX; i = i+lfTransIncrementValue){
 							double RadThetaX = (double)(thetaX * CIRCLEDEGREE);
							double RadThetaY = (double)(thetaY * CIRCLEDEGREE);
							double RadThetaZ = (double)(thetaZ * CIRCLEDEGREE);

//							nCorrelation = ComputeCorrelation(m_pRefer->pVol, m_ptFloatTrace, i, j, k, m_lfGlobalScaleX, m_lfGlobalScaleY, m_lfGlobalScaleZ, RadThetaX, RadThetaY,RadThetaZ);
//							nCorrelation = ComputeCorrelation(pMap, m_pFltTraceInfo, i, j, k, m_lfGlobalScaleX, m_lfGlobalScaleY, m_lfGlobalScaleZ, RadThetaX, RadThetaY,RadThetaZ);
							lfCorrelation = ComputeCorrelation(pMap, m_pFltTraceInfo, i, j, k, m_lfGlobalScaleX, m_lfGlobalScaleY, m_lfGlobalScaleZ, RadThetaX, RadThetaY, RadThetaZ);
							TRACE(_T("Scale=(%lf,%lf,%lf), Translate=(%lf,%lf,%lf), Theta=(%lf,%lf,%lf), Correlation=%lf\n"),m_lfGlobalScaleX, m_lfGlobalScaleY, m_lfGlobalScaleZ,i,j,k, thetaX, thetaY, thetaZ, lfCorrelation);

 							if(lfMax<lfCorrelation){
 		//						nMax = nCorrelation;
								lfMax = lfCorrelation;
// 								KeepParameter(m_lfGlobalScaleX, m_lfGlobalScaleY, m_lfGlobalScaleZ, i, j, k,RadThetaX, RadThetaY, RadThetaZ);
								KeepParameter(m_lfGlobalScaleX, m_lfGlobalScaleY, m_lfGlobalScaleZ, i, j, k,thetaX, thetaY, thetaZ);
 							}
 						}
 					}
				}
			}
		}
	}
	delete[] pMap; 

	DWORD ElapsedTime = GetTickCount() - CurrentTime;	
//	QueryPerformanceCounter();		

	double lfTransXMM = m_lfTranslateXofFlt*m_pFloat->lfVoxelX;
	double lfTransYMM = m_lfTranslateYofFlt*m_pFloat->lfVoxelY;
	double lfTransZMM = m_lfTranslateZofFlt*m_pFloat->lfVoxelZ;	


	CString str;
	str.Format(_T("S = (%0.3lf, %0.3lf, %0.3lf)\r\nT = (%0.3lf, %0.3lf, %0.3lf)\r\nT = (%0.3lfmm, %0.3lfmm, %0.3lfmm)\r\nR = (%0.3lf, %0.3lf, %0.3lf)\r\n\r\nMaxCC = %lf\r\n\r\nTotal Processing Time = %dms"), 
				m_lfOptimizedScaleX, m_lfOptimizedScaleY, m_lfOptimizedScaleZ, 
				m_lfTranslateXofFlt, m_lfTranslateYofFlt, m_lfTranslateZofFlt,
				lfTransXMM, lfTransYMM, lfTransZMM,
				m_lfThetaX, m_lfThetaY, m_lfThetaZ, lfMax, ElapsedTime);				
	
	AfxMessageBox(str);
}

double RxxCrossCorrelationAgent::ComputeCorrelation(unsigned short *pTemp, _FPoint *ptBndyTrace, double lfTransX, double lfTransY, double lfTransZ,
										 double lfScaleX, double lfScaleY, double lfScaleZ, double fThetaX, double fThetaY, double fThetaZ)
{	
	int nBndyptCnt = m_nFloatTraceCnt;

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

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

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

	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;

	double lfCorrelation=0.0;
	
	for(int i=0;i<nBndyptCnt;i++){
		double lfShiftOriginX = (double)((ptBndyTrace[i].x - lfFltCenterX)*lfScaleX);
		double lfShiftOriginY = (double)((ptBndyTrace[i].y - lfFltCenterY)*lfScaleY);
		double lfShiftOriginZ = (double)((ptBndyTrace[i].z - 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;

		lfCorrelation += PartialVolumeInterplation(lfTransformedX, lfTransformedY, lfTransformedZ, (unsigned char*)pTemp);

		/*
		int nTransformedX = (int)(lfTransformedX+0.5);
		int nTransformedY = (int)(lfTransformedY+0.5);
		int nTransformedZ = (int)(lfTransformedZ+0.5);

		int nLoc = (int)(m_pRefer->nVolX*m_pRefer->nVolY*nTransformedZ + m_pRefer->lfVoxelX*nTransformedY + nTransformedX);
//		if(pTemp[nLoc]!=0)		nCnt++;
//		if(pTemp[nLoc]==0)	nCnt++;
		nCnt+=pTemp[nLoc];
		*/
	}		
	return lfCorrelation;
}

double RxxCrossCorrelationAgent::ComputeCorrelation(unsigned char *pRefMap, _TraceInfo *pFltTraceInfo, double lfTransX, double lfTransY, double lfTransZ,
										 double lfScaleX, double lfScaleY, double lfScaleZ, double fThetaX, double fThetaY, double fThetaZ)
{	
	double lfCosThetaX = cos(fThetaX);
	double lfSinThetaX = sin(fThetaX);

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

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

	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;

	double lfCorrelation=0.0;	
	
	for(int z = 0; z<m_pFloat->nVolZ ; z++){
		if(m_nFltContourAve > pFltTraceInfo[z].nTraceCnt)	continue;
		for(int i=0;i<pFltTraceInfo[z].nTraceCnt;i+=1){
			double lfShiftOriginX = (double)((pFltTraceInfo[z].ptTraceContour[i].x - lfFltCenterX)*lfScaleX);
			double lfShiftOriginY = (double)((pFltTraceInfo[z].ptTraceContour[i].y - lfFltCenterY)*lfScaleY);