autosample3d.cpp

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

										  VoxelOffset(x,y-1,z-1)-VoxelOffset(x,y+1,z-1)))+
								         (VoxelOffset(x,y,z+1)-VoxelOffset(x,y,z-1)));

				// calculate the magnitude of gradient
				double	absgrad = (double)sqrt((double)(gradx*gradx+grady*grady+gradz*gradz));

				// calculate opacity
				if(absgrad == 0.0 && val == m_Threshold)
					opacity = 1.0;
				else if(absgrad == 0.0)
					opacity = 0.0;
				else {
					double	upper = abs(m_Threshold-val);
					double	lower = absgrad;
					double	temp1 = (1/m_Tolerance)*(upper/lower);
					if(temp1 <= 1.0)	opacity = 1.0;
					else				opacity = 0.0;
				}

				// Assign opacity value to output memory
				if(opacity > 0.4) {
					unsigned short temp2 = (unsigned short)(opacity*255);
					if(temp2 > m_Threshold)
						(*levoy)[((z*m_pResol->Rv)+y)*m_pResol->Ru+x] = 255;
				} // end of if
			} //end of x for
		} // end of y for
	} // end of z for

	return TRUE;
}

//*********************************************************************
//	Sobel3D(): Edge detection with sobel operator
//			   (using Intel Image Processing Library)
//*********************************************************************
BOOL RxAutoSample3D::Sobel3D(BYTE **sobel)
{
	// Convolve Sx, Sy, Sz mask with original volume dataset
	//
	//	-1	-2	-1			0	0	0			1	2	1
	//	-2	-1	-2			0	0	0			2	1	2
	//	-1	-2	-1			0	0	0			1	2	1
	//
	for (int z = 1; z < m_pResol->Rw-1; z++) {
		// Allocate x, y, z magnitude buffer
		AllocateImageData(&m_MagXBuf, *m_pResol);
		AllocateImageData(&m_MagYBuf, *m_pResol);
		AllocateImageData(&m_MagZBuf, *m_pResol);

		// Calculate x, y, z gradient
		MakeDX(z, &m_MagXBuf);
		MakeDY(z, &m_MagYBuf);
		MakeDZ(z, &m_MagZBuf);

		// Calculate the magnitude of gradient in each slice
		for (int y = 1; y < m_pResol->Rv-1; y++) {
			for (int x = 1; x < m_pResol->Ru-1; x++) {
				int	index2D = y * m_pResol->Ru + x;
				int index3D = ((z*m_pResol->Rv)+y)*m_pResol->Ru+x;

				unsigned short gradx = m_MagXBuf[index2D];
				unsigned short grady = m_MagYBuf[index2D];
				unsigned short gradz = m_MagZBuf[index2D];

				double absgrad = (double)sqrt((double)(gradx*gradx+grady*grady+gradz*gradz));

				if(absgrad > m_Threshold) {
					(*sobel)[index3D] = 255;
				}
			} // end of x for
		} // end of y for

		// Deallocate x, y, z magnitude buffer
		ReleaseImageData(m_MagXBuf, *m_pResol);
		ReleaseImageData(m_MagYBuf, *m_pResol);
		ReleaseImageData(m_MagZBuf, *m_pResol);
	} // end of z for

	return TRUE;
}

