spectrum.cpp

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

#include "StdAfx.h"
#include <xmmintrin.h>
#include <malloc.h>
#include <math.h>
#include "Spectrum.h"

static const int MIN_DENSITY=0;
static const int MAX_DENSITY=4095;
static const int CACHESIZE=32;

RxSpectrum::RxSpectrum()
{
	m_pTable = (__m128*)_mm_malloc(DENSITY_RANGE*16, CACHESIZE);
	m_pMPRTable = (__m128*)_mm_malloc(DENSITY_RANGE*16, CACHESIZE);

	m_nVolumeNum = 0;
}

RxSpectrum::~RxSpectrum()
{
	_mm_free(m_pTable);
	_mm_free(m_pMPRTable);
}



void RxSpectrum::Initialize(int nStart, int nEnd)
{
	// do nothing
	ASSERT(nEnd == nStart+DENSITY_RANGE-1);
	m_nVolumeNum = 0;
}

// fcolor : b, g, r, a
void RxSpectrum::SetVolumeInfo(short nDensityValue[4], float fColor[13])
{
	ASSERT(IsInputValid(nDensityValue, fColor));
	for(int i=0; i<4; i++) 
		m_pVolumeInfo[m_nVolumeNum].m_nDensityValue[i] = nDensityValue[i];

	for( i=0; i<13; i++) 
		m_pVolumeInfo[m_nVolumeNum].m_fColor[i] = fColor[i];
	
	m_nVolumeNum++;
}


// 2001.11.17. density狼 裹困啊 MIN_DENSITY,MAX_DENSITY俊 乐瘤 臼酒档 登档废
BOOL RxSpectrum::IsInputValid(short nDensityValue[4], float fColor[4])
{
	int i;
	for(i=0; i<4; i++) {
//		if(MIN_DENSITY <= nDensityValue[i] && nDensityValue[i] <= MAX_DENSITY) {
//		} else {
//			return FALSE;
//		};
		if(0.0f <= fColor[i] && fColor[i] <= 1.0f) {
		} else {
			return FALSE;
		}
	}
	for(i=0; i<3; i++) {
		if(nDensityValue[i] <= nDensityValue[i+1]) {
		} else {
			return FALSE;
		}
	}
	return TRUE;
}


BOOL inline RxSpectrum::InRange(int density)
{
	return (MIN_DENSITY <= density && density <= MAX_DENSITY);
}
// alpha甫 拌魂秦辑 table阑 积己茄饶 color甫 饶俊 蚌茄促.
/*               ........       
 *              /        \      
 *         ..../          \.....
 *            p1 p2     p3 p4   
 */


