vrraycastingagent.cpp

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

	int *bix = m_piBlockIdx[0], *biy = m_piBlockIdx[1], *biz = m_piBlockIdx[2];
	int *six = m_piSkipIdx[0], *siy = m_piSkipIdx[1], *siz = m_piSkipIdx[2];
	int *rix = m_piIdx[0], *riy = m_piIdx[1], *riz = m_piIdx[2];

	for (iy = iStartY; iy < iEndY; iy+=iSkip) {

		if (::WaitForSingleObject(g_hVRThreadKillEvent, 0)==WAIT_OBJECT_0) {
			return FALSE;
		}

		for (ix = iMinX; ix < iMaxX; ix+=iSkip) {
//			if (pIsoEndoInfo->m_bReuse&& ix%(iSkip*2)==0 && iy%(iSkip*2)==0) {
//				continue;
//			}

			r = 0;
			g = 0;
			b = 0;
			opacity = 0;
			iCurLocalizer = 0;
			m_iLocalizerPlane[iPart] = 0;
			int iPrevSampleDensity = INIT_VALUE;

			float d1[4] = {IntToFloat(ix), IntToFloat(iy), 0 /*m_fWMin*/, 1};
			m4v4mulf(afRayOrigin, m_afInvViewingM, d1);

			float afRefRayOrigin[4];
			m4v4mulf(afRefRayOrigin, afRefInvViewing, d1);


			int idx = 4*((ix/iSkip)+(iy/iSkip)*(m_iImgX/iSkip)), iType;

			if(GetRayBoundInReference(afRefRayOrigin, min_t, max_t)==FALSE)
				continue;

			
		switch ((iType = GetRaySampleBound(iPart, afRayOrigin, &min_t, &max_t, &r, &g, &b, &opacity))) {
			case 0:
				continue;
			case 1:
				iSample = FloatToInt(GetTLSample(m_afCropMin[0], 
					__min(__max(m_afCropMin[1], afRayOrigin[1]+rv[1]*min_t), m_afCropMax[1]-0.01f), 
					__min(__max(m_afCropMin[2], afRayOrigin[2]+rv[2]*min_t), m_afCropMax[2]-0.01f)));
				break;
			case 4:
				iSample = FloatToInt(GetTLSample(m_afCropMax[0]-0.01f, 
					__min(__max(m_afCropMin[1], afRayOrigin[1]+rv[1]*min_t), m_afCropMax[1]-0.01f), 
					__min(__max(m_afCropMin[2], afRayOrigin[2]+rv[2]*min_t), m_afCropMax[2]-0.01f)));
				break;
			case 2:
				iSample = FloatToInt(GetTLSample(__min(__max(m_afCropMin[0], afRayOrigin[0]+rv[0]*min_t), m_afCropMax[0]-0.01f), 
					m_afCropMin[1], 
					__min(__max(m_afCropMin[2], afRayOrigin[2]+rv[2]*min_t), m_afCropMax[2]-0.01f)));
				break;
			case 5:
				iSample = FloatToInt(GetTLSample(__min(__max(m_afCropMin[0], afRayOrigin[0]+rv[0]*min_t), m_afCropMax[0]-0.01f),
					m_afCropMax[1]-0.01f, 
					__min(__max(m_afCropMin[2], afRayOrigin[2]+rv[2]*min_t), m_afCropMax[2]-0.01f)));
				break;
			case 3:
				iSample = FloatToInt(GetTLSample(__min(__max(m_afCropMin[0], afRayOrigin[0]+rv[0]*min_t), m_afCropMax[0]-0.01f),
					__min(__max(m_afCropMin[1], afRayOrigin[1]+rv[1]*min_t), m_afCropMax[1]-0.01f), 
					m_afCropMin[2]));
				break;
			case 6:
				iSample = FloatToInt(GetTLSample(__min(__max(m_afCropMin[0], afRayOrigin[0]+rv[0]*min_t), m_afCropMax[0]-0.01f),
					__min(__max(m_afCropMin[1], afRayOrigin[1]+rv[1]*min_t), m_afCropMax[1]-0.01f), 
					m_afCropMax[2]-0.01f));
				break;
			case 7:
				iSample = FloatToInt(GetTLSample(afRayOrigin[0]+rv[0]*min_t,
					afRayOrigin[1]+rv[1]*min_t,
					afRayOrigin[2]+rv[2]*min_t));
				break;
			case 8:
				pImage[idx] = (unsigned char)FloatToInt(b * 255.0f);
				pImage[idx+1] = (unsigned char)FloatToInt(g * 255.0f);
				pImage[idx+2] = (unsigned char)FloatToInt(r * 255.0f);
				pImage[idx+3] = (unsigned char)FloatToInt(opacity * 255.0f);
				continue;
			}


		float fColor[4];
		if(iType == 9)
			fColor[0]=fColor[1]=fColor[2]=fColor[3]=0;
		else
			m_pShader->GetColor(iSample, fColor);

			if (fColor[0]!=0 /* && iSample>=m_iMin && iSample<=m_iMax */ ) {
				r = fColor[3];
				g = fColor[2];
				b = fColor[1];
				opacity = fColor[0];
				min_t++;
			}


			for (t = min_t; t < max_t; t += blk_sample_step) {
				ipx = FloatToInt(afRayOrigin[0] + rv[0] * t);
				ipy = FloatToInt(afRayOrigin[1] + rv[1] * t);
				ipz = FloatToInt(afRayOrigin[2] + rv[2] * t);

				if( !m_pShader->IsTransparent(m_pmmBlock[bix[ipx+aiBlockAdjust[0]] + biy[ipy+aiBlockAdjust[1]] + biz[ipz+aiBlockAdjust[2]]].nMax,
					m_pmmBlock[bix[ipx+aiBlockAdjust[0]] + biy[ipy+aiBlockAdjust[1]] + biz[ipz+aiBlockAdjust[2]]].nMin) ) {

					limit_t = t + blk_sample_step;
					if (limit_t > max_t)
						limit_t = max_t;
					for (float tt = t; tt < limit_t; tt += m_fSampleDistance) {
						ipx = FloatToInt(afRayOrigin[0] + rv[0] * tt);
						ipy = FloatToInt(afRayOrigin[1] + rv[1] * tt);
						ipz = FloatToInt(afRayOrigin[2] + rv[2] * tt);
						
						if ( !m_pShader->IsTransparent(m_pmmDen[six[ipx+aiBlockAdjust[0]] + siy[ipy+aiBlockAdjust[1]] + siz[ipz+aiBlockAdjust[2]]].nMax,
							m_pmmDen[six[ipx+aiBlockAdjust[0]] + siy[ipy+aiBlockAdjust[1]] + siz[ipz+aiBlockAdjust[2]]].nMin) ) {

							GetRaySample(iPart, afRayOrigin, tt, __min(max_t, tt+1), &cur_intensity, afVoxel,
								r, g, b, opacity, iCurLocalizer, iPrevSampleDensity);
							if (opacity >= OPACITY_THRESHOLD) {
								t = max_t;
								break;
							}
						}
					}
				}
			}
			
			r = m_pIsoEndoInfo->m_fBkR * (1-opacity) + r;
			g = m_pIsoEndoInfo->m_fBkG * (1-opacity) + g;
			b = m_pIsoEndoInfo->m_fBkB * (1-opacity) + b;

			if (r >= 1.0f) r = 1.0f;  // Oh! these case should not happen
			if (g >= 1.0f) g = 1.0f;  // because this may cause color shift
			if (b >= 1.0f) b = 1.0f;

			pImage[idx] = (unsigned char)FloatToInt(b * 255.0f);
			pImage[idx+1] = (unsigned char)FloatToInt(g * 255.0f);
			pImage[idx+2] = (unsigned char)FloatToInt(r * 255.0f);
			pImage[idx+3] = (unsigned char)FloatToInt(opacity*255.0f);
		}
	}
	
	pIsoEndoInfo->m_iElapsedTime = int(timer.GetElapsedTime() * 1000000);
	return TRUE;
}

