jjskeleton.cpp

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

#include "stdafx.h"
#include "JJSkeleton.h"

#include <math.h>

#define BAND		0
#define INSIDE		1
#define KNOWN		2
#define BACKGROUND	3
#define MAX_VALUE	10000

#define TXT_WRITE			0
#define RAW_WRITE			1
#define LUNG_EDGE_WRITE		2
#define LUNG_VOLUME_WRITE	3

#define LEVEL_SET_F_1		0
#define LEVEL_SET_F_D		1
#define LEVEL_SET_F_D_FAST	2
#define CLUSTER_GRAPH		3

#define BKValue 0	// background value

JJSkeleton::JJSkeleton()
{
	m_pLongDistanceMap = NULL;
	m_pPointType = NULL;
	m_pDistanceMapFromPS = NULL;
	m_pDistanceMapFromPS_FAST_WAVE = NULL;
	m_pClusterGraph = NULL;
	m_pLabel = NULL;

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

JJSkeleton::~JJSkeleton()
{
	if (m_pLongDistanceMap != NULL)
		free(m_pLongDistanceMap);

	if (m_pPointType != NULL)
		free(m_pPointType);

	if (m_pDistanceMapFromPS != NULL)
		free(m_pDistanceMapFromPS);

	if (m_pDistanceMapFromPS_FAST_WAVE != NULL)
		free(m_pDistanceMapFromPS_FAST_WAVE);
	
	if (m_pClusterGraph != NULL)
		free(m_pClusterGraph);

	if (m_pLabel != NULL)
		free(m_pLabel);
}

void JJSkeleton::Initialization(void)
{
	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_pVesselEdge->GetMask(i + m_OffsetX, j + m_OffsetY, k + m_OffsetZ))
				{
					m_pPointType[lIndex] = BAND;
					m_pLongDistanceMap[lIndex] = 1.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_pPointType[lIndex] = INSIDE;
					m_pLongDistanceMap[lIndex] = MAX_VALUE;					
				} else
				{
					m_pPointType[lIndex] = BACKGROUND;
					m_pLongDistanceMap[lIndex] = 0.0f;										
				}
			}
}

Node JJSkeleton::HeapRemoveMin(void)
{
	if (m_lHeapLength == 0)
	{
		Node Null;
		Null.time = -1;		

		return Null;
	}
	
	Node Min = m_arrHeap[1];
	Node Last = m_arrHeap[m_lHeapLength];
	m_lHeapLength = m_lHeapLength - 1;

	unsigned long lIndex = 1, lChildIndex = 2;

	while (lChildIndex <= m_lHeapLength) 
	{
		if ((lChildIndex < m_lHeapLength) && 
			(m_arrHeap[lChildIndex].time >= m_arrHeap[lChildIndex + 1].time))
			lChildIndex = lChildIndex + 1;

		if (Last.time <= m_arrHeap[lChildIndex].time)
			break; // check whether the desired operation is performed or not
		else
		{
			m_arrHeap[lIndex] = m_arrHeap[lChildIndex];
			lIndex = lChildIndex;
			lChildIndex = 2*lChildIndex;
		}
	}

	m_arrHeap[lIndex] = Last;

	return Min;
}

void JJSkeleton::HeapInsert(Node NewNode)
{
	m_lHeapLength++;

	if (m_lHeapLength > 1000000)
		AfxMessageBox(_T("Heap Overflow"));

	unsigned long lIndex = m_lHeapLength;

	while (lIndex > 1 && m_arrHeap[lIndex/2].time > NewNode.time)
	{
		m_arrHeap[lIndex] = m_arrHeap[lIndex/2];
		lIndex = lIndex/2;
	}

	m_arrHeap[lIndex] = NewNode;
}

BOOL JJSkeleton::HeapSearch(Node SearchNode)
{
	unsigned long lIndex = 1;

	for (int i = 1; i <= m_lHeapLength; i++)
	{
		if (m_arrHeap[i].x == SearchNode.x && m_arrHeap[i].y == SearchNode.y && m_arrHeap[i].z == SearchNode.z)
		{
			return TRUE;
		}
	}

	return FALSE;
}

