[请注意，这个程序中的求逆的那个部分的程序的健壮性不强，无法处理a[0][0] = 0 的情形，当遇到这种情形的时候，会发生溢出的情形，具体的改进，请看下一篇博文]

几个特殊点的测试是成功。求逆用的是一个网友的一部分程序。

很粗糙，没有整理，讲究看吧。

其实，程序的结构可以大大的改进。有时间在说了

http://s4/middle/9151e730xbfc3581471f3&690—C++实现" TITLE="世界坐标系到局部坐标系的转换程序 —C++实现" />



http://s12/middle/9151e730xbfc35b1936bb&690—C++实现" TITLE="世界坐标系到局部坐标系的转换程序 —C++实现" />

http://s3/middle/9151e730xbfc35af40302&690—C++实现" TITLE="世界坐标系到局部坐标系的转换程序 —C++实现" />

#include <iostream>
#include <math.h>

using namespace std;

#define n 4

class node
{
 public :
  node(double c,double b,double w):x(c),y(b),z(w){}
  node()
  {
   x = 0;

  }
  node(node &p)
  {
   x = p.x;
   y = p.y;
   z = p.z;
  }
  void print()
  {
   cout << x << ends << y << ends << z << endl;
  }
  double x;
  double y;
  double z;
};

node caculate(node &a,node &b)
{
 node c;
 c.x = a.y*b.z-a.z*b.y;
 c.y = a.z*b.x-a.x*b.z;
 c.z = a.x*b.y-a.y*b.x;
 return c;
}

double total(node &a)
{
 return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);
};

node nodeave(node a,double m)
{
 a.x = a.x/m;
 a.y = a.y/m;
 a.z = a.z/m;
 return a;
}

int main()
{
 //待测试的坐标点
 node test(0,0,0);
 //假设平面的坐标点，假设p1为局部坐标系的原点
 node p1(0,0,1);
 node p2(1,0,1);
 node p3(0,1,1); 
 //求局部X坐标
 node p12(p2.x-p1.x,p2.y-p1.y,p2.z-p1.z);
 node p13(p3.x-p1.x,p3.y-p1.y,p3.z-p1.z);
 node X = nodeave(p12,total(p12));
 //求局部Y坐标
 node Z = nodeave(caculate(p12,p13),total(caculate(p12,p13)));
 //求局部Y坐标
 node Y = caculate(X,Z);

 //一下的代码用于求逆矩阵
 int i,j,k;double m;
 double a[n][n],E[n][n];

 a[0][0]=X.x; a[0][1] = X.y; a[0][2] = X.z; a[0][3] = 0;
 a[1][0]=Y.x; a[1][1] = Y.y; a[1][2] = Y.z; a[1][3] = 0;
 a[2][0]=Z.x; a[2][1] = Z.y; a[2][2] = Z.z; a[2][3] = 0;
 a[3][0]=p1.x; a[3][1] = p1.y; a[3][2] = p1.z; a[3][3] = 1;

 for(int i = 0 ; i < n ; ++i)
 {
  for(int j = 0; j < n;++j)
   cout << a[i][j] << ends;
  cout <<  endl;
 }

 //单位矩阵E[n][n]
 for(i=0;i<n;i++){
  for(j=0;j<n;j++){
   if(i==j)
    E[i][j]=1;
   else
    E[i][j]=0;
  }
 }

 //上三角变换
 for(k=0;k<n-1;k++){
  for(i=k+1;i<n;i++){
   m=a[i][k]/a[k][k];
   for(j=0;j<n;j++){
    a[i][j]=a[i][j]-m*a[k][j];
    E[i][j]=E[i][j]-m*E[k][j];
   }
  }
 }
 //下三角变换
 for(k=n-1;k>0;k--){
  for(i=k-1;i>=0;i--){
   m=a[i][k]/a[k][k];
   for(j=0;j<n;j++){
    a[i][j]=a[i][j]-m*a[k][j];
    E[i][j]=E[i][j]-m*E[k][j];
   }
  }
 }
 //单位矩阵……
 for(i=0;i<n;i++)
  for(j=0;j<n;j++){
   E[i][j]=E[i][j]/a[i][i];
  }

 node result;
 result.x = E[0][0]*test.x + E[1][0]*test.y+E[2][0]*test.z+E[3][0];
 result.y = E[0][1]*test.x + E[1][1]*test.y+E[2][1]*test.z+E[3][1];
 result.z = E[0][2]*(test.x) + E[1][2]*(test.y)+E[2][2]*(test.z)+E[3][2];

 cout << result.x  << ends <<result.y << ends <<result.z << endl;
 return 0;
}