C++ DotProd函數代碼示例

void InterpolatorKNN::CalculateBaryzentricCoordinates(NNCandidate candidates[3], Vertex v) {
	D3DXVECTOR3 v_01 = Normalize(candidates[1].vertex.pos-candidates[0].vertex.pos);
	D3DXVECTOR3 v_02 = Normalize(candidates[2].vertex.pos-candidates[0].vertex.pos);	
	
	// project v onto plane <v_01, v_02> with hesse normal form
	D3DXVECTOR3 normal = Normalize(CrossProd(v_01, v_02));
	D3DXVECTOR3 pointOnPlane = candidates[0].vertex.pos;
	float d = DotProd(normal, pointOnPlane);
	float distance = DotProd(normal, v.pos) - d;
	v.pos = v.pos + distance * (-normal);
	
	float translate = 0.0f;
	if(v.pos.x == 0 && v.pos.y == 0 && v.pos.z == 0) {
		translate = 10;
	}
		
	candidates[0].vertex.pos.x += translate;
	candidates[1].vertex.pos.x += translate;
	candidates[2].vertex.pos.x += translate;
	v.pos.x += translate;
	
	D3DXVECTOR3 res = SolveLES( candidates[0].vertex.pos, 
															candidates[1].vertex.pos, 
															candidates[2].vertex.pos, 
															v.pos);
	
	candidates[0].weight = res.x;
	candidates[1].weight = res.y;
	candidates[2].weight = res.z;

	D3DXVECTOR3 check = candidates[0].weight * candidates[0].vertex.pos + 
											candidates[1].weight * candidates[1].vertex.pos +
											candidates[2].weight * candidates[2].vertex.pos -
											v.pos;

	float error = abs(check.x) + abs(check.y) + abs(check.z);
	if(error > 0.00001) {
		PD(L"big error solving lgs: ", error);
	}
}

inline void computeStiffnessMat(const double3x3& strainDeriv, const double& K, double3x3 *ppJacobian[3])
{
    const Vector3d *a = (const Vector3d *)&strainDeriv.x[0];
    const Vector3d *b = (const Vector3d *)&strainDeriv.x[3];
    const Vector3d *c = (const Vector3d *)&strainDeriv.x[6];
    const double aa = DotProd(*a, *a);
    const double ab = DotProd(*a, *b);
    const double ac = DotProd(*a, *c);
    const double bb = DotProd(*b, *b);
    const double bc = DotProd(*b, *c);
    const double cc = DotProd(*c, *c);
    const double S11 = aa*K;
    const double S12 = (ab+ac)*K;
    const double S22 = (bb+cc+bc+bc)*K;

    ppJacobian[0]->ZeroMatrix();
    ppJacobian[1]->ZeroMatrix();
    ppJacobian[2]->ZeroMatrix();
    double *x;
    x=ppJacobian[0]->x;
    x[0]=x[4]=x[8]=S11;
    x=ppJacobian[1]->x;
    x[0]=x[4]=x[8]=S12;
    x=ppJacobian[2]->x;
    x[0]=x[4]=x[8]=S22;
}

static Point3 RefractVector(ShadeContext &sc, Point3 N, Point3 V, float ior) { 
	float VN,nur,k1;
	VN = DotProd(-V,N);
	if (sc.backFace) nur = ior;
	else nur = (ior!=0.0f) ? 1.0f/ior: 1.0f;
	k1 = 1.0f-nur*nur*(1.0f-VN*VN);
	if (k1<=0.0f) {
		// Total internal reflection: 
		return FNormalize(2.0f*VN*N + V);
		}
	else 
		return (nur*VN-(float)sqrt(k1))*N + nur*V;
	}

//if damping is required, we compute additional damping force in the 
//returned force vector
void CSpringElement::computeNodalForce( 
	const Vector3d &p0,			//node 0's world position
	const Vector3d &p1,			//node 1's world position
	const Vector3d &v0,			//node 0's velocity
	const Vector3d &v1,			//node 1's velocity
	Vector3d &f0,				//elastic force for node 0
	double3x3 *pstiffness)
{
	Vector3d Xij;
	_computeElasticForce(p0, p1, f0, Xij);

	if (m_Kd > 1e-10){
		Vector3d Vij = v0 - v1;
		const double dot1 = DotProd(Vij, Xij);
		const double dot2 = DotProd(Xij, Xij);
		Vector3d Fvisco = (m_Kd*dot1/dot2)*Xij;
		f0+=Fvisco;
	}
	if (pstiffness){

	}
}