//*********************************************************************
//	SobelLevoy3D(): Edge detection with sobel operator and levoy filtering
//*********************************************************************
BOOL RxAutoSample3D::SobelLevoy3D(BYTE **levoy)
{
	// Convolve Sx, Sy, Sz mask with original volume dataset
	//
	//	-1	-2	-1			0	0	0			1	2	1
	//	-2	-1	-2			0	0	0			2	1	2
	//	-1	-2	-1			0	0	0			1	2	1
	//
	unsigned short	val;
	double			opacity = 0.0;
	for (int z = 1; z < m_pResol->Rw-1; z++) {
		// Allocate x, y, z magnitude buffer
		AllocateImageData(&m_MagXBuf, *m_pResol);
		AllocateImageData(&m_MagYBuf, *m_pResol);
		AllocateImageData(&m_MagZBuf, *m_pResol);

		// Calculate gradient
		MakeDX(z, &m_MagXBuf);
		MakeDY(z, &m_MagYBuf);
		MakeDZ(z, &m_MagZBuf);

		// Calculate the magnitude of gradient
		for (int y = 1; y < m_pResol->Rv-1; y++) {
			for (int x = 1; x < m_pResol->Ru-1; x++) {
				int	index2D = y * m_pResol->Ru + x;
				int index3D = ((z*m_pResol->Rv)+y)*m_pResol->Ru+x;

				if((val=m_pnVolume[index3D]) < m_Threshold)
					opacity = 0.0;

				unsigned short gradx = m_MagXBuf[index2D];
				unsigned short grady = m_MagYBuf[index2D];
				unsigned short gradz = m_MagZBuf[index2D];

				double absgrad = (double)sqrt((double)(gradx*gradx+grady*grady+gradz*gradz));
				if(absgrad == 0.0 && val == m_Threshold)
					opacity = 1.0;
				else if(absgrad == 0.0)
					opacity = 0.0;
				else {
					double	upper = abs(m_Threshold-val);
					double	lower = absgrad;
					double	temp1 = (1/m_Tolerance) * (upper/lower);
					if(temp1 <= 1.0)	opacity = 1.0;
					else				opacity = 0.0;
				}

				// Assign opacity value to output memory
				if(opacity > 0.4) {
					unsigned short temp2 = (unsigned short)(opacity * 255);
					if(temp2 > m_Threshold)
						(*levoy)[index3D] = 255;
				} // end of if
			} // end of x for
		} // end of y for

		// Deallocate x, y, z magnitude buffer
		ReleaseImageData(m_MagXBuf, *m_pResol);
		ReleaseImageData(m_MagYBuf, *m_pResol);
		ReleaseImageData(m_MagZBuf, *m_pResol);
	} // end of z for

	return TRUE;
}

//*********************************************************************
//	MakeDX(): Calculate X gradient
//*********************************************************************
BOOL RxAutoSample3D::MakeDX(int z, unsigned short **pnDstSlice)
{
	IplROI *srcROI1 = iplCreateROI(0, 2, 0, m_pResol->Ru-2, m_pResol->Rv);
	IplROI *srcROI2 = iplCreateROI(0, 0, 0, m_pResol->Ru-2, m_pResol->Rv);
	IplROI *tmpROI  = iplCreateROI(0, 1, 0, m_pResol->Ru-2, m_pResol->Rv);

	IplImage *srcImage1 = iplCreateImageHeader(
		1, 0, IPL_DEPTH_16U,		// data of byte type : unsigned short
		"Gray", "Gray",				// color order
		IPL_DATA_ORDER_PIXEL,		// channel arrangement
		IPL_ORIGIN_TL,				// top-left orientation
		IPL_ALIGN_QWORD,			// 8-byte align
		m_pResol->Ru, m_pResol->Rv,	// image width and height
		NULL, NULL, NULL, NULL);	// not tiled
	if(srcImage1 == NULL)
		return FALSE;

	int index = z * m_pResol->Rv * m_pResol->Ru;
	srcImage1->imageData = (char *)(m_pnVolume + index);
	srcImage1->roi = srcROI1;

	IplImage *srcImage2 = iplCreateImageHeader(
		1, 0, IPL_DEPTH_16U,		// data of byte type : unsigned short
		"Gray", "Gray",				// color order
		IPL_DATA_ORDER_PIXEL,		// channel arrangement
		IPL_ORIGIN_TL,				// top-left orientation
		IPL_ALIGN_QWORD,			// 8-byte align
		m_pResol->Ru, m_pResol->Rv,	// image width and height
		NULL, NULL, NULL, NULL);	// not tiled
	if(srcImage2 == NULL)
		return FALSE;

	srcImage2->imageData = (char *)(m_pnVolume + index);
	srcImage2->roi = srcROI2;

	IplImage *tmpImage = iplCreateImageHeader(
		1, 0, IPL_DEPTH_16U,		// data of byte type
		"Gray", "Gray",				// color order
		IPL_DATA_ORDER_PIXEL,		// channel arrangement
		IPL_ORIGIN_TL,				// top-left orientation
		IPL_ALIGN_QWORD,			// 8-byte align
		m_pResol->Ru, m_pResol->Rv,	// image width and height
		NULL, NULL, NULL, NULL);	// not tiled
	if(tmpImage == NULL)
		return FALSE;

	iplAllocateImage(tmpImage, 1, 0);
	tmpImage->roi = tmpROI;
	iplSubtract(srcImage1, srcImage2, tmpImage);
	tmpImage->roi = NULL;

	IplImage *dstImage = iplCreateImageHeader(
		1, 0, IPL_DEPTH_16U,		// data of byte type
		"Gray", "Gray",				// color order
		IPL_DATA_ORDER_PIXEL,		// channel arrangement
		IPL_ORIGIN_TL,				// top-left orientation
		IPL_ALIGN_QWORD,			// 8-byte align
		m_pResol->Ru, m_pResol->Rv,	// image width and height
		NULL, NULL, NULL, NULL);	// not tiled
	if(dstImage == NULL)
		return FALSE;
	dstImage->imageData = (char *)*pnDstSlice;

	iplCopy(tmpImage, dstImage);

	iplDeallocate(srcImage1, IPL_IMAGE_HEADER);
	iplDeallocate(srcImage2, IPL_IMAGE_HEADER);
	iplDeallocate(tmpImage, IPL_IMAGE_ALL);
	iplDeallocate(dstImage, IPL_IMAGE_HEADER);

	iplDeleteROI(srcROI1);
	iplDeleteROI(srcROI2);
	iplDeleteROI(tmpROI);

	return TRUE;
}

