vrraycastingagent.cpp

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

// VRRayCastingAgent.cpp : implementation file
//
//////////////////////////////////////////////////////////////////////
//
//	Title: Main class to render volume
//
//////////////////////////////////////////////////////////////////////
//
//	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 <math.h>
#include "matrix4d.h"
#include "FastInlineFuncs.h"
#include "BlockData.h"
#include "VRUtil.h"
#include "VoxelShader.h"
#include "VRRayCastingAgent.h"
#include "RayCastingCommonInfo.h"
#include "list.h"
#include "evaltime.h"
#include "VolumeData.h"
#include "fusionglobal.h"

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

#define AMBIENT_BASE 0.2f
#define DIFFUSE_BASE 0.7f
#define SPECULAR_BASE 0.1f
#define SAMPLING_RATE 1.0f
const static int INIT_VALUE = -1;
extern HANDLE g_hVRThreadKillEvent;
#define SETZERO(buf) (memset(buf, 0, sizeof(buf)))

//////////////////////////////////////////////////////////////////////
// RxVRRayCastingAgent

RxVRRayCastingAgent::RxVRRayCastingAgent(int iLSN) : RxRayCastingCommonAgent(iLSN)
{
	SETZERO(m_afLight);
	SETZERO(m_afRayDirection);
	SETZERO(m_afInvRayDirection);
	SETZERO(m_afRefRayDirection);
	SETZERO(m_afRefInvRayDirection);
	SETZERO(m_fLocalizerPos);
	SETZERO(m_iPlaneOrder);
	SETZERO(m_afCropMin);
	SETZERO(m_afCropMax);
	SETZERO(m_afByteTable);
	SETZERO(m_afRayTable);
	SETZERO(m_fLocalizerPos);
	SETZERO(m_fLocalizerPos);
	SETZERO(m_aiCropMin);
	SETZERO(m_aiCropMax);
	SETZERO(m_iLocalizerPlane);

	m_fWMin = m_fWMax = 0;
	m_fInvGradFactor = 0;
	m_fGradFactor = 0;
	m_pShader = NULL;
	m_fA = m_fB = m_fC = m_fD1 = m_fD2 = m_fDenom = 0;
	m_iMin = 0;
	m_iMax = 0;
	m_pObject2x2 = m_pObject8x8 = NULL;
	m_fBrightnessFactor = 0;
	m_fSampleDistance = 0;
}

RxVRRayCastingAgent::~RxVRRayCastingAgent()
{
}