// color绰 mem俊 alpha, bb, gg, rr狼 鉴栏肺 涝仿 : note曼绊.
void RxSpectrum::MakeDensityTable()
{
	int		i,j;
	float	fLengthInverse;
	float color[12];
	float deltaC[2][3];
	float fAlpha;

	__m128	vDeltaWeight;
	__m128	vWeight;
	__m128	vColor;

	__m128  *pTable = NULL;
	__m128	vMaxAlpha[DENSITY_RANGE];

	float bb,gg,rr;

	
	for(int count=0; count<2; count++) {

		if(count==0) { // mpr table父甸扁
			pTable = m_pMPRTable;
		} else { // rendering table父甸扁
			pTable = m_pTable;
		}

		memset(pTable, 0, DENSITY_RANGE*16);
		memset(vMaxAlpha, 0 , DENSITY_RANGE*16);

		// for each volume
		for(i=0; i<m_nVolumeNum; i++) {

			fAlpha = (count == 0)?m_pVolumeInfo[i].m_fColor[12]: m_pVolumeInfo[i].m_fColor[12]*m_pVolumeInfo[i].m_fColor[12];

			for(j=0; j<12; j++) 
				color[j]= m_pVolumeInfo[i].m_fColor[j] * fAlpha ; 
			// p1....p2
			if(m_pVolumeInfo[i].m_nDensityValue[1] != m_pVolumeInfo[i].m_nDensityValue[0])
				fLengthInverse = 1.0f / (float)(m_pVolumeInfo[i].m_nDensityValue[1] - m_pVolumeInfo[i].m_nDensityValue[0]);
			else
				fLengthInverse = 0.0f;

			vDeltaWeight	= _mm_set_ps1(fLengthInverse);
			vWeight			= _mm_setzero_ps();

			int pos=0;
			for(j=0; j<3; j++) 
				deltaC[0][j]=color[pos++]*fLengthInverse;
			
			for(j=0; j<3; j++) 
				deltaC[1][j]=color[pos++]*fLengthInverse;
			
			bb = color[0];
			gg = color[1];
			rr = color[2];


			for(j=m_pVolumeInfo[i].m_nDensityValue[0]; j<m_pVolumeInfo[i].m_nDensityValue[1]; j++) {

				vWeight = _mm_add_ps(vWeight, vDeltaWeight);
				
				//tykim
				vColor = _mm_set_ps(rr, gg, bb, fAlpha);		// r,g,b,1
				bb=bb-deltaC[0][0]+deltaC[1][0];
				gg=gg-deltaC[0][1]+deltaC[1][1];
				rr=rr-deltaC[0][2]+deltaC[1][2];

				__m128 vResultColor = _mm_mul_ps(vColor, vWeight);

				if(InRange(j)) {
					pTable[j] = _mm_add_ps(pTable[j], vResultColor);
					vMaxAlpha[j] = _mm_max_ss(vMaxAlpha[j], vResultColor); // ss档 惑包 绝促.
				}

			}


			if(m_pVolumeInfo[i].m_nDensityValue[2] != m_pVolumeInfo[i].m_nDensityValue[1])
				fLengthInverse = 1.0f / (float)(m_pVolumeInfo[i].m_nDensityValue[2] - m_pVolumeInfo[i].m_nDensityValue[1]);
			else
				fLengthInverse = 0.0f;

			//tykim
			pos=3;
			for(j=0; j<3; j++) 
				deltaC[0][j]=color[pos++]*fLengthInverse;
			
			for(j=0; j<3; j++) 
				deltaC[1][j]=color[pos++]*fLengthInverse;
			
			
			bb = color[3];
			gg = color[4];
			rr = color[5];

			// p2....p3
			for(j=m_pVolumeInfo[i].m_nDensityValue[1]; j<m_pVolumeInfo[i].m_nDensityValue[2]; j++) {
				__m128 vResultColor = vColor;

				//tykim
				vResultColor = _mm_set_ps(rr, gg, bb, fAlpha);		// r,g,b,1

				bb=bb-deltaC[0][0]+deltaC[1][0];
				gg=gg-deltaC[0][1]+deltaC[1][1];
				rr=rr-deltaC[0][2]+deltaC[1][2];

				if(InRange(j)) {
					pTable[j] = _mm_add_ps(pTable[j], vResultColor);
					vMaxAlpha[j] = _mm_max_ss(vMaxAlpha[j], vResultColor); // ss档 惑包 绝促
				}
			}

			// p3....p4
			if(m_pVolumeInfo[i].m_nDensityValue[2] != m_pVolumeInfo[i].m_nDensityValue[3])
				fLengthInverse = 1.0f / (float)(m_pVolumeInfo[i].m_nDensityValue[3] - m_pVolumeInfo[i].m_nDensityValue[2]);
			else
				fLengthInverse = 0.0f;
			vDeltaWeight	= _mm_set_ps1(fLengthInverse);
			vWeight			= _mm_set_ps1(1.0f);

			pos=6;
			for(j=0; j<3; j++) 
				deltaC[0][j]=color[pos++]*fLengthInverse;
			
			for(j=0; j<3; j++) 
				deltaC[1][j]=color[pos++]*fLengthInverse;

			bb = color[6];
			gg = color[7];
			rr = color[8];

			for(j=m_pVolumeInfo[i].m_nDensityValue[2]; j<m_pVolumeInfo[i].m_nDensityValue[3]; j++) {
			
				vColor = _mm_set_ps(rr, gg, bb, fAlpha);		// r,g,b,1
				vWeight = _mm_sub_ps(vWeight, vDeltaWeight);
				__m128 vResultColor = _mm_mul_ps(vColor, vWeight);

				bb=bb-deltaC[0][0]+deltaC[1][0];
				gg=gg-deltaC[0][1]+deltaC[1][1];
				rr=rr-deltaC[0][2]+deltaC[1][2];

				if(InRange(j)) {
					pTable[j] = _mm_add_ps(pTable[j], vResultColor);
					vMaxAlpha[j] = _mm_max_ss(vMaxAlpha[j], vResultColor); // ss档 惑包 绝促
				}
			}
		}


		// prefetch nta
		for(i=0; i<=MAX_DENSITY-MIN_DENSITY; i++) {

			int bEq = _mm_comigt_ss(pTable[i], _mm_setzero_ps());
			if(bEq) {

				__m128 vInvAlpha = _mm_rcp_ss(pTable[i]);
				vInvAlpha = _mm_mul_ss(vInvAlpha, vMaxAlpha[i]);
				vInvAlpha = _mm_shuffle_ps(vInvAlpha, vInvAlpha, 0);
				pTable[i] = _mm_mul_ps(pTable[i], vInvAlpha);
			}
		}

	} // end of two table
}



