jjskeleton.cpp

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

					}
					if (!bFind)
					{
						m_BranchPoint[m_nBranchNum].x = (short)x;
						m_BranchPoint[m_nBranchNum].y = (short)y;
						m_BranchPoint[m_nBranchNum].z = (short)z;
						m_BranchPoint[m_nBranchNum].time = m_pDistanceMapFromPS[(short)x + (short)y*m_VolX + (short)z*m_VolX*m_VolY];
						m_nBranchNum++;
					}
					
					bCont = FALSE;
				}
				else if (m_pPointType[(short)x + (short)y*m_VolX + (short)z*m_VolX*m_VolY] == BACKGROUND)
					bCont = FALSE;
				else
					m_pSkeleton->SetMask((short)x + m_OffsetX, (short)y + m_OffsetY, (short)z + m_OffsetZ);

				min = 10000.0f;
				short min_nr;

				for (int nr = 0; nr < 26; nr++)
				{
					next_x = x + nr_x[nr];
					next_y = y + nr_y[nr];
					next_z = z + nr_z[nr];
					
					if (m_pPointType[next_x + next_y*m_VolX + next_z*m_VolX*m_VolY] == KNOWN &&
						m_pDistanceMapFromPS_FAST_WAVE[next_x + next_y*m_VolX + next_z*m_VolX*m_VolY] < min)
					{
						min = m_pDistanceMapFromPS_FAST_WAVE[next_x + next_y*m_VolX + next_z*m_VolX*m_VolY];
						min_nr = nr;
					}
				}

				x = x + nr_x[min_nr];
				y = y + nr_y[min_nr];
				z = z + nr_z[min_nr];

			} while (bCont);
		}
	}
}

void JJSkeleton::FindConnectedComponent(void)
{
	int i, j, k;
	unsigned long lIndex;

	memset(m_pLabel, 0, sizeof(unsigned short)*m_VolX*m_VolY*m_VolZ);

	for (k = 0; k < m_VolZ; k++)
		for (j = 0; j < m_VolY; j++)
			for (i = 0; i < m_VolX; i++)
			{
				lIndex = i + j*m_VolX + k*m_VolX*m_VolY;

				if (m_pVesselVolume->GetMask(i + m_OffsetX, j + m_OffsetY, k + m_OffsetZ))
					m_pLabel[lIndex] = 0xff;
				else
					m_pLabel[lIndex] = 0;
			}

	GetConnectedComponent(m_VolX, m_VolY, m_VolZ);
}

int JJSkeleton::GetConnectedComponent(int nWidth, int nHeight, int nDepth)
{
	unsigned int iVolumeSize = nWidth*nHeight*nDepth;

	memset(pTable, 0, iVolumeSize*sizeof(unsigned int));

	for(int i = 0; i < iVolumeSize; i++)
		pTable[i] = i; 

	int nLabelCnt = 1;
	for(int z = 0; z < nDepth; z++) 
		for(int y = 0; y < nHeight; y++) 
			for(int x = 0; x < nWidth; x++) 
			{
				unsigned long nIndex = x + y*nWidth + z*nWidth*nHeight;
				
				if (m_pLabel[nIndex] == BKValue)	
					continue;
							
				int nState = GetState(m_pLabel, nWidth, nHeight, nDepth, x, y, z);

				SetLabelByState(m_pLabel, nWidth, nHeight, nDepth, x, y, z, &nLabelCnt, nState, pTable);
			}

	for(i = 0;i < iVolumeSize; i++) 
		pTable[i] = FindFinalLabel(pTable, i);
	
	int nMaxLabel = MapFillingWithSortedLabel(pTable, m_pLabel, iVolumeSize);

	memset(pTable, 0, iVolumeSize*sizeof(unsigned int));
	for(i = 0; i < iVolumeSize; i++)
		if (m_pLabel[i] != 0) 
			pTable[m_pLabel[i]]++;

	int nMax = 0, nMin = 0;


	if (g_nPatient <= 5 || g_nPatient >= 21)
	{
		for(i = 0; i < nMaxLabel; i++)
		{
			if (pTable[i] > 100) 
			{
				m_nObjectID[m_nObjectNum] = i;
				m_nObjectNum++;
				nMax++;
			}
			else
				nMin++;
		}
	} else
	{
		for(i = 0; i < nMaxLabel; i++)
		{
			if (pTable[i] > 100) 
			{
				m_nObjectID[m_nObjectNum] = i;
				m_nObjectNum++;
				nMax++;
			}
			else
				nMin++;
		}
	}

	if (pTable)
	{
		free(pTable);
		pTable = NULL;
	}

	return m_nObjectNum;
}

int JJSkeleton::FindFinalLabel(unsigned int *pTable, unsigned short nNowLabel)
{
	if (pTable[nNowLabel] == nNowLabel) 
		return (nNowLabel);
	else 
		return (FindFinalLabel(pTable, pTable[nNowLabel]));
}