BOOL RxVRRayCastingAgent::InitRendering(RxRayCastingCommonInfo* pIsoEndoInfo)
{
	m_pIsoEndoInfo = pIsoEndoInfo;

	m_pnVolume = pIsoEndoInfo->m_pnVolume;
	// info肺 哗具 窃. RxGet...葛滴.
	m_piSliceTable = pIsoEndoInfo->m_piSliceTable; // RxGetVolumeData(m_iLoadedSeriesNum)->m_piInterpolatedMap;

	m_iVolX = pIsoEndoInfo->m_iVolX;
	m_iVolY = pIsoEndoInfo->m_iVolY;
	m_iVolZ = pIsoEndoInfo->m_iVolZ;

	m_afMaxVol[0] = IntToFloat(m_iVolX-1);
	m_afMaxVol[1] = IntToFloat(m_iVolY-1);
	m_afMaxVol[2] = IntToFloat(m_iVolZ-1);

	RxBlockData* pBlock = (RxBlockData*)pIsoEndoInfo->m_pBlockData;

	m_pmmDen = pBlock->m_pmmBuffer;
	m_pmmBlock = pBlock->m_pmmBuffer8x8;

	m_piIdx[0] = pBlock->m_piIdx[0];
	m_piIdx[1] = pBlock->m_piIdx[1];
	m_piIdx[2] = pBlock->m_piIdx[2];

	m_piBlockIdx[0] = pBlock->m_piBlockIdx[0];
	m_piBlockIdx[1] = pBlock->m_piBlockIdx[1];
	m_piBlockIdx[2] = pBlock->m_piBlockIdx[2];

	m_piSkipIdx[0] = pBlock->m_piSkipIdx[0];
	m_piSkipIdx[1] = pBlock->m_piSkipIdx[1];
	m_piSkipIdx[2] = pBlock->m_piSkipIdx[2];

	m_bCreated = TRUE;

	if(pIsoEndoInfo->m_szImage.cx == 0 || pIsoEndoInfo->m_szImage.cy == 0)
		return FALSE;

	if (m_iImgX!=pIsoEndoInfo->m_szImage.cx || m_iImgY!=pIsoEndoInfo->m_szImage.cy) {
		m_iImgX = pIsoEndoInfo->m_szImage.cx;
		m_iImgY = pIsoEndoInfo->m_szImage.cy;

		TRACE(_T("\nimage deleted"));
		if (m_pImage)
			delete[] m_pImage;
		m_pImage = new unsigned char [m_iImgX*m_iImgY*4];
		if (m_pImageOnce)
			delete[] m_pImageOnce;
		m_pImageOnce = new unsigned char [((m_iImgX+1)/2)*((m_iImgY+1)/2)*4];
		if (m_pImageTwice)
			delete[] m_pImageTwice;
		m_pImageTwice = new unsigned char [((m_iImgX+3)/4)*((m_iImgY+3)/4)*4];
	}

	for (int i=0;i<4;i++) {
		for (int j=0;j<4;j++) {
			m_afViewingM[i][j] = float(pIsoEndoInfo->m_mxViewing.m[i][j]);
			m_afInvViewingM[i][j] = float(pIsoEndoInfo->m_mxInvViewing.m[i][j]);
		}
	}
	m_fRatio = pIsoEndoInfo->m_fRatioZ;

	m_fAmbient = AMBIENT_BASE-IntToFloat(pIsoEndoInfo->m_iContrast-30)*0.01f*0.8f;
	if (m_fAmbient < 0)
		m_fAmbient = 0;
	m_fDiffuse = DIFFUSE_BASE+IntToFloat(pIsoEndoInfo->m_iContrast-30)*0.01f;

	float fBright = IntToFloat(pIsoEndoInfo->m_iBrightness-30)*0.02f;
	m_fBrightnessFactor = 1.f + fBright*fBright*fBright;

	m_pShader = m_pIsoEndoInfo->GetShader();
	m_pShader->InitRendering(SAMPLING_RATE);

	if(m_pIsoEndoInfo->IsRef()) {
		m_afCropMin[0] = IntToFloat(m_pIsoEndoInfo->m_iCrop_XMin);
		m_afCropMin[1] = IntToFloat(m_pIsoEndoInfo->m_iCrop_YMin);
		m_afCropMin[2] = IntToFloat(m_pIsoEndoInfo->m_iCrop_ZMin);

		m_afCropMax[0] = IntToFloat(m_pIsoEndoInfo->m_iCrop_XMax);
		m_afCropMax[1] = IntToFloat(m_pIsoEndoInfo->m_iCrop_YMax);
		m_afCropMax[2] = IntToFloat(m_pIsoEndoInfo->m_iCrop_ZMax);

		m_aiCropMin[0] = m_pIsoEndoInfo->m_iCrop_XMin;
		m_aiCropMin[1] = m_pIsoEndoInfo->m_iCrop_YMin;
		m_aiCropMin[2] = m_pIsoEndoInfo->m_iCrop_ZMin;

		m_aiCropMax[0] = m_pIsoEndoInfo->m_iCrop_XMax;
		m_aiCropMax[1] = m_pIsoEndoInfo->m_iCrop_YMax;
		m_aiCropMax[2] = m_pIsoEndoInfo->m_iCrop_ZMax;
	} else {

		for(i=0; i<3; i++) {
			m_afCropMin[i] = m_pIsoEndoInfo->m_afBoundingMin[i];
			m_afCropMax[i] = m_pIsoEndoInfo->m_afBoundingMax[i];

			m_aiCropMin[i] = FloatToInt(m_afCropMin[i]);
			m_aiCropMax[i] = FloatToInt(m_afCropMax[i]);
		}
	}


	m_fGradFactor = IntToFloat(m_pIsoEndoInfo->m_iGradientMag);
	m_fInvGradFactor = 1.f/(m_fGradFactor*m_fGradFactor);

	for (i=0;i<256;i++)
		m_afByteTable[i] = IntToFloat(i)/255.f;


	float afTemp[4] = {0, 0, 1, 0};
	m4v4mulf(m_afRayDirection, m_afInvViewingM, afTemp);

	// | m_afRayDirection | != 1
	m_afRayDirection[0] /= m_pIsoEndoInfo->m_fRayMagnitude;
	m_afRayDirection[1] /= m_pIsoEndoInfo->m_fRayMagnitude;
	m_afRayDirection[2] /= m_pIsoEndoInfo->m_fRayMagnitude;
	m_afRayDirection[3] /= m_pIsoEndoInfo->m_fRayMagnitude;

	for(i=0; i<4; i++) {
		m_afRayDirection[i] *= SAMPLING_RATE;
	}

	// | m_afLight | == 1
	memcpy(m_afLight, m_afRayDirection, 12);
	v3nrmf(m_afLight);

	m_afInvRayDirection[0] = 1.f/m_afRayDirection[0];
	m_afInvRayDirection[1] = 1.f/m_afRayDirection[1];
	m_afInvRayDirection[2] = 1.f/m_afRayDirection[2];

	m_fWMin = IntToFloat(m_pIsoEndoInfo->m_iWMin);
	m_fWMax = IntToFloat(m_pIsoEndoInfo->m_iWMax);
	m_hThreadKillEvent = pIsoEndoInfo->m_hThreadKillEvent;

	// m_fSampleDistance * | m_afRayDirection | == 1
	m_fSampleDistance = 1.f/(float)(sqrt(m_afRayDirection[0]*m_afRayDirection[0]
		+ m_afRayDirection[1]*m_afRayDirection[1]
		+ m_afRayDirection[2]*m_afRayDirection[2]));

//	m_fSampleDistance = 1;

	TRACE(_T("\nInit Rendering... volume size=%d"),m_iVolX);
	return TRUE;
}