int DGCubeMap::_getCloestAnchorVertex(const Vector3d &pos, const AnchorVertexInfo *pVertex, const int nv) const
{
	const Vector3d n = Normalize(pos); 
	double maxdotval = -100000.0;
	int index=-1;
	for (int i=0; i<nv; i++){
		double dotval = DotProd(n, pVertex[i].m_dir);
		if (maxdotval<dotval){
			maxdotval = dotval, index = i;
		}
	}
	return index;
}

void InterpolatorKNN::KollinearWithDistinctPositions(NNCandidate candidates[3], Vertex v) {
	// project v onto line <v_01, v_02> with hesse normal form	
	D3DXVECTOR3 v_01 = candidates[1].vertex.pos-candidates[0].vertex.pos;
	D3DXVECTOR3 v_02 = candidates[2].vertex.pos-candidates[0].vertex.pos;
	D3DXVECTOR3 v_0v = v.pos-candidates[0].vertex.pos;
	D3DXVECTOR3 v_1v = v.pos-candidates[1].vertex.pos;
	D3DXVECTOR3 v_2v = v.pos-candidates[2].vertex.pos;
	D3DXVECTOR3 v_2 = CrossProd(v_01, v_0v); 
	if(v_2.x != 0.0f || v_2.y != 0.0f || v_2.z != 0.0f) {
		D3DXVECTOR3 normal = Normalize(CrossProd(v_01, v_2));
		D3DXVECTOR3 pointOnPlane = candidates[0].vertex.pos;
		float d = DotProd(normal, pointOnPlane);
		float distance = DotProd(normal, v.pos) - d;
		v.pos = v.pos + distance * (-normal);
	}

	if(DotProd(v_01, v_02) < -0.9f) {
		// 0 is in the middle
		if(DotProd(v_01, v_0v) > 0.9f){
			// v on the same side as 1
			InterpolateLinear(candidates[0], candidates[1], v);
			candidates[2].weight = 0.0f;
		}
		else{
			// v on the same side as 2
			InterpolateLinear(candidates[0], candidates[2], v);
			candidates[1].weight = 0.0f;
		}
	}
	if(Length(v_01) < Length(v_02)) {
		// 1 is in the middle
		if(DotProd(-v_01, v_1v) > 0.9f){
			// v on the same side as 0
			InterpolateLinear(candidates[1], candidates[0], v);
			candidates[2].weight = 0.0f;
		}
		else{
			// v on the same side as 2
			InterpolateLinear(candidates[1], candidates[2], v);
			candidates[0].weight = 0.0f;
		}
	}
	else{
		// 2 is in the middle
		if(DotProd(-v_02, v_2v) > 0.9f){
			// v on the same side as 0
			InterpolateLinear(candidates[2], candidates[0], v);
			candidates[1].weight = 0.0f;
		}
		else{
			// v on the same side as 1
			InterpolateLinear(candidates[2], candidates[1], v);
			candidates[0].weight = 0.0f;
		}
	}
}

void InterpolatorKNN::AllSamePosition(NNCandidate candidates[3], Vertex v) {
	float weights[3];
	float sum = 0.0f;
	for(int i=0; i < 3; ++i) {
	weights[i] = max(0, DotProd(candidates[i].vertex.normal, v.normal));
		sum += weights[i];
	}
	if(sum == 0){
		sum = 3.0f;
		weights[0] = 1.0f; weights[1] = 1.0f; weights[2] = 1.0f;  
	}
	candidates[0].weight = weights[0] / sum;
	candidates[1].weight = weights[1] / sum;
	candidates[2].weight = weights[2] / sum;
}