int JJSkeleton::GetState(unsigned short *pMap, int nWidth, int nHeight, int nDepth, int x, int y, int z)
{
	int nIdx = x + y*nWidth + z*nWidth*nHeight;
	int nIdxXbefore = (x - 1) + y*nWidth + z*nWidth*nHeight;
	int nIdxYbefore = x + (y - 1)*nWidth + z*nWidth*nHeight;
	int nIdxZbefore = x + y*nWidth + (z - 1)*nWidth*nHeight;

	if (x == 0 && y == 0 && z == 0)
		return (JJSkeleton::COMP_NEW);

	else if (y == 0 && z == 0 && pMap[nIdxXbefore] != BKValue)
		return (JJSkeleton::COMP_X);
	
	else if (y == 0 && z == 0 && pMap[nIdxXbefore] == BKValue)
		return (JJSkeleton::COMP_NEW);

	else if (x == 0 && z == 0 && pMap[nIdxYbefore] != BKValue)
		return (JJSkeleton::COMP_Y);
	
	else if (x == 0 && z == 0 && pMap[nIdxYbefore] == BKValue)
		return (JJSkeleton::COMP_NEW);
	
	else if (y == 0 && x == 0 && pMap[nIdxZbefore] != BKValue) 
		return (JJSkeleton::COMP_Z);
	
	else if (y == 0 && x == 0 && pMap[nIdxZbefore] == BKValue) 
		return (JJSkeleton::COMP_NEW);

	else if (z == 0 && pMap[nIdxXbefore] != BKValue)
		return (JJSkeleton::COMP_X);

	else if (z == 0 && pMap[nIdxXbefore] == BKValue)
		return (JJSkeleton::COMP_NEW);

	else if (z == 0 && pMap[nIdxYbefore] != BKValue)
		return (JJSkeleton::COMP_Y);
	
	else if (z == 0 && pMap[nIdxYbefore] == BKValue)
		return (JJSkeleton::COMP_NEW);

	else if (pMap[nIdxXbefore] != BKValue && pMap[nIdxYbefore] != BKValue && pMap[nIdxZbefore] != BKValue) 
	{
		if (pMap[nIdxXbefore] == pMap[nIdxYbefore] && pMap[nIdxXbefore] == pMap[nIdxZbefore])
			return (JJSkeleton::COMP_SAME);
	
		else if (pMap[nIdxXbefore] != pMap[nIdxYbefore] && pMap[nIdxXbefore] != pMap[nIdxZbefore])
			return (JJSkeleton::COMP_DIFF_XYZ);
		
		else if (pMap[nIdxXbefore] != pMap[nIdxYbefore])
			return (JJSkeleton::COMP_DIFF_XY);
		
		else if (pMap[nIdxYbefore] != pMap[nIdxZbefore])
			return (JJSkeleton::COMP_DIFF_YZ);

		else if (pMap[nIdxXbefore] != pMap[nIdxZbefore])
			return (JJSkeleton::COMP_DIFF_ZX);
	}
	
	else if (pMap[nIdxXbefore] != BKValue)
		return (JJSkeleton::COMP_X);

	else if (pMap[nIdxYbefore] != BKValue)
		return (JJSkeleton::COMP_Y);
	
	else if (pMap[nIdxZbefore] != BKValue)
		return (JJSkeleton::COMP_Z);

	else return (JJSkeleton::COMP_NEW);

	return 0;
}

