crosscorrelationagent.cpp

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

			double lfShiftOriginZ = (double)((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;

			
			if(lfTransformedX>=0 && lfTransformedX<m_pRefer->nVolX-1 &&
				lfTransformedY>=0 && lfTransformedY<m_pRefer->nVolY-1 &&
					lfTransformedZ>=0 && lfTransformedZ<m_pRefer->nVolZ-1)
				lfCorrelation += PartialVolumeInterplation(lfTransformedX, lfTransformedY, lfTransformedZ, pRefMap);
//			lfCorrelation += NearestInterpolation(lfTransformedX, lfTransformedY, lfTransformedZ, pRefMap);	

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

			if(nTransformedZ >=0 && nTransformedZ<m_pRefer->nVolZ &&
				nTransformedY >=0 && nTransformedY <m_pRefer->nVolY &&
				nTransformedX >=0 && nTransformedX <m_pRefer->nVolX){

				nLoc = (int)(m_pRefer->nVolX*m_pRefer->nVolY*nTransformedZ + m_pRefer->nVolX*nTransformedY + nTransformedX);
				nCnt+=pRefMap[nLoc];
			}
			*/
		}	
	}
	return lfCorrelation;
}

void RxxCrossCorrelationAgent::KeepParameter(double lfScaleVectorX, double lfScaleVectorY, double lfScaleVectorZ, double lfTransX, double lfTransY, double lfTransZ, double thetaX, double thetaY, double thetaZ)
{	
	m_lfTranslateXofFlt = lfTransX;
	m_lfTranslateYofFlt = lfTransY;
	m_lfTranslateZofFlt = lfTransZ;

	m_lfOptimizedScaleX = lfScaleVectorX;
	m_lfOptimizedScaleY = lfScaleVectorY;
	m_lfOptimizedScaleZ = lfScaleVectorZ;

	/*
	m_lfThetaX = thetaX;
	m_lfThetaY = thetaY;
	m_lfThetaZ = thetaZ;

	m_lfRotX = m_lfThetaX/CIRCLEDEGREE;
	m_lfRotY = m_lfThetaY/CIRCLEDEGREE;
	m_lfRotZ = m_lfThetaZ/CIRCLEDEGREE;
	*/

	//degree
	m_lfRotX = thetaX;
	m_lfRotY = thetaY;
	m_lfRotZ = thetaZ;

	//radian
	m_lfThetaX = m_lfRotX*CIRCLEDEGREE;
	m_lfThetaY = m_lfRotY*CIRCLEDEGREE;
	m_lfThetaZ = m_lfRotZ*CIRCLEDEGREE;
}

BOOL RxxCrossCorrelationAgent::GetOutBoundary()
{
	int nTh = 130;
	int nFltVolSize = m_pFloat->nVolX*m_pFloat->nVolY*m_pFloat->nVolZ;
//	unsigned short* pFloatVol = new unsigned short[nFltVolSize];
//	memset(pFloatVol, 0x00, sizeof(unsigned short)*nFltVolSize);
//	for(int nIdx=0 ; nIdx<nFltVolSize ; nIdx++){
//		if(m_pFloat->pVol[nIdx]>=nTh){
//			pFloatVol[nIdx] = m_pFloat->pVol[nIdx];			
//		}
//	}
	
//	int nPixelCnt=0;
	int nCenterIdx = 13;
	int nMaskSize = 3;
	int nMargin = nMaskSize/2;
	int nLength = 0;
	_FPoint* pTempBndy = new _FPoint[nFltVolSize];
	memset(pTempBndy,0x00, sizeof(_FPoint)*nFltVolSize);

	for(int z=1 ; z<m_pFloat->nVolZ-1 ; z++){
		for(int y=1 ; y<m_pFloat->nVolY-1 ; y++){
			for(int x=1 ; x<m_pFloat->nVolX-1 ; x++){
				int nIdx = z*m_pFloat->nVolX*m_pFloat->nVolY+y*m_pFloat->nVolX+x;				
				BOOL bStatus= FALSE;

				int nPreZIdx = (z-1)*m_pFloat->nVolX*m_pFloat->nVolY+y*m_pFloat->nVolX+x;
				int nNextZIdx = (z+1)*m_pFloat->nVolX*m_pFloat->nVolY+y*m_pFloat->nVolX+x;
				int nPreYIdx = m_pFloat->nVolX*m_pFloat->nVolY+(y-1)*m_pFloat->nVolX+x;
				int nNextYIdx = m_pFloat->nVolX*m_pFloat->nVolY+(y+1)*m_pFloat->nVolX+x;
				int nPreXIdx = m_pFloat->nVolX*m_pFloat->nVolY+y*m_pFloat->nVolX+x-1;
				int nNextXIdx = m_pFloat->nVolX*m_pFloat->nVolY+y*m_pFloat->nVolX+x+1;

				if(m_pFloat->pVol[nIdx]<nTh){
					if(m_pFloat->pVol[nPreZIdx]>=nTh){	bStatus = TRUE; }
					else if(m_pFloat->pVol[nNextZIdx]>=nTh){ bStatus = TRUE; }
					else if(m_pFloat->pVol[nPreYIdx]>=nTh){ bStatus = TRUE; }
					else if(m_pFloat->pVol[nNextYIdx]>=nTh){ bStatus = TRUE; }
					else if(m_pFloat->pVol[nPreXIdx]>=nTh){ bStatus = TRUE; }
					else if(m_pFloat->pVol[nNextXIdx]>=nTh){ bStatus = TRUE; }
				}

				if(bStatus){
					pTempBndy[nLength].x = x;
					pTempBndy[nLength].y = y;
					pTempBndy[nLength].z = z;
					nLength++;				
				}				
			}
		}
//		TRACE(_T("%d\n"), nPixelCnt);		
//		nPixelCnt =0;
	}
	if(m_ptFloatTrace){ delete[] m_ptFloatTrace;	m_ptFloatTrace=NULL; }
	m_ptFloatTrace = new _FPoint[nLength];
	memcpy(m_ptFloatTrace, pTempBndy, sizeof(_FPoint)*nLength);
	if(pTempBndy)	delete[] pTempBndy;
	m_nFloatTraceCnt = nLength;

	return TRUE;
}