BOOL RxVRRayCastingAgent::GetRayDirectionInRef(float afRefInvViewing[4][4])
{
	for(int i=0; i<4; i++)
		for(int j=0; j<4; j++)
			afRefInvViewing[i][j] = (float)m_pIsoEndoInfo->m_mxInvRefViewing.m[i][j];

	float afTemp[4] = {0, 0, 1, 0};
	m4v4mulf(m_afRefRayDirection, afRefInvViewing, afTemp);
	// | m_afRefRayDirection | == 1
	m_afRefRayDirection[0] /= m_pIsoEndoInfo->m_fRayMagnitude;
	m_afRefRayDirection[1] /= m_pIsoEndoInfo->m_fRayMagnitude;
	m_afRefRayDirection[2] /= m_pIsoEndoInfo->m_fRayMagnitude;
	m_afRefRayDirection[3] /= m_pIsoEndoInfo->m_fRayMagnitude;

	m_afRefInvRayDirection[0] = 1.f/m_afRefRayDirection[0];
	m_afRefInvRayDirection[1] = 1.f/m_afRefRayDirection[1];
	m_afRefInvRayDirection[2] = 1.f/m_afRefRayDirection[2];

	return TRUE;
}

BOOL RxVRRayCastingAgent::GetRayBoundInReference(float *afRefRayOrigin, float &min_t, float &max_t)
{
	min_t = -MAX_T;
	max_t = MAX_T;

	float tmin,  tmax, tmp;

	// ref mode牢 版快 skip
	if(m_pIsoEndoInfo->IsRef()) {
		return TRUE;
	}

	for(int i=0; i<3; i++) {
		if (__abs(m_afRefRayDirection[i])>0.00001f) {
			tmax = (m_pIsoEndoInfo->m_iRefCropMax[i]-0.1f-afRefRayOrigin[i])*m_afRefInvRayDirection[i];
			tmin = (m_pIsoEndoInfo->m_iRefCropMin[i]+0.1f-afRefRayOrigin[i])*m_afRefInvRayDirection[i];

			if (tmax < tmin) {
				tmp = tmax;
				tmax = tmin;
				tmin = tmp;
			}
			if (tmin > min_t) {
				min_t = tmin;
			}
			if (tmax < max_t) {
				max_t = tmax;
			}
		}
	}
	for(i=0; i<3; i++) {
		float x = afRefRayOrigin[i]+(min_t+1)*m_afRefRayDirection[i];
		if(!((m_pIsoEndoInfo->m_iRefCropMin[i] <= x) && (x <= m_pIsoEndoInfo->m_iRefCropMax[i])))
			return FALSE;
		x = afRefRayOrigin[i]+(max_t-1)*m_afRefRayDirection[i];
		if(!((m_pIsoEndoInfo->m_iRefCropMin[i] <= x) && (x <= m_pIsoEndoInfo->m_iRefCropMax[i])))
			return FALSE;
	}
	return TRUE;
}

