jjdistancemap.cpp

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

						iCurIndex = (i + x) + (j + y)*m_VolX + (k + z)*nPlane;
						if (fabs(m_pLongDistanceMap[iCurIndex] - 500.0f) < 0.001f)
						{
							if (z == 0)
								m_pLongDistanceMap[iCurIndex] = nCenterDist + 1.0;
							else
								m_pLongDistanceMap[iCurIndex] = nCenterDist + m_lfZRatioXY;
							xArr2[nNumEdge2] = i + x;
							yArr2[nNumEdge2] = j + y;
							zArr2[nNumEdge2] = k + z;
							nNumEdge2++;
						}
					}
				}
	}

	// distance 4
	nNumEdge1 = 0;
	for(w = 0; w < nNumEdge2; w++)
	{			
		i = xArr2[w];
		j = yArr2[w];
		k = zArr2[w];
		iCurIndex = i + j*m_VolX + k*nPlane;
		nCenterDist = m_pLongDistanceMap[iCurIndex];

		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 (fabs(m_pLongDistanceMap[iCurIndex] - 500.0f) < 0.001f)
						{
							if (z == 0)
								m_pLongDistanceMap[iCurIndex] = nCenterDist + 1.0;
							else
								m_pLongDistanceMap[iCurIndex] = nCenterDist + m_lfZRatioXY;
							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];
		iCurIndex = i + j*m_VolX + k*nPlane;
		nCenterDist = m_pLongDistanceMap[iCurIndex];

		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 (fabs(m_pLongDistanceMap[iCurIndex] - 500.0f) < 0.001f)
						{
							if (z == 0)
								m_pLongDistanceMap[iCurIndex] = nCenterDist + 1.0;
							else
								m_pLongDistanceMap[iCurIndex] = nCenterDist + m_lfZRatioXY;
							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];
		iCurIndex = i + j*m_VolX + k*nPlane;
		nCenterDist = m_pLongDistanceMap[iCurIndex];

		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 (fabs(m_pLongDistanceMap[iCurIndex] - 500.0f) < 0.001f)
						{
							if (z == 0)
								m_pLongDistanceMap[iCurIndex] = nCenterDist + 1.0;
							else
								m_pLongDistanceMap[iCurIndex] = nCenterDist + m_lfZRatioXY;
							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];
		iCurIndex = i + j*m_VolX + k*nPlane;
		nCenterDist = m_pLongDistanceMap[iCurIndex];

		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 (fabs(m_pLongDistanceMap[iCurIndex] - 500.0f) < 0.001f)
						{
							if (z == 0)
								m_pLongDistanceMap[iCurIndex] = nCenterDist + 1.0;
							else
								m_pLongDistanceMap[iCurIndex] = nCenterDist + m_lfZRatioXY;
							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];
		iCurIndex = i + j*m_VolX + k*nPlane;
		nCenterDist = m_pLongDistanceMap[iCurIndex];

		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 (fabs(m_pLongDistanceMap[iCurIndex] - 500.0f) < 0.001f)
						{
							if (z == 0)
								m_pLongDistanceMap[iCurIndex] = nCenterDist + 1.0;
							else
								m_pLongDistanceMap[iCurIndex] = nCenterDist + m_lfZRatioXY;
							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];
		iCurIndex = i + j*m_VolX + k*nPlane;
		nCenterDist = m_pLongDistanceMap[iCurIndex];

		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 (fabs(m_pLongDistanceMap[iCurIndex] - 500.0f) < 0.001f)
						{
							if (z == 0)
								m_pLongDistanceMap[iCurIndex] = nCenterDist + 1.0;
							else
								m_pLongDistanceMap[iCurIndex] = nCenterDist + m_lfZRatioXY;
							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];
		iCurIndex = i + j*m_VolX + k*nPlane;
		nCenterDist = m_pLongDistanceMap[iCurIndex];

		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 (fabs(m_pLongDistanceMap[iCurIndex] - 500.0f) < 0.001f)
						{
							if (z == 0)
								m_pLongDistanceMap[iCurIndex] = nCenterDist + 1.0;
							else
								m_pLongDistanceMap[iCurIndex] = nCenterDist + m_lfZRatioXY;
							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_pLongDistanceMap[(x) + (y)*m_VolX + (z)*m_VolX*m_VolY];
				
				fprintf(tfp, "%c", dist);
			}
		}
		fclose(tfp);
	}*/
}

void JJDistanceMap::CalPropagationDistance(void)
{
	int i, j, k, x, y, z, w, iCurIndex, nNumEdge1 = 0, nNumEdge2 = 0;
	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_pDistanceMap = (unsigned char *)malloc(m_VolX*m_VolY*m_VolZ*sizeof(unsigned char));
	memset((void *)m_pDistanceMap, 0, m_VolX*m_VolY*m_VolZ*sizeof(unsigned char));

	short* xArr1 = new short[4000000]; 
	short* yArr1 = new short[4000000];
	short* zArr1 = new short[4000000];	

	short* xArr2 = new short[4000000]; 
	short* yArr2 = new short[4000000];
	short* zArr2 = new short[4000000];

	// distance = 0
	nNumEdge1 = 0;
	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_pDistanceMap[iCurIndex] = 0; 
					xArr1[nNumEdge1] = i;
					yArr1[nNumEdge1] = j;
					zArr1[nNumEdge1] = k;
					nNumEdge1++;
				}
				else
					m_pDistanceMap[iCurIndex] = m_iMaxDistance; 
			}
		}
	}
	
	// distance 1
	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] = 1;
							xArr2[nNumEdge2] = i + x;
							yArr2[nNumEdge2] = j + y;
							zArr2[nNumEdge2] = k + z;
							nNumEdge2++;
						}
					}
				}
	}

	// distance 2
	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] = 2;
							xArr1[nNumEdge1] = i + x;
							yArr1[nNumEdge1] = j + y;
							zArr1[nNumEdge1] = k + z;
							nNumEdge1++;
						}
					}
				}
	}

	// distance 3
	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] = 3;
							xArr2[nNumEdge2] = i + x;
							yArr2[nNumEdge2] = j + y;
							zArr2[nNumEdge2] = k + z;
							nNumEdge2++;
						}
					}
				}
	}

	// distance 4
	nNumEdge1 = 0;