#include "2DTraceContour.h"
BOOL RxxCrossCorrelationAgent::GetZaxisSampling2DTraceContour(int nFltModality, int nRefModality)
{	
	Rx2DTraceContour TraceContour;
	m_pFltTraceInfo = new _TraceInfo[m_pFloat->nVolZ];
	memset(m_pFltTraceInfo, 0x00, sizeof(_TraceInfo)*m_pFloat->nVolZ);
	int nVolIdx=0;
	int nTh;
	
	// PET = 40 SKIN
	// PET = 100 Brain
	if(nFltModality == PET)		nTh = 40;	
	else if(nFltModality == MR)		nTh = 70;
	else if(nFltModality == CT)		nTh = 150; // FOR CTA
	int nFltSliceSize = m_pFloat->nVolX*m_pFloat->nVolY;
	
	for(int z=0 ; z<m_pFloat->nVolZ ; z++){

		BOOL bContour = FALSE;
		int nZIdx = z*m_pFloat->nVolX*m_pFloat->nVolY;
		BYTE* pMap = new BYTE[nFltSliceSize];
		memset(pMap, 0x00, nFltSliceSize);

		for(int nSliceIdx = 0 ; nSliceIdx < nFltSliceSize ; nSliceIdx++){
			nVolIdx = nZIdx+nSliceIdx;
			if(m_pFloat->pVol[nVolIdx]>=nTh){
				pMap[nSliceIdx] = 255;
				bContour = TRUE;
			}			
		}

		if(bContour){
			m_pFltTraceInfo[z].ptTraceContour = TraceContour.TraceContour(pMap, &m_pFltTraceInfo[z].nTraceCnt, m_pFloat->nVolX, m_pFloat->nVolY);
			m_nFltContourAve+=m_pFltTraceInfo[z].nTraceCnt;			
		}
		TRACE(_T("(Row=%d, TraceCnt=%d)\n"),z, m_pFltTraceInfo[z].nTraceCnt);
		delete[] pMap;
	}
	m_nFltContourAve/=z;

	m_pRefTraceInfo = new _TraceInfo[m_pRefer->nVolZ];
	memset(m_pRefTraceInfo, 0x00, sizeof(_TraceInfo)*m_pRefer->nVolZ);

	// MR = 70;
	if(nFltModality == PET)		nTh = 40;	
	else if(nFltModality == MR)		nTh = 70;
	else if(nFltModality == CT)		nTh = 150;  //For CTA
	
	int nRefSliceSize = m_pRefer->nVolX*m_pRefer->nVolY;

	for(z=0 ; z<m_pRefer->nVolZ ; z++){
		BOOL bContour = FALSE;
		int nZIdx = z*m_pRefer->nVolX*m_pRefer->nVolY;
		BYTE* pMap = new BYTE[nRefSliceSize];
		memset(pMap, 0x00, nRefSliceSize);
		for(int nSliceIdx = 0 ; nSliceIdx < nRefSliceSize ; nSliceIdx++){
			nVolIdx = nZIdx+nSliceIdx;
			if(m_pRefer->pVol[nVolIdx]>=nTh){
				pMap[nSliceIdx] = 255;
				bContour = TRUE;
			}			
		}
		if(bContour){
			m_pRefTraceInfo[z].ptTraceContour = TraceContour.TraceContour(pMap, &m_pRefTraceInfo[z].nTraceCnt, m_pRefer->nVolX, m_pRefer->nVolY);
			m_nRefContourAve+=m_pRefTraceInfo[z].nTraceCnt;
		}
		TRACE(_T("(Row=%d, TraceCnt=%d)\n"),z, m_pRefTraceInfo[z].nTraceCnt);
		delete[] pMap;		
	}
	m_nRefContourAve/=z;

	return TRUE;
}