void JJSkeleton::SetLabelByState(unsigned short *pMap, int nWidth, int nHeight, int nDepth, int x, int y, int z, int *nLabel, int nState,unsigned int *pTable)
{
	int nIdx = x+y*nWidth+z*nWidth*nHeight;
	int nMinTemp, nMaxTemp;
	int Min_TempFinalLabel, Max_TempFinalLabel;

	switch(nState){
		case COMP_NONE	:
		case COMP_NEW	:
			pMap[nIdx] = (*nLabel)++;
			break;			
		case COMP_X		:
			pMap[nIdx] = pMap[nIdx-1];
			break;
		case COMP_Y		:
			pMap[nIdx] = pMap[nIdx-nWidth];
			break;
		case COMP_Z     :
		case COMP_SAME	:
			pMap[nIdx] = pMap[nIdx-nWidth*nHeight];
			break;
		case COMP_DIFF_XY	:
			nMinTemp = min(pMap[nIdx-1], pMap[nIdx-nWidth]);
			nMaxTemp = max(pMap[nIdx-1], pMap[nIdx-nWidth]);
			Max_TempFinalLabel = FindFinalLabel(pTable, nMaxTemp);
			Min_TempFinalLabel = FindFinalLabel(pTable, nMinTemp);
			pTable[Max_TempFinalLabel] = pMap[nIdx] = Min_TempFinalLabel;

			break;
		case COMP_DIFF_YZ	:
			nMinTemp = min(pMap[nIdx-nWidth], pMap[nIdx-nWidth*nHeight]);
			nMaxTemp = max(pMap[nIdx-nWidth], pMap[nIdx-nWidth*nHeight]);
			Max_TempFinalLabel = FindFinalLabel(pTable, nMaxTemp);
			Min_TempFinalLabel = FindFinalLabel(pTable, nMinTemp);
			pTable[Max_TempFinalLabel] = pMap[nIdx] = Min_TempFinalLabel;

			break;
		case COMP_DIFF_ZX	:
			nMinTemp = min(pMap[nIdx-nWidth*nHeight], pMap[nIdx-1]);
			nMaxTemp = max(pMap[nIdx-nWidth*nHeight], pMap[nIdx-1]);
			Max_TempFinalLabel = FindFinalLabel(pTable, nMaxTemp);
			Min_TempFinalLabel = FindFinalLabel(pTable, nMinTemp);
			pTable[Max_TempFinalLabel] = pMap[nIdx] = Min_TempFinalLabel;

			break;
		case COMP_DIFF_XYZ	:
			nMinTemp = min(min(pMap[nIdx-1], pMap[nIdx-nWidth]), pMap[nIdx-nWidth*nHeight]);
			nMaxTemp = max(max(pMap[nIdx-1], pMap[nIdx-nWidth]), pMap[nIdx-nWidth*nHeight]);
			Max_TempFinalLabel = FindFinalLabel(pTable, nMaxTemp);
			Min_TempFinalLabel = FindFinalLabel(pTable, nMinTemp);
			pTable[Max_TempFinalLabel] = pMap[nIdx] = Min_TempFinalLabel;

			break;
	}
}

int JJSkeleton::MapFillingWithSortedLabel(unsigned int *pTable, unsigned short *pMap, int iVolumeSize)
{
	unsigned short *pSortTable = (unsigned short *)GlobalAlloc(GMEM_FIXED, iVolumeSize*sizeof(unsigned short));
	ZeroMemory(pSortTable , iVolumeSize*sizeof(unsigned short));

	int nSortedCnt = 1, old = 0;
	for(int i=0; i<iVolumeSize; i++) {
		int nUnsortedLabel = pTable[pMap[i]]; 
		pMap[i] = nUnsortedLabel;

		if( nUnsortedLabel!=0 && pSortTable[nUnsortedLabel] == 0) {
			//int diff = abs(old-nUnsortedLabel);
			//if(diff > 100) 
				pSortTable[nUnsortedLabel] = nSortedCnt++;	
			//old = nUnsortedLabel;
		}
	}

	for(i=0; i<iVolumeSize; i++) {
		pMap[i] = pSortTable[pMap[i]];
	}
	TRACE(_T("nSortedCnt=%d\n"), nSortedCnt);

	GlobalFree((unsigned short *)pSortTable);

	return(nSortedCnt);
}

int JJSkeleton::RemoveSmallComponent(unsigned char *pMap, int nWidth, int nHeight, int nDepth)
{
	unsigned int iVolumeSize = nWidth*nHeight*nDepth;
	unsigned short *pTempMap = (unsigned short *)GlobalAlloc(GMEM_FIXED, iVolumeSize*sizeof(unsigned short));
	unsigned int *pTable   = (unsigned int *)GlobalAlloc(GMEM_FIXED, iVolumeSize*sizeof(unsigned int));
//	ZeroMemory(pTempMap , iVolumeSize*sizeof(unsigned short));
//	ZeroMemory(pTable   , iVolumeSize*sizeof(unsigned int));
	memset(pTempMap, 0, iVolumeSize*sizeof(unsigned short));
	memset(pTable, 0, iVolumeSize*sizeof(unsigned int));
	
	for(int i=0;i<iVolumeSize;i++) // i > 512*512老 锭 巩力 积辨 荐 乐澜. 
	{
		pTempMap[i] = pMap[i];
		pTable[i] = i; 
	}

	int nLabelCnt=1;
	for(int z=0; z<nDepth; z++) {
		int nSliceOffset = z*nWidth*nHeight;
		for(int y=0; y<nHeight; y++) {			
			int nRowOffset = y*nWidth;
			for(int x=0; x<nWidth; x++) {
				int nIdx = x + nRowOffset + nSliceOffset;
				if(pMap[nIdx]==BKValue )	continue;
							
				int nState = GetState(pTempMap, nWidth, nHeight, nDepth, x, y, z);

				SetLabelByState(pTempMap, nWidth, nHeight, nDepth, x, y, z, 
						&nLabelCnt, nState, pTable);
			}
		}
	}

	for(i=0;i<iVolumeSize;i++) 
		pTable[i] = FindFinalLabel(pTable, i);
	
	int nMaxLabel = MapFillingWithSortedLabel(pTable, pTempMap, iVolumeSize);
	
	///////////////////////////////////////////////////////////////////////////////////////////////
	memset(pTable, 0, iVolumeSize*sizeof(unsigned int));
	for( i=0; i < iVolumeSize; i++){
		if(pTempMap[i]!=0 ) 
			pTable[pTempMap[i]]++;
	}		

	int *pBigCompIndex = new int[nMaxLabel];
	int count = 0, j;

	// small component甫 吧矾辰促.
	// pBigCompIndex绰 器窃等 汗伎捞 8000俺 捞惑牢 巴甸狼 牢郸胶
	for(i=0; i<nMaxLabel; i++) {
		//if( pTable[i] > 700 && pTable[i] < 1000 ) {
		if( pTable[i] > 10000 ) { // test 
			pBigCompIndex[count++] = i;
		}
	}
	for( i=0;i<iVolumeSize;i++) {
		for( j=0; j<count; j++) {
			pMap[i] = 0x00;
			if( pTempMap[i] == pBigCompIndex[j] ) {
				pMap[i] = 0xff;
				j = count;
			}
		}
	} 
	delete []pBigCompIndex;

	GlobalFree((unsigned short *)pTempMap);
	GlobalFree((unsigned int *)pTable);

	return(nMaxLabel);
}