#define SET_MASK_OFFSET(iOFFSET) {iOffset=iOFFSET; iQuotient = iOffset >> 3; iRemainder = iOffset & 0x07;}
#define GET_MASK(pbyMASK) ((pbyMASK[iQuotient] >> (7 - iRemainder)) & 0x01)
#define SET_MASK(pbyMASK) pbyMASK[iQuotient] = pbyMASK[iQuotient] | (0x01 << (7 - iRemainder))

BOOL RxVRRayCastingAgent::IsInsideInCrop(int x, int y, int z)
{
	if(
		(m_aiCropMin[0] < x && x < m_aiCropMax[0]) &&
		(m_aiCropMin[1] < y && y < m_aiCropMax[1]) &&
		(m_aiCropMin[2] < z && z < m_aiCropMax[2])
		)
		return TRUE;
	else
		return FALSE;

}
int RxVRRayCastingAgent::Growing3D(int x, int y, int z)
{
	int iOffset, iQuotient, iRemainder;
	int nVoxelNum=0;

	ListItem liVoxel;
	RxList cList(65536);
	cList.ClearAll();
	LPBYTE pMask=(LPBYTE)VirtualAlloc(NULL, m_iVolX*m_iVolY*m_iVolZ/8,MEM_RESERVE | MEM_TOP_DOWN | MEM_COMMIT, PAGE_READWRITE);
	ASSERT(pMask);
	memset(pMask, 0, m_iVolX*m_iVolY*m_iVolZ/8);

	memset(pMask+m_iVolX*m_iVolY*(m_aiCropMin[2])/8, 0xff, m_iVolX*m_iVolY/8);
	memset(pMask+m_iVolX*m_iVolY*(m_aiCropMax[2])/8, 0xff, m_iVolX*m_iVolY/8);
	for(int iSlice=m_aiCropMin[2]+1; iSlice<m_aiCropMax[2]; iSlice++) {
		memset(pMask+m_iVolX*m_iVolY*iSlice/8 + m_iVolX*m_aiCropMin[1]/8, 0xff, m_iVolX/8);
		memset(pMask+m_iVolX*m_iVolY*iSlice/8 + m_iVolX*m_aiCropMax[1]/8, 0xff, m_iVolX/8);
	}

	cList.AddTail(x,y,z);
	while(!cList.IsEmpty()) {
		liVoxel = cList.RemoveHead();

		int iFrom, iTo, iStep;
		for(int c=0; c<2; c++) {
			if(c==0) {
				iFrom=liVoxel.x;
				iTo=m_aiCropMin[0];
				iStep=-1;
			} else {
				iFrom=liVoxel.x+1;
				iTo=m_aiCropMax[0];
				iStep=1;
			}

			BOOL bUpQ=FALSE, bDownQ=FALSE, bLeftQ=FALSE, bRightQ=FALSE;
			for(int i=iFrom; i!=iTo; i+=iStep) {


				SET_MASK_OFFSET(m_piIdx[0][i] + m_piIdx[1][liVoxel.y] + m_piIdx[2][liVoxel.z]);
				if(!m_pShader->IsTransparent(m_pnVolume[iOffset]) && !GET_MASK(pMask)) {
					SET_MASK(pMask);
					nVoxelNum++;
				} else
					break;

				SET_MASK_OFFSET(m_piIdx[0][i] + m_piIdx[1][liVoxel.y] + m_piIdx[2][liVoxel.z+1]);
				if(!m_pShader->IsTransparent(m_pnVolume[iOffset]) && !GET_MASK(pMask)) {
					if(!bUpQ) {
						bUpQ = TRUE;
						cList.AddTail(i, liVoxel.y, liVoxel.z + 1);
					}
				} else
					bUpQ = FALSE;

				// bDownQ (z - 1)
				SET_MASK_OFFSET(m_piIdx[0][i] + m_piIdx[1][liVoxel.y] + m_piIdx[2][liVoxel.z-1]);
				if(!m_pShader->IsTransparent(m_pnVolume[iOffset]) && !GET_MASK(pMask)) {
					if (!bDownQ) {
						bDownQ = TRUE;
						cList.AddTail(i, liVoxel.y, liVoxel.z - 1);
					}
				} else
					bDownQ = FALSE;

				// bLeftQ (y - 1)
				SET_MASK_OFFSET(m_piIdx[0][i] + m_piIdx[1][liVoxel.y-1] + m_piIdx[2][liVoxel.z]);
				if(!m_pShader->IsTransparent(m_pnVolume[iOffset]) && !GET_MASK(pMask)) {
					if (!bLeftQ) {
						bLeftQ = TRUE;
						cList.AddTail(i, liVoxel.y - 1, liVoxel.z);
					}
				} else
					bLeftQ = FALSE;

				// bRightQ (y + 1)
				SET_MASK_OFFSET(m_piIdx[0][i] + m_piIdx[1][liVoxel.y+1] + m_piIdx[2][liVoxel.z]);
				if(!m_pShader->IsTransparent(m_pnVolume[iOffset]) && !GET_MASK(pMask)) {
					if (!bRightQ) {
						bRightQ = TRUE;
						cList.AddTail(i, liVoxel.y + 1, liVoxel.z);
					}
				} else
					bRightQ = FALSE;
				
			} // end of a span
		} // end of a seed
	}

	VirtualFree(pMask, m_iVolX*m_iVolY*m_iVolZ/8, MEM_DECOMMIT);
	VirtualFree(pMask, 0, MEM_RELEASE);
	return nVoxelNum;
}