void JJSkeleton::HeapModify(Node NewNode)
{
	unsigned long lIndex = 1;

	for (int i = 1; i <= m_lHeapLength; i++)
	{
		if (m_arrHeap[i].x == NewNode.x && m_arrHeap[i].y == NewNode.y && m_arrHeap[i].z == NewNode.z)
		{
			lIndex = i;
			break;
		}
	}

	if (m_arrHeap[lIndex].time > NewNode.time)
	{
		while (lIndex > 1 && m_arrHeap[lIndex/2].time > NewNode.time)
		{
			m_arrHeap[lIndex] = m_arrHeap[lIndex/2];
			lIndex = lIndex/2;
		}
	} else
	{
		while ((lIndex*2 + 1 <= m_lHeapLength) && 
			((m_arrHeap[lIndex*2].time < NewNode.time) || (m_arrHeap[lIndex*2 + 1].time < NewNode.time)))
		{
			if (m_arrHeap[lIndex*2].time < m_arrHeap[lIndex*2 + 1].time)
			{
				m_arrHeap[lIndex] = m_arrHeap[lIndex*2];
				lIndex = lIndex*2;
			} else
			{
				m_arrHeap[lIndex] = m_arrHeap[lIndex*2 + 1];
				lIndex = lIndex*2 + 1;
			}
		}
	}

	m_arrHeap[lIndex] = NewNode;
}

void JJSkeleton::HeapInitialization(void)
{
	m_lHeapLength = 0;
}

BOOL JJSkeleton::HeapEmpty(void)
{
	if (m_lHeapLength == 0)
		return TRUE;
	else
		return FALSE;
}

void JJSkeleton::Iteration_ORD_1(short nDistanceMapMode)
{
	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};

	// 26 neighbor
/*	short nr_x[26];
	short nr_y[26];
	short nr_z[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))
				{
					nr_x[count] = x;
					nr_y[count] = y;
					nr_z[count] = z;
					count++;
				}
			}*/

	float sol;
	BOOL bInHeap;

	while (!HeapEmpty())
	{
		Head = HeapRemoveMin(); /* STEP 1 */
		
		lIndex = Head.x + Head.y*m_VolX + Head.z*m_VolX*m_VolY;
		m_pPointType[lIndex] = KNOWN;

		if (Head.x > 1 && Head.y > 1 && Head.z > 1 && Head.z < m_VolZ - 1)
		{

		for (int nr = 0; nr < 6; nr++)
		{
			k = Head.x + nr_x[nr];
			l = Head.y + nr_y[nr];
			m = Head.z + nr_z[nr];
			lIndex = k + l*m_VolX + m*m_VolX*m_VolY;

			if (m_pPointType[lIndex] != KNOWN && m_pPointType[lIndex] != BACKGROUND)
			{
				if (m_pPointType[lIndex] == BAND) 
					bInHeap = TRUE;
				else if (m_pPointType[lIndex] == INSIDE) 
				{
					bInHeap = FALSE;
					m_pPointType[lIndex] = BAND; /* STEP 2 */
				}
				
				sol = MAX_VALUE; /* STEP 3 */
				
				switch (nDistanceMapMode)
				{
				case LEVEL_SET_F_1:
					m_fVel = 1.0f;

					sol = solve_ORD_1(k - 1, l, m, k, l - 1, m, k, l, m - 1, sol, LEVEL_SET_F_1);
					sol = solve_ORD_1(k + 1, l, m, k, l - 1, m, k, l, m - 1, sol, LEVEL_SET_F_1);
					sol = solve_ORD_1(k - 1, l, m, k, l + 1, m, k, l, m - 1, sol, LEVEL_SET_F_1);
					sol = solve_ORD_1(k + 1, l, m, k, l + 1, m, k, l, m - 1, sol, LEVEL_SET_F_1);
					
					sol = solve_ORD_1(k - 1, l, m, k, l - 1, m, k, l, m + 1, sol, LEVEL_SET_F_1);
					sol = solve_ORD_1(k + 1, l, m, k, l - 1, m, k, l, m + 1, sol, LEVEL_SET_F_1);
					sol = solve_ORD_1(k - 1, l, m, k, l + 1, m, k, l, m + 1, sol, LEVEL_SET_F_1);
					sol = solve_ORD_1(k + 1, l, m, k, l + 1, m, k, l, m + 1, sol, LEVEL_SET_F_1);
					break;
				case LEVEL_SET_F_D:
					m_fVel = pow(m_pLongDistanceMap[lIndex], 0.5f);

					sol = solve_ORD_1(k - 1, l, m, k, l - 1, m, k, l, m - 1, sol, LEVEL_SET_F_D);
					sol = solve_ORD_1(k + 1, l, m, k, l - 1, m, k, l, m - 1, sol, LEVEL_SET_F_D);
					sol = solve_ORD_1(k - 1, l, m, k, l + 1, m, k, l, m - 1, sol, LEVEL_SET_F_D);
					sol = solve_ORD_1(k + 1, l, m, k, l + 1, m, k, l, m - 1, sol, LEVEL_SET_F_D);
					
					sol = solve_ORD_1(k - 1, l, m, k, l - 1, m, k, l, m + 1, sol, LEVEL_SET_F_D);
					sol = solve_ORD_1(k + 1, l, m, k, l - 1, m, k, l, m + 1, sol, LEVEL_SET_F_D);
					sol = solve_ORD_1(k - 1, l, m, k, l + 1, m, k, l, m + 1, sol, LEVEL_SET_F_D);
					sol = solve_ORD_1(k + 1, l, m, k, l + 1, m, k, l, m + 1, sol, LEVEL_SET_F_D);
					break;
				case LEVEL_SET_F_D_FAST:
					m_fVel = exp(m_pLongDistanceMap[lIndex]/m_fMaxDistanceMapFromPS*25.0f);

					sol = solve_ORD_1(k - 1, l, m, k, l - 1, m, k, l, m - 1, sol, LEVEL_SET_F_D_FAST);
					sol = solve_ORD_1(k + 1, l, m, k, l - 1, m, k, l, m - 1, sol, LEVEL_SET_F_D_FAST);
					sol = solve_ORD_1(k - 1, l, m, k, l + 1, m, k, l, m - 1, sol, LEVEL_SET_F_D_FAST);
					sol = solve_ORD_1(k + 1, l, m, k, l + 1, m, k, l, m - 1, sol, LEVEL_SET_F_D_FAST);
					
					sol = solve_ORD_1(k - 1, l, m, k, l - 1, m, k, l, m + 1, sol, LEVEL_SET_F_D_FAST);
					sol = solve_ORD_1(k + 1, l, m, k, l - 1, m, k, l, m + 1, sol, LEVEL_SET_F_D_FAST);
					sol = solve_ORD_1(k - 1, l, m, k, l + 1, m, k, l, m + 1, sol, LEVEL_SET_F_D_FAST);
					sol = solve_ORD_1(k + 1, l, m, k, l + 1, m, k, l, m + 1, sol, LEVEL_SET_F_D_FAST);
					break;
				}

				if (bInHeap)
				{
				//	if (m_pLongDistanceMap[lIndex] > sol)
					Node New;

					switch (nDistanceMapMode)
					{
					case LEVEL_SET_F_1:
						if (m_pLongDistanceMap[lIndex] > 1.0f)
						{
							m_pLongDistanceMap[lIndex] = sol;
		
							New.time = sol;
							New.x = k;
							New.y = l;
							New.z = m;
							HeapModify(New); /* STEP 4 */
						}
						break;
					case LEVEL_SET_F_D:
						m_pDistanceMapFromPS[lIndex] = sol;

						New.time = sol;
						New.x = k;
						New.y = l;
						New.z = m;
						HeapModify(New); /* STEP 4 */
						break;
					case LEVEL_SET_F_D_FAST:
						m_pDistanceMapFromPS_FAST_WAVE[lIndex] = sol;

						New.time = sol;
						New.x = k;
						New.y = l;
						New.z = m;
						HeapModify(New); /* STEP 4 */
						break;
					}
				} else
				{
					switch (nDistanceMapMode)
					{
					case LEVEL_SET_F_1:
						m_pLongDistanceMap[lIndex] = sol;
						break;
					case LEVEL_SET_F_D:
						m_pDistanceMapFromPS[lIndex] = sol;
						break;
					case LEVEL_SET_F_D_FAST:
						m_pDistanceMapFromPS_FAST_WAVE[lIndex] = sol;
						break;
					}

					Node New;
					New.time = sol;
					New.x = k;
					New.y = l;
					New.z = m;
					HeapInsert(New); /* STEP 4 */
				}
			}
		}

		}
	}
}

