voxelshader.cpp

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

	int px, py, pz;
	int shadeindex=0;
	double pxd, pyd, pzd, plength;
	int pxlimit;
	double nx, ny, nz;
	double n_dot_l_xy;
	double n_dot_l_z;
	double n_dot_l;
	double n_dot_h_xy;
	double n_dot_h_z;
	double n_dot_h;


	RxMatrix4D	mtxViewInverse = m_mtxView.Inverse();

	RxVect4D		initial_vec(0.0, 0.0, -1.0, 0.0);

	RxVect4D		obj_light = mtxViewInverse*initial_vec;
	RxVect4D		obj_viewpoint = mtxViewInverse*initial_vec;

	obj_light.m[2] *= m_fVoxelSizeZ / m_fVoxelSizeX;
	obj_viewpoint.m[2] *= m_fVoxelSizeZ / m_fVoxelSizeX;

	obj_light.Normalize();
	obj_viewpoint.Normalize();

	RxVect4D		obj_highlight = obj_light + obj_viewpoint;
	obj_highlight.Normalize();


	// default
	__m128 vDefaultIntensity = _mm_set_ps(0.2f, 0.7f, 0.1f, 1);
	//__m128 vDefaultIntensity = _mm_set_ps(0.3f, 0.4f, 0.6f, 1); wet-look
	//__m128 vDefaultIntensity = _mm_set_ps(0.4f, 0.6f, 0.0f, 1);


	// contrast
	float fContrast = 0.01f*(m_iContrast+20-50); // -0.5 ~ +0.5
	__m128 vContrast = _mm_set_ps(-fContrast*0.8f, fContrast, 0.0f, 0.0f);
	vDefaultIntensity = _mm_add_ps(vDefaultIntensity, vContrast);


	// brightness
	float fBright = 0.02f*(m_iBright+20-50);
	fBright = 1.0f + fBright * fBright * fBright;
	__m128 vBright = _mm_set_ps1(fBright);
	vDefaultIntensity = _mm_mul_ps(vDefaultIntensity, vBright);

	__m128 vCurrentIntensity = _mm_setzero_ps();


	// floating point version
	m_pfIntensity[0] = m_pfIntensity[1] = 0.2f;

	shadeindex=2;

	for ( py = -NORM_N; py <= NORM_N; py++ ) {
		pxlimit = (py < 0) ? (NORM_N+py) : (NORM_N-py);
		pz = -1;
		pxd = (double)(-pxlimit-1);
		pyd = (double)py;
		pzd = (double)-1;
		for ( px = -pxlimit; px <= pxlimit; px++ ) {
			pxd += 1.;
			if ( px <= 0 ) {
				pz++;
				pzd += 1.;
			} else {
				pz--;
				pzd -= 1.;
			}
			plength = 1./sqrt(pxd*pxd+pyd*pyd+pzd*pzd);

			nx = pxd*plength;
			ny = pyd*plength;
			nz = pzd*plength;

			/* store diffuse and specular terms */
			n_dot_l_xy = nx*obj_light.m[0] + ny*obj_light.m[1];
			n_dot_l_z  = nz*obj_light.m[2];

			n_dot_h_xy = nx*obj_highlight.m[0] + ny*obj_highlight.m[1];
			n_dot_h_z  = nz*obj_highlight.m[2];


			// plus z
			n_dot_l = n_dot_l_xy + n_dot_l_z;
			n_dot_h = n_dot_h_xy + n_dot_h_z;
/*
			if(n_dot_l > 0)
				n_dot_l = 0;
			else
				n_dot_l = -n_dot_l;

			if(n_dot_h > 0)
				n_dot_h = 0;
			else
				n_dot_h = -n_dot_h;
*/
			n_dot_l = fabs(n_dot_l);
			n_dot_h = fabs(n_dot_h);


			vCurrentIntensity = _mm_set_ps(1.0f, (float)n_dot_l, (float)pow(n_dot_h, MATL_SHINY), 1.0f);
			__m128 vResult = _mm_mul_ps(vCurrentIntensity, vDefaultIntensity);

			float fResult[4];
			_mm_storeu_ps(fResult, vResult);
			fResult[0] = fResult[1] + fResult[2] + fResult[3];
			if(fResult[0] < 0)
				fResult[0] = 0;


			// floating point version
			m_pfIntensity[shadeindex] = fResult[0];

			shadeindex++;			


			// minus z
			n_dot_l = n_dot_l_xy - n_dot_l_z;
			n_dot_h = n_dot_h_xy - n_dot_h_z;
/*
			if(n_dot_l > 0)
				n_dot_l = 0;
			else
				n_dot_l = -n_dot_l;

			if(n_dot_h > 0)
				n_dot_h = 0;
			else
				n_dot_h = -n_dot_h;
*/
			n_dot_l = fabs(n_dot_l);
			n_dot_h = fabs(n_dot_h);


			vCurrentIntensity = _mm_set_ps(1.0f, (float)n_dot_l, (float)pow(n_dot_h, MATL_SHINY), 1.0f);


			vResult = _mm_mul_ps(vCurrentIntensity, vDefaultIntensity);
			_mm_storeu_ps(fResult, vResult);
			fResult[0] = fResult[1] + fResult[2] + fResult[3];
			if(fResult[0] < 0)
				fResult[0] = 0;
			
			// floating point version
			m_pfIntensity[shadeindex] = fResult[0];


			shadeindex++;

		}
	}

	return TRUE;
}

//////////////////////////////////////////////////////////////////////
// get operator interface

