自学教程：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 vectorvoid 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 vorticitiesintInitializeTurbModes(	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;}

void SymmetryMod::ModifyPolyObject (TimeValue t, ModContext &mc, PolyObject *pobj, INode *inode) {	MNMesh &mesh = pobj->GetMesh();	if (mUseRampageWeldMath)		mesh.SetFlag(MN_MESH_USE_MAX2012_WELD_MATH,TRUE);	// Luna task 747	// We cannot support specified normals in Symmetry at this time.	mesh.ClearSpecifiedNormals();	Interval iv = FOREVER;	int axis, slice, weld, flip;	float thresh;	mp_pblock->GetValue (kSymAxis, t, axis, iv);	mp_pblock->GetValue (kSymFlip, t, flip, iv);	mp_pblock->GetValue (kSymSlice, t, slice, iv);	mp_pblock->GetValue (kSymWeld, t, weld, iv);	mp_pblock->GetValue (kSymThreshold, t, thresh, iv);	if (thresh<0) thresh=0;	// Get transform from mp_mirror controller:	Matrix3 tm  = CompMatrix (t, NULL, &mc, &iv);	Matrix3 itm = Inverse (tm);	// Get DotProd(N,x)=off plane definition from transform	Point3 Axis(0,0,0);	Axis[axis] = flip ? -1.0f : 1.0f;	Point3 or = tm.GetTrans();	Point3 N = Normalize(tm*Axis - or);	float off = DotProd (N, or);	if (slice) SlicePolyObject (mesh, N, off);	MirrorPolyObject (mesh, axis, tm, itm);	if (weld) {		WeldPolyObject (mesh, N, off, thresh);		RemovePolySpurs (mesh);	}	pobj->UpdateValidity (GEOM_CHAN_NUM, iv);	pobj->UpdateValidity (TOPO_CHAN_NUM, iv);	pobj->UpdateValidity (VERT_COLOR_CHAN_NUM, iv);	pobj->UpdateValidity (TEXMAP_CHAN_NUM, iv);	pobj->UpdateValidity (SELECT_CHAN_NUM, iv);}

static void DisplayPlaneAngle(HWND hDlg, int id, Point3 &dirPt, Point3 &planePt){    float len = 0.0f, cosAng = 0.0f;    TCHAR buf[32];    len = Length(planePt);    if(len != 0) {        cosAng = DotProd(planePt, dirPt) / len;        if(cosAng > 0.99999f)   // beyond float accuracy!            SetDlgItemText(hDlg, id, _T("0"));        else {            _stprintf(buf, _T("%g"), acos((double)cosAng) * RadToDegDbl);            SetDlgItemText(hDlg, id, buf);        }    }    else        SetDlgItemText(hDlg, id, _T("90"));}

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;  #pragma omp parallel for default(shared) private(i)   for (i = 0; i < NUM_SPECIES; i++)    ratesxy[i] = cxy[i] * ( bcoef[i] * fac + ratesxy[i] );  }

void CVolume3D::createSphereVolume(const int nx, const float r){	int i, j, k;	//realeas old space	_free();	m_nx = m_ny = m_nz = nx;	m_nDensityVolumeType = GL_FLOAT;	//density volume allocation	const int nsize = nx*nx*nx;	m_pDensityVolume = new float [nsize];	assert(m_pDensityVolume!=NULL);	float *p = m_pDensityVolume;	m_pGradientVolume = new Vector4f[nsize];	assert(m_pGradientVolume!=NULL);	Vector4f *pgrad = m_pGradientVolume;	//assign dist val using non-binary voxelization	const Vector3f h(0.5f, 0.5f, 0.5f);	const float cc = (nx-1)*0.5f;	const Vector3f center(cc, cc, cc);	for (k=0; k<nx; k++){		for (j=0; j<nx; j++){			for (i=0; i<nx; i++, p++, pgrad++){				Vector3f v(i, j, k);				Vector3f dist=v-center;				float d2 = DotProd(dist, dist);				const float d = (float)sqrt(d2);				const float d1 = 0.5f/(d+1e-6f);    //make sure none-zero				float u=dist2Density(d, r);				*p  = u;				const Vector3f N=dist*d1+h;				pgrad->x = N.x;				pgrad->y = N.y;				pgrad->z = N.z;				pgrad->w=1;			}		}	}}

Point3 UVW_ChannelClass::UVFaceNormal(int index){	if (index < 0) return Point3(0.0f,0.0f,1.0f);	if (index >= f.Count()) return Point3(0.0f,0.0f,1.0f);	if (f[index]->count < 3) 		return Point3(0.0f,0.0f,1.0f);	Point3 vec1,vec2;	if (f[index]->count == 3)	{		int a = f[index]->t[0];		int b = f[index]->t[1];		int c = f[index]->t[2];		vec1 = Normalize(v[b].GetP()-v[a].GetP());		vec2 = Normalize(v[c].GetP()-v[a].GetP());	}	else	{		int a = f[index]->t[0];		int b = f[index]->t[1];				vec1 = Normalize(v[b].GetP()-v[a].GetP());		for (int i = 2; i < f[index]->count; i++)		{			b = f[index]->t[i];			vec2 = Normalize(v[b].GetP()-v[a].GetP());			float dot = DotProd(vec1,vec2);			if (fabs(dot) != 1.0f) 				i = f[index]->count;		}	}	Point3 norm = CrossProd(vec1,vec2);	return Normalize(norm);}

static void FindVertexAngles(PatchMesh &pm, float *vang) {	int i;	for (i=0; i<pm.numVerts + pm.numVecs; i++) vang[i] = 0.0f;	for (i=0; i<pm.numPatches; i++) {		Patch &p = pm.patches[i];		for (int j=0; j<p.type; j++) {			Point3 d1 = pm.vecs[p.vec[j*2]].p - pm.verts[p.v[j]].p;			Point3 d2 = pm.vecs[p.vec[((j+p.type-1)%p.type)*2+1]].p - pm.verts[p.v[j]].p;			float len = LengthSquared(d1);			if (len == 0) continue;			d1 /= Sqrt(len);			len = LengthSquared (d2);			if (len==0) continue;			d2 /= Sqrt(len);			float cs = DotProd (d1, d2);			if (cs>=1) continue;	// angle of 0			if (cs<=-1) vang[p.v[j]] += PI;			else vang[p.v[j]] += (float) acos (cs);		}	}}

inline void _computeReflectedLeastSquareMatrix(	const int nv, 	const double *_P, const int strideP, 	const double *_Q, const int strideQ, 	const double *_Weight, const int strideW, 	const Vector3d &n,	double3x3 &A){	A.setZeroMatrix();	for (int i=0; i<nv; i++){		Vector3d p = _getVertex3dUsingStride(_P, i, strideP);		Vector3d q = _getVertex3dUsingStride(_Q, i, strideQ);		const double w = _getDoubleUsingStride(_Weight, i, strideW);		Vector3d &s = q;		const double dist = DotProd(s, n)*2.0;		s -= dist*n;		double3x3 r; 		VectorTensorProduct(p*w, q, r);		A += r;	}}

// If |projGridPt - GridPt| < gridCubeRadius// and probBitmap->GetPixel(src->UVW(bary)) < RandomZeroToOne()// Also a generic random factor.bool plDistributor::IProbablyDoIt(int iFace, Point3& del, const Point3& bary) const{    if( fRand.RandZeroToOne() >= fOverallProb )    {        return false;    }    if( (kIsoNone == fIsolation) || (kIsoLow == fIsolation) )    {        if( LengthSquared(del) >= fSpacing*fSpacing )        {            return false;        }        Point3 faceNorm = fSurfMesh->FaceNormal(iFace, FALSE);        if( DotProd(del, faceNorm) < 0 )        {            return false;        }    }    if( IFailsAngProb(iFace, bary) )    {        return false;    }    if( IFailsAltProb(iFace, bary) )    {        return false;    }    if( IFailsProbBitmap(iFace, bary) )    {        return false;    }    return true;}

//=================================================================// Methods for DumpNodesTEP//int DumpNodesTEP::callback(INode *pnode){	ASSERT_MBOX(!(pnode)->IsRootNode(), "Encountered a root node!");	if (::FNodeMarkedToSkip(pnode))		return TREE_CONTINUE;	// Get node's parent	INode *pnodeParent;	pnodeParent = pnode->GetParentNode();		// The model's root is a child of the real "scene root"	TSTR strNodeName(pnode->GetName());	BOOL fNodeIsRoot = pnodeParent->IsRootNode( );		int iNode = ::GetIndexOfINode(pnode);	int iNodeParent = ::GetIndexOfINode(pnodeParent, !fNodeIsRoot/*fAssertPropExists*/);	// Convenient time to cache this	m_phec->m_rgmaxnode[iNode].imaxnodeParent = fNodeIsRoot ? SmdExportClass::UNDESIRABLE_NODE_MARKER : iNodeParent;	// Root node has no parent, thus no translation	if (fNodeIsRoot)		iNodeParent = -1;			// check to see if the matrix isn't right handed	m_phec->m_rgmaxnode[iNode].isMirrored = DotProd( CrossProd( m_phec->m_rgmaxnode[iNode].mat3ObjectTM.GetRow(0).Normalize(), m_phec->m_rgmaxnode[iNode].mat3ObjectTM.GetRow(1).Normalize() ).Normalize(), m_phec->m_rgmaxnode[iNode].mat3ObjectTM.GetRow(2).Normalize() ) < 0;	// Dump node description	fprintf(m_pfile, "%3d /"%s/" %3d/n", 		iNode, 		strNodeName, 		iNodeParent );	return TREE_CONTINUE;}

void LUFactorSymm(const IX neq, R8 **a){  IX i, j;  R8 dot, tmp;  a[1][1] = 1.0 / a[1][1];  for(i=2; i<=neq; i++) { /* process column i */    for(j=2; j<i; j++) {      a[i][j] -= DotProd(j-1, a[i], a[j]);    }    /* process diagonal i */    for(dot=0.0, j=1; j<i; j++) {      tmp = a[i][j] * a[j][j];      dot += a[i][j] * tmp;      a[i][j] = tmp;    }    a[i][i] -= dot;    if(a[i][i] == 0.0) {      error(3, __FILE__, __LINE__, "Zero on the diagonal, row %d", i);    }    a[i][i] = 1.0 / a[i][i];  }  /*  end of i loop  */}   /*  end of LUFactorSymm  */

// Determine is the node has negative scaling.// This is used for mirrored objects for example. They have a negative scale factor// so when calculating the normal we should take the vertices counter clockwise.// If we don't compensate for this the objects will be 'inverted'.BOOL Exporter::TMNegParity(Matrix3 &m){	return (DotProd(CrossProd(m.GetRow(0),m.GetRow(1)),m.GetRow(2))<0.0)?1:0;}