rxysegment_temp.cpp

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

	// to replace the value of background to 65535
	// Background means outside the chest.
	// SRG Inside the bounding box
	if(bSrgSuccess = SRG3D_LineMinMax(0, 0, int((iSliceNum-1)*0.5),
		pTempVolume, m_pLabel, iVolX, iVolY, iSliceNum, 
		thresholdmin, thresholdmax)) {	
		ControlProgress(pWndProgress);
	}

	// SRG threshold for chest
	thresholdmin = THRESHOLD, thresholdmax = 4095;

	// 4. Second Seeded Region Growing to get lung area
	// SRG region is chest outside the lung. After SRG, we invert the region.
	// Initialize seed points of 2D SRG 
	m_ptSeed = new CPoint[iSliceNum];
	for( i=0; i<iSliceNum; i++) {
		m_ptSeed[i].x = -1;
		m_ptSeed[i].y = -1;
	}
	int offset = OFFSET;
	// code for finding seed point of each slice
	for( i = 0; i < iVolumeSize; i++) {
		if( pTempVolume[i] > THRESHOLD ) // for chest
		{
			iIndexX = (i%iSliceSize)%iVolX;
			iIndexY = (i%iSliceSize)/iVolX;
			iIndexZ = (int)(i/iSliceSize);
			iSeedIndex = i;
			if(m_ptSeed[iIndexZ].x == -1) {
				// code for removing seed point in isolated area
				if( i > offset*iVolX ) {
					if( (pTempVolume[i+offset] > THRESHOLD) && (pTempVolume[i-offset] > THRESHOLD) && 
						(pTempVolume[i+offset*iVolX] > THRESHOLD) && (pTempVolume[i-offset*iVolX] > THRESHOLD))
						
						m_ptSeed[iIndexZ] = CPoint( iIndexX, iIndexY );
				}
				else {
					if( (pTempVolume[i+offset] > THRESHOLD) && (pTempVolume[i-offset] > THRESHOLD) && 
						(pTempVolume[i+offset*iVolX] > THRESHOLD) )
						
						m_ptSeed[iIndexZ] = CPoint( iIndexX, iIndexY );
				}

			}
			else i += iSliceSize;
		}
	}
	for( i=0; i<iVolumeSize; i++) {
		iIndexX = (i%iSliceSize)%iVolX;
		iIndexY = (i%iSliceSize)/iVolX;
		iIndexZ = i/iSliceSize;

		if(m_pLabel[i] != 0) // inside of the bounding box
			pTempVolume[i] = 65535;	
	}
	ControlProgress(pWndProgress);

	memset(m_pLabel, 0x00, sizeof(unsigned char)*iVolumeSize);
	thresholdmin = THRESHOLD, thresholdmax = 4095;
	
	for(i=0 ; i<iSliceNum ; i++) 
		SRG2D_MinMax(m_ptSeed[i].x, m_ptSeed[i].y, 
			&(pTempVolume[i*iSliceSize]), &(m_pLabel[i*iSliceSize]), iVolX, iVolY, 
			thresholdmin, thresholdmax );			
	ControlProgress(pWndProgress);

	for( i=0; i<iVolumeSize; i++) {
		int y = (i%(iSliceSize))/iVolX;
		if(m_pLabel[i] != 0) // chest
			pTempVolume[i] = 65535;
		else if( ptMaxBound.y - y < (ptMaxBound.y-ptMinBound.y)*REMOVEBACK ) // 60 : bounding box狼 10%
			pTempVolume[i] = 65535;
	}	
	// 企 何盒父 255肺 汲沥.(企客 扁包瘤 + alpha)
	for( i=0; i<iVolumeSize; i++) {	
		if(pTempVolume[i]==65535) { // 企 官冰 康开 
			pTempVolume[i] = 0;
			m_pLabel[i] = 0;
		}
		else  {
			pTempVolume[i] = 0x0fff;
			m_pLabel[i] = 255;
		}
	}
	// 5. Multi-resolution connected component labeling(企客 扁包瘤父 吧矾晨)
	// reduce resolution of x and y axis in 2 times
	// but not z axis because restoration of volume is not easy.
	int iLowCount = 0;
	RxyConnectedComponent *pConnected;
	int iSliceOffset, iLowLabelSize, iLowSliceSize;

