autosample3d.cpp

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

// AutoSample3D.cpp: implementation of the RxAutoSample3D class.
//
//////////////////////////////////////////////////////////////////////
//
//	Title: Automatic Sampling
//
//	1. Random sampling
//	2. Uniform sampling
//	3. Sobel-Levoy sampling
//		- Zucker operator with Levoy filtering
//		- Sobel operator with Levoy filtering (Intel Image Processing Library)
//
//////////////////////////////////////////////////////////////////////
//
//	Author: Helen Hong, 3DMed co. LTD
//	138-dong 417-ho Seoul National Univ.
//	San 56-1 Shinlim-dong Kwanak-gu, Seoul, Korea
//	Email. [e-mail]
//
//	Date	: 2002. 9. 10.
//	Update	: 2003. 1. 22.
//
////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////
//	include
////////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "AutoSample3D.h"

#include "VolumeInfo.h"
#include "AutoSample3DInfo.h"

#include "ipl.h"

#include <time.h>
#include <cmath>

////////////////////////////////////////////////////////////////////////
//	define
////////////////////////////////////////////////////////////////////////

#define	DESTROY_CLASS_INSTANCE(ci)	{if(ci) {delete ci; ci = NULL;}}

#define	SOBEL_LEVOY_OPERATOR
#define	SOBEL_OPERATOR

////////////////////////////////////////////////////////////////////////
//	declaration
////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

RxAutoSample3D::RxAutoSample3D()
{
	m_pnVolume		= NULL;

	m_pResol		= new RxVolumeInfo;
}

RxAutoSample3D::~RxAutoSample3D()
{
	DESTROY_CLASS_INSTANCE(m_pResol);
}

////////////////////////////////////////////////////////////////////////

//********************************************************************
//	SetVolumeData(): set volume data
//********************************************************************
BOOL RxAutoSample3D::SetVolumeData(unsigned short *Buf)
{
	m_pnVolume		= Buf;

	return TRUE;
}

//********************************************************************
//	SetVolumeResolution(): set volume resolution
//********************************************************************
BOOL RxAutoSample3D::SetVolumeResolution(RxVolumeInfo Resol)
{
	m_pResol->Ru	= Resol.Ru;
	m_pResol->Rv	= Resol.Rv;
	m_pResol->Rw	= Resol.Rw;

	return TRUE;
}

//********************************************************************
//	SetSampleInfo(): set sampling information
//********************************************************************
BOOL RxAutoSample3D::SetSampleInfo(int number)
{
	m_SampleNum		= number;			// number of sample points
	m_Threshold		= 200;				// threshold for levoy filtering
	m_Tolerance		= 1.0;				// tolerance for levoy filtering

	return TRUE;
}

//********************************************************************
//	AllocVolumeData(): allocate volume data
//********************************************************************
BOOL RxAutoSample3D::AllocateVolumeData(BYTE **Buf,RxVolumeInfo Resol)
{
	unsigned Size = Resol.Ru * Resol.Rv * Resol.Rw;
	*Buf = (BYTE *)VirtualAlloc(NULL, Size, 
		MEM_RESERVE | MEM_TOP_DOWN | MEM_COMMIT, PAGE_READWRITE);

	return TRUE;
}

//********************************************************************
//	ReleaseVolumeData(): release volume data
//********************************************************************
BOOL RxAutoSample3D::ReleaseVolumeData(BYTE *Buf,RxVolumeInfo Resol)
{
	unsigned Size = Resol.Ru * Resol.Rv * Resol.Rw;
	VirtualFree(Buf, Size, MEM_DECOMMIT);
	VirtualFree(Buf, 0, MEM_RELEASE);
	Buf = NULL;

	return TRUE;
}

//********************************************************************
//	AllocateImageDAta(): allocate image data
//********************************************************************
BOOL RxAutoSample3D::AllocateImageData(unsigned short **Buf,RxVolumeInfo Resol)
{
	unsigned Size = Resol.Ru * Resol.Rv * sizeof(unsigned short);
	*Buf = (unsigned short *)VirtualAlloc(NULL, Size, 
		MEM_RESERVE | MEM_TOP_DOWN | MEM_COMMIT, PAGE_READWRITE);

	return TRUE;
}

//********************************************************************
//	ReleaseImageData(): release image data
//********************************************************************
BOOL RxAutoSample3D::ReleaseImageData(unsigned short *Buf,RxVolumeInfo Resol)
{
	unsigned Size = Resol.Ru * Resol.Rv * sizeof(unsigned short);
	VirtualFree(Buf, Size, MEM_DECOMMIT);
	VirtualFree(Buf, 0, MEM_RELEASE);
	Buf = NULL;

	return TRUE;
}

