jjdistancemap.cpp

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

	for(w = 0; w < nNumEdge2; w++)
	{			
		i = xArr2[w];
		j = yArr2[w];
		k = zArr2[w];
		for (x = -1; x <= 1; x++)
			for (y = -1; y <= 1; y++)
				for (z = -1; z <= 1; z++)
				{
					if ((0 <= (i + x)) && ((i + x) < m_VolX) &&
						(0 <= (j + y)) && ((j + y) < m_VolY) &&
						(0 <= (k + z)) && ((k + z) < m_VolZ))
					{
						iCurIndex = (i + x) + (j + y)*m_VolX + (k + z)*nPlane;
						if (m_pDistanceMap[iCurIndex] == m_iMaxDistance)
						{
							m_pDistanceMap[iCurIndex] = 4;
							xArr1[nNumEdge1] = i + x;
							yArr1[nNumEdge1] = j + y;
							zArr1[nNumEdge1] = k + z;
							nNumEdge1++;
						}
					}
				}
	}

	// distance 5
	nNumEdge2 = 0;
	for(w = 0; w < nNumEdge1; w++)
	{			
		i = xArr1[w];
		j = yArr1[w];
		k = zArr1[w];
		for (x = -1; x <= 1; x++)
			for (y = -1; y <= 1; y++)
				for (z = -1; z <= 1; z++)
				{
					if ((0 <= (i + x)) && ((i + x) < m_VolX) &&
						(0 <= (j + y)) && ((j + y) < m_VolY) &&
						(0 <= (k + z)) && ((k + z) < m_VolZ))
					{
						iCurIndex = (i + x) + (j + y)*m_VolX + (k + z)*nPlane;
						if (m_pDistanceMap[iCurIndex] == m_iMaxDistance)
						{
							m_pDistanceMap[iCurIndex] = 5;
							xArr2[nNumEdge2] = i + x;
							yArr2[nNumEdge2] = j + y;
							zArr2[nNumEdge2] = k + z;
							nNumEdge2++;
						}
					}
				}
	}

	// distance 6
	nNumEdge1 = 0;
	for(w = 0; w < nNumEdge2; w++)
	{			
		i = xArr2[w];
		j = yArr2[w];
		k = zArr2[w];
		for (x = -1; x <= 1; x++)
			for (y = -1; y <= 1; y++)
				for (z = -1; z <= 1; z++)
				{
					if ((0 <= (i + x)) && ((i + x) < m_VolX) &&
						(0 <= (j + y)) && ((j + y) < m_VolY) &&
						(0 <= (k + z)) && ((k + z) < m_VolZ))
					{
						iCurIndex = (i + x) + (j + y)*m_VolX + (k + z)*nPlane;
						if (m_pDistanceMap[iCurIndex] == m_iMaxDistance)
						{
							m_pDistanceMap[iCurIndex] = 6;
							xArr1[nNumEdge1] = i + x;
							yArr1[nNumEdge1] = j + y;
							zArr1[nNumEdge1] = k + z;
							nNumEdge1++;
						}
					}
				}
	}

