rigidreg3dagent.cpp

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

//********************************************************************
//	SameTransformFactor(): Check same resolution and voxel size
//********************************************************************
BOOL RxRigidReg3DAgent::SameTransformFactor()
{
	if(m_pTransform->xmove == 0.0 && m_pTransform->ymove == 0.0 && m_pTransform->zmove == 0.0) 
		if(m_pTransform->xscale == 1.0 && m_pTransform->yscale == 1.0 && m_pTransform->zscale == 1.0) 
			return TRUE;

	return FALSE;
}

//********************************************************************
//	AutomaticConverge(): automatic convergence
//********************************************************************
BOOL RxRigidReg3DAgent::AutomaticConverge()
{
	int		bCONVERGE = FALSE, repeat = 0;
	double	vector[4], transform[4], gdm[6];

	// Set parameters related to mutual information
	m_pDerivateMI3D->SetVolumeData(m_pnReferVolume);
	m_pDerivateMI3D->SetReferVolInfo(m_pRigidReg3DInfo->m_ReferVolInfo);
	m_pDerivateMI3D->SetSampleNumber(m_pRigidReg3DInfo->m_SampleNum);
	m_pDerivateMI3D->SetSampleA(m_sampleA);
	m_pDerivateMI3D->SetTransform(m_pTransform);
	m_pDerivateMI3D->SetMIParameters(m_pRigidReg3DInfo->m_varFloat,m_pRigidReg3DInfo->m_covarRefer,m_pRigidReg3DInfo->m_covarFloat);

	// Automatic convergence
	do {
		gdm[0] = gdm[1] = gdm[2] = gdm[3] = gdm[4] = gdm[5] = 0.0;

		// Interpolate and Generate corresponding sample points in float volume
		m_sampleB = new RxAutoSample3DInfo[m_pRigidReg3DInfo->m_SampleNum];
		for (int row = 0; row < m_pRigidReg3DInfo->m_SampleNum; row++) {
			vector[0] = (double)m_sampleA[row].x;
			vector[1] = (double)m_sampleA[row].y;
			vector[2] = (double)m_sampleA[row].z;
			vector[3] = 1.0;

			m_pMatrix3D->matrix3D_vector4(m_inverse_obj2eye,vector,transform);
			//m_pMatrix3D->matrix3D_vector4(m_obj2eye,vector,transform);
	
#ifdef	NEAREST_NEIGHBOR
			NearestNeighbor(transform,&m_sampleB[row].x,&m_sampleB[row].y,&m_sampleB[row].z,&m_sampleB[row].val);
#endif
#ifdef	TRILINEAR_INTERPOLATION
			TrilinearInterpol(transform,&m_sampleB[row].x,&m_sampleB[row].y,&m_sampleB[row].z,&m_sampleB[row].val);
#endif
#ifdef	TRILINEAR_INTERPOLATION_ORIGIN
			TrilinearInterpolOrigin(transform,&m_sampleB[row].x,&m_sampleB[row].y,&m_sampleB[row].z,&m_sampleB[row].val);
#endif
		} 

		// Calculate the derivative of mutual information
		double	totalfactor[6];
		totalfactor[0] = totalfactor[1] = totalfactor[2] = totalfactor[3] = totalfactor[4] = totalfactor[5] = 0.0;
		m_pDerivateMI3D->Calculate(m_sampleB,&totalfactor);

		// Stochastic steepest descending method
		gdm[0]	= m_pRigidReg3DInfo->m_etaTrans * totalfactor[0];
		gdm[1]	= m_pRigidReg3DInfo->m_etaTrans * totalfactor[1];
		gdm[2]	= m_pRigidReg3DInfo->m_etaTrans * totalfactor[2];
		gdm[3]	= m_pRigidReg3DInfo->m_etaRot   * totalfactor[3];
		gdm[4]	= m_pRigidReg3DInfo->m_etaRot   * totalfactor[4];
		gdm[5]	= m_pRigidReg3DInfo->m_etaRot	* totalfactor[5];

		m_pTransform->xmove	+= gdm[0];
		m_pTransform->ymove	+= gdm[1];
		m_pTransform->zmove	+= gdm[2];
		m_pTransform->xrot	+= gdm[3];
		m_pTransform->yrot	+= gdm[4];
		m_pTransform->zrot	+= gdm[5];

		// Generate updated view matrix on updated transformation vector
		m_pMatrix3D->make_viewmat3D(m_obj2eye,m_pTransform);
		m_pMatrix3D->invert_mat3D(m_obj2eye,m_inverse_obj2eye);

		delete [] m_sampleB;

		if(fabs(gdm[0]) < m_pRigidReg3DInfo->m_TransLimit && fabs(gdm[1]) < m_pRigidReg3DInfo->m_TransLimit && fabs(gdm[2]) < m_pRigidReg3DInfo->m_TransLimit &&
			fabs(gdm[3]) < m_pRigidReg3DInfo->m_RotLimit && fabs(gdm[4]) < m_pRigidReg3DInfo->m_RotLimit && fabs(gdm[5]) < m_pRigidReg3DInfo->m_RotLimit) {
			bCONVERGE	= TRUE;
			m_Repeat	= repeat+1;
		}
		repeat++;
	} while (repeat < m_pRigidReg3DInfo->m_maxIter && !bCONVERGE);

	return TRUE;
}