//********************************************************************
//	UniformSampling(): Generate sample points using uniform sampling
//********************************************************************
BOOL RxAutoSample3D::UniformSampling(RxAutoSample3DInfo **sample)
{
	*sample = new RxAutoSample3DInfo[m_SampleNum];

	int iskipX = (int)((m_pResol->Ru/m_SampleNum)-1);
	int	iskipY = (int)((m_pResol->Rv/m_SampleNum)-1);
	int	iskipZ = (int)((m_pResol->Rw/m_SampleNum)-1);

	for (int count = 0; count < m_SampleNum; count++) {
		(*sample)[count].x	= (count+1)*iskipX;
		(*sample)[count].y	= (count+1)*iskipY;
		(*sample)[count].z	= (count+1)*iskipZ;
		(*sample)[count].val = VoxelOffset((*sample)[count].x,(*sample)[count].y,(*sample)[count].z);
	}

	return TRUE;
}

//********************************************************************
//	RandomSampling(): generate sample points using random sampling
//********************************************************************
BOOL RxAutoSample3D::RandomSampling(RxAutoSample3DInfo **sample)
{
	*sample = new RxAutoSample3DInfo[m_SampleNum];

	// seed the random-number generator with current time
	// so that the number will be different every time we run
	srand((unsigned)time(NULL));
	for (int count = 0; count < m_SampleNum; count++) {
		(*sample)[count].x	= (int)(((double)rand()/(double)RAND_MAX)*m_pResol->Ru);
		(*sample)[count].y	= (int)(((double)rand()/(double)RAND_MAX)*m_pResol->Rv);
		(*sample)[count].z	= (int)(((double)rand()/(double)RAND_MAX)*m_pResol->Rw);
		(*sample)[count].val = VoxelOffset((*sample)[count].x,(*sample)[count].y,(*sample)[count].z);
	}


	return TRUE;
}

//********************************************************************
//	SobelSampling(): Generate sample points using sobel-levoy sampling
//********************************************************************
BOOL RxAutoSample3D::SobelSampling(RxAutoSample3DInfo **sample)
{
	// Allocate memory for magnitude volume
	AllocateVolumeData(&m_MagBuf,*m_pResol);

#ifdef	ZUCKER_OPERATOR
	Zucker3D(&m_MagBuf);
#endif
#ifdef	SOBEL_OPERATOR
	Sobel3D(&m_MagBuf);
#endif

	// Allocate temporary memory
	RxAutoSample3DInfo	*temp;
	int					count=0, increment=0;
	temp = new RxAutoSample3DInfo[5000000];
	for (int z = 0; z < m_pResol->Rw; z++) {
		for (int y = 0; y < m_pResol->Rv; y++) {
			for (int x = 0; x < m_pResol->Ru; x++) {
				if(m_MagBuf[((z*m_pResol->Rv)+y)*m_pResol->Ru+x] == 255) {
					temp[count].x = x;	temp[count].y = y;
					temp[count].z = z;	temp[count].val = 0;
					count++;
				} // end of if
			} // end of x for
		} // end of y for
	} // end of z for
	ReleaseVolumeData(m_MagBuf,*m_pResol);

	// Allocate sample points A memory
	*sample = new RxAutoSample3DInfo[m_SampleNum];
	if(count <= m_SampleNum) {
		increment = 1;	m_SampleNum = count;
	}
	else {
		increment = count / m_SampleNum;
	}

	int i, j;
	for (i = 0, j = 0; j < m_SampleNum; i += increment, j++) {
		(*sample)[j].x	= temp[i].x;
		(*sample)[j].y	= temp[i].y;
		(*sample)[j].z	= temp[i].z;
		(*sample)[j].val = VoxelOffset(temp[i].x,temp[i].y,temp[i].z);
	}
	delete []temp;

	return TRUE;
}