// 蜡公磊
	
	// distance 7
	nNumEdge2 = 0;
	for(w = 0; w < nNumEdge1; w++)
	{			
		i = xArr1[w];
		j = yArr1[w];
		k = zArr1[w];
		for (x = -1; x <= 1; x++)
			for (y = -1; y <= 1; y++)
				for (z = -1; z <= 1; z++)
				{
					if ((0 <= (i + x)) && ((i + x) < m_VolX) &&
						(0 <= (j + y)) && ((j + y) < m_VolY) &&
						(0 <= (k + z)) && ((k + z) < m_VolZ))
					{
						iCurIndex = (i + x) + (j + y)*m_VolX + (k + z)*nPlane;
						if (m_pDistanceMap[iCurIndex] == m_iMaxDistance)
						{
							m_pDistanceMap[iCurIndex] = 7;
							xArr2[nNumEdge2] = i + x;
							yArr2[nNumEdge2] = j + y;
							zArr2[nNumEdge2] = k + z;
							nNumEdge2++;
						}
					}
				}
	}

	// distance 8
	nNumEdge1 = 0;
	for(w = 0; w < nNumEdge2; w++)
	{			
		i = xArr2[w];
		j = yArr2[w];
		k = zArr2[w];
		for (x = -1; x <= 1; x++)
			for (y = -1; y <= 1; y++)
				for (z = -1; z <= 1; z++)
				{
					if ((0 <= (i + x)) && ((i + x) < m_VolX) &&
						(0 <= (j + y)) && ((j + y) < m_VolY) &&
						(0 <= (k + z)) && ((k + z) < m_VolZ))
					{
						iCurIndex = (i + x) + (j + y)*m_VolX + (k + z)*nPlane;
						if (m_pDistanceMap[iCurIndex] == m_iMaxDistance)
						{
							m_pDistanceMap[iCurIndex] = 8;
							xArr1[nNumEdge1] = i + x;
							yArr1[nNumEdge1] = j + y;
							zArr1[nNumEdge1] = k + z;
							nNumEdge1++;
						}
					}
				}
	}

	// distance 9
	nNumEdge2 = 0;
	for(w = 0; w < nNumEdge1; w++)
	{			
		i = xArr1[w];
		j = yArr1[w];
		k = zArr1[w];
		for (x = -1; x <= 1; x++)
			for (y = -1; y <= 1; y++)
				for (z = -1; z <= 1; z++)
				{
					if ((0 <= (i + x)) && ((i + x) < m_VolX) &&
						(0 <= (j + y)) && ((j + y) < m_VolY) &&
						(0 <= (k + z)) && ((k + z) < m_VolZ))
					{
						iCurIndex = (i + x) + (j + y)*m_VolX + (k + z)*nPlane;
						if (m_pDistanceMap[iCurIndex] == m_iMaxDistance)
						{
							m_pDistanceMap[iCurIndex] = 9;
							xArr2[nNumEdge2] = i + x;
							yArr2[nNumEdge2] = j + y;
							zArr2[nNumEdge2] = k + z;
							nNumEdge2++;
						}
					}
				}
	}

	// distance 10
	nNumEdge1 = 0;
	for(w = 0; w < nNumEdge2; w++)
	{			
		i = xArr2[w];
		j = yArr2[w];
		k = zArr2[w];
		for (x = -1; x <= 1; x++)
			for (y = -1; y <= 1; y++)
				for (z = -1; z <= 1; z++)
				{
					if ((0 <= (i + x)) && ((i + x) < m_VolX) &&
						(0 <= (j + y)) && ((j + y) < m_VolY) &&
						(0 <= (k + z)) && ((k + z) < m_VolZ))
					{
						iCurIndex = (i + x) + (j + y)*m_VolX + (k + z)*nPlane;
						if (m_pDistanceMap[iCurIndex] == m_iMaxDistance)
						{
							m_pDistanceMap[iCurIndex] = 10;
							xArr1[nNumEdge1] = i + x;
							yArr1[nNumEdge1] = j + y;
							zArr1[nNumEdge1] = k + z;
							nNumEdge1++;
						}
					}
				}
	}

	delete [] xArr1;
	delete [] yArr1;
	delete [] zArr1;

	delete [] xArr2;
	delete [] yArr2;
	delete [] zArr2;

/*	char a, b, c;

	for (z = 0; z < m_VolZ; z++) {
		strcpy(tfile, "temp");
		a = z % 10 + '0';
		b = (z % 100) / 10 + '0';
		c = (z - z % 100) / 100 + '0';
		tfile[4] = c;
		tfile[5] = b;
		tfile[6] = a;
		tfile[7] = '\0';
		strcat(tfile, ".raw");
		if ((tfp = fopen(tfile, "wb"))==NULL) {
			printf("\n Cannot open file: %s", tfile);
			exit(0);
		}

// raw file print
		for (y = 0; y < m_VolY; y++) {
			for (x = 0; x < m_VolX; x++) {
				unsigned char dist;

				dist = (unsigned char)m_pDistanceMap[(x) + (y)*m_VolX + (z)*m_VolX*m_VolY];
				
				fprintf(tfp, "%c", dist);
			}
		}
		fclose(tfp);
	}*/
}