BOOL RxVRRayCastingAgent::Get3DPoint(int i,int j, int &x, int &y, int &z)
{
	ASSERT(m_pShader);

	x=y=z=-1;
	float afRayOrigin[4];
	float d1[4] = {(float)i, (float)j, 0, 1.};
	m4v4mulf(afRayOrigin, m_afInvViewingM, d1);
	float fStep[4]={m_afRayDirection[0],m_afRayDirection[1],m_afRayDirection[2],m_afRayDirection[3]};
	v3nrmf(fStep);

	float min_t = -MAX_T;
	float max_t = MAX_T;
	int iType = GetRaySampleBound(0, afRayOrigin, &min_t, &max_t, NULL, NULL, NULL, NULL);
	if(iType == 0)
		return FALSE;

	for (float t = min_t+1; t < max_t; t += 1.) {
		int ipx = FloatToInt(afRayOrigin[0] + m_afRayDirection[0] * t);
		int ipy = FloatToInt(afRayOrigin[1] + m_afRayDirection[1] * t);
		int ipz = FloatToInt(afRayOrigin[2] + m_afRayDirection[2] * t);

		unsigned short nDensity = m_pnVolume[m_piIdx[0][ipx] + m_piIdx[1][ipy] + m_piIdx[2][ipz]];
		if(!m_pShader->IsTransparent(nDensity)) {
			x = ipx;
			y = ipy;
			z = ipz;
			return TRUE;;
		}
	}
	return FALSE;
}