// return value: RGBQuad(bb,gg,rr,reserved)
unsigned int RxSpectrum::GetDIBColor(int nDensity)
{
	// 0-1
	__m128 vColor = m_pMPRTable[nDensity];

	// 0-255
	vColor = _mm_mul_ps(vColor, _mm_set_ps1(255.0f));

	// convert to __m64
	__m64 v64Color = _mm_packs_pu16(_mm_cvtps_pi16(vColor), _mm_setzero_si64());

	// convert to unsigned int (not signed int)
	unsigned int iABGRColor = (unsigned int)_mm_cvtsi64_si32(v64Color);

	unsigned int iBGRAColor = (iABGRColor << 24) | (iABGRColor >> 8);	// don't care alpha
	_mm_empty();									// don't remove (255.0f)
	return iBGRAColor;
}

/*
// 0 <= nWindowingValue <= 255
unsigned int RxSpectrum::GetDIBColor(int nDensity, int nWindowingValue)
{
	// 0-1
	__m128 vColor = m_pMPRTable[nDensity];

	// 0-255
	vColor = _mm_mul_ps(vColor, _mm_set_ps1((float)nWindowingValue));

	// convert to __m64
	__m64 v64Color = _mm_packs_pu16(_mm_cvtps_pi16(vColor), _mm_setzero_si64());

	// convert to unsigned int (not signed int)
	unsigned int iABGRColor = (unsigned int)_mm_cvtsi64_si32(v64Color);

	unsigned int iBGRAColor = (iABGRColor << 24) | (iABGRColor >> 8);	// don't care alpha
	_mm_empty();									// don't remove (255.0f)
	return iBGRAColor;
}
*/

unsigned int RxSpectrum::GetInterpolatedDIBColor(int nDensity1, int nDensity2, float fWeight1, float fWeight2)
{
	__m128 vColor1 = m_pMPRTable[nDensity1];
	__m128 vColor2 = m_pMPRTable[nDensity2];

	__m128 vWeight1 = _mm_set_ss(fWeight1);
	__m128 vWeight2 = _mm_set_ss(fWeight2);
	__m128 vWeight = _mm_add_ss(vWeight1, vWeight2);

	vWeight1 = _mm_shuffle_ps(vWeight1, vWeight1, 0);
	vWeight2 = _mm_shuffle_ps(vWeight2, vWeight2, 0);

	__m128 vWeightInverse = _mm_rcp_ss(vWeight);
	vWeightInverse = _mm_mul_ss(vWeightInverse, _mm_set_ss(255.0f));
	vWeightInverse = _mm_shuffle_ps(vWeightInverse, vWeightInverse, 0);

	vColor1 = _mm_mul_ps(vColor1, vWeight2);
	vColor2 = _mm_mul_ps(vColor2, vWeight1);

	// weighted sum
	__m128 vColor = _mm_add_ps(vColor1, vColor2);

	// 0 - 255.0f
	vColor = _mm_mul_ps(vColor, vWeightInverse);
	vColor = _mm_shuffle_ps(vColor, vColor, _MM_SHUFFLE(0, 1, 2, 3));	// color啊 捞惑窍搁 捞巴阑 官层杭 巴

	__m64 v32Color = _mm_packs_pu16 (_mm_cvtps_pi16(vColor), _mm_setzero_si64());

	unsigned int nRetval = _mm_cvtsi64_si32(v32Color);
	_mm_empty();
	return nRetval;
}