float RxVoxelShader::GetAlpha(int iI, int iJ, int iK)
{
	int nOffset = ((iK*m_iVoxY) + iJ)*m_iVoxX + iI;
	unsigned short unDensity = m_pVolume[nOffset];
	__m128 vColor = GetColor(unDensity);

	float fRetval;
	_mm_store_ss(&fRetval, vColor);
	return fRetval;
}





//////////////////////////////////////////////////////////////////////
// Shading interface


__m128 RxVoxelShader::GetShadedVector(__m128 &vColor, int iOffset)
{
	int iShadeTableIndex = GetNormalIndexFromPosition(iOffset);
	__m128 vShade = _mm_set_ps1(m_pfIntensity[iShadeTableIndex]);
	__m128 vShadedColor = _mm_mul_ps(vColor, vShade);
	__m128 vResultColor = _mm_move_ss(vShadedColor, vColor);
	return vResultColor;
}

__m128 RxVoxelShader::GetShadedVector(__m128 &vColor, int iOffset, int iInternalNum)
{
	int iShadeTableIndex = GetNormalIndexFromPosition(iOffset, iInternalNum);
	__m128 vShade = _mm_set_ps1(m_pfIntensity[iShadeTableIndex]);
	__m128 vShadedColor = _mm_mul_ps(vColor, vShade);
	__m128 vResultColor = _mm_move_ss(vShadedColor, vColor);
	return vResultColor;
}


//////////////////////////////////////////////////////////////////////
// windowing

__m128 RxVoxelShader::GetRawWindowing(int iI, int iJ, int iK, float fAlpha)
{
	int nOffset = ((iK*m_iVoxY) + iJ)*m_iVoxX + iI;
	unsigned short unDensity = m_pVolume[nOffset];
	__m128 vAlpha = _mm_set1_ps(fAlpha);
	__m128 vWindow = m_aByte2M128Table[m_pbyUnshadedWindowing[unDensity]];

	return _mm_mul_ps(vAlpha, vWindow);
}


__m128 RxVoxelShader::GetRawWindowing(int iI, int iJ, int iK, int iInternalNum, float fAlpha)
{
	unsigned short unDensity = GetInterpolatedDensity(iI, iJ, iK, iInternalNum);
	__m128 vAlpha = _mm_set1_ps(fAlpha);
	__m128 vWindow = m_aByte2M128Table[m_pbyUnshadedWindowing[unDensity]];

	return _mm_mul_ps(vAlpha, vWindow);
}



BOOL RxVoxelShader::SetShadingMode(BOOL bMode)
{
	m_bShade = bMode;
	return m_bShade;
}


void RxVoxelShader::operator=(const RxVoxelShader& info)
{
	m_pVolume = info.m_pVolume;
	m_pnOriginalVolume = info.m_pnOriginalVolume;
	
	memcpy(m_pfIntensity, info.m_pfIntensity, sizeof(float)*MAX_NORMAL);


	m_pnNormalPy = info.m_pnNormalPy;
	memcpy(m_anNormalPyStorage, info.m_anNormalPyStorage, sizeof(short)*(1 + 2 * NORM_N));
	m_pnNormalXLimit = info.m_pnNormalXLimit;							/* max abs(x) allowed for a given y */
	memcpy(m_anNormalXLimitStorage, info.m_anNormalXLimitStorage, sizeof(short)*(1 + 2 * NORM_N));

	m_bShade = info.m_bShade;


	m_iBright = info.m_iBright;
	m_iContrast = info.m_iContrast;
	m_mtxView = info.m_mtxView;

	// volume size
	m_iVoxX = info.m_iVoxX;
	m_iVoxY = info.m_iVoxY;
	m_iVoxZ = info.m_iVoxZ;


	// voi
	m_iVOIXm=info.m_iVOIXm, m_iVOIXM=info.m_iVOIXM;
	m_iVOIYm=info.m_iVOIYm, m_iVOIYM=info.m_iVOIYM;
	m_iVOIZm=info.m_iVOIZm, m_iVOIZM=info.m_iVOIZM;


	// voxel size
	m_fVoxelSizeX = info.m_fVoxelSizeX;
	m_fVoxelSizeY = info.m_fVoxelSizeY;
	m_fVoxelSizeZ = info.m_fVoxelSizeZ;
	m_iRatioZ = info.m_iRatioZ;

	// inverse
	m_iInverseVoxelSizeX = info.m_iInverseVoxelSizeX;
	m_iInverseVoxelSizeY = info.m_iInverseVoxelSizeY;
	m_iInverseVoxelSizeZ = info.m_iInverseVoxelSizeZ;

	m_iInterpolatedZ = info.m_iInterpolatedZ;


	// offset
	m_iOffsetX = info.m_iOffsetX;
	m_iOffsetY = info.m_iOffsetY;
	m_iOffsetZ = info.m_iOffsetZ;


	// view
	m_iPrincipalAxis = info.m_iPrincipalAxis;
	memcpy(m_iDistance, info.m_iDistance, sizeof(int)*2*3);
	m_fSamplingDistance = info.m_fSamplingDistance;


	m_iShadeBufferSize = info.m_iShadeBufferSize;

	m_pbyShadedWindowing = info.m_pbyShadedWindowing;	// 寇何俊辑 且寸
	m_pbyUnshadedWindowing = info.m_pbyUnshadedWindowing;	// 寇何俊辑 且寸
	m_pLocalWindowing = info.m_pLocalWindowing;			// 寇何俊辑 且寸

	memcpy(m_aByte2M128Table, info.m_aByte2M128Table, 256*16);

	RxDensityTable::operator =(info);
}