int RxVRRayCastingAgent::GetRaySampleBound(int iPart, float *afRayOrigin, float *min_t, float *max_t, float* r, float* g, float* b, float* opacity)
{
	float tmin,  tmax, tmp;
	int i;
	int iFirst = 7, iPos;


	for (i=0;i<3;i++) {
		iPos = i+1;
		if (__abs(m_afRayDirection[i])>0.00001f) {
			tmax = (m_afCropMax[i]-0.1f-afRayOrigin[i])*m_afInvRayDirection[i];
			tmin = (m_afCropMin[i]+0.1f-afRayOrigin[i])*m_afInvRayDirection[i];

			if (tmax < tmin) {
				tmp = tmax;
				tmax = tmin;
				tmin = tmp;
				iPos+=3;
			}
			if (tmin > *min_t) {
				*min_t = tmin;
				iFirst = iPos;
			}
			if (tmax < *max_t)
				*max_t = tmax;
		}
	}


	// for crop
	if (!IS_VOXEL_IN_CROPVR(afRayOrigin[0]+(*min_t+1)*m_afRayDirection[0],
			afRayOrigin[1]+(*min_t+1)*m_afRayDirection[1],
			afRayOrigin[2]+(*min_t+1)*m_afRayDirection[2])) {
		return 0;
	}

	if (!IS_VOXEL_IN_CROPVR(afRayOrigin[0]+(*max_t-1)*m_afRayDirection[0],
			afRayOrigin[1]+(*max_t-1)*m_afRayDirection[1],
			afRayOrigin[2]+(*max_t-1)*m_afRayDirection[2])) {
		return 0;
	}

	BOOL bFirst = FALSE;

	bFirst = FALSE;
	if (*min_t > *max_t)
		return 0;

	RxVolumeData* pVolume = (m_pIsoEndoInfo->IsRef()) ? RxGetVolumeData(0) : RxGetVolumeData(1);
	RxMaskVolume *pMask = pVolume->GetEditVolume();
	unsigned char cMask = pMask->GetMask(int(afRayOrigin[0]+(*min_t+1)*m_afRayDirection[0]),
		int(afRayOrigin[1]+(*min_t+1)*m_afRayDirection[1]),
		int(afRayOrigin[2]+(*min_t+1)*m_afRayDirection[2]/m_fRatio));
	if(cMask == 0)
		return 9;

	return iFirst;
}