BOOL RxVRRayCastingAgent::PartialRender(int iPart)
{
	TRACE(_T("\nlevel %d starting..."),m_iLevel);
	BOOL bRetval = RenderNormal(iPart);
	if(bRetval) {
		TRACE(_T("\nlevel %d end..."),m_iLevel);
	} else {
		TRACE(_T("\nlevel %d interrupted..."),m_iLevel);
	}
	return bRetval;
}

BOOL RxVRRayCastingAgent::RenderNormal(int iPart)
{
	RxEvalTime timer;

	register int ix, iy;
	register float t;
	float min_t, max_t, limit_t, blk_last_t = -1;
	float opacity, cur_intensity;
	float r, g, b;
	int iSample;
	float blk_sample_step = 4.f; // * m_fSampleDistance;	// t俊 措茄 概俺函荐扁霖
	int ipx, ipy, ipz, iCurLocalizer;
	float afVoxel[3], afRayOrigin[4], rv[3];

	memcpy(rv, m_afRayDirection, 12);
	int aiBlockAdjust[3];
	aiBlockAdjust[0] = (m_afRayDirection[0] > 0)?1:0;
	aiBlockAdjust[1] = (m_afRayDirection[1] > 0)?1:0;
	aiBlockAdjust[2] = (m_afRayDirection[2] > 0)?1:0;

	float afRefInvViewing[4][4];
	GetRayDirectionInRef(afRefInvViewing);

	RxRayCastingCommonInfo* pIsoEndoInfo = m_pIsoEndoInfo;
	COLORREF rgbLocalizer = 0; //RxGetOptionInfo()->Get3DDisplayLocalPlaneColor();

	unsigned char *pImage;
	switch (m_iLevel) {
	case 0 :
		pImage = m_pImage;
		break;
	case 1 :
		pImage = m_pImageOnce;
		break;
	case 2 :
		pImage = m_pImageTwice;
		break;
	}
	unsigned char *pOrgImage = pImage;
	int iStartY, iEndY;

	int iMinX=0, iMinY=0, iMaxX=m_iImgX-1, iMaxY=m_iImgY-1;
	GetRayImageBound(&iMinX, &iMinY, &iMaxX, &iMaxY);

	int iSkip = 1<<m_iLevel;
	iStartY = __max(iMinY, m_iImgY/(int)m_dwNumberOfCPU*(iPart-1));
	iEndY = __min(iMaxY, m_iImgY/(int)m_dwNumberOfCPU*(iPart));

#ifdef MY_DEBUG
	float tmpColor[4];
	for(int i=0; i<4096; i++) {
		m_pShader->GetColor(i, tmpColor);
		if(tmpColor[0] > 0)
			exit(0);
	}
#endif