voxelshader.cpp

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

// VoxelShader.cpp : implementation file
//
//////////////////////////////////////////////////////////////////////
//
//	Title: This class manage shading voxel
//
//////////////////////////////////////////////////////////////////////
//
//	Author: H.W.Kye
//	138-dong 417-ho Seoul National University
//	San 56-1 Shinlim-dong Kwanak-gu Seoul, Korea
//	Email. 
//
//	Date	:
//	Update	:
//
//////////////////////////////////////////////////////////////////////

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

#include "StdAfx.h"
#include <xmmintrin.h>
#include <math.h>
#include "Matrix4D.h"
#include "VoxelShader.h"

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

static const int NORM_BITS=13;	/*   number of bits to encode a normal:(1+ceil(log2(2*(N*N+N+1)))) */
static const float STANDARDSAMPLES_IN_MM=1;
static const int MATL_SHINY=10;

//////////////////////////////////////////////////////////////////////
// RxVoxelShader

RxVoxelShader::RxVoxelShader()
{
	int iSseArraySize = ((MAX_NORMAL)+3)/4;

	m_iVoxX = m_iVoxY = m_iVoxZ = 0;
	m_iVOIXm = m_iVOIXM = 0;
	m_iVOIYm = m_iVOIYM = 0;
	m_iVOIZm = m_iVOIZM = 0;
	m_iBright = 30;
	m_iContrast = 30;
	m_iRatioZ = 1;
	m_fSamplingDistance = 1.0f;
	InitNormalTables();
	SetShadingMode(TRUE);


	m_aByte2M128Table = (__m128*)_mm_malloc(256*16, 16);
	for(int i=0; i<256; i++) {
		float f = (float) i / 255.0f;
		m_aByte2M128Table[i] = _mm_set_ps(f,f,f,1.0f); // 父甸扁
	}
}

RxVoxelShader::~RxVoxelShader()
{
	_mm_free(m_aByte2M128Table);
	m_aByte2M128Table = NULL;
}


//////////////////////////////////////////////////////////////////////
// Initialize

BOOL RxVoxelShader::InitNormalTables()
{
	// initialize NormalPy
	int	iXcount = 1;
	int	iCodePoint = 2;
	m_pnNormalPy = &m_anNormalPyStorage[NORM_N];
	m_pnNormalXLimit = &m_anNormalXLimitStorage[NORM_N];

	for (int y = -NORM_N; y <= NORM_N; y++) {
		m_pnNormalPy[y] = short(iCodePoint + (((iXcount - 1) / 2) << 1));
		iCodePoint += (iXcount << 1);
		m_pnNormalXLimit[y] = short((iXcount - 1) / 2);
		iXcount += (y < 0)? 2 : -2;
	}

	return TRUE;
}


//////////////////////////////////////////////////////////////////////
// shading

int RxVoxelShader::GetNormalIndexFromPosition(int iVoxelOffset)
{
	int dx,dy,dz;
	dx = int(m_pnOriginalVolume[iVoxelOffset+m_iOffsetX]) - int(m_pnOriginalVolume[iVoxelOffset-m_iOffsetX]);
	dx *= m_iInverseVoxelSizeX;

	dy = int(m_pnOriginalVolume[iVoxelOffset+m_iOffsetY]) - int(m_pnOriginalVolume[iVoxelOffset-m_iOffsetY]);
	dy *= m_iInverseVoxelSizeY;

	dz = int(m_pnOriginalVolume[iVoxelOffset+m_iOffsetZ]) - int(m_pnOriginalVolume[iVoxelOffset-m_iOffsetZ]);
	dz *= m_iInverseVoxelSizeZ;

	short normal;
	if ( dx==0 && dy==0 && dz==0) {
		normal = 0;
	} else {
		normal = GetNormalIndexFromNormalVector(dx,dy,dz);
	}

	return normal;
}