//********************************************************************
//	SobelLevoySampling(): Generate sample points using sobel-levoy sampling
//********************************************************************
BOOL RxAutoSample3D::SobelLevoySampling(RxAutoSample3DInfo **sample)
{
	// Allocate memory for magnitude volume
	AllocateVolumeData(&m_MagBuf,*m_pResol);

#ifdef	ZUCKER_LEVOY_OPERATOR
	ZuckerLevoy3D(&m_MagBuf);
#endif
#ifdef	SOBEL_LEVOY_OPERATOR
	SobelLevoy3D(&m_MagBuf);
#endif

	// Allocate temporary memory
	RxAutoSample3DInfo	*temp;
	int					count=0, increment=0;
	temp = new RxAutoSample3DInfo[5000000];
	for (int z = 0; z < m_pResol->Rw; z++) {
		for (int y = 0; y < m_pResol->Rv; y++) {
			for (int x = 0; x < m_pResol->Ru; x++) {
				if(m_MagBuf[((z*m_pResol->Rv)+y)*m_pResol->Ru+x] == 255) {
					temp[count].x = x;	temp[count].y = y;
					temp[count].z = z;	temp[count].val = 0;
					count++;
				} // end of if
			} // end of x for
		} // end of y for
	} // end of z for
	ReleaseVolumeData(m_MagBuf,*m_pResol);

	// Allocate sample points A memory
	*sample = new RxAutoSample3DInfo[m_SampleNum];
	if(count <= m_SampleNum) {
		increment = 1;	m_SampleNum = count;
	}
	else {
		increment = count / m_SampleNum;
	}

	int i, j;
	for (i = 0, j = 0; j < m_SampleNum; i += increment, j++) {
		(*sample)[j].x	= temp[i].x;
		(*sample)[j].y	= temp[i].y;
		(*sample)[j].z	= temp[i].z;
		(*sample)[j].val = VoxelOffset(temp[i].x,temp[i].y,temp[i].z);
	}
	delete []temp;

	return TRUE;
}