float BerconGradient::getGradientValueDist(ShadeContext& sc) {
	switch (p_normalType) {	 
		case 0: { // View			 
			return -sc.P().z; //Length(sc.OrigView()); //(sc.PointTo(sc.P(), REF_CAMERA)).z;
		}
		case 1: { // Local X
			return (sc.PointTo(sc.P(), REF_OBJECT)).x;
		}
		case 2: { // Local Y
			return (sc.PointTo(sc.P(), REF_OBJECT)).y;
		}
		case 3: { // Local Z
			return (sc.PointTo(sc.P(), REF_OBJECT)).z;
		}
		case 4: { // World X
			return (sc.PointTo(sc.P(), REF_WORLD)).x;
		}
		case 5: { // World Y
			return (sc.PointTo(sc.P(), REF_WORLD)).y;
		}
		case 6: { // World Z
			return (sc.PointTo(sc.P(), REF_WORLD)).z;
		}
		case 7: { // Camera X
			return sc.P().x; //(sc.PointTo(sc.P(), REF_CAMERA)).x;
		}
		case 8: { // Camera Y
			return sc.P().y; //(sc.PointTo(sc.P(), REF_CAMERA)).y;
		}
		case 9: { // Camera Z
			return -sc.P().z; //-(sc.PointTo(sc.P(), REF_CAMERA)).z;
		}
		case 10: { // To Object
			if (sc.InMtlEditor() || !p_node)
				return -sc.P().z; //(sc.PointTo(sc.P(), REF_CAMERA)).z;			
			return Length((p_node->GetNodeTM(sc.CurTime())).GetTrans() - sc.PointTo(sc.P(), REF_WORLD));
		}
		case 11: { // Object Z			
			if (sc.InMtlEditor() || !p_node)
				return -sc.P().z; //(sc.PointTo(sc.P(), REF_CAMERA)).z;			
			Matrix3 tm = p_node->GetNodeTM(sc.CurTime());
			Point3 a = tm.GetTrans() - sc.PointTo(sc.P(), REF_WORLD);
			Point3 b = FNormalize(tm.GetRow(2));
			return (-DotProd(b, a) / Length(b));
		}
	}
	return 0.f;
}

void SymmetryMod::WeldTriObject (Mesh & mesh, Point3 & N, float offset, float threshold) {
	// Find vertices in target zone of mirror plane:
	BitArray targetVerts;
	targetVerts.SetSize (mesh.numVerts, true);
	targetVerts.ClearAll ();
	for (int i=0; i<mesh.numVerts; i++) {
		float dist = DotProd (N, mesh.verts[i]) - offset;
		if (fabsf(dist) > threshold) continue;
		targetVerts.Set (i);
	}

	// Weld the suitable border vertices:
	MeshDelta tmd(mesh);
	BOOL found = tmd.WeldByThreshold (mesh, targetVerts, threshold);
	tmd.Apply (mesh);
}

// Foley & vanDam: Computer Graphics: Principles and Practice, 
//     2nd Ed. pp 756ff.
Point3 SContext::RefractVector(float ior)
{
	Point3 N = Normal();
	float VN,nur,k;
	VN = DotProd(-viewDir,N);
	if (backFace) nur = ior;
	else nur = (ior!=0.0f) ? 1.0f/ior: 1.0f;
	k = 1.0f-nur*nur*(1.0f-VN*VN);
	if (k<=0.0f) {
		// Total internal reflection: 
		return ReflectVector();
	}
	else {
		return (nur*VN-(float)sqrt(k))*N + nur*viewDir;
	}
}

/*
====================
NormalAngle
====================
*/
static float NormalAngle( Point3 v1, Point3 v2 )
{
	float len = Length( v1 );
	if (len == 0) return 0;
	v1 /= len;

	len = Length( v2 );
	if (len == 0) return 0;
	v2 /= len;

	float normal_angle = DotProd( v1, v2 );

	normal_angle = min( 1.f, max( normal_angle, -1.f ) );

	return acosf( normal_angle );
}

