derivatemi3d.cpp

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

					
						factor1 *= factor4;

						result[0] = factor1 * factor5[0];
						result[1] = factor1 * factor5[1];
						result[2] = factor1 * factor5[2];
						result[3] = factor1 * factor5[3];
						result[4] = factor1 * factor5[4];
						result[5] = factor1 * factor5[5];

						SumofFactor[0] += result[0];
						SumofFactor[1] += result[1];
						SumofFactor[2] += result[2];
						SumofFactor[3] += result[3];
						SumofFactor[4] += result[4];
						SumofFactor[5] += result[5];
					} // end of if : factor4 != 0.0
				} // end of factor1 != 0.0
			} // end of if : j != i
		} // end of j for
	} // end of i for

	(*totalfactor)[0] = SumofFactor[0] / (double)m_SampleNum;
	(*totalfactor)[1] = SumofFactor[1] / (double)m_SampleNum;
	(*totalfactor)[2] = SumofFactor[2] / (double)m_SampleNum;
	(*totalfactor)[3] = SumofFactor[3] / (double)m_SampleNum;
	(*totalfactor)[4] = SumofFactor[4] / (double)m_SampleNum;
	(*totalfactor)[5] = SumofFactor[5] / (double)m_SampleNum;

	/*
	(*totalfactor)[0] = -(SumofFactor[0] / (double)m_SampleNum);
	(*totalfactor)[1] = -(SumofFactor[1] / (double)m_SampleNum);
	(*totalfactor)[2] = -(SumofFactor[2] / (double)m_SampleNum);
	(*totalfactor)[3] = -(SumofFactor[3] / (double)m_SampleNum);
	(*totalfactor)[4] = -(SumofFactor[4] / (double)m_SampleNum);
	(*totalfactor)[5] = -(SumofFactor[5] / (double)m_SampleNum);
	*/

	return TRUE;
}