//*********************************************************************
//	Zucker3D(): Edge detection with zucker 3D operator
//*********************************************************************
BOOL RxAutoSample3D::Zucker3D(BYTE **sobel)
{
	// Convolve Sx, Sy, Sz mask with original volume dataset
	//
	// -0.577  -0.707  -0.577	  0  0  0	  0.577  0.707  0.577
	// -0.707  -1.0    -0.707     0  0  0     0.707  1.0    0.707
	// -0.577  -0.707  -0.577     0  0  0     0.577  0.707  0.577
	//
	for (int z = 1; z < m_pResol->Rw-1; z++) {
		for (int y = 1; y < m_pResol->Rv-1; y++) {
			for (int x = 1; x < m_pResol->Ru-1; x++) {
				int gradx = (int)((0.577*(VoxelOffset(x+1,y-1,z-1)+VoxelOffset(x+1,y+1,z-1)+VoxelOffset(x+1,y-1,z+1)+
					                      VoxelOffset(x+1,y+1,z+1)-VoxelOffset(x-1,y-1,z-1)-VoxelOffset(x-1,y+1,z-1)-
										  VoxelOffset(x-1,y-1,z+1)-VoxelOffset(x-1,y+1,z+1)))+
								  (0.707*(VoxelOffset(x+1,y-1,z)+VoxelOffset(x+1,y+1,z)+VoxelOffset(x+1,y,z-1)+
								          VoxelOffset(x+1,y,z+1)-VoxelOffset(x-1,y-1,z)-VoxelOffset(x-1,y+1,z)-
										  VoxelOffset(x-1,y,z-1)-VoxelOffset(x-1,y,z+1)))+
										 (VoxelOffset(x+1,y,z)-VoxelOffset(x-1,y,z)));
				int grady = (int)((0.577*(VoxelOffset(x-1,y+1,z-1)+VoxelOffset(x+1,y+1,z-1)+VoxelOffset(x-1,y-1,z+1)+
					                      VoxelOffset(x+1,y+1,z+1)-VoxelOffset(x-1,y-1,z-1)-VoxelOffset(x+1,y-1,z-1)-
										  VoxelOffset(x-1,y-1,z+1)-VoxelOffset(x+1,y-1,z+1)))+
								  (0.707*(VoxelOffset(x-1,y+1,z)+VoxelOffset(x+1,y+1,z)+VoxelOffset(x,y+1,z-1)+
								          VoxelOffset(x,y+1,z+1)-VoxelOffset(x-1,y-1,z)-VoxelOffset(x+1,y-1,z)-
										  VoxelOffset(x,y-1,z-1)-VoxelOffset(x,y-1,z+1)))+
										 (VoxelOffset(x,y+1,z)-VoxelOffset(x,y-1,z)));
				int	gradz = (int)((0.577*(VoxelOffset(x-1,y-1,z+1)+VoxelOffset(x+1,y-1,z+1)+VoxelOffset(x-1,y+1,z+1)+
					                      VoxelOffset(x+1,y+1,z+1)-VoxelOffset(x-1,y-1,z-1)-VoxelOffset(x+1,y-1,z-1)-
										  VoxelOffset(x-1,y+1,z-1)-VoxelOffset(x+1,y+1,z-1)))+
								  (0.707*(VoxelOffset(x-1,y,z-1)+VoxelOffset(x+1,y,z-1)+VoxelOffset(x,y-1,z-1)+
								          VoxelOffset(x,y+1,z-1)-VoxelOffset(x-1,y,z-1)-VoxelOffset(x+1,y,z-1)-
										  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));
				if(absgrad > m_Threshold)
					(*sobel)[((z*m_pResol->Rv)+y)*m_pResol->Ru+x] = 255;
			} // end of x for
		} // end of y for
	} // end of z for

	return TRUE;
}

//*********************************************************************
//	ZuckerLevoy3D(): Edge detection with zucker 3D operator and levoy filtering
//*********************************************************************
BOOL RxAutoSample3D::ZuckerLevoy3D(BYTE **levoy)
{
	unsigned short	val;
	double			opacity = 0.0;

	// Convolve Sx, Sy, Sz mask with original volume dataset
	//
	// -0.577  -0.707  -0.577	  0  0  0	  0.577  0.707  0.577
	// -0.707  -1.0    -0.707     0  0  0     0.707  1.0    0.707
	// -0.577  -0.707  -0.577     0  0  0     0.577  0.707  0.577
	//
	for (int z = 1; z < m_pResol->Rw-1; z++) {
		for (int y = 1; y < m_pResol->Rv-1; y++) {
			for (int x = 1; x < m_pResol->Ru-1; x++) {
				if((val=VoxelOffset(x,y,z)) < m_Threshold)
					opacity = 0.0;

				int gradx = (int)((0.577*(VoxelOffset(x+1,y-1,z-1)+VoxelOffset(x+1,y+1,z-1)+VoxelOffset(x+1,y-1,z+1)+
					                      VoxelOffset(x+1,y+1,z+1)-VoxelOffset(x-1,y-1,z-1)-VoxelOffset(x-1,y+1,z-1)-
										  VoxelOffset(x-1,y-1,z+1)-VoxelOffset(x-1,y-1,z+1)))+
								  (0.707*(VoxelOffset(x+1,y-1,z)+VoxelOffset(x+1,y+1,z)+VoxelOffset(x+1,y,z-1)+
								          VoxelOffset(x+1,y,z+1)-VoxelOffset(x-1,y-1,z)-VoxelOffset(x-1,y+1,z)-
										  VoxelOffset(x-1,y,z-1)-VoxelOffset(x-1,y,z+1)))+
										 (VoxelOffset(x+1,y,z)-VoxelOffset(x-1,y,z)));

				int grady = (int)((0.577*(VoxelOffset(x-1,y+1,z-1)+VoxelOffset(x+1,y+1,z-1)+VoxelOffset(x-1,y-1,z+1)+
					                      VoxelOffset(x+1,y+1,z+1)-VoxelOffset(x-1,y-1,z-1)-VoxelOffset(x+1,y-1,z-1)-
										  VoxelOffset(x-1,y-1,z+1)-VoxelOffset(x+1,y-1,z+1)))+
								  (0.707*(VoxelOffset(x-1,y+1,z)+VoxelOffset(x+1,y+1,z)+VoxelOffset(x,y+1,z-1)+
								          VoxelOffset(x,y+1,z+1)-VoxelOffset(x-1,y-1,z)-VoxelOffset(x+1,y-1,z)-
										  VoxelOffset(x,y-1,z-1)-VoxelOffset(x,y-1,z+1)))+
										 (VoxelOffset(x,y+1,z)-VoxelOffset(x,y-1,z)));

				int	gradz = (int)((0.577*(VoxelOffset(x-1,y-1,z+1)+VoxelOffset(x+1,y-1,z+1)+VoxelOffset(x-1,y+1,z+1)+
					                      VoxelOffset(x+1,y+1,z+1)-VoxelOffset(x-1,y-1,z-1)-VoxelOffset(x+1,y-1,z-1)-
										  VoxelOffset(x-1,y+1,z-1)-VoxelOffset(x+1,y+1,z-1)))+
								  (0.707*(VoxelOffset(x-1,y,z-1)+VoxelOffset(x+1,y,z-1)+VoxelOffset(x,y-1,z-1)+
								          VoxelOffset(x,y+1,z-1)-VoxelOffset(x-1,y,z-1)-VoxelOffset(x+1,y,z-1)-