jjskeleton.cpp

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

				r = sqrt((2.0f/m_fVel/m_fVel - (T1 - T3)*(T1 - T3)));
				s = (T1 + T3 - r)/2.0f;

				if (s >= T1 && s >= T3) 
					sol = minimum(sol, s);
				else
				{
					s += r;
					if (s >= T1 && s >= T3) 
						sol = minimum(sol, s);
				}
			} else // (T1, X, X)
			{
				sol = minimum(sol, 1.0f/m_fVel + T1);
			}
		}
	} else // (X, ?, ?)
	{
		if (m_pPointType[lIndex2] == KNOWN) // (X, T2, ?)
			if (m_pPointType[lIndex3] == KNOWN) // (X, T2, T3)
			{
				r = sqrt((2.0f/m_fVel/m_fVel - (T2 - T3)*(T2 - T3)));
				s = (T2 + T3 - r)/2.0f;

				if (s >= T2 && s >= T3) 
					sol = minimum(sol, s);
				else
				{
					s += r;
					if (s >= T2 && s >= T3) 
						sol = minimum(sol, s);
				}
			} else // (X, T2, X)
				sol = minimum(sol, 1.0f/m_fVel + T2);
		else if (m_pPointType[lIndex3] == KNOWN) // (X, X, T3)
			sol = minimum(sol, 1.0f/m_fVel + T3);
	}

	return sol;
}

inline float JJSkeleton::minimum(float num1, float num2) 
{	
	if (num1 < num2) 
		return num1; 
	else 
		return num2;
}

void JJSkeleton::Skeletonization(void)
{
	DWORD CurrentTime, ElapsedTime;	
	double fTotalTime = 0.0;

	FindConnectedComponent();

	Initialization();

	Iteration_ORD_1(LEVEL_SET_F_1);
	
	for (int c = 0; c < m_nObjectNum; c++)
	{
		FindPointSource(c);
		PropagateFromPS(c);
		BuildClusterGraph();
		PropagateFromPS_FAST_WAVE(c);
		ExtractCenterline();
	}

	VisualizeSkeleton();
}

void JJSkeleton::FindPointSource(short nObjectLabel)
{
	unsigned long lIndex;
	int i, j, k;
	m_fMaxDistanceMapFromPS = 0.0f;

	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_pLabel[lIndex] == m_nObjectID[nObjectLabel] && m_pVesselVolume->GetMask(i + m_OffsetX, j + m_OffsetY, k + m_OffsetZ))
				{
					if (m_pLongDistanceMap[lIndex] < MAX_VALUE && m_pLongDistanceMap[lIndex] > m_fMaxDistanceMapFromPS)
					{
						m_fMaxDistanceMapFromPS = m_pLongDistanceMap[lIndex];

						m_PS_X = i;
						m_PS_Y = j;
						m_PS_Z = k;
					}
				}
			}
}

void JJSkeleton::PropagateFromPS(short nObjectLabel)
{
	int i, j, k;
	unsigned long lIndex;

	HeapInitialization();

	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_PS_X == i && m_PS_Y == j && m_PS_Z == k)
				{
					m_pPointType[lIndex] = BAND;
					m_pDistanceMapFromPS[lIndex] = 0.0f;

					Node NewNode;
					NewNode.time = 0.0f;
					NewNode.x = i;
					NewNode.y = j;
					NewNode.z = k;

					HeapInsert(NewNode);
				// include vessel edge or not when making cluster graph
				} else if (m_pLabel[lIndex] == m_nObjectID[nObjectLabel] 
					&& m_pVesselVolume->GetMask(i + m_OffsetX, j + m_OffsetY, k + m_OffsetZ) 
					&& !m_pVesselEdge->GetMask(i + m_OffsetX, j + m_OffsetY, k + m_OffsetZ))