int RxVoxelShader::GetNormalIndexFromPosition(int iVoxelOffset, int iInternalNum)
{
	int dx1,dx2;
	int dy1,dy2;

	int dx,dy,dz;
	dx1 = int(m_pnOriginalVolume[iVoxelOffset+m_iOffsetX]) - int(m_pnOriginalVolume[iVoxelOffset-m_iOffsetX]);
	dx2 = int(m_pnOriginalVolume[iVoxelOffset+m_iOffsetZ+m_iOffsetX]) - int(m_pnOriginalVolume[iVoxelOffset+m_iOffsetZ-m_iOffsetX]);

	dx = iInternalNum*dx2 + ((1 << PIXEDPOINT_ACC) - iInternalNum)*dx1;
	dx = dx >> PIXEDPOINT_ACC;
	dx *= m_iInverseVoxelSizeX;

	dy1 = int(m_pnOriginalVolume[iVoxelOffset+m_iOffsetY]) - int(m_pnOriginalVolume[iVoxelOffset-m_iOffsetY]);
	dy2 = int(m_pnOriginalVolume[iVoxelOffset+m_iOffsetZ+m_iOffsetY]) - int(m_pnOriginalVolume[iVoxelOffset+m_iOffsetZ-m_iOffsetY]);

	dy = iInternalNum*dy2 + ((1 << PIXEDPOINT_ACC) - iInternalNum)*dy1;
	dy = dy >> PIXEDPOINT_ACC;
	dy *= m_iInverseVoxelSizeY;

	dz = int(m_pnOriginalVolume[iVoxelOffset+m_iOffsetZ]) - int(m_pnOriginalVolume[iVoxelOffset]);
	dz += dz;
	dz *= m_iInverseVoxelSizeZ;

	short normal;
	if ( dx==0 && dy==0 && dz==0) {
		normal = 0;
	} else {
		normal = GetNormalIndexFromNormalVector(dx,dy,dz);
	}

	return normal;
}


int	RxVoxelShader::GetGradientMagnitude(int nOffset)
{
	int volume = m_iVoxX*m_iVoxY*m_iVoxZ-1;

	int X2 = (nOffset+m_iOffsetX < volume)?nOffset+m_iOffsetX:nOffset;
	int X1 = (nOffset-m_iOffsetX >= 0)?	nOffset-m_iOffsetX:nOffset;

	int Y2 = (nOffset+m_iOffsetY < volume)?nOffset+m_iOffsetY:nOffset;
	int Y1 = (nOffset-m_iOffsetY >= 0)?	nOffset-m_iOffsetY:nOffset;

	int Z2 = (nOffset+m_iOffsetZ < volume)?nOffset+m_iOffsetZ:nOffset;
	int Z1 = (nOffset-m_iOffsetZ >= 0)?	nOffset-m_iOffsetZ:nOffset;


	int dx = m_pnOriginalVolume[X2] - m_pnOriginalVolume[X1];
	int dy = m_pnOriginalVolume[Y2] - m_pnOriginalVolume[Y1];
	int dz = m_pnOriginalVolume[Z2] - m_pnOriginalVolume[Z1];

	float sqr = (float)(dx*dx + dy*dy + dz*dz);
	int igm = _mm_cvttss_si32(_mm_sqrt_ss(_mm_load_ss(&sqr)));

	return igm>>1;
}



//////////////////////////////////////////////////////////////////////
//	Interface functions

BOOL RxVoxelShader::SetVolumeParameters(unsigned short *p3DVolume, unsigned short *pnOriginalVolume,
						const int iVolumeSizeX, const int iVolumeSizeY, const int iVolumeSizeZ,
						const double fVoxelSizeX, const double fVoxelSizeY, const double fVoxelSizeZ, int iRatioZ)
{
	// p3DVolume data
	m_pVolume = p3DVolume;
	m_pnOriginalVolume = pnOriginalVolume;

	m_iVoxX = iVolumeSizeX;
	m_iVoxY = iVolumeSizeY;
	m_iVoxZ = iVolumeSizeZ;

	m_iOffsetX = 1;
	m_iOffsetY = m_iVoxX;
	m_iOffsetZ = m_iVoxX * m_iVoxY;

	m_fVoxelSizeX = fVoxelSizeX;
	m_fVoxelSizeY = fVoxelSizeY;
	m_fVoxelSizeZ = fVoxelSizeZ;
	m_iRatioZ = iRatioZ;

	m_iInverseVoxelSizeX = int(1000 * fVoxelSizeX / m_fVoxelSizeX);	// prevent variables from overflow
	m_iInverseVoxelSizeY = int(1000 * fVoxelSizeX / m_fVoxelSizeY);
	m_iInverseVoxelSizeZ = int(1000 * fVoxelSizeX / m_fVoxelSizeZ);

	return TRUE;
}