//*********************************************************************
//	MakeDY(): Calculate Y gradient
//*********************************************************************
BOOL RxAutoSample3D::MakeDY(int z, unsigned short **pnDstSlice)
{
	IplROI *srcROI1 = iplCreateROI(0, 0, 2, m_pResol->Ru, m_pResol->Rv-2);
	IplROI *srcROI2 = iplCreateROI(0, 0, 0, m_pResol->Ru, m_pResol->Rv-2);
	IplROI *tmpROI  = iplCreateROI(0, 0, 1, m_pResol->Ru, m_pResol->Rv-2);

	IplImage *srcImage1 = iplCreateImageHeader(
		1, 0, IPL_DEPTH_16U,			// data of byte type : unsigned short
		"Gray", "Gray",					// color order
		IPL_DATA_ORDER_PIXEL,			// channel arrangement
		IPL_ORIGIN_TL,					// top-left orientation
		IPL_ALIGN_QWORD,				// 8-byte align
		m_pResol->Ru, m_pResol->Rv,		// image width and height
		NULL, NULL, NULL, NULL);		// not tiled
	if(srcImage1 == NULL)
		return FALSE;
	
	int index = z * m_pResol->Rv * m_pResol->Ru;
	srcImage1->imageData = (char *)(m_pnVolume + index);
	srcImage1->roi = srcROI1;

	IplImage *srcImage2 = iplCreateImageHeader(
		1, 0, IPL_DEPTH_16U,			// data of byte type
		"Gray", "Gray",					// color order
		IPL_DATA_ORDER_PIXEL,			// channel arrangement
		IPL_ORIGIN_TL,					// top-left orientation
		IPL_ALIGN_QWORD,				// 8-byte align
		m_pResol->Ru, m_pResol->Rv,		// image width and height
		NULL, NULL, NULL, NULL);		// not tiled
	if(srcImage2 == NULL)
		return FALSE;
	
	srcImage2->imageData = (char *)(m_pnVolume + index);
	srcImage2->roi = srcROI2;

	IplImage *tmpImage = iplCreateImageHeader(
		1, 0, IPL_DEPTH_16U,			// data of byte type : unsigned short
		"Gray", "Gray",					// color order
		IPL_DATA_ORDER_PIXEL,			// channel arrangement
		IPL_ORIGIN_TL,					// top-left orientation
		IPL_ALIGN_QWORD,				// 8-byte align
		m_pResol->Ru, m_pResol->Rv,		// image width and height
		NULL, NULL, NULL, NULL);		// not tiled
	if(tmpImage == NULL)
		return FALSE;

	iplAllocateImage(tmpImage, 1, 0);
	tmpImage->roi = tmpROI;
	iplSubtract(srcImage1, srcImage2, tmpImage);
	tmpImage->roi = NULL;

	IplImage *dstImage = iplCreateImageHeader(
		1, 0, IPL_DEPTH_16U,			// data of byte type : unsigned short
		"Gray", "Gray",					// color order
		IPL_DATA_ORDER_PIXEL,			// channel orientation
		IPL_ORIGIN_TL,					// top-left orientation
		IPL_ALIGN_QWORD,				// 8-byte align
		m_pResol->Ru, m_pResol->Rv,		// image width and height
		NULL, NULL, NULL, NULL);		// not tiled
	if(dstImage == NULL)
		return FALSE;

	dstImage->imageData = (char *)*pnDstSlice;

	iplCopy(tmpImage, dstImage);

	iplDeallocate(srcImage1, IPL_IMAGE_HEADER);
	iplDeallocate(srcImage2, IPL_IMAGE_HEADER);
	iplDeallocate(tmpImage, IPL_IMAGE_ALL);
	iplDeallocate(dstImage, IPL_IMAGE_HEADER);

	iplDeleteROI(srcROI1);
	iplDeleteROI(srcROI2);
	iplDeleteROI(tmpROI);

	return TRUE;
}