//				} else if (m_pVesselVolume->GetMask(i + m_OffsetX, j + m_OffsetY, k + m_OffsetZ))
				{
					m_pPointType[lIndex] = INSIDE;
					m_pDistanceMapFromPS[lIndex] = MAX_VALUE;					
				} else 
//					if (m_pVesselEdge->GetMask(i + m_OffsetX, j + m_OffsetY, k + m_OffsetZ) 
//					|| !m_pVesselVolume->GetMask(i + m_OffsetX, j + m_OffsetY, k + m_OffsetZ))
				{
					m_pPointType[lIndex] = BACKGROUND;
					m_pDistanceMapFromPS[lIndex] = 0.0f;										
				}
			}
		
	Iteration_ORD_1(LEVEL_SET_F_D);

	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_pLabel[lIndex] == m_nObjectID[nObjectLabel] && m_pVesselVolume->GetMask(i + m_OffsetX, j + m_OffsetY, k + m_OffsetZ) && !m_pVesselEdge->GetMask(i + m_OffsetX, j + m_OffsetY, k + m_OffsetZ))
				{
					if (m_pDistanceMapFromPS[lIndex] != MAX_VALUE)
					{
						m_pDistanceMapFromPS[lIndex] *= 0.5;
					}
				}
			}
}

void JJSkeleton::PropagateFromPS_FAST_WAVE(short nObjectLabel)
{
	int i, j, k;
	unsigned long lIndex;

	HeapInitialization();

	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_PS_X == i && m_PS_Y == j && m_PS_Z == k)
				{
					m_pPointType[lIndex] = BAND;
					m_pDistanceMapFromPS_FAST_WAVE[lIndex] = 0.0f;

					Node NewNode;
					NewNode.time = 0.0f;
					NewNode.x = i;
					NewNode.y = j;
					NewNode.z = k;

					HeapInsert(NewNode);
//				} else if (m_pVesselVolume->GetMask(i + m_OffsetX, j + m_OffsetY, k + m_OffsetZ) && !m_pVesselEdge->GetMask(i + m_OffsetX, j + m_OffsetY, k + m_OffsetZ))
				} else if (m_pLabel[lIndex] == m_nObjectID[nObjectLabel] && m_pVesselVolume->GetMask(i + m_OffsetX, j + m_OffsetY, k + m_OffsetZ))
				{
					m_pPointType[lIndex] = INSIDE;
					m_pDistanceMapFromPS_FAST_WAVE[lIndex] = MAX_VALUE;					
				} else
				{
					m_pPointType[lIndex] = BACKGROUND;
					m_pDistanceMapFromPS_FAST_WAVE[lIndex] = 0.0f;										
				}
			}
		
	Iteration_ORD_1(LEVEL_SET_F_D_FAST);
}