float JJSkeleton::solve_ORD_1(int i1, int j1, int k1, int i2, int j2, int k2, int i3, int j3, int k3, float sol, short nDistanceMapMode)
{
	float r, s;
	unsigned long lIndex1 = i1 + j1*m_VolX + k1*m_VolX*m_VolY;
	unsigned long lIndex2 = i2 + j2*m_VolX + k2*m_VolX*m_VolY; 
	unsigned long lIndex3 = i3 + j3*m_VolX + k3*m_VolX*m_VolY; 

	float T1;
	float T2;
	float T3;

	switch (nDistanceMapMode)
	{
	case LEVEL_SET_F_1:
		T1 = m_pLongDistanceMap[lIndex1];
		T2 = m_pLongDistanceMap[lIndex2];
		T3 = m_pLongDistanceMap[lIndex3];	
		break;
	case LEVEL_SET_F_D:
		T1 = m_pDistanceMapFromPS[lIndex1];
		T2 = m_pDistanceMapFromPS[lIndex2];
		T3 = m_pDistanceMapFromPS[lIndex3];	
		break;
	case LEVEL_SET_F_D_FAST:
		T1 = m_pDistanceMapFromPS_FAST_WAVE[lIndex1];
		T2 = m_pDistanceMapFromPS_FAST_WAVE[lIndex2];
		T3 = m_pDistanceMapFromPS_FAST_WAVE[lIndex3];	
		break;
	}

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

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