//*********************************************************************
//	DerivateTrans3D(): Involve the volume gradient and the derivative
//					   of transformed coordinate with respect to
//				       transformation
//*********************************************************************
BOOL RxDerivateMI3D::DerivateTrans3D(RxAutoSample3DInfo vi,RxAutoSample3DInfo vj,double (*result)[6])
{
	// Generate gradient of vi and vj vectors
	double	igradient[4];
	igradient[0] = (double)GetIntensity3D(vi.x+1,vi.y,vi.z)-(double)GetIntensity3D(vi.x-1,vi.y,vi.z);
	igradient[1] = (double)GetIntensity3D(vi.x,vi.y+1,vi.z)-(double)GetIntensity3D(vi.x,vi.y-1,vi.z);
	igradient[2] = (double)GetIntensity3D(vi.x,vi.y,vi.z+1)-(double)GetIntensity3D(vi.x,vi.y,vi.z-1);
	igradient[3] = 0.0;

	double	jgradient[4];
	jgradient[0] = (double)GetIntensity3D(vj.x+1,vj.y,vj.z)-(double)GetIntensity3D(vj.x-1,vj.y,vj.z);
	jgradient[1] = (double)GetIntensity3D(vj.x,vj.y+1,vj.z)-(double)GetIntensity3D(vj.x,vj.y-1,vj.z);
	jgradient[2] = (double)GetIntensity3D(vj.x,vj.y,vj.z+1)-(double)GetIntensity3D(vj.x,vj.y,vj.z-1);
	jgradient[3] = 0.0;

	// Generate the derivative of transformed coordinate with respect to the transformation
	double	iderivate[9];
	iderivate[0] = (cos(m_Trans.zrot*pi_180)*sin(m_Trans.yrot*pi_180)*cos(m_Trans.xrot*pi_180)*vi.y) +
		           (sin(m_Trans.zrot*pi_180)*sin(m_Trans.xrot*pi_180)*vi.y) -
				   (cos(m_Trans.zrot*pi_180)*sin(m_Trans.yrot*pi_180)*sin(m_Trans.xrot*pi_180)*vi.z) +
				   (sin(m_Trans.zrot*pi_180)*cos(m_Trans.xrot*pi_180)*vi.z);
	iderivate[1] = (-cos(m_Trans.zrot*pi_180)*sin(m_Trans.yrot*pi_180)*vi.x)+
		           (cos(m_Trans.zrot*pi_180)*cos(m_Trans.yrot*pi_180)*sin(m_Trans.xrot*pi_180)*vi.y)+
				   (cos(m_Trans.zrot*pi_180)*cos(m_Trans.yrot*pi_180)*cos(m_Trans.xrot*pi_180)*vi.z);
	iderivate[2] = (-sin(m_Trans.zrot*pi_180)*cos(m_Trans.yrot*pi_180)*vi.x)-
		           (sin(m_Trans.zrot*pi_180)*sin(m_Trans.yrot*pi_180)*sin(m_Trans.xrot*pi_180)*vi.y)-
				   (cos(m_Trans.zrot*pi_180)*cos(m_Trans.xrot*pi_180)*vi.y)-
				   (sin(m_Trans.zrot*pi_180)*sin(m_Trans.yrot*pi_180)*cos(m_Trans.xrot*pi_180)*vi.z)+
				   (cos(m_Trans.zrot*pi_180)*sin(m_Trans.xrot*pi_180)*vi.z);
	iderivate[3] = (sin(m_Trans.zrot*pi_180)*sin(m_Trans.yrot*pi_180)*cos(m_Trans.xrot*pi_180)*vi.x)-
		           (cos(m_Trans.zrot*pi_180)*sin(m_Trans.xrot*pi_180)*vi.y)-
				   (sin(m_Trans.zrot*pi_180)*sin(m_Trans.yrot*pi_180)*sin(m_Trans.xrot*pi_180)*vi.z)-
				   (cos(m_Trans.zrot*pi_180)*cos(m_Trans.xrot*pi_180)*vi.z);
	iderivate[4] = (-sin(m_Trans.zrot*pi_180)*sin(m_Trans.yrot*pi_180)*vi.x)+
				   (sin(m_Trans.zrot*pi_180)*cos(m_Trans.yrot*pi_180)*sin(m_Trans.xrot*pi_180)*vi.y)+
				   (sin(m_Trans.zrot*pi_180)*cos(m_Trans.yrot*pi_180)*cos(m_Trans.xrot*pi_180)*vi.z);
	iderivate[5] = (cos(m_Trans.zrot*pi_180)*cos(m_Trans.yrot*pi_180)*vi.x)+
				   (cos(m_Trans.zrot*pi_180)*sin(m_Trans.yrot*pi_180)*sin(m_Trans.xrot*pi_180)*vi.y)-
				   (sin(m_Trans.zrot*pi_180)*cos(m_Trans.xrot*pi_180)*vi.y)+
				   (cos(m_Trans.zrot*pi_180)*sin(m_Trans.yrot*pi_180)*cos(m_Trans.xrot*pi_180)*vi.z)+
				   (sin(m_Trans.zrot*pi_180)*sin(m_Trans.xrot*pi_180)*vi.z);
	iderivate[6] = (cos(m_Trans.yrot*pi_180)*cos(m_Trans.xrot*pi_180)*vi.y)-
		           (cos(m_Trans.yrot*pi_180)*sin(m_Trans.xrot*pi_180)*vi.z);
	iderivate[7] = (-cos(m_Trans.yrot*pi_180)*vi.x)-
		           (sin(m_Trans.yrot*pi_180)*sin(m_Trans.xrot*pi_180)*vi.y)-
				   (sin(m_Trans.yrot*pi_180)*cos(m_Trans.xrot*pi_180)*vi.z);
	iderivate[8] = 0.0;

	double	jderivate[9];
	jderivate[0] = (cos(m_Trans.zrot*pi_180)*sin(m_Trans.yrot*pi_180)*cos(m_Trans.xrot*pi_180)*vj.y) +
		           (sin(m_Trans.zrot*pi_180)*sin(m_Trans.xrot*pi_180)*vj.y) -
				   (cos(m_Trans.zrot*pi_180)*sin(m_Trans.yrot*pi_180)*sin(m_Trans.xrot*pi_180)*vj.z) +
				   (sin(m_Trans.zrot*pi_180)*cos(m_Trans.xrot*pi_180)*vj.z);
	jderivate[1] = (-cos(m_Trans.zrot*pi_180)*sin(m_Trans.yrot*pi_180)*vj.x)+
		           (cos(m_Trans.zrot*pi_180)*cos(m_Trans.yrot*pi_180)*sin(m_Trans.xrot*pi_180)*vj.y)+
				   (cos(m_Trans.zrot*pi_180)*cos(m_Trans.yrot*pi_180)*cos(m_Trans.xrot*pi_180)*vj.z);
	jderivate[2] = (-sin(m_Trans.zrot*pi_180)*cos(m_Trans.yrot*pi_180)*vj.x)-
		           (sin(m_Trans.zrot*pi_180)*sin(m_Trans.yrot*pi_180)*sin(m_Trans.xrot*pi_180)*vj.y)-
				   (cos(m_Trans.zrot*pi_180)*cos(m_Trans.xrot*pi_180)*vj.y)-
				   (sin(m_Trans.zrot*pi_180)*sin(m_Trans.yrot*pi_180)*cos(m_Trans.xrot*pi_180)*vj.z)+
				   (cos(m_Trans.zrot*pi_180)*sin(m_Trans.xrot*pi_180)*vj.z);
	jderivate[3] = (sin(m_Trans.zrot*pi_180)*sin(m_Trans.yrot*pi_180)*cos(m_Trans.xrot*pi_180)*vj.x)-
		           (cos(m_Trans.zrot*pi_180)*sin(m_Trans.xrot*pi_180)*vj.y)-
				   (sin(m_Trans.zrot*pi_180)*sin(m_Trans.yrot*pi_180)*sin(m_Trans.xrot*pi_180)*vj.z)-
				   (cos(m_Trans.zrot*pi_180)*cos(m_Trans.xrot*pi_180)*vj.z);
	jderivate[4] = (-sin(m_Trans.zrot*pi_180)*sin(m_Trans.yrot*pi_180)*vj.x)+
				   (sin(m_Trans.zrot*pi_180)*cos(m_Trans.yrot*pi_180)*sin(m_Trans.xrot*pi_180)*vj.y)+
				   (sin(m_Trans.zrot*pi_180)*cos(m_Trans.yrot*pi_180)*cos(m_Trans.xrot*pi_180)*vj.z);
	jderivate[5] = (cos(m_Trans.zrot*pi_180)*cos(m_Trans.yrot*pi_180)*vj.x)+
				   (cos(m_Trans.zrot*pi_180)*sin(m_Trans.yrot*pi_180)*sin(m_Trans.xrot*pi_180)*vj.y)-
				   (sin(m_Trans.zrot*pi_180)*cos(m_Trans.xrot*pi_180)*vj.y)+
				   (cos(m_Trans.zrot*pi_180)*sin(m_Trans.yrot*pi_180)*cos(m_Trans.xrot*pi_180)*vj.z)+
				   (sin(m_Trans.zrot*pi_180)*sin(m_Trans.xrot*pi_180)*vj.z);
	jderivate[6] = (cos(m_Trans.yrot*pi_180)*cos(m_Trans.xrot*pi_180)*vj.y)-
		           (cos(m_Trans.yrot*pi_180)*sin(m_Trans.xrot*pi_180)*vj.z);
	jderivate[7] = (-cos(m_Trans.yrot*pi_180)*vj.x)-
		           (sin(m_Trans.yrot*pi_180)*sin(m_Trans.xrot*pi_180)*vj.y)-
				   (sin(m_Trans.yrot*pi_180)*cos(m_Trans.xrot*pi_180)*vj.z);
	jderivate[8] = 0.0;

	(*result)[0] = igradient[0] - jgradient[0];
	(*result)[1] = igradient[1] - jgradient[1];
	(*result)[2] = igradient[2] - jgradient[2];
	(*result)[3] = (igradient[0]*iderivate[0]+igradient[1]*iderivate[3]+igradient[2]*iderivate[6])-
		           (jgradient[0]*jderivate[0]+jgradient[1]*jderivate[3]+jgradient[2]*jderivate[6]);
	(*result)[4] = (igradient[0]*iderivate[1]+igradient[1]*iderivate[4]+igradient[2]*iderivate[7])-
		           (jgradient[0]*jderivate[1]+jgradient[1]*jderivate[4]+jgradient[2]*jderivate[7]);
	(*result)[5] = (igradient[0]*iderivate[2]+igradient[1]*iderivate[5]+igradient[2]*iderivate[8])-
		           (jgradient[0]*jderivate[2]+jgradient[1]*jderivate[5]+jgradient[2]*jderivate[8]);

	return TRUE;
}