void JJDistanceMap::CalEuclideanDistance(void)
{
	int i, j, k, n, iCurIndex;
	int L = m_VolX, M = m_VolY, N = m_VolZ;
	int df, db, d, w;
	int buff[2025];
	int rMax, rStart, rEnd;
	int nPlane = m_VolX*m_VolY;

	char		t;		/* temporary character variable */
	char		tfile[20];	/* temporary name file */
	FILE		*tfp;		/* temporary file description */

	// 512*512*512*0.5*2B = 125 MB
	m_pDBF = (unsigned int *)malloc(m_VolX*m_VolY*m_VolZ*sizeof(unsigned int));
	memset((void *)m_pDBF, 0, m_VolX*m_VolY*m_VolZ*sizeof(unsigned int));
	
	// New Algorithms for Euclidean Distance Transformation of An n-Dimensional Digitized
	// Picture with Applications
	// initialization m_pDBF
	for (k = 0; k < m_VolZ; k++)
	{
		for (j = 0; j < m_VolY; j++)
		{
			for (i = 0; i < m_VolX; i++)
			{
				iCurIndex = i + j*m_VolX + k*nPlane;
				if (m_pBinaryVolume->IsEdge(i, j, k))
					m_pDBF[iCurIndex] = 0; 
				else
					m_pDBF[iCurIndex] = 1; 
			}
		}
	}

	// (Step 1)
	// forward scan
	for (k = 1; k <= N; k++)
	{
		for (j = 1; j <= M; j++)
		{
			df = L;
			for (i = 1; i <= L; i++)
			{
				iCurIndex = (i - 1) + (j - 1)*m_VolX + (k - 1)*nPlane;
				if (m_pDBF[iCurIndex])
					df = df + 1;
				else
					df = 0;
				m_pDBF[iCurIndex] = df*df;
			}
		}
	}

	// backward scan
	for (k = 1; k <= N; k++)
	{
		for (j = 1; j <= M; j++)
		{
			db = L;
			for (i = L; i >= 1; i--)
			{
				iCurIndex = (i - 1) + (j - 1)*m_VolX + (k - 1)*nPlane;
				if (m_pDBF[iCurIndex])
					db = db + 1;
				else
					db = 0;
				m_pDBF[iCurIndex] = minimum(m_pDBF[iCurIndex], db*db);
			}
		}
	}

	// (Step 2)
	for (k = 1; k <= N; k++)
	{
		for (i = 1; i <= L; i++)
		{
			for (j = 1; j <= M; j++)
			{
				iCurIndex = (i - 1) + (j - 1)*m_VolX + (k - 1)*nPlane;
				buff[j] = m_pDBF[iCurIndex];
			}
			for (j = 1; j <= M; j++)
			{
				d = buff[j];
				if (d)
				{
					rMax = (int)sqrt(d) + 1;
					rStart = Min(rMax, (j - 1));
					rEnd = Min(rMax, (M - j));
					for (n = -rStart; n <= rEnd; n++)
					{
						w = buff[j + n] + n*n;
						if (w < d)
							d = w;
					}
				}
				iCurIndex = (i - 1) + (j - 1)*m_VolX + (k - 1)*nPlane;
				m_pDBF[iCurIndex] = d;
			}
		}
	}

	// (Step 3)
	for (j = 1; j <= M; j++)
	{
		for (i = 1; i <= L; i++)
		{
			for (k = 1; k <= N; k++)
			{
				iCurIndex = (i - 1) + (j - 1)*m_VolX + (k - 1)*nPlane;
				buff[k] = m_pDBF[iCurIndex];
			}
			for (k = 1; k <= N; k++)
			{
				d = buff[k];
				if (d)
				{
					rMax = (int)sqrt(d) + 1;
					rStart = Min(rMax, (k - 1));
					rEnd = Min(rMax, (N - k));
					for (n = -rStart; n <= rEnd; n++)
					{
						w = buff[k + n] + n*n;
						if (w < d)
							d = w;
					}
				}
				iCurIndex = (i - 1) + (j - 1)*m_VolX + (k - 1)*nPlane;
				m_pDBF[iCurIndex] = d;
			}
		}
	}

	// square root of the distance, m_pDBF
	for (k = 0; k < m_VolZ; k++)
		for (j = 0; j < m_VolY; j++)
			for (i = 0; i < m_VolX; i++)
			{
				iCurIndex = i + j*m_VolX + k*nPlane;
				m_pDBF[iCurIndex] = sqrt(m_pDBF[iCurIndex]);
			}

/*	char a, b, c;
	int x, y, z;

	for (z = 0; z < m_VolZ; z++) {
		strcpy(tfile, "temp");
		a = z % 10 + '0';
		b = (z % 100) / 10 + '0';
		c = (z - z % 100) / 100 + '0';
		tfile[4] = c;
		tfile[5] = b;
		tfile[6] = a;
		tfile[7] = '\0';
		strcat(tfile, ".raw");
		if ((tfp = fopen(tfile, "wb"))==NULL) {
			printf("\n Cannot open file: %s", tfile);
			exit(0);
		}

// raw file print
		for (y = 0; y < m_VolY; y++) {
			for (x = 0; x < m_VolX; x++) {
				unsigned char dist;

				dist = (unsigned char)m_pDBF[(x) + (y)*m_VolX + (z)*m_VolX*m_VolY];
				
				fprintf(tfp, "%c", dist);
			}
		}
// text file print
//		for (y = 0; y < m_VolY; y++) {
//			for (x = 0; x < m_VolX; x++) {
//				fprintf(tfp, "%3d", (unsigned short)m_pDBF[(x) + (y)*m_VolX + (z)*m_VolX*m_VolY]);
//			}
//			fprintf(tfp, "\n\r");
//		}
		fclose(tfp);
	}*/
}

short JJDistanceMap::Min(short param1, short param2)
{
	if (param1 < param2)
		return param1;
	else
		return param2;
}