近日在写一个Parser,用到了Matrix变换,一直以来都是使用库函数,然此次要求不得使用任何函数库(当然了,基本库还是可以用的),而且必须用C语言写~ 难啊,小弟不才,Google,baidu,翻了个遍,没有太合适的,而且把线形代数的原理忘得一干二净...
感谢Google,感谢Graphics Gems I,感谢所有在网络上瞎发评论的同仁们... Just a temporary version updated from Craphics Gems.
static
double
det2x2(
double
a,
double
b,
double
c,
double
d)
... {
double ans = a * d - b * c;
return ans;
}
static double det3x3( double * ctm) ... {
double ans;
ans = ctm[0] * det2x2(ctm[4], ctm[7], ctm[5], ctm[8])
- ctm[1] * det2x2(ctm[3], ctm[6], ctm[5], ctm[8]);
+ ctm[2] * det2x2(ctm[3], ctm[6], ctm[4], ctm[7]);
return ans;
}
static void adjoint( double * ctm1, double * ctm2)
... {
ctm2[0] = det2x2(ctm1[4], ctm1[7], ctm1[5], ctm1[8]);
ctm2[1] = det2x2(ctm1[3], ctm1[6], ctm1[5], ctm1[8]);
ctm2[2] = det2x2(ctm1[3], ctm1[6], ctm1[4], ctm1[7]);
ctm2[3] = det2x2(ctm1[1], ctm1[7], ctm1[2], ctm1[8]);
ctm2[4] = det2x2(ctm1[0], ctm1[6], ctm1[2], ctm1[8]);
ctm2[5] = det2x2(ctm1[0], ctm1[6], ctm1[1], ctm1[7]);
ctm2[6] = det2x2(ctm1[1], ctm1[4], ctm1[2], ctm1[5]);
ctm2[7] = det2x2(ctm1[0], ctm1[3], ctm1[2], ctm1[5]);
ctm2[8] = det2x2(ctm1[0], ctm1[3], ctm1[1], ctm1[4]);
//Need to add sign later...
}
static void matrixToPhalanx( double * ctm, double * dPhalanx)
... {
dPhalanx[0] = ctm[0];
dPhalanx[1] = ctm[1];
dPhalanx[2] = 0;
dPhalanx[3] = ctm[2];
dPhalanx[4] = ctm[3];
dPhalanx[5] = 0;
dPhalanx[6] = ctm[4];
dPhalanx[7] = ctm[5];
dPhalanx[8] = 1;
}
static void phalanxToMatrix( double * dPhalanx, double * ctm)
... {
double dtmp1 = dPhalanx[8];
double dtmp2 = dPhalanx[2];
dPhalanx[0] = dPhalanx[0] * dtmp1 - dPhalanx[6] * dtmp2;
dPhalanx[1] = dPhalanx[1] * dtmp1 - dPhalanx[7] * dtmp2;
dPhalanx[2] = 0;
//dPhalanx[2] = dPhalanx[2] * dtmp1 - dPhalanx[8] * dtmp2;
dtmp2 = dPhalanx[5];
dPhalanx[3] = dPhalanx[3] * dtmp1 - dPhalanx[6] * dtmp2;
dPhalanx[4] = dPhalanx[4] * dtmp1 - dPhalanx[7] * dtmp2;
dPhalanx[5] = 0;
dPhalanx[6] = dPhalanx[6] / dtmp1;
dPhalanx[7] = dPhalanx[7] / dtmp1;
dPhalanx[8] = 1;
//dPhalanx[8] = dPhalanx[8] / dtmp;
ctm[0] = dPhalanx[0];
ctm[1] = dPhalanx[1];
ctm[2] = dPhalanx[3];
ctm[3] = dPhalanx[4];
ctm[4] = dPhalanx[6];
ctm[5] = dPhalanx[7];
}
static void matrixInversion( double * ctm1, double * ctm2) ... {
//
// inverse matrix
int i = 0;
double det;
double dPhalanx1[9], dPhalanx2[9];
matrixToPhalanx(ctm1, dPhalanx1);
adjoint(dPhalanx1, dPhalanx2);
det = det3x3(dPhalanx2);
if( fabs( det ) < 0.0000001)
...{
//error...
}
for( i=0; i<6; i++)
dPhalanx2[i] = dPhalanx2[i] / det;
phalanxToMatrix(dPhalanx2,ctm2);
}
... {
double ans = a * d - b * c;
return ans;
}
static double det3x3( double * ctm) ... {
double ans;
ans = ctm[0] * det2x2(ctm[4], ctm[7], ctm[5], ctm[8])
- ctm[1] * det2x2(ctm[3], ctm[6], ctm[5], ctm[8]);
+ ctm[2] * det2x2(ctm[3], ctm[6], ctm[4], ctm[7]);
return ans;
}
static void adjoint( double * ctm1, double * ctm2)
... {
ctm2[0] = det2x2(ctm1[4], ctm1[7], ctm1[5], ctm1[8]);
ctm2[1] = det2x2(ctm1[3], ctm1[6], ctm1[5], ctm1[8]);
ctm2[2] = det2x2(ctm1[3], ctm1[6], ctm1[4], ctm1[7]);
ctm2[3] = det2x2(ctm1[1], ctm1[7], ctm1[2], ctm1[8]);
ctm2[4] = det2x2(ctm1[0], ctm1[6], ctm1[2], ctm1[8]);
ctm2[5] = det2x2(ctm1[0], ctm1[6], ctm1[1], ctm1[7]);
ctm2[6] = det2x2(ctm1[1], ctm1[4], ctm1[2], ctm1[5]);
ctm2[7] = det2x2(ctm1[0], ctm1[3], ctm1[2], ctm1[5]);
ctm2[8] = det2x2(ctm1[0], ctm1[3], ctm1[1], ctm1[4]);
//Need to add sign later...
}
static void matrixToPhalanx( double * ctm, double * dPhalanx)
... {
dPhalanx[0] = ctm[0];
dPhalanx[1] = ctm[1];
dPhalanx[2] = 0;
dPhalanx[3] = ctm[2];
dPhalanx[4] = ctm[3];
dPhalanx[5] = 0;
dPhalanx[6] = ctm[4];
dPhalanx[7] = ctm[5];
dPhalanx[8] = 1;
}
static void phalanxToMatrix( double * dPhalanx, double * ctm)
... {
double dtmp1 = dPhalanx[8];
double dtmp2 = dPhalanx[2];
dPhalanx[0] = dPhalanx[0] * dtmp1 - dPhalanx[6] * dtmp2;
dPhalanx[1] = dPhalanx[1] * dtmp1 - dPhalanx[7] * dtmp2;
dPhalanx[2] = 0;
//dPhalanx[2] = dPhalanx[2] * dtmp1 - dPhalanx[8] * dtmp2;
dtmp2 = dPhalanx[5];
dPhalanx[3] = dPhalanx[3] * dtmp1 - dPhalanx[6] * dtmp2;
dPhalanx[4] = dPhalanx[4] * dtmp1 - dPhalanx[7] * dtmp2;
dPhalanx[5] = 0;
dPhalanx[6] = dPhalanx[6] / dtmp1;
dPhalanx[7] = dPhalanx[7] / dtmp1;
dPhalanx[8] = 1;
//dPhalanx[8] = dPhalanx[8] / dtmp;
ctm[0] = dPhalanx[0];
ctm[1] = dPhalanx[1];
ctm[2] = dPhalanx[3];
ctm[3] = dPhalanx[4];
ctm[4] = dPhalanx[6];
ctm[5] = dPhalanx[7];
}
static void matrixInversion( double * ctm1, double * ctm2) ... {
//
// inverse matrix
int i = 0;
double det;
double dPhalanx1[9], dPhalanx2[9];
matrixToPhalanx(ctm1, dPhalanx1);
adjoint(dPhalanx1, dPhalanx2);
det = det3x3(dPhalanx2);
if( fabs( det ) < 0.0000001)
...{
//error...
}
for( i=0; i<6; i++)
dPhalanx2[i] = dPhalanx2[i] / det;
phalanxToMatrix(dPhalanx2,ctm2);
}
还没完善,只是一个3*3矩阵的逆矩阵求法....... 还要改进,太慢了~ 不过似乎没有其它方法了~ 高斯消元?也许挺好,再说啦,此番先留个记。