//*********************************************************************
//	Gauss_point(): calculate univariate gaussian density function
//
//		v1: intensity of transformed coordinate v1
//		v2: intensity of transformed coordinate v2
//		variance: variance of gaussian density function
//*********************************************************************
double RxDerivateMI3D::Gauss_point(unsigned short v1, unsigned short v2, double variance)
{
	double	density = 0.0, dist = 0.0, factor = 1.0;
	factor	*= PI_2;
	factor	= 1.0/((double)sqrt((double)(factor*variance)));
	dist	= (double)((v1-v2)*(v1-v2))/variance;
	density	= (double)(factor * (double)exp((double)(-0.5)*dist));

	if(density < 0.0000005 && density > -0.0000005)
		density = 0.0;

	return density;
}

//*********************************************************************
//	Gauss_pair(): calculate bivariate gaussian density function
//
//		u1, v1: intensity of u1, v1
//		u2, v2: intensity of u2, v2
//		cov1: covariance of reference volume
//		cov2: covariance of test volume
//*********************************************************************
double RxDerivateMI3D::Gauss_pair(unsigned short u1,unsigned short v1,unsigned short u2,unsigned short v2,double cov1,double cov2)
{
	// Calculate first part
	double	determ	= cov1 * cov2;
	double	factor	= 1.0;
	factor	*= PI_2;
	factor	= (double)(1.0/((double)sqrt((double)(factor*determ))));

	// Generate transpose vector variable of (wi-wj)
	double	vector1[2], vector2[2], vector[2];
	double	matrix[2][2];
	double	density = 0.0, dist = 0.0;
	vector1[0]		= vector2[0] = (double)(u1-v1);
	vector1[1]		= vector2[1] = (double)(u2-v2);
	matrix[0][0]	= 1.0 / cov1;
	matrix[0][1]	= 0.0;
	matrix[1][0]	= 0.0;
	matrix[1][1]	= 1.0 / cov2;

	m_pMatrix3D->matrix22_vector21(matrix,vector2,vector);
	m_pMatrix3D->vector12_vector21(vector1,vector,&dist);

	// Calculate gaussian density
	density = factor * (double)exp((double)((-0.5)*dist));

	if(density < 0.0000005 && density > -0.0000005) 
		density = 0.0;

	return density;
}

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