void JJSkeleton::BuildClusterGraph(void)
{	
	Node Head;
	unsigned long lIndex;
	unsigned short k, l, m;

	// 6 neighbor
	short nr_x[6] = {-1, 0, 0, 1, 0, 0};
	short nr_y[6] = {0, -1, 0, 0, 1, 0};
	short nr_z[6] = {0, 0, -1, 0, 0, 1};	

	memset(m_pClusterGraph, 0x00, m_VolX*m_VolY*m_VolZ*sizeof(unsigned short));
	m_nClusterIndex = 0;
	HeapInitialization();

	for (int i = 0; i < 65535; i++)
	{
		m_GL_Index[i] = 0;
		for (int j = 0; j < 20; j++)
			m_GraphLink[i][j] = 0;
		
		m_ExtremePoint[i].time = 0.0;
	}

	for (i = 0; i < 65535; i++)
		m_nClusterMapping[i] = 0;

	Node RootNode;
	RootNode.time = m_pDistanceMapFromPS[m_PS_X + m_PS_Y*m_VolX + m_PS_Z*m_VolX*m_VolY];
	RootNode.x = m_PS_X;
	RootNode.y = m_PS_Y;
	RootNode.z = m_PS_Z;
	RootNode.pointer = 1;

	HeapInsert(RootNode);

	while (!HeapEmpty())
	{
		Head = HeapRemoveMin();
		lIndex = Head.x + Head.y*m_VolX + Head.z*m_VolX*m_VolY;

		if (m_pClusterGraph[lIndex] == 0)
		{
			m_nClusterIndex++;

			if (m_nClusterMapping[Head.pointer] < (short)Head.time)
				BuildCGLink(Head.pointer, m_nClusterIndex);
			else
				BuildCGLink(m_nClusterIndex, Head.pointer);

			m_pClusterGraph[lIndex] = m_nClusterIndex;

			m_nClusterMapping[m_nClusterIndex] = (short)Head.time;
			m_ExtremePoint[m_nClusterIndex] = Head;

			SetNode(Head.x, Head.y, Head.z, (short)Head.time, m_nClusterIndex);
		} else
		{
			if (m_nClusterMapping[Head.pointer] < (short)Head.time)
				BuildCGLink(Head.pointer, m_pClusterGraph[lIndex]);
			else
				BuildCGLink(m_pClusterGraph[lIndex], Head.pointer);			
		}
	}

	for (i = 1; i <= m_nClusterIndex; i++)
	{
		if (m_GL_Index[i] == 0)
		{
			m_ControlPoint[m_nExtremeNum].x = m_ExtremePoint[i].x;
			m_ControlPoint[m_nExtremeNum].y = m_ExtremePoint[i].y;
			m_ControlPoint[m_nExtremeNum].z = m_ExtremePoint[i].z;
			m_nExtremeNum++;
		}
	}
}

void JJSkeleton::SetNode(short nSeed_X, short nSeed_Y, short nSeed_Z, short nDistance, short nClusterIndex)
{
	Node NewNode;
	unsigned short k, l, m;
	unsigned long lIndex;

	// 6 neighbor
/*	short nr_x[6] = {-1, 0, 0, 1, 0, 0};
	short nr_y[6] = {0, -1, 0, 0, 1, 0};
	short nr_z[6] = {0, 0, -1, 0, 0, 1};
	BOOL bGoDeeper[6] = {FALSE, FALSE, FALSE, FALSE, FALSE, FALSE};*/

	// 26 neighbor
	short nr_x[26];
	short nr_y[26];
	short nr_z[26];
	BOOL bGoDeeper[26];
	short count = 0;

	for (int z = -1; z <= 1; z++)
		for (int y = -1; y <= 1; y++)
			for (int x = -1; x <= 1; x++)
			{
				if (!(x == 0 && y == 0 && z == 0))
				{
					bGoDeeper[count] = FALSE;
					nr_x[count] = x;
					nr_y[count] = y;
					nr_z[count] = z;
					count++;
				}
			}

	for (int nr = 0; nr < 26; nr++)
	{
		k = nSeed_X + nr_x[nr];
		l = nSeed_Y + nr_y[nr];
		m = nSeed_Z + nr_z[nr];
		lIndex = k + l*m_VolX + m*m_VolX*m_VolY;

		if (m_pPointType[lIndex] == KNOWN && m_pClusterGraph[lIndex] == 0)
		{
			if ((short)m_pDistanceMapFromPS[lIndex] == nDistance)
			{
				m_pClusterGraph[lIndex] = nClusterIndex;
				bGoDeeper[nr] = TRUE;
				
				if (m_ExtremePoint[nClusterIndex].time > m_pDistanceMapFromPS[lIndex])
				{
					m_ExtremePoint[nClusterIndex].time = m_pDistanceMapFromPS[lIndex];
					m_ExtremePoint[nClusterIndex].x = k;
					m_ExtremePoint[nClusterIndex].y = l;
					m_ExtremePoint[nClusterIndex].z = m;
				}
			} else
			{
				NewNode.time = m_pDistanceMapFromPS[lIndex];
				NewNode.x = k;
				NewNode.y = l;
				NewNode.z = m;
				NewNode.pointer = nClusterIndex;

				if (!HeapSearch(NewNode))
					HeapInsert(NewNode);
			}
		}
	}

	for (nr = 0; nr < 26; nr++)
	{
		if (bGoDeeper[nr])
		{
			k = nSeed_X + nr_x[nr];
			l = nSeed_Y + nr_y[nr];
			m = nSeed_Z + nr_z[nr];

			SetNode(k, l, m, nDistance, nClusterIndex);
		}
	}
}