//********************************************************************
//	NearestNeighbor(): nearest neighbor 
//********************************************************************
BOOL RxRigidReg3DAgent::NearestNeighbor(double transform[4],int *x,int *y,int *z,unsigned short *val)
{
	*x		= (int)(transform[0]+0.5);
	*y		= (int)(transform[1]+0.5);
	*z		= (int)(transform[2]+0.5);
	*val	= GetIntensity3D(*x, *y, *z);

	return TRUE;
}

//********************************************************************
//	TrilinearInterpol(): trilinear interpolation
//********************************************************************
BOOL RxRigidReg3DAgent::TrilinearInterpol(double transform[4],int *x,int *y,int *z,unsigned short *val)
{
	// Trilinear interpolation
	int	xlow	= (int)transform[0];
	int	ylow	= (int)transform[1];
	int	zlow	= (int)transform[2];

	// Boundary check and defined high value
	int	xhigh, yhigh, zhigh;
	BoundaryCheck(&xlow,&xhigh,m_pRigidReg3DInfo->m_FloatVolInfo.Ru-1);
	BoundaryCheck(&ylow,&yhigh,m_pRigidReg3DInfo->m_FloatVolInfo.Rv-1);
	BoundaryCheck(&zlow,&zhigh,m_pRigidReg3DInfo->m_FloatVolInfo.Rw-1);

	// Interpolation fractions
	double	xx = transform[0] - (double)xlow;
	double	xX = 1.0 - xx;
	double	yy = transform[1] - (double)ylow;
	double	zz = transform[2] - (double)zlow;

	// Calculate weighting factor for all 8-neighbors
	double xxyy = xx * yy;
	double xxyY = xx - xxyy;
	double xXyy = yy - xxyy;
	double xXyY = xX - xXyy;

	// Calculate weighting factors
	double w2	= xXyY * zz;
	double w4	= xXyy * zz;
	double w6	= xxyY * zz;
	double w8	= xxyy * zz;

	double w1	= xXyY - w2;
	double w3	= xXyy - w4;
	double w5	= xxyY - w6;
	double w7	= xxyy - w8;

	// Trilinear grayvalue interpolation
	unsigned short value = (unsigned short)
		((GetIntensity3D(xlow,ylow,zlow)*w1	  + GetIntensity3D(xlow,ylow,zhigh)*w2 +
		  GetIntensity3D(xlow,yhigh,zlow)*w3  + GetIntensity3D(xlow,yhigh,zhigh)*w4 +
		  GetIntensity3D(xhigh,ylow,zlow)*w5  + GetIntensity3D(xhigh,ylow,zhigh)*w6 +
		  GetIntensity3D(xhigh,yhigh,zlow)*w7 + GetIntensity3D(xhigh,yhigh,zhigh)*w8)+0.5);

	// Assign transformed coordinate and grayvalue to sampleB array
	*x = xlow;	*y = ylow;	*z = zlow;
	*val = value;
	
	return TRUE;
}