void RxxCrossCorrelationAgent::FindOptimizedCoodinate(int nRelation)
{
	DWORD CurrentTime = GetTickCount();	

	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;
								}
							}
						}
					}
				}
			}
		}
	}

	
	int nCorrelation = ComputeCorrelation(pMap, m_pFltTraceInfo, 0, 0, 0, m_lfGlobalScaleX, m_lfGlobalScaleY, m_lfGlobalScaleZ, 0.0, 0.0,0.0);
	KeepParameter(m_lfGlobalScaleX, m_lfGlobalScaleY, m_lfGlobalScaleZ, 0, 0, 0, 0.0, 0.0, 0.0);
	
	double lfIncS = 0.0;
	double lfIncT = 1.0;
	double lfIncR = 0.15*CIRCLEDEGREE;

	int nIteration = 3;	
	
	
	for(int i=0 ; i<nIteration ; i++){
		nCorrelation = IncreaseTransXvector(nCorrelation, lfIncT, pMap);
		nCorrelation = DecreaseTransXvector(nCorrelation, lfIncT, pMap);

		nCorrelation = IncreaseTransYvector(nCorrelation, lfIncT, pMap);
		nCorrelation = DecreaseTransYvector(nCorrelation, lfIncT, pMap);		

		nCorrelation = CWRotZvector(nCorrelation, lfIncR, pMap);
		nCorrelation = CCWRotZvector(nCorrelation, lfIncR, pMap);		

		nCorrelation = CWRotXvector(nCorrelation, lfIncR, pMap);
		nCorrelation = CCWRotXvector(nCorrelation, lfIncR, pMap);

		nCorrelation = CWRotYvector(nCorrelation, lfIncR, pMap);
		nCorrelation = CCWRotYvector(nCorrelation, lfIncR, pMap);

		nCorrelation = IncreaseTransZvector(nCorrelation, lfIncT, pMap);
		nCorrelation = DecreaseTransZvector(nCorrelation, lfIncT, pMap);		
	}	 	
	
	delete[] pMap; 

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

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

	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(Rad) = (%0.3lf, %0.3lf, %0.3lf)\r\nR(Degree) = (%0.3lf, %0.3lf, %0.3lf)\r\n\r\nMaxCC = %d\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, 
				m_lfRotX, m_lfRotY, m_lfRotZ, nCorrelation, ElapsedTime);				
	
	AfxMessageBox(str);
}

double RxxCrossCorrelationAgent::IncreaseTransXvector(int lfPresentCC, double lfIncT, BYTE* pMap)
{	
	double CC =0.0;
	int nIsTranslateCOUNTERState = EMPTYCOUNTER;
	for(double lfTransX=m_lfTranslateXofFlt+lfIncT ; ; lfTransX +=lfIncT){		
		CC = ComputeCorrelation(pMap, m_pFltTraceInfo, lfTransX, m_lfTranslateYofFlt, m_lfTranslateZofFlt, m_lfGlobalScaleX, m_lfGlobalScaleY, m_lfGlobalScaleZ, m_lfThetaX, m_lfThetaY, m_lfThetaZ);

//		TRACE(_T("Scale=(%lf,%lf,%lf), Translate=(%lf,%lf,%lf), Theta=(%lf,%lf,%lf), MI=%d\n"),m_lfGlobalScaleX, m_lfGlobalScaleY, m_lfGlobalScaleZ,lfTransX,m_lfTranslateYofFlt,m_lfTranslateZofFlt, m_lfThetaX, m_lfThetaY, m_lfThetaZ, MI);
		if(CC>lfPresentCC){
			KeepParameter(m_lfGlobalScaleX, m_lfGlobalScaleY, m_lfGlobalScaleZ, lfTransX, m_lfTranslateYofFlt, m_lfTranslateZofFlt, m_lfThetaX, m_lfThetaY, m_lfThetaZ);
			lfPresentCC = CC;
			nIsTranslateCOUNTERState = EMPTYCOUNTER;
		}
		else if(nIsTranslateCOUNTERState!=FULLCOUNTER)	nIsTranslateCOUNTERState++;		
		else	break;		
	}
	return lfPresentCC;
}