void JJSkeleton::SetParameter(RxMaskVolume* pVesselEdge, RxMaskVolume* pVesselVolume, RxMaskVolume* pSkeleton, int minX, int maxX, int minY, int maxY, int minZ, int maxZ, BOOL bRef)
{
	m_pVesselEdge = pVesselEdge;
	m_pVesselVolume = pVesselVolume;
	m_pSkeleton = pSkeleton;

	m_VolX = maxX - minX + 1 + 4;
	m_VolY = maxY - minY + 1 + 4;
	m_VolZ = maxZ - minZ + 1 + 4;
	m_OffsetX = minX - 2;
	m_OffsetY = minY - 2;
	m_OffsetZ = minZ - 2;

	if (m_OffsetZ <= 0)
	{
		m_VolZ = maxZ - minZ + 1 + 2 + 1 + m_OffsetZ;
		m_OffsetZ = 1;
		if (m_VolZ >= m_pVesselVolume->m_iVolZ)
			m_VolZ = m_pVesselVolume->m_iVolZ - 1;
	}

	// memory allocation
	if (m_pLabel != NULL)
	{
		free(m_pLabel);
		m_pLabel = NULL;
	}
	if (pTable != NULL)
	{
		free(pTable);
		pTable = NULL;
	}
	if (m_pLongDistanceMap != NULL)
	{
		free(m_pLongDistanceMap);
		m_pLongDistanceMap = NULL;
	}
	if (m_pPointType != NULL)
	{
		free(m_pPointType);
		m_pPointType = NULL;
	}
	if (m_pDistanceMapFromPS != NULL)
	{
		free(m_pDistanceMapFromPS);
		m_pDistanceMapFromPS = NULL;
	}
	if (m_pClusterGraph != NULL)
	{
		free(m_pClusterGraph);
		m_pClusterGraph = NULL;
	}
	if (m_pDistanceMapFromPS_FAST_WAVE != NULL)
	{
		free(m_pDistanceMapFromPS_FAST_WAVE);
		m_pDistanceMapFromPS_FAST_WAVE = NULL;
	}

	m_pLabel = (unsigned short *)malloc(m_VolX*m_VolY*m_VolZ*sizeof(unsigned short));	
	pTable  = (unsigned int *)malloc(m_VolX*m_VolY*m_VolZ*sizeof(unsigned int));
	m_pLongDistanceMap = (float *)malloc(m_VolX*m_VolY*m_VolZ*sizeof(float));
	m_pPointType = (unsigned char*)malloc(m_VolX*m_VolY*m_VolZ*sizeof(unsigned char));
	m_pDistanceMapFromPS = (float *)malloc(m_VolX*m_VolY*m_VolZ*sizeof(float));
	m_pClusterGraph = (unsigned short*)malloc(m_VolX*m_VolY*m_VolZ*sizeof(unsigned short));
	m_pDistanceMapFromPS_FAST_WAVE = (float *)malloc(m_VolX*m_VolY*m_VolZ*sizeof(float));

	m_nObjectNum = 0;
	m_nExtremeNum = 0;
	m_nBranchNum = 0;

	int i, j, k;

	for (k = 0; k < m_pVesselEdge->m_iVolZ; k++)
		for (j = 0; j < m_pVesselEdge->m_iVolY; j++)
			for (i = 0; i < m_pVesselEdge->m_iVolX; i++)
				m_pSkeleton->ClearMask(i, j, k);
}