/* Apply the method developed by Matthew Brown (see BMVC 02 paper) to
   fit a 3D quadratic function through the DOG function values around
   the location (s,r,c), i.e., (scale,row,col), at which a peak has
   been detected.  Return the interpolated peak position by setting
   the vector "offset", which gives offset from position (s,r,c).  The
   returned function value is the interpolated DOG magnitude at this peak.
*/
float FitQuadratic(float offset[3], Image *dogs, int s, int r, int c)
{
    float g[3], **dog0, **dog1, **dog2;
    static float **H = NULL;

    /* First time through, allocate space for Hessian matrix, H. */
    if (H == NULL)
      H = AllocMatrix(3, 3, PERM_POOL);

    /* Select the dog images at peak scale, dog1, as well as the scale
       below, dog0, and scale above, dog2.
    */
    dog0 = dogs[s-1]->pixels;
    dog1 = dogs[s]->pixels;
    dog2 = dogs[s+1]->pixels;

    /* Fill in the values of the gradient from pixel differences. */
    g[0] = (dog2[r][c] - dog0[r][c]) / 2.0;
    g[1] = (dog1[r+1][c] - dog1[r-1][c]) / 2.0;
    g[2] = (dog1[r][c+1] - dog1[r][c-1]) / 2.0;

    /* Fill in the values of the Hessian from pixel differences. */
    H[0][0] = dog0[r][c] - 2.0 * dog1[r][c] + dog2[r][c];
    H[1][1] = dog1[r-1][c] - 2.0 * dog1[r][c] + dog1[r+1][c];
    H[2][2] = dog1[r][c-1] - 2.0 * dog1[r][c] + dog1[r][c+1];
    H[0][1] = H[1][0] = ((dog2[r+1][c] - dog2[r-1][c]) -
			 (dog0[r+1][c] - dog0[r-1][c])) / 4.0;
    H[0][2] = H[2][0] = ((dog2[r][c+1] - dog2[r][c-1]) -
			 (dog0[r][c+1] - dog0[r][c-1])) / 4.0;
    H[1][2] = H[2][1] = ((dog1[r+1][c+1] - dog1[r+1][c-1]) -
			 (dog1[r-1][c+1] - dog1[r-1][c-1])) / 4.0;

    /* Solve the 3x3 linear sytem, Hx = -g.  Result gives peak offset.
       Note that SolveLinearSystem destroys contents of H.
    */
    offset[0] = - g[0];
    offset[1] = - g[1];
    offset[2] = - g[2];
    SolveLinearSystem(offset, H, 3); 

    /* Also return value of DOG at peak location using initial value plus
       0.5 times linear interpolation with gradient to peak position
       (this is correct for a quadratic approximation).  
    */
    return (dog1[r][c] + 0.5 * DotProd(offset, g, 3));
}

static void WebRate(realtype xx, realtype yy, realtype *cxy, realtype *ratesxy, 
                    void *user_data)
{
  long int i;
  realtype fac;
  UserData data;
  
  data = (UserData)user_data;
  
  for (i = 0; i<NUM_SPECIES; i++)
    ratesxy[i] = DotProd(NUM_SPECIES, cxy, acoef[i]);
  
  fac = ONE + ALPHA * xx * yy;
  
  for (i = 0; i < NUM_SPECIES; i++)
    ratesxy[i] = cxy[i] * ( bcoef[i] * fac + ratesxy[i] );  
}

/* const removal: const R8 **a */
void LUSolveSymm(const IX neq, R8 **a, R8 *b)
{
  IX i;

  for(i=2; i<=neq; i++) {    /*  forward substitution  */
    b[i] -= DotProd(i-1, a[i], b);
  }

  for(i=neq; i; i--) {
    b[i] *= a[i][i];
  }

  for(i=neq; i>1; i--) {     /*  back substitution  */
    DAXpY(i-1, -b[i], a[i], b);
  }

}  /*  end of LUSolveSymm  */

//TODO MDInitialize vorticities
int
InitializeTurbModes
(
	const MDConstants K,
	const Molecule* molecule,
	const TurbConstVecs* turb_vecs,
	const TurbConsts* turb,
	KraichnanMode** kraich_modes,
	const unsigned t
)
{

	//TODO profile & parallelise
	//TODO check valid input/output 
/*
#pragma omp parallel for default (none) \
shared (K.PartNum, molecule, turb_vecs, turb, K.Lx, K.Ly, K.Lz, t, delta_t, NF)\
schedule(dynamic, 5000)
*/
	for(unsigned i = 0; i < K.PartNum; ++i)
	{
		for(unsigned modeIndex = 0; modeIndex < K.NF; ++modeIndex)
		{
			double pos_vec[kDIM] = {0};

			NormalizeVector (K.L, molecule[i].position, pos_vec);
			double kn_dot_x = DotProd ( turb_vecs[modeIndex].kn, pos_vec);

			unsigned real_time = t * K.delta_t;                

			//in LaTeX would be:  \Omega_n = |\kappa| (\vec k \cdot \vec x) + |\omega| t
			double Omega_n = turb[modeIndex].k * kn_dot_x +
				turb[modeIndex].omega * real_time;

			//cos and sin
			kraich_modes[i][modeIndex].sin = sin (Omega_n);
			kraich_modes[i][modeIndex].cos = cos (Omega_n);
                       assert (kraich_modes[i][modeIndex].sin < 1);
                       assert (kraich_modes[i][modeIndex].sin > -1);
                       assert (kraich_modes[i][modeIndex].cos < 1);
                       assert (kraich_modes[i][modeIndex].cos > -1);
		}
	}
	return 0;
}