//	if(iSliceNum%4==0) iSliceOffset = iSliceNum/4;
//	else 
		iSliceOffset = iSliceNum/4 + iSliceNum%4;
	
	iLowSliceSize = iVolX*iVolY/16;
	iLowLabelSize = iLowSliceSize*iSliceOffset;
	RxPoint3D<float> ptRestoreSeed;

	if(iSliceNum > 300) {
		m_pLowLabel = new unsigned char[iLowLabelSize];
		for( i=0; i<iVolumeSize; i++) {
			iIndexX = (i%(iSliceSize))%iVolX;
			iIndexY = (i%(iSliceSize))/iVolX;
			iIndexZ = i/iSliceSize;

			if( iIndexX%4==0 && iIndexY%4==0 && iIndexZ%4==0) 
				m_pLowLabel[iLowCount++] = m_pLabel[i];
		}
		if(m_pLabel)			delete []m_pLabel;
		ControlProgress(pWndProgress);

		// Select the largest connected component 
		pConnected = new RxyConnectedComponent;
		label = pConnected->GetConnectedComponent( m_pLowLabel, iVolX/4, iVolY/4, iSliceOffset, 1, iVolumeSize/(100*64));
		delete pConnected;
		ControlProgress(pWndProgress);
		
		// restoration of m_pLabel from m_pLowLabel
		iLowCount = 0;

		for(i=0; i<iVolumeSize; i++) {
			iIndexX = (i%(iSliceSize))%iVolX;
			iIndexY = (i%(iSliceSize))/iVolX;
			iIndexZ = i/iSliceSize;

			if( iIndexX%4==0 && iIndexY%4==0 && iIndexZ%4==0 ) { 
				if( m_pLowLabel[iLowCount] && m_pLowLabel[iLowCount+offset] && m_pLowLabel[iLowCount+iVolX/4*offset] 
					&& m_pLowLabel[iLowCount+iLowSliceSize*offset] && m_pLowLabel[iLowCount+offset*2] && m_pLowLabel[iLowCount+iVolX/4*offset*2] 
					&& m_pLowLabel[iLowCount+iLowSliceSize*offset*2]) {
					ptRestoreSeed = RxPoint3D<float>(iIndexX, iIndexY, iIndexZ);
					break;
				}
				iLowCount++;
			}
		}
		thresholdmin = 4000, thresholdmax = 4095;
		m_pLabel = new unsigned char[iVolumeSize];
		memset(m_pLabel, 0x00, sizeof(unsigned char)*iVolumeSize);

		SRG3D_LineMinMax(ptRestoreSeed.x, ptRestoreSeed.y, ptRestoreSeed.z,
			pTempVolume, m_pLabel, iVolX, iVolY, iSliceNum, 
			thresholdmin, thresholdmax );

		delete []m_pLowLabel;
	}
	else {
			pConnected = new RxyConnectedComponent;
			label = pConnected->GetConnectedComponent(m_pLabel, iVolX, iVolY, iSliceNum, 1, iVolumeSize/(100));
			delete pConnected;
			ControlProgress(pWndProgress);
	}
	// 企+扁包瘤 何盒父 吧矾晨 
	for( i=0; i<iVolumeSize; i++) {	
		if( m_pLabel[i] )	pTempVolume[i] = 0x0fff;
		else	pTempVolume[i] = 0;
	}
	ControlProgress(pWndProgress);

	SRGTime = GetTickCount () - CurrentTime;
	g_outFile << "SRG time: " << SRGTime << '\n';
	CurrentTime = GetTickCount ();

	// 6. Erosion using circle mask by yeni
	int nErose =  NUM_OF_EROSION; 	
	unsigned char* pSrc = new unsigned char[iSliceSize];					
	unsigned char* pTmp = new unsigned char[iSliceSize];

	for(i=0 ; i<iSliceNum ; i++) 
	{
		memcpy(pSrc, &(m_pLabel[i*iSliceSize]), sizeof(unsigned char)*iSliceSize);
		Erose_CircleMask(pSrc, pTmp, iVolX, iVolY, nErose);//, pErosionMask);
		//Erose(pSrc, pTmp, iVolX, iVolY, nErose);
		memcpy( &(m_pLabel[i*iSliceSize]), pTmp, sizeof(unsigned char)*iSliceSize);		
		if(i%2==1)
			ControlProgress(pWndProgress);
	}
	delete[] pSrc;
	delete[] pTmp;

	ErosionTime = GetTickCount () - CurrentTime;