BOOL RxVRRayCastingAgent::GetRaySample(int iPart, float *afRayOrigin, float tl, float tu, 
									   float *intensity, float *afVoxel, 
									   float &r, float &g, float &b, float& cur_opacity, int& iCurLocalizer,
									   int &iPrevSampleDensity)
{
	float diffuse;
	float afLightColor[3] = {1, 1, 1};
	float ocr, ocg, ocb;
	float object_opacity;
	BOOL bRef;

	if (m_pIsoEndoInfo->IsRef())
		bRef = TRUE;
	else
		bRef = FALSE;

	float t=tl;
	afVoxel[0] = afRayOrigin[0] + m_afRayDirection[0] * t;
	afVoxel[1] = afRayOrigin[1] + m_afRayDirection[1] * t;
	afVoxel[2] = afRayOrigin[2] + m_afRayDirection[2] * t;

	RxVolumeData* pVolume = (m_pIsoEndoInfo->IsRef()) ? RxGetVolumeData(0) : RxGetVolumeData(1);
	RxMaskVolume *pSegmentMask;
	RxMaskVolume *pEdgeMask;

	pSegmentMask = pVolume->GetEditVolume();
	pEdgeMask = pVolume->GetEdgeVolume();

//	RxMaskVolume *pNoduleMask = pVolume->GetNoduleVolume();
	RxMaskVolume *pVesselMask = pVolume->GetVesselVolume();

	unsigned char cMask = pSegmentMask->GetMask(int(afVoxel[0]), int(afVoxel[1]), int(afVoxel[2]/m_fRatio));
//	unsigned char cNoduleMask = pNoduleMask->GetMask(int(afVoxel[0]), int(afVoxel[1]), int(afVoxel[2]/m_fRatio));
	unsigned char cVesselMask = FALSE;//pVesselMask->GetMask(int(afVoxel[0]), int(afVoxel[1]), int(afVoxel[2]/m_fRatio));
//	unsigned char cEdgeMask = pEdgeMask->GetMask(int(afVoxel[0]), int(afVoxel[1]), int(afVoxel[2]/m_fRatio));

	if(cMask == 0)// && cNoduleMask == 0 /*&& cEdgeMask == 0*/) // for segmentation ?
		return TRUE;

	float sample_d = GetTLSample(afVoxel[0], afVoxel[1], afVoxel[2]);
	float fColor[4];

	RxVoxelShader* m_pShader = m_pIsoEndoInfo->GetShader();

	int iSample = (int)(sample_d);
	if(iPrevSampleDensity == INIT_VALUE) {
		iPrevSampleDensity = iSample;
	}

/*	if (cNoduleMask) {
		if (pVolume->m_nMode == 0)
		{
			fColor[0] = 1;	// alpha
			fColor[1] = 0;
			fColor[2] = 0;
			fColor[3] = 1;	// color
		} else
		{
			fColor[0] = 1;	// alpha
			fColor[1] = 0;
			fColor[2] = 1;
			fColor[3] = 0;	// color
		}
	} else*/ if (cVesselMask)
	{
		if (bRef)
		{
			fColor[0] = 1;	// alpha
			fColor[1] = 0;
			fColor[2] = 0;
			fColor[3] = 1;	// color
		} else {
			fColor[0] = 1;	// alpha
			fColor[1] = 0;
			fColor[2] = 1;
			fColor[3] = 0;	// color
		}
	} else {
//		m_pShader->GetColor(iSample, fColor);
		m_pShader->GetPreIntegratedColor(iPrevSampleDensity, iSample, fColor);
	}
	iPrevSampleDensity = iSample;
	object_opacity = fColor[0];

	if (object_opacity < 0.00001)
		return TRUE;
	float fInvOpacity = 1.f/fColor[0];

	ocr = fColor[3]*fInvOpacity;
	ocg = fColor[2]*fInvOpacity;
	ocb = fColor[1]*fInvOpacity;

	float nrVec[3];
	GetNormal(afVoxel[0], afVoxel[1], afVoxel[2], &nrVec[0], &nrVec[1], &nrVec[2]);
	float fGrad = v3nrmf(nrVec);

	if (m_fGradFactor > 0) {
		if (fGrad < m_fGradFactor)
			object_opacity *= fGrad*fGrad*m_fInvGradFactor;
	}

	diffuse = v3dot(nrVec, m_afLight);
	if (diffuse < 0)
		diffuse = - diffuse; // * 0.3;

	float hwVec[3];
	hwVec[0] = m_afLight[0];
	hwVec[1] = m_afLight[1];
	hwVec[2] = m_afLight[2];
	v3nrmf(hwVec);

	float specular = v3dot(nrVec, m_afLight);
	if (specular < 0) {
		specular = -specular;
	}
	float specular2 = specular*specular;
	float specular4 = specular2*specular2;
	specular = specular4*specular4*specular2;

	float fatt = 1;
	float fSpecularColor[3] = {1, 1, 1};

/*	if (cNoduleMask) 
	{
		object_opacity = 1.0;
		cur_opacity = 0.0;

		if (pVolume->m_nMode == 0)
		{
			ocr = 1.0;
			r = (m_fAmbient*ocr+fatt*afLightColor[0]*(m_fDiffuse*ocr*diffuse+SPECULAR_BASE*specular))*(1-cur_opacity)*object_opacity*m_fBrightnessFactor;
			g = 0;
			b = 0;
			cur_opacity = 1.0;	// color
		} else
		{
			ocb = 1.0;
			r = 0;
			g = (m_fAmbient*ocg+fatt*afLightColor[1]*(m_fDiffuse*ocg*diffuse+SPECULAR_BASE*specular))*(1-cur_opacity)*object_opacity*m_fBrightnessFactor;
			b = 0;
			cur_opacity = 1.0;	// color
		}
	} else*/ if (cVesselMask)
	{
		object_opacity = 1.0;
		cur_opacity = 0.0;

		if (bRef)
		{
			ocr = 1.0;
			r = (m_fAmbient*ocr+fatt*afLightColor[0]*(m_fDiffuse*ocr*diffuse+SPECULAR_BASE*specular))*(1-cur_opacity)*object_opacity*m_fBrightnessFactor;
			g = 0;
			b = 0;
			cur_opacity = 1.0;	// color
		} else
		{
			ocb = 1.0;
			r = 0;
			g = (m_fAmbient*ocg+fatt*afLightColor[1]*(m_fDiffuse*ocg*diffuse+SPECULAR_BASE*specular))*(1-cur_opacity)*object_opacity*m_fBrightnessFactor;
			b = 0;
			cur_opacity = 1.0;	// color
		}
	} else {
		r += (m_fAmbient*ocr+fatt*afLightColor[0]*(m_fDiffuse*ocr*diffuse+SPECULAR_BASE*specular))*(1-cur_opacity)*object_opacity*m_fBrightnessFactor;
		g += (m_fAmbient*ocg+fatt*afLightColor[1]*(m_fDiffuse*ocg*diffuse+SPECULAR_BASE*specular))*(1-cur_opacity)*object_opacity*m_fBrightnessFactor;
		b += (m_fAmbient*ocb+fatt*afLightColor[2]*(m_fDiffuse*ocb*diffuse+SPECULAR_BASE*specular))*(1-cur_opacity)*object_opacity*m_fBrightnessFactor;
		cur_opacity	+= (1 - cur_opacity) * object_opacity;
	}

	return TRUE;
}