void JJSkeleton::BuildCGLink(unsigned short nParent, unsigned short nChild)
{
	for (int i = 0; i < m_GL_Index[nParent]; i++)
		if (m_GraphLink[nParent][i] == nChild)
			return;

	m_GraphLink[nParent][m_GL_Index[nParent]] = nChild;
	m_GL_Index[nParent]++;

	return;
}

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

	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;

//				m_pVesselEdge->ClearMask(i + m_OffsetX, j + m_OffsetY, k + m_OffsetZ);
//				m_pVesselVolume->ClearMask(i + m_OffsetX, j + m_OffsetY, k + m_OffsetZ);
/*				if (m_pDistanceMapFromPS_FAST_WAVE[lIndex] != MAX_VALUE && m_pDistanceMapFromPS_FAST_WAVE[lIndex] > 0.0)
					m_pVesselVolume->SetMask(i, j, k + m_OffsetZ);
				else
					m_pVesselVolume->ClearMask(i, j, k + m_OffsetZ);*/
			}

/*	for (c = 0; c < m_nObjectNum; c++)
	{
		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_pLabel[lIndex] == m_nObjectID[c])
						m_pVesselVolume->SetMask(i, j, k + m_OffsetZ);
	/*				if (m_pDistanceMapFromPS_FAST_WAVE[lIndex] != MAX_VALUE && m_pDistanceMapFromPS_FAST_WAVE[lIndex] > 0.0)
						m_pVesselVolume->SetMask(i, j, k + m_OffsetZ);
					else
						m_pVesselVolume->ClearMask(i, j, k + m_OffsetZ);*/
//				}
//	}
}

void JJSkeleton::ExtractCenterline(void)
{
	double k1_x, k1_y, k1_z, h = 1;
	double f_x, f_y, f_z, mag;
	short x, y, z, x2, y2, z2;
	short next_x, next_y, next_z;
	int i, j, k;

	// 26 neighbor
	short nr_x[26];
	short nr_y[26];
	short nr_z[26];
	short count = 0;
	double min;

	for (k = -1; k <= 1; k++)
		for (j = -1; j <= 1; j++)
			for (i = -1; i <= 1; i++)
			{
				if (!(i == 0 && j == 0 && k == 0))
				{
					nr_x[count] = i;
					nr_y[count] = j;
					nr_z[count] = k;
					count++;
				}
			}

	BOOL bCont;
	
	for (i = 1; i <= m_nClusterIndex; i++)
	{
		if (m_GL_Index[i] == 0)
		{
			x = m_ExtremePoint[i].x;
			y = m_ExtremePoint[i].y;
			z = m_ExtremePoint[i].z;
			
			do
			{
				bCont = TRUE;
				if (x == m_PS_X && y == m_PS_Y && z == m_PS_Z)
					bCont = FALSE;
				else if (m_pSkeleton->GetMask((short)x + m_OffsetX, (short)y + m_OffsetY, (short)z + m_OffsetZ))
				{
					BOOL bFind = FALSE;
					for (int tr = 0; tr < m_nBranchNum; tr++)
					{
						if ((short)x == m_BranchPoint[tr].x && (short)y == m_BranchPoint[tr].y	&& (short)z == m_BranchPoint[tr].z)
							bFind = TRUE;