BOOL RxVoxelShader::SetViewing(const RxMatrix4D &mtrxView)
{
	// opacity correction
	m_mtxView = mtrxView;
	RxVect4D vecView(m_mtxView.ViewVectX(), m_mtxView.ViewVectY(), m_mtxView.ViewVectZ(), 1.0);
	vecView.Normalize();


	double fViewingX = fabs(vecView.m[0]);
	double fViewingY = fabs(vecView.m[1]);
	double fViewingZ = fabs(vecView.m[2]);

	double dKx = (fViewingX==0)? 10000000.0 : m_fVoxelSizeX/fViewingX;
	double dKy = (fViewingY==0)? 10000000.0 : m_fVoxelSizeY/fViewingY;
	double dKz = (fViewingZ==0)? 10000000.0 : m_fVoxelSizeZ/m_iRatioZ/fViewingZ;


	m_fSamplingDistance = (float)__min(__min(dKx, dKy), dKz);

	RxDensityTable::InitRendering(m_fSamplingDistance * STANDARDSAMPLES_IN_MM);
	MakeShadeTable();
	return TRUE;
}

BOOL RxVoxelShader::SetBrightContrast(int iBright, int iContrast)
{
	ASSERT(iBright >= 0 && iBright <= 100);
	ASSERT(iContrast >= 0 && iContrast <= 100);

	m_iBright = iBright;
	m_iContrast = iContrast;

	MakeShadeTable();

	return TRUE;
}


BOOL RxVoxelShader::SetWindowing(unsigned char *pShadedWindowing, unsigned char *pUnshadedWindowing, unsigned char *pLocalWindowing)
{
	m_pbyShadedWindowing = pShadedWindowing;
	m_pbyUnshadedWindowing = pUnshadedWindowing;
	m_pLocalWindowing = pLocalWindowing;

	m_aiSATUnshadedWindowing[0] = 0;
	m_aiSATShadedWindowing[0] = 0;

	for(int i=0; i<DENSITY_RANGE; i++) {
		m_aiSATUnshadedWindowing[i+1] = m_aiSATUnshadedWindowing[i] + m_pbyUnshadedWindowing[i];
		m_aiSATShadedWindowing[i+1] = m_aiSATShadedWindowing[i] + m_pbyShadedWindowing[i];
	}

	return TRUE;
}

BOOL RxVoxelShader::SetWindowingTable(unsigned char *pShadedWindowing, unsigned char *pUnshadedWindowing)
{
	m_pbyShadedWindowing = pShadedWindowing;
	m_pbyUnshadedWindowing = pUnshadedWindowing;

	m_aiSATUnshadedWindowing[0] = 0;
	m_aiSATShadedWindowing[0] = 0;

	for(int i=0; i<DENSITY_RANGE; i++) {
		m_aiSATUnshadedWindowing[i+1] = m_aiSATUnshadedWindowing[i] + m_pbyUnshadedWindowing[i];
		m_aiSATShadedWindowing[i+1] = m_aiSATShadedWindowing[i] + m_pbyShadedWindowing[i];
	}

	return TRUE;
}

BOOL RxVoxelShader::Set3DVOI(const int iVOIXm, const int iVOIXM, const int iVOIYm, const int iVOIYM, const int iVOIZm, const int iVOIZM)
{
	m_iVOIXm = iVOIXm;
	m_iVOIXM = iVOIXM;

	m_iVOIYm = iVOIYm;
	m_iVOIYM = iVOIYM;

	m_iVOIZm = iVOIZm, 
	m_iVOIZM = iVOIZM;

	return TRUE;
}



// table order
// ambient, diffuse, specular, 1
BOOL RxVoxelShader::MakeShadeTable()
{