//********************************************************************
//	TrilienarInterpolOrigin(): Original trilinear interpolation
//********************************************************************
BOOL RxRigidReg3DAgent::TrilinearInterpolOrigin(double transform[4],int *x,int *y,int *z,unsigned short *val)
{
	// Trilinear interpolation
	int	xlow = (int)transform[0];
	int	ylow = (int)transform[1];
	int	zlow = (int)transform[2];

	// Boundary check and defined high value
	int	xhigh, yhigh, zhigh;
	BoundaryCheck(&xlow,&xhigh,m_pRigidReg3DInfo->m_FloatVolInfo.Ru-1);
	BoundaryCheck(&ylow,&yhigh,m_pRigidReg3DInfo->m_FloatVolInfo.Rv-1);
	BoundaryCheck(&zlow,&zhigh,m_pRigidReg3DInfo->m_FloatVolInfo.Rw-1);

	// Interpolation fractions
	double	fWX = transform[0] - xlow;
	double	fWY = transform[1] - ylow;
	double	fWZ = transform[2] - zlow;

	double	fInvWX = 1.0 - fWX;
	double	fInvWY = 1.0 - fWY;
	double	fInvWZ = 1.0 - fWZ;

	double	w1 = fInvWX*fInvWY*fInvWZ;		// (1-x)(1-y)(1-z)
	double	w2 = fWX*fInvWY*fInvWZ;			// x(1-y)(1-z)
	double	w3 = fInvWX*fWY*fInvWZ;			// (1-x)y(1-z)
	double	w4 = fWX*fWY*fInvWZ;			// xy(1-z)
	double	w5 = fInvWX*fInvWY*fWZ;			// (1-x)(1-y)z
	double	w6 = fWX*fInvWY*fWZ;			// x(1-y)z
	double	w7 = fInvWX*fWY*fWZ;			// (1-x)yz
	double	w8 = fWX*fWY*fWZ;				// xyz

	// Trilinear grayvalue interpolation
	unsigned short value = (unsigned short)
		((GetIntensity3D(xlow,ylow,zlow)*w1	  + GetIntensity3D(xlow,ylow,zhigh)*w2 +
		  GetIntensity3D(xlow,yhigh,zlow)*w3  + GetIntensity3D(xlow,yhigh,zhigh)*w4 +
		  GetIntensity3D(xhigh,ylow,zlow)*w5  + GetIntensity3D(xhigh,ylow,zhigh)*w6 +
		  GetIntensity3D(xhigh,yhigh,zlow)*w7 + GetIntensity3D(xhigh,yhigh,zhigh)*w8)+0.5);

	// Assign transformed coordinate and grayvalue to sampleB array
	*x = xlow;	*y = ylow;	*z = zlow;
	*val = value;
	
	return TRUE;
}

//********************************************************************
//	GetIntensity3D(): Get intensity value from volume memory
//********************************************************************
unsigned short RxRigidReg3DAgent::GetIntensity3D(int x,int y,int z)
{
	if(x < 0 || y < 0 || z < 0 ||
		x > m_pRigidReg3DInfo->m_FloatVolInfo.Ru-1 || y > m_pRigidReg3DInfo->m_FloatVolInfo.Rv-1 || z > m_pRigidReg3DInfo->m_FloatVolInfo.Rw-1)
		return 0;
	return m_pnFloatVolume[((z*m_pRigidReg3DInfo->m_FloatVolInfo.Rv)+y)*m_pRigidReg3DInfo->m_FloatVolInfo.Ru+x];
}

//********************************************************************
//	BoundaryCheck(): boundary check
//********************************************************************
BOOL RxRigidReg3DAgent::BoundaryCheck(int *low,int *high,int boundary)
{
	if(*low >= boundary) {
		*low = *high = boundary;
	}
	else if(*low < 0) {
		*low = 0;
		*high = 1;
	}
	else {
		*high = *low + 1;
	}

	return TRUE;
}

//********************************************************************
//	IsValid_MIAnalysis(): Check whether MI analysis perform or not
//********************************************************************
BOOL RxRigidReg3DAgent::IsValid_MIAnalysis(int valid)
{
	if(valid == 1) 
		return TRUE;
	return FALSE;
}