//	g_outFile << "Erosion time: " << ErosionTime << '\n';
	CurrentTime = GetTickCount ();

	// 7. Multi-resolution connected component labeling(企父 急琶)
	
	// reduce resolution of x and y axis in 2 times
	// but not z axis because restoration of volume is not easy.
	if(iSliceNum > 300) {
		iLowCount = 0;
		m_pLowLabel = new unsigned char[iLowLabelSize];

		for(i=0; i<iVolumeSize; i++) {
			iIndexX = (i%iSliceSize)%iVolX;
			iIndexY = (i%iSliceSize)/iVolX;
			iIndexZ = i/iSliceSize;

			if( iIndexX%4==0 && iIndexY%4==0 && iIndexZ%4==0 )
				m_pLowLabel[iLowCount++] = m_pLabel[i];
		}
		ControlProgress(pWndProgress);

		// Select the largest connected component 
		pConnected = new RxyConnectedComponent;
		label = pConnected->GetConnectedComponent(m_pLowLabel, iVolX/4, iVolY/4, iSliceOffset, 2, iVolumeSize/(100*64));
		delete pConnected;
		ControlProgress(pWndProgress);

		// restoration of m_pLabel from m_pLowLabel
		// use the value of nearest neighbors
		iLowCount = 0;
		RxPoint3D<float> ptSeedLeft, ptSeedRight;
		ptSeedLeft.x = -1, ptSeedRight.x = -1;
		thresholdmin = 4000, thresholdmax = 4095;

		for(i=0; i<iVolumeSize; i++) {
			iIndexX = (i%(iSliceSize))%iVolX;
			iIndexY = (i%(iSliceSize))/iVolX;
			iIndexZ = i/iSliceSize;

			if( iIndexX%4==0 && iIndexY%4==0 && iIndexZ%4==0 ) {
				if(ptSeedLeft.x !=-1 && ptSeedRight.x !=-1)
					break;
				if(ptSeedLeft.x==-1 && m_pLowLabel[iLowCount]==0xff ) { 
					ptSeedLeft = RxPoint3D<float>(iIndexX, iIndexY, iIndexZ);
				}
				else if(ptSeedRight.x==-1 && m_pLowLabel[iLowCount]==0xee ) {
					ptSeedRight = RxPoint3D<float>(iIndexX, iIndexY, iIndexZ);
				}
				iLowCount++;
			}
		}
		delete []m_pLowLabel;

		unsigned short *pTempVolume2 = new unsigned short[iVolumeSize];
		for(i=0; i<iVolumeSize; i++) {
			if(m_pLabel[i]) pTempVolume2[i] = 0x0fff; // yeni0616
			else pTempVolume2[i] = 0x0000;
		}
		memset(m_pLabel, 0x00, sizeof(unsigned char)*iVolumeSize);

		SRG3D_LineMinMax(ptSeedLeft.x, ptSeedLeft.y, ptSeedLeft.z,
			pTempVolume2, m_pLabel, iVolX, iVolY, iSliceNum, 
			thresholdmin, thresholdmax);

		unsigned char *pTempLabel = new unsigned char[iVolumeSize];
		memset(pTempLabel, 0, sizeof(unsigned char)*iVolumeSize);

		SRG3D_LineMinMax(ptSeedRight.x, ptSeedRight.y, ptSeedRight.z,
			pTempVolume2, pTempLabel, iVolX, iVolY, iSliceNum, 
			thresholdmin, thresholdmax );

		delete []pTempVolume2;
		for(i=0; i<iVolumeSize; i++) {
			if(pTempLabel[i])
				m_pLabel[i] = 0xff;
		}
		delete []pTempLabel;
		
		ControlProgress(pWndProgress);
	}
	else {
		pConnected = new RxyConnectedComponent;
		label = pConnected->GetConnectedComponent(m_pLabel, iVolX, iVolY, iSliceNum, 2, iVolumeSize/(100));
		delete pConnected;
		ControlProgress(pWndProgress);
	}

	CCTime = GetTickCount () - CurrentTime;