//*********************************************************************
//	MakeDY(): Calculate Z gradient
//*********************************************************************
BOOL RxAutoSample3D::MakeDZ(int z, unsigned short **pnDstSlice)
{
	IplImage *srcImage1 = iplCreateImageHeader(
		1, 0, IPL_DEPTH_16U,			// data of byte type : unsigned short
		"Gray", "Gray",					// color order
		IPL_DATA_ORDER_PIXEL,			// channel arrangement
		IPL_ORIGIN_TL,					// top-left orientation
		IPL_ALIGN_QWORD,				// 8-byte align
		m_pResol->Ru, m_pResol->Rv,		// image width and height
		NULL, NULL, NULL, NULL);		// not tiled
	if(srcImage1 == NULL)
		return FALSE;

	// z+1 slice
	int index1 = (z+1) * m_pResol->Rv * m_pResol->Ru;
	srcImage1->imageData = (char *)(m_pnVolume + index1);

	IplImage *srcImage2 = iplCreateImageHeader(
		1, 0, IPL_DEPTH_16U,			// data of byte type
		"Gray", "Gray",					// color order
		IPL_DATA_ORDER_PIXEL,			// channel arrangement
		IPL_ORIGIN_TL,					// top-left orientation
		IPL_ALIGN_QWORD,				// 8 byte align
		m_pResol->Ru, m_pResol->Rv,		// image width and height
		NULL, NULL, NULL, NULL);		// not tiled
	if(srcImage2 == NULL)
		return FALSE;

	// z-1 slice
	int index2 = (z-1) * m_pResol->Rv * m_pResol->Ru;
	srcImage2->imageData = (char *)(m_pnVolume + index2);

	IplImage *tmpImage = iplCreateImageHeader(
		1, 0, IPL_DEPTH_16U,			// data of byte type : unsigned short
		"Gray", "Gray",					// color order
		IPL_DATA_ORDER_PIXEL,			// channel arrangement
		IPL_ORIGIN_TL,					// top-left orientation
		IPL_ALIGN_QWORD,				// 8-byte align
		m_pResol->Ru, m_pResol->Rv,		// image width and height
		NULL, NULL, NULL, NULL);		// not tiled
	if(tmpImage == NULL)
		return FALSE;

	iplAllocateImage(tmpImage, 1, 0);
	iplSubtract(srcImage1, srcImage2, tmpImage);

	IplImage *dstImage = iplCreateImageHeader(
		1, 0, IPL_DEPTH_16U,			// data of byte type : unsigned short
		"Gray", "Gray",					// color order
		IPL_DATA_ORDER_PIXEL,			// channel arrangement
		IPL_ORIGIN_TL,					// top-left orientation
		IPL_ALIGN_QWORD,				// 8-byte align
		m_pResol->Ru, m_pResol->Rv,		// image width and height
		NULL, NULL, NULL, NULL);		// not tiled
	if(dstImage == NULL)
		return FALSE;

	dstImage->imageData = (char *)*pnDstSlice;
	iplCopy(tmpImage, dstImage);

	iplDeallocate(srcImage1, IPL_IMAGE_HEADER);
	iplDeallocate(srcImage2, IPL_IMAGE_HEADER);
	iplDeallocate(tmpImage, IPL_IMAGE_ALL);
	iplDeallocate(dstImage, IPL_IMAGE_HEADER);

	return TRUE;
}

//*********************************************************************
//	VoxelOffset(): calculate voxel offset for indexing
//*********************************************************************
inline unsigned short RxAutoSample3D::VoxelOffset(int x, int y, int z)
{
	return m_pnVolume[((z*m_pResol->Rv)+y)*m_pResol->Ru+x];
}