//	g_outFile << "Multi-resolution CC time: " << CCTime << '\n';
	CurrentTime = GetTickCount ();

	// 8. Dilation using circle mask by yeni
	pSrc = new unsigned char[iSliceSize];
	pTmp = new unsigned char[iSliceSize];

	for(i=0 ; i<iSliceNum ; i++) 
	{
		memcpy(pSrc, &(m_pLabel[i*iSliceSize]), sizeof(unsigned char)*iSliceSize);
		Dilate_CircleMask(pSrc, pTmp, iVolX, iVolY, nErose);
		//Dilate(pSrc, pTmp, iVolX, iVolY, nErose);
		memcpy( &(m_pLabel[i*iSliceSize]), pTmp, sizeof(unsigned char)*iSliceSize);		
		if(i%2==1)
			ControlProgress(pWndProgress);
	}

	delete[] pSrc;
	delete[] pTmp;

	DilationTime = GetTickCount () - CurrentTime;
//	g_outFile << "Dilation time: " << DilationTime << '\n';
	g_outFile << "Airway Separation time: " << ErosionTime+DilationTime+CCTime << '\n';
	CurrentTime = GetTickCount ();

	// 9. clear boundary
	// m_pLabel == 0牢 汗伎甸俊 措秦 connected component
	unsigned short *pTempVolume2 = new unsigned short[iVolumeSize];
	for( i=0; i<iVolumeSize; i++) {
		// erosion-cc-dilation 搬苞 
		if( pTempVolume[i]==0x0fff && m_pLabel[i]==0 ) {
			m_pLabel[i] = 255;
			pTempVolume2[i] = 0x0fff;
		}
		else { 
			m_pLabel[i] = 0;
			pTempVolume2[i] = 0;
		}
	}		

	// Multi-resolution connected component labeling(魔侥等 康开 吝 扁包瘤父 急琶)
	// reduce resolution of x and y axis in 2 times
	// but not z axis because restoration of volume is not easy.
	// if slicenum is large, execute low-resolution connected component.
	// else, execute original connected component.
	if(iSliceNum > 300) {
		iLowCount = 0;
		m_pLowLabel = new unsigned char[iLowLabelSize];

		for(i=0; i<iVolumeSize; i++) {
			iIndexX = (i%(iSliceSize))%iVolX;
			iIndexY = (i%(iSliceSize))/iVolX;
			iIndexZ = i/iSliceSize;

			if( iIndexX%4==0 && iIndexY%4==0 && iIndexZ%4==0 )
				m_pLowLabel[iLowCount++] = m_pLabel[i];
		}
		if(m_pLabel)			delete []m_pLabel;
		ControlProgress(pWndProgress);
		
		// Select the largest connected component 
		pConnected = new RxyConnectedComponent;
		label = pConnected->GetConnectedComponent(m_pLowLabel, iVolX/4, iVolY/4, iSliceOffset, 1, iVolumeSize/(64*10000));
		delete pConnected;
		ControlProgress(pWndProgress);

		// restoration of m_pLabel from m_pLowLabel
		iLowCount = 0;
		offset = 2;

		for(i=0; i<iVolumeSize; i++) {
			iIndexX = (i%(iSliceSize))%iVolX;
			iIndexY = (i%(iSliceSize))/iVolX;
			iIndexZ = i/iSliceSize;

			if( iIndexX%4==0 && iIndexY%4==0 && iIndexZ%4==0) { 
				if( m_pLowLabel[iLowCount] && m_pLowLabel[iLowCount+offset] && m_pLowLabel[iLowCount+iVolX/4*offset] 
					&& m_pLowLabel[iLowCount+iLowSliceSize*offset] && m_pLowLabel[iLowCount+offset*2] && m_pLowLabel[iLowCount+iVolX/4*offset*2] 
					&& m_pLowLabel[iLowCount+iLowSliceSize*offset*2]) {
					ptRestoreSeed = RxPoint3D<float>(iIndexX, iIndexY, iIndexZ);
					break;
				}
				iLowCount++;
			}
		}
		delete []m_pLowLabel;

		thresholdmin = 4000, thresholdmax = 4096;
		m_pLabel = new unsigned char[iVolumeSize];
		memset(m_pLabel, 0x00, sizeof(unsigned char)*iVolumeSize);

		SRG3D_LineMinMax(ptRestoreSeed.x, ptRestoreSeed.y, ptRestoreSeed.z,
			pTempVolume2, m_pLabel, iVolX, iVolY, iSliceNum, 
			thresholdmin, thresholdmax );
	}
	else {
		pConnected = new RxyConnectedComponent;
		label = pConnected->GetConnectedComponent(m_pLabel, iVolX, iVolY, iSliceNum, 1, iVolumeSize/(10000));
		delete pConnected;
		ControlProgress(pWndProgress);
	}
	

	// m_pLabel == 1牢 何盒阑 盔贰 杭俘俊辑 力芭 
	for( i=0; i<iVolumeSize; i++ ) {
		// SRG 搬苞俊辑 big component(ex. 扁包瘤)父 猾 搬苞 
		if( pTempVolume[i]==0x0fff && m_pLabel[i]) 
			pTempVolume[i] = 0x0000;
	}
	ControlProgress(pWndProgress);

	delete []pTempVolume2;	
/*
	for(i=0; i<iVolumeSize; i++) {
		if(m_pLabel[i])
			pTempVolume[i] = 0x0fff;
		else
			pTempVolume[i] = 0x0000;
	}
*/
	// Edge Detection
	// store segmentation results to mask volume of pVolume
	if(m_pLabel) delete []m_pLabel;

	EdgeDetection(pTempVolume, iVolX, iVolY, iSliceNum, nMode, pWndProgress, ptMaxBound, ptMinBound); 

	ControlProgress(pWndProgress);

	delete []m_ptSeed;
	delete []pTempVolume;

	if (pWndProgress) {	delete pWndProgress;	pWndProgress = NULL;	}	

	CorrectionTime = GetTickCount () - CurrentTime;
	g_outFile << "False Negative Correction time: " << CorrectionTime << '\n';

	ElapsedTime = SRGTime + ErosionTime + CCTime + DilationTime + CorrectionTime;
	g_outFile << "Volume Segmentation Time: " << ElapsedTime << '\n';
/*	g_outFile << "x : " << iVolX << '\n';
	g_outFile << "y : " << iVolY << '\n';
	g_outFile << "z : " << iSliceNum << '\n';*/
	g_outFile << "Bounding box size: " << iVolumeSize << '\n';

	g_outFile.close();

	(RxGetFrameMain()->m_pFMWndVR[nMode])->RedrawWnd(TRUE);
	(RxGetFrameMain()->m_pFMWnd3DLocalA[nMode])->Redraw();
	(RxGetFrameMain()->m_pFMWnd3DLocalC[nMode])->Redraw();
	(RxGetFrameMain()->m_pFMWnd3DLocalS[nMode])->Redraw();	
}

void RxySegment::ControlProgress(RxProgressWnd *pProWnd)
{
	if(pProWnd && pProWnd->Cancelled()){
		delete pProWnd;
		pProWnd = NULL;
		return;
	}			
	else if(pProWnd){
		pProWnd->StepIt();
		pProWnd->PeekAndPump();
	}
}

void RxySegment::NoduleSegmentation(int nMode)
{
	// for calculating elapsed time
	DWORD CurrentTime, ElapsedTime;
	char TimeStr[255];
	int i, j, k;
	int iVolX, iVolY, iSliceNum;
	RxVolumeData *pVolumeData;
	long nImageSize, nSliceSize;
	int tempX, tempY, tempZ, seedIndex;
	RxPoint3D<double> ptNudulePos;
	RxPoint3D<float> ptMaxBound(-100, -100, -100), ptMinBound(5000, 5000, 5000);
	int label = 0;
	int thresholdmin, thresholdmax;
	int NoduleSizeXY = 50, NoduleSizeZ = 10;
	int index;
	unsigned char* pLabel;
	unsigned short* pTempVol;

	int nNumNodule;
	int threshold[100];
	int nErose[100];
	int RefNodulePos[100];
	int FloatNodulePos[100];
//		1 全扁后
//		2 蜡公磊
//		3 捞犁帮
//		4 辫荐沥
//		5 颊措牢, Float 2锅掳 辆剧 92锅掳 浇扼捞胶
//	    6 捞荐磊
//		7 决沥巢
//	    8 捞篮雀, 2锅掳 巴 肋 救凳
//      9 烙篮惑, 2锅掳 巴 肋 救凳
//      10 捞荐康, 4锅掳 巴 肋 救凳
	int nPatient;

	if (PATIENT_1)
		nPatient = 1;
	if (PATIENT_2)
		nPatient = 